openGPMP
Open Source Mathematics Package
Public Member Functions | List of all members
gpmp::PrimalityTest Class Reference

#include <prime_test.hpp>

Public Member Functions

bool is_prime (uint64_t n)
 Determine if an integer is prime. More...
 
bool compute_miller_rabin (uint64_t d, uint64_t n)
 Algorithm determining liklihood a number is prime. More...
 
bool miller_rabin_prime (uint64_t n, uint64_t iters)
 Modified primes algorithm. More...
 
void miller_rabin (uint64_t iters, uint64_t min_val, uint64_t max_val)
 Miller-Rabin driver, prints values that satisfy conditions. More...
 
bool witness (uint64_t n, uint64_t d, uint64_t a, uint64_t s)
 
bool AKS (uint64_t n)
 
uint64_t jacobian_number (uint64_t a, uint64_t n)
 
bool solovoy_strassen (uint64_t p, uint64_t iters)
 
uint64_t mod_mul (uint64_t a, uint64_t b, uint64_t m)
 
uint64_t mod_pow (uint64_t a, uint64_t b, uint64_t m)
 
bool carmichael_num (uint64_t n)
 
uint64_t ETF (uint64_t n)
 

Detailed Description

Primality Class dealing with prime number manipulations

Examples
primes.cpp.

Definition at line 50 of file prime_test.hpp.

Member Function Documentation

◆ AKS()

bool gpmp::PrimalityTest::AKS ( uint64_t  n)

Definition at line 226 of file prime_test.cpp.

226  {
227  // Check if n is a perfect power
228  for (uint64_t b = 2; b <= std::sqrt(n); b++) {
229  uint64_t a = std::round(std::pow(n, 1.0 / b));
230  // if (std::pow(a, b) == n) {
231  if (std::pow(a, b) < std::numeric_limits<double>::epsilon()) {
232  return false;
233  }
234  }
235 
236  // Find the smallest r such that ord_r(n) > log2(n)^2
237  uint64_t r = 2;
238  while (r < n) {
239  bool is_coprime = true;
240  for (uint64_t i = 2; i <= std::sqrt(r); i++) {
241  if (r % i == 0 && n % i == 0) {
242  is_coprime = false;
243  break;
244  }
245  }
246  if (is_coprime) {
247  uint64_t t = 1;
248  for (uint64_t i = 0; i < std::floor(2 * std::log2(n)); i++) {
249  t = (t * r) % n;
250  }
251  if (t == 1) {
252  r++;
253  continue;
254  }
255  bool is_prime = true;
256  for (uint64_t i = 0; i < std::floor(std::log2(n)); i++) {
257  t = (t * t) % n;
258  if (t == 1) {
259  is_prime = false;
260  break;
261  }
262  }
263  if (is_prime) {
264  return true;
265  }
266  }
267  r++;
268  }
269  return false;
270 }
gpmp::PrimalityTest::is_prime
bool is_prime(uint64_t n)
Determine if an integer is prime.
Definition: prime_test.cpp:102

◆ carmichael_num()

bool gpmp::PrimalityTest::carmichael_num ( uint64_t  n)
Examples
primes.cpp.

Definition at line 349 of file prime_test.cpp.

349  {
350  for (uint64_t b = 2; b < n; b++) {
351  // If "b" is relatively prime to n
352  if (ba.op_gcd(b, n) == 1)
353  // And pow(b, n-1)%n is not 1, return false.
354  if (mod_pow(b, n - 1, n) != 1)
355  return false;
356  }
357  return true;
358 }
gpmp::Basics::op_gcd
int64_t op_gcd(int64_t x, int64_t y)
Find Greatest Common Divisor of 2 integers.
Definition: arith.cpp:50
gpmp::PrimalityTest::mod_pow
uint64_t mod_pow(uint64_t a, uint64_t b, uint64_t m)
Definition: prime_test.cpp:72
ba
gpmp::Basics ba
Definition: prime_test.cpp:55

References ba, and gpmp::Basics::op_gcd().

Referenced by main().

◆ compute_miller_rabin()

bool gpmp::PrimalityTest::compute_miller_rabin ( uint64_t  d,
uint64_t  n 
)

Algorithm determining liklihood a number is prime.

Parameters
[in]d: target number (uint64_t)
[in]n: target - 1 (uint64_t)
Precondition
miller_rabin_prime()
Returns
true/false (bool)

Definition at line 119 of file prime_test.cpp.

