1 /*
2 Copyright (C) 2009 William Hart
3
4 This file is part of FLINT.
5
6 FLINT is free software: you can redistribute it and/or modify it under
7 the terms of the GNU Lesser General Public License (LGPL) as published
8 by the Free Software Foundation; either version 2.1 of the License, or
9 (at your option) any later version. See <http://www.gnu.org/licenses/>.
10 */
11
12 #define ulong ulongxx /* interferes with system includes */
13 #include <stdlib.h>
14 #include <stdio.h>
15 #undef ulong
16 #include <gmp.h>
17 #include "flint.h"
18 #include "ulong_extras.h"
19
20 mp_limb_t flint_pseudosquares[] = {17, 73, 241, 1009, 2641, 8089, 18001,
21 53881, 87481, 117049, 515761, 1083289, 3206641, 3818929, 9257329,
22 22000801, 48473881, 48473881, 175244281, 427733329, 427733329,
23 898716289u, 2805544681u, 2805544681u, 2805544681u
24 #ifndef FLINT64
25 };
26 #else
27 , 10310263441u, 23616331489u, 85157610409u, 85157610409u,
28 196265095009u, 196265095009u, 2871842842801u, 2871842842801u,
29 2871842842801u, 26250887023729u, 26250887023729u, 112434732901969u,
30 112434732901969u, 112434732901969u, 178936222537081u,
31 178936222537081u, 696161110209049u, 696161110209049u,
32 2854909648103881u, 6450045516630769u, 6450045516630769u,
33 11641399247947921u, 11641399247947921u, 190621428905186449u,
34 196640248121928601u, 712624335095093521u, 1773855791877850321u };
35 #endif
36
37 #if FLINT64
38 #define FLINT_NUM_PSEUDOSQUARES 52
39 #else
40 #define FLINT_NUM_PSEUDOSQUARES 25
41 #endif
42
n_is_prime_pseudosquare(mp_limb_t n)43 int n_is_prime_pseudosquare(mp_limb_t n)
44 {
45 unsigned int i, j, m1;
46 mp_limb_t p, B, NB, exp, mod8;
47 const mp_limb_t * primes;
48 const double * inverses;
49
50 if (n < UWORD(2)) return 0;
51
52 if ((n & UWORD(1)) == UWORD(0))
53 {
54 return (n == UWORD(2));
55 }
56
57 primes = n_primes_arr_readonly(FLINT_PSEUDOSQUARES_CUTOFF+1);
58 inverses = n_prime_inverses_arr_readonly(FLINT_PSEUDOSQUARES_CUTOFF+1);
59
60 for (i = 0; i < FLINT_PSEUDOSQUARES_CUTOFF; i++)
61 {
62 double ppre;
63 p = primes[i];
64 if (p*p > n) return 1;
65 ppre = inverses[i];
66 if (!n_mod2_precomp(n, p, ppre)) return 0;
67 }
68
69 B = primes[FLINT_PSEUDOSQUARES_CUTOFF];
70 NB = (n - 1)/B + 1;
71 m1 = 0;
72
73 for (i = 0; i < FLINT_NUM_PSEUDOSQUARES; i++)
74 {
75 if (flint_pseudosquares[i] > NB) break;
76 }
77
78 exp = (n - 1)/2;
79
80 for (j = 0; j <= i; j++)
81 {
82 mp_limb_t mod = n_powmod2(primes[j], exp, n);
83 if ((mod != UWORD(1)) && (mod != n - 1)) return 0;
84 if (mod == n - 1) m1 = 1;
85 }
86
87 mod8 = n % 8;
88
89 if ((mod8 == 3) || (mod8 == 7)) return 1;
90
91 if (mod8 == 5)
92 {
93 mp_limb_t mod = n_powmod2(UWORD(2), exp, n);
94 if (mod == n - 1) return 1;
95 flint_printf("Whoah, %wu is a probable prime, but not prime, please report!!\n", n);
96 flint_abort();
97 }
98 else
99 {
100 if (m1) return 1;
101 for (j = i + 1; j < FLINT_NUM_PSEUDOSQUARES + 1; j++)
102 {
103 mp_limb_t mod = n_powmod2(primes[j], exp, n);
104 if (mod == n - 1) return 1;
105 if (mod != 1)
106 {
107 flint_printf("Whoah, %wu is a probable prime, but not prime, please report!!\n", n);
108 flint_abort();
109 }
110 }
111 flint_printf("Whoah, %wu is a probable prime, but not prime, please report!!\n", n);
112 flint_abort();
113 }
114
115 return 0; /* not reached, but silence the compiler */
116 }
117