Line data Source code
1 : /*************************************************************************
2 : *
3 : * Project
4 : * _____ _____ __ __ _____
5 : * / ____| __ \| \/ | __ \
6 : * ___ _ __ ___ _ __ | | __| |__) | \ / | |__) |
7 : * / _ \| '_ \ / _ \ '_ \| | |_ | ___/| |\/| | ___/
8 : *| (_) | |_) | __/ | | | |__| | | | | | | |
9 : * \___/| .__/ \___|_| |_|\_____|_| |_| |_|_|
10 : * | |
11 : * |_|
12 : *
13 : * Copyright (C) Akiel Aries, <akiel@akiel.org>, et al.
14 : *
15 : * This software is licensed as described in the file LICENSE, which
16 : * you should have received as part of this distribution. The terms
17 : * among other details are referenced in the official documentation
18 : * seen here : https://akielaries.github.io/openGPMP/ along with
19 : * important files seen in this project.
20 : *
21 : * You may opt to use, copy, modify, merge, publish, distribute
22 : * and/or sell copies of the Software, and permit persons to whom
23 : * the Software is furnished to do so, under the terms of the
24 : * LICENSE file.
25 : *
26 : *
27 : *
28 : * This software is distributed on an AS IS basis, WITHOUT
29 : * WARRANTY OF ANY KIND, either express or implied.
30 : *
31 : ************************************************************************/
32 :
33 : /*
34 : * Some core concepts seen in Number Theory primarily prime number
35 : * based equations, formulas, algorithms.
36 : */
37 : #include <openGPMP/arithmetic.hpp>
38 : #include <openGPMP/nt/prime_test.hpp>
39 :
40 : #include <algorithm>
41 : #include <cmath>
42 : #include <cstdlib>
43 :
44 : #ifdef __APPLE__
45 : #include <stdlib.h>
46 : #endif
47 :
48 : #include <cstring>
49 : #include <iostream>
50 : #include <random>
51 : #include <stdio.h>
52 : #include <string>
53 :
54 : // declare Basics and Primality class objects
55 : gpmp::Basics ba;
56 : gpmp::PrimalityTest prim;
57 :
58 : // Calculates (a * b) mod m
59 4533 : uint64_t gpmp::PrimalityTest::mod_mul(uint64_t a, uint64_t b, uint64_t m) {
60 4533 : uint64_t res = 0;
61 37265 : while (b > 0) {
62 32732 : if (b & 1) {
63 18625 : res = (res + a) % m;
64 : }
65 32732 : a = (2 * a) % m;
66 32732 : b >>= 1;
67 : }
68 4533 : return res;
69 : }
70 :
71 : // Calculates (a ^ b) mod m
72 394 : uint64_t gpmp::PrimalityTest::mod_pow(uint64_t a, uint64_t b, uint64_t m) {
73 394 : uint64_t res = 1;
74 394 : a %= m;
75 3808 : while (b > 0) {
76 3414 : if (b & 1) {
77 1112 : res = mod_mul(res, a, m);
78 : }
79 3414 : a = mod_mul(a, a, m);
80 3414 : b >>= 1;
81 : }
82 394 : return res;
83 : }
84 :
85 : /*
86 : uint64_t gpmp::PrimalityTest::mod_pow(uint64_t base, uint64_t exponent,
87 : uint64_t mod) {
88 : uint64_t x = 1;
89 : uint64_t y = base;
90 :
91 : while (exponent > 0) {
92 : if (exponent & 1) {
93 : x = (x * y) % mod;
94 : }
95 : exponent = exponent >> 1;
96 : y = (y * y) % mod;
97 : }
98 :
99 : return x;
100 : }*/
101 :
102 8 : bool gpmp::PrimalityTest::is_prime(uint64_t n) {
103 8 : if (n <= 1) {
104 1 : return false;
105 : }
106 :
107 : // Check from 2 to n-1
108 532483 : for (uint64_t iter = 2; iter < n; iter++) {
109 532479 : if (n % iter == 0) {
110 3 : return false;
111 : }
112 : }
113 4 : return true;
114 : }
115 :
116 : /*
117 : * determining if a given number is likely to be prime
118 : */
119 0 : bool gpmp::PrimalityTest::compute_miller_rabin(uint64_t d, uint64_t n) {
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 0 : 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 0 : uint64_t x = mod_pow(a, d, n);
140 :
141 0 : if (x == 1 || x == n - 1) {
142 0 : 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 0 : while (d != n - 1) {
151 0 : x = (x * x) % n;
152 0 : d *= 2;
153 :
154 0 : if (x == 1) {
155 0 : return false;
156 : }
157 0 : if (x == n - 1) {
158 0 : return true;
159 : }
160 : }
161 :
162 : // Return composite
163 0 : return false;
164 : }
165 :
166 12 : bool gpmp::PrimalityTest::witness(uint64_t n,
167 : uint64_t d,
168 : uint64_t a,
169 : uint64_t s) {
170 12 : uint64_t x = mod_pow(a, d, n);
171 12 : if (x == 1 || x == n - 1) {
172 6 : return false;
173 : }
174 7 : for (uint64_t r = 1; r < s; r++) {
175 7 : x = mod_mul(x, x, n);
176 7 : if (x == n - 1) {
177 6 : return false;
178 : }
179 : }
180 0 : return true;
181 : }
182 :
183 7 : bool gpmp::PrimalityTest::miller_rabin_prime(uint64_t n, uint64_t iters = 10) {
184 7 : if (n < 2) {
185 0 : return false;
186 : }
187 7 : if (n == 2 || n == 3) {
188 1 : return true;
189 : }
190 6 : if (n % 2 == 0) {
191 3 : return false;
192 : }
193 3 : uint64_t d = n - 1, s = 0;
194 9 : while (d % 2 == 0) {
195 6 : d /= 2;
196 6 : s++;
197 : }
198 15 : for (uint64_t i = 0; i < iters; i++) {
199 12 : uint64_t a = rand() % (n - 3) + 2;
200 12 : if (witness(n, d, a, s)) {
201 0 : return false;
202 : }
203 : }
204 3 : return true;
205 : }
206 :
207 0 : void gpmp::PrimalityTest::miller_rabin(uint64_t iters,
208 : uint64_t min_val,
209 : uint64_t max_val) {
210 0 : std::cout << "Primes between " << min_val << " and " << max_val
211 0 : << 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 0 : for (; min_val < max_val; min_val++) {
218 : // while (min_val < max_val) {
219 0 : if (miller_rabin_prime(min_val, iters)) {
220 0 : std::cout << min_val << " ";
221 : }
222 : }
223 0 : std::cout << "\n";
224 0 : }
225 :
226 0 : bool gpmp::PrimalityTest::AKS(uint64_t n) {
227 : // Check if n is a perfect power
228 0 : for (uint64_t b = 2; b <= std::sqrt(n); b++) {
229 0 : uint64_t a = std::round(std::pow(n, 1.0 / b));
230 : // if (std::pow(a, b) == n) {
231 0 : if (std::pow(a, b) < std::numeric_limits<double>::epsilon()) {
232 0 : return false;
233 : }
234 : }
235 :
236 : // Find the smallest r such that ord_r(n) > log2(n)^2
237 0 : uint64_t r = 2;
238 0 : while (r < n) {
239 0 : bool is_coprime = true;
240 0 : for (uint64_t i = 2; i <= std::sqrt(r); i++) {
241 0 : if (r % i == 0 && n % i == 0) {
242 0 : is_coprime = false;
243 0 : break;
244 : }
245 : }
246 0 : if (is_coprime) {
247 0 : uint64_t t = 1;
248 0 : for (uint64_t i = 0; i < std::floor(2 * std::log2(n)); i++) {
249 0 : t = (t * r) % n;
250 : }
251 0 : if (t == 1) {
252 0 : r++;
253 0 : continue;
254 : }
255 0 : bool is_prime = true;
256 0 : for (uint64_t i = 0; i < std::floor(std::log2(n)); i++) {
257 0 : t = (t * t) % n;
258 0 : if (t == 1) {
259 0 : is_prime = false;
260 0 : break;
261 : }
262 : }
263 0 : if (is_prime) {
264 0 : return true;
265 : }
266 : }
267 0 : r++;
268 : }
269 0 : return false;
270 : }
271 :
272 : /*
273 : * another algorithm capable of finding primes
274 : */
275 58 : uint64_t gpmp::PrimalityTest::jacobian_number(uint64_t a, uint64_t n) {
276 58 : if (!a)
277 0 : return 0; // (0/n) = 0
278 :
279 58 : 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 58 : if (a == 1)
288 25 : return ans; // (1/n) = 1
289 :
290 102 : 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 116 : while (a % 2 == 0) {
298 47 : a = a / 2;
299 47 : if (n % 8 == 3 || n % 8 == 5)
300 33 : ans = -ans;
301 : }
302 69 : std::swap(a, n);
303 :
304 69 : if (a % 4 == 3 && n % 4 == 3)
305 12 : ans = -ans;
306 :
307 69 : a = a % n;
308 :
309 69 : if (a > n / 2)
310 10 : a = a - n;
311 : }
312 :
313 33 : if (n == 1)
314 30 : return ans;
315 :
316 3 : return 0;
317 : }
318 :
319 : /*
320 : * primality test determining if a number is composite or probably
321 : * prime.
322 : *
323 : * a^(p-1)/2 = (a/p) (mod p)
324 : */
325 6 : bool gpmp::PrimalityTest::solovoy_strassen(uint64_t p, uint64_t iters) {
326 6 : if (p < 2)
327 0 : return false;
328 :
329 6 : if (p != 2 && p % 2 == 0)
330 1 : return false;
331 :
332 59 : for (uint64_t i = 0; i < iters; i++) {
333 : // Generate a random number a
334 58 : uint64_t a = rand() % (p - 1) + 1;
335 58 : uint64_t jacobian = (p + jacobian_number(a, p)) % p;
336 58 : uint64_t mod = mod_pow(a, (p - 1) / 2, p);
337 :
338 58 : if (!jacobian || mod != jacobian)
339 4 : return false;
340 : }
341 1 : return true;
342 : }
343 :
344 : /*
345 : * Absolute Fermat Pseudoprimes referred to as Carmichael Numbers
346 : * are composite integers satisfying the congruence forumla below
347 : * b^n - 1 = b (mod n)
348 : */
349 2 : bool gpmp::PrimalityTest::carmichael_num(uint64_t n) {
350 562 : for (uint64_t b = 2; b < n; b++) {
351 : // If "b" is relatively prime to n
352 561 : if (ba.op_gcd(b, n) == 1)
353 : // And pow(b, n-1)%n is not 1, return false.
354 320 : if (mod_pow(b, n - 1, n) != 1)
355 1 : return false;
356 : }
357 1 : return true;
358 : }
359 :
360 : /*
361 : * since this module deals with a great deal of number theory and
362 : * various logical arithmetic operations, the greek algorithm sieve of
363 : * Eratosthenes is able to capture all prime numbers to any given
364 : * limit
365 : */
366 : /*
367 : void gpmp::PrimalityTest::sieve_of_eratosthenes(uint64_t n) {
368 : // Create a boolean array "prime[0..n]" and initialize
369 : // all entries it as true. A value in prime[i] will
370 : // finally be false if i is Not a prime, else true.
371 : bool prime[n + 1];
372 : memset(prime, true, sizeof(prime));
373 :
374 : for (uint64_t p = 2; p * p <= n; p++) {
375 : // If prime[p] is not changed, then it is a prime
376 : if (prime[p] == true) {
377 : // Update all multiples of p greater than or
378 : // equal to the square of it numbers which are
379 : // multiple of p and are less than p^2 are
380 : // already been marked.
381 : for (uint64_t i = p * p; i <= n; i += p)
382 : prime[i] = false;
383 : }
384 : }
385 :
386 : // Print all prime numbers
387 : for (uint64_t p = 2; p <= n; p++) {
388 : if (prime[p]) {
389 : std::cout << p << " " << std::endl;
390 : }
391 : }
392 : }*/
393 :
394 : /*
395 : * Euler's Totient Function counts the positive numbers up until a
396 : * given integer
397 : *
398 : * Definition: For any positive integer n, Φ(n) is the number of all
399 : * positive integers less than or equal to n that are relatively prime
400 : * to n.
401 : *
402 : * Example: Φ(1) = 1
403 : * Φ(2) = 1
404 : * Φ(3) = 2
405 : * Φ(4) = 2
406 : * Φ(5) = 4
407 : */
408 0 : uint64_t gpmp::PrimalityTest::ETF(uint64_t n) {
409 0 : uint64_t result = 1;
410 :
411 0 : for (uint64_t index = 2; uint64_t(index) < n; index++) {
412 0 : if (ba.op_gcd(index, n) == 1) {
413 0 : result++;
414 : }
415 : }
416 :
417 0 : return result;
418 : }
|