119  {
120  // Pick a random number in [2..n-2] Corner cases make sure that n
121  // > 4
122 
123 #ifdef __APPLE__
124  std::random_device rd;
125  std::mt19937 gen(rd());
126  std::uniform_int_distribution<uint64_t> dist(2, n - 3);
127  uint64_t a = dist(gen);
128 #else
129  uint64_t a = 2 + rand() % (n - 4);
130 #endif
131  /*
132  std::random_device rd;
133  std::mt19937_64 gen(rd());
134  std::uniform_int_distribution<uint64_t> dis(2, n - 4);
135  uint64_t a = dis(gen);
136  */
137 
138  // Compute a^d % n
139  uint64_t x = mod_pow(a, d, n);
140 
141  if (x == 1 || x == n - 1) {
142  return true;
143  }
144 
145  // Keep squaring x while one of the following doesn't
146  // happen
147  // (I) d does not reach n-1
148  // (II) (x^2) % n is not 1
149  // (III) (x^2) % n is not n-1
150  while (d != n - 1) {
151  x = (x * x) % n;
152  d *= 2;
153 
154  if (x == 1) {
155  return false;
156  }
157  if (x == n - 1) {
158  return true;
159  }
160  }
161 
162  // Return composite
163  return false;
164 }

◆ ETF()

uint64_t gpmp::PrimalityTest::ETF ( uint64_t  n)

Definition at line 408 of file prime_test.cpp.

408  {
409  uint64_t result = 1;
410 
411  for (uint64_t index = 2; uint64_t(index) < n; index++) {
412  if (ba.op_gcd(index, n) == 1) {
413  result++;
414  }
415  }
416 
417  return result;
418 }

References ba, and gpmp::Basics::op_gcd().

◆ is_prime()

bool gpmp::PrimalityTest::is_prime ( uint64_t  n)

Determine if an integer is prime.

Parameters
[in]n: any number (uint64_t)
Returns
: true/false (bool)
Examples
primes.cpp.

Definition at line 102 of file prime_test.cpp.

102  {
103  if (n <= 1) {
104  return false;
105  }
106 
107  // Check from 2 to n-1
108  for (uint64_t iter = 2; iter < n; iter++) {
109  if (n % iter == 0) {
110  return false;
111  }
112  }
113  return true;
114 }

Referenced by main().

◆ jacobian_number()

uint64_t gpmp::PrimalityTest::jacobian_number ( uint64_t  a,
uint64_t  n 
)

Definition at line 275 of file prime_test.cpp.

275  {
276  if (!a)
277  return 0; // (0/n) = 0
278 
279  uint64_t ans = 1;
280 
281  /*if (a < 0) {
282  a = -a; // (a/n) = (-a/n)*(-1/n)
283  if (n % 4 == 3)
284  ans = -ans; // (-1/n) = -1 if n = 3 (mod 4)
285  }*/
286 
287  if (a == 1)
288  return ans; // (1/n) = 1
289 
290  while (a) {
291  /*if (a < 0) {
292  a = -a; // (a/n) = (-a/n)*(-1/n)
293  if (n % 4 == 3)
294  ans = -ans; // (-1/n) = -1 if n = 3 (mod 4)
295  }*/
296 
297  while (a % 2 == 0) {
298  a = a / 2;
299  if (n % 8 == 3 || n % 8 == 5)
300  ans = -ans;
301  }
302  std::swap(a, n);
303 
304  if (a % 4 == 3 && n % 4 == 3)
305  ans = -ans;
306 
307  a = a % n;
308 
309  if (a > n / 2)
310  a = a - n;
311  }
312 
313  if (n == 1)
314  return ans;
315 
316  return 0;
317 }

◆ miller_rabin()

void gpmp::PrimalityTest::miller_rabin ( uint64_t  iters,
uint64_t  min_val,
uint64_t  max_val 
)

Miller-Rabin driver, prints values that satisfy conditions.

Note
Finds the primes in a given range
Parameters
[in]iters: iterations determine accuracy (uint64_t)
[in]min_val: bottom end of range (uint64_t)
[in]max_val: top end of range (uint64_t)
[out]result: values within range that satisfy
Returns
Void
Examples
primes.cpp.

Definition at line 207 of file prime_test.cpp.

