prime.h

#pragma once
 
#include <cmath>    // For std::sqrt()
#include <cstdlib>  // For randomness.
#include <cstdint>  // For fixed width integers.
 
namespace prime {
 
  inline bool is(const uint64_t& num) {
    if (num == 1) return false;
 
    // We only need to check up to the square root of the number in order to know if it's prime.
    auto root = static_cast<uint64_t>(std::sqrt(num)) + 1;
 
    // Start at 2, so we don't get a 0/1 false positive.
    for (size_t x = 2; x <= root; ++x) {
      if (num % x == 0)
        return false;
    }
    return true;
  }
 
 
  template <typename T = uint64_t> inline void next(T& num) {
    // Get to an odd number.
    if (num % 2 == 0) num++;
 
    // Loop until we find one.
    for (; !prime::is(num); num += 2) {}
  }
 
 
  inline uint64_t raise(uint64_t value, uint64_t exp, const uint64_t& mod) {
    uint64_t ret = 1;
 
    // Ensure it's bounded by the mod.
    value = value % mod;
 
    // March down the exponent until it's been 0d.
    while (exp > 0) {
 
      // If the current bit is 1, multiply ret by our value, and mod it.
      if (exp & 1) ret = (ret*value) % mod;
 
      // Shift exp down.
      exp = exp >> 1;
      value = (value*value) % mod;
    }
    return ret;
  }
 
 
  inline std::pair<uint64_t, uint64_t> generate() {
    auto q = static_cast<uint32_t>(std::rand());
    prime::next<uint32_t>(q);
 
    // Sometime's q will not be prime, due to casting.
    // We can just re-roll the number if that's the case.
    auto p = (static_cast<uint64_t>(q) * 2) + 1;
    if (!prime::is(p)) {
      return prime::generate();
    }
    return {p, q};
  }
}