209  {
210  std::cout << "Primes between " << min_val << " and " << max_val
211  << std::endl;
212  // int d = iters - 1;
213  // int min_val = minimum;
214  // int max_val = maximum;
215 
216  // traverse from the min to max
217  for (; min_val < max_val; min_val++) {
218  // while (min_val < max_val) {
219  if (miller_rabin_prime(min_val, iters)) {
220  std::cout << min_val << " ";
221  }
222  }
223  std::cout << "\n";
224 }
gpmp::PrimalityTest::miller_rabin_prime
bool miller_rabin_prime(uint64_t n, uint64_t iters)
Modified primes algorithm.
Definition: prime_test.cpp:183

Referenced by main().

◆ miller_rabin_prime()

bool gpmp::PrimalityTest::miller_rabin_prime ( uint64_t  n,
uint64_t  iters = 10 
)

Modified primes algorithm.

Parameters
[in]n: target number (uint64_t)
[in]iters: iterations determine accuracy (uint64_t)
Precondition
miller_rabin()

return true/false (bool)

Examples
primes.cpp.

Definition at line 183 of file prime_test.cpp.

183  {
184  if (n < 2) {
185  return false;
186  }
187  if (n == 2 || n == 3) {
188  return true;
189  }
190  if (n % 2 == 0) {
191  return false;
192  }
193  uint64_t d = n - 1, s = 0;
194  while (d % 2 == 0) {
195  d /= 2;
196  s++;
197  }
198  for (uint64_t i = 0; i < iters; i++) {
199  uint64_t a = rand() % (n - 3) + 2;
200  if (witness(n, d, a, s)) {
201  return false;
202  }
203  }
204  return true;
205 }
gpmp::PrimalityTest::witness
bool witness(uint64_t n, uint64_t d, uint64_t a, uint64_t s)
Definition: prime_test.cpp:166

Referenced by main(), testing_miller(), and testing_new_miller().

◆ mod_mul()

uint64_t gpmp::PrimalityTest::mod_mul ( uint64_t  a,
uint64_t  b,
uint64_t  m 
)

Definition at line 59 of file prime_test.cpp.

59  {
60  uint64_t res = 0;
61  while (b > 0) {
62  if (b & 1) {
63  res = (res + a) % m;
64  }
65  a = (2 * a) % m;
66  b >>= 1;
67  }
68  return res;
69 }

References python.linalg::res.

◆ mod_pow()

uint64_t gpmp::PrimalityTest::mod_pow ( uint64_t  a,
uint64_t  b,
uint64_t  m 
)
Examples
primes.cpp.

Definition at line 72 of file prime_test.cpp.

72  {
73  uint64_t res = 1;
74  a %= m;
75  while (b > 0) {
76  if (b & 1) {
77  res = mod_mul(res, a, m);
78  }
79  a = mod_mul(a, a, m);
80  b >>= 1;
81  }
82  return res;
83 }
gpmp::PrimalityTest::mod_mul
uint64_t mod_mul(uint64_t a, uint64_t b, uint64_t m)
Definition: prime_test.cpp:59

References python.linalg::res.

Referenced by gpmp::Logarithms::BSGS(), main(), and gpmp::Factorization::pollard_rho().

◆ solovoy_strassen()

bool gpmp::PrimalityTest::solovoy_strassen ( uint64_t  p,
uint64_t  iters 
)
Examples
primes.cpp.

Definition at line 325 of file prime_test.cpp.

325  {
326  if (p < 2)
327  return false;
328 
329  if (p != 2 && p % 2 == 0)
330  return false;
331 
332  for (uint64_t i = 0; i < iters; i++) {
333  // Generate a random number a
334  uint64_t a = rand() % (p - 1) + 1;
335  uint64_t jacobian = (p + jacobian_number(a, p)) % p;
336  uint64_t mod = mod_pow(a, (p - 1) / 2, p);
337 
338  if (!jacobian || mod != jacobian)
339  return false;
340  }
341  return true;
342 }
gpmp::PrimalityTest::jacobian_number
uint64_t jacobian_number(uint64_t a, uint64_t n)
Definition: prime_test.cpp:275

Referenced by main().

◆ witness()

bool gpmp::PrimalityTest::witness ( uint64_t  n,
uint64_t  d,
uint64_t  a,
uint64_t  s 
)

Definition at line 166 of file prime_test.cpp.

169  {
170  uint64_t x = mod_pow(a, d, n);
171  if (x == 1 || x == n - 1) {
172  return false;
173  }
174  for (uint64_t r = 1; r < s; r++) {
175  x = mod_mul(x, x, n);
176  if (x == n - 1) {
177  return false;
178  }
179  }
180  return true;
181 }
The documentation for this class was generated from the following files:
Generated by doxygen 1.9.1