1 /* mpz_millerrabin(n,reps) -- An implementation of the probabilistic primality
2 test found in Knuth's Seminumerical Algorithms book. If the function
3 mpz_millerrabin() returns 0 then n is not prime. If it returns 1, then n is
4 'probably' prime. The probability of a false positive is (1/4)**reps, where
5 reps is the number of internal passes of the probabilistic algorithm. Knuth
6 indicates that 25 passes are reasonable.
7
8 With the current implementation, the first 24 MR-tests are substituted by a
9 Baillie-PSW probable prime test.
10
11 This implementation of the Baillie-PSW test was checked up to 2463*10^12,
12 for smaller values no MR-test is performed, regardless of reps, and
13 2 ("surely prime") is returned if the number was not proved composite.
14
15 If GMP_BPSW_NOFALSEPOSITIVES_UPTO_64BITS is defined as non-zero,
16 the code assumes that the Baillie-PSW test was checked up to 2^64.
17
18 THE FUNCTIONS IN THIS FILE ARE FOR INTERNAL USE ONLY. THEY'RE ALMOST
19 CERTAIN TO BE SUBJECT TO INCOMPATIBLE CHANGES OR DISAPPEAR COMPLETELY IN
20 FUTURE GNU MP RELEASES.
21
22 Copyright 1991, 1993, 1994, 1996-2002, 2005, 2014, 2018-2022 Free
23 Software Foundation, Inc.
24
25 Contributed by John Amanatides.
26 Changed to "BPSW, then Miller Rabin if required" by Marco Bodrato.
27
28 This file is part of the GNU MP Library.
29
30 The GNU MP Library is free software; you can redistribute it and/or modify
31 it under the terms of either:
32
33 * the GNU Lesser General Public License as published by the Free
34 Software Foundation; either version 3 of the License, or (at your
35 option) any later version.
36
37 or
38
39 * the GNU General Public License as published by the Free Software
40 Foundation; either version 2 of the License, or (at your option) any
41 later version.
42
43 or both in parallel, as here.
44
45 The GNU MP Library is distributed in the hope that it will be useful, but
46 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
47 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
48 for more details.
49
50 You should have received copies of the GNU General Public License and the
51 GNU Lesser General Public License along with the GNU MP Library. If not,
52 see https://www.gnu.org/licenses/. */
53
54 #include "gmp-impl.h"
55
56 #ifndef GMP_BPSW_NOFALSEPOSITIVES_UPTO_64BITS
57 #define GMP_BPSW_NOFALSEPOSITIVES_UPTO_64BITS 0
58 #endif
59
60 static int millerrabin (mpz_srcptr,
61 mpz_ptr, mpz_ptr,
62 mpz_srcptr, unsigned long int);
63
64 int
65 mpz_millerrabin (mpz_srcptr n, int reps)
66 {
67 mpz_t nm, x, y, q;
68 mp_bitcnt_t k;
69 int is_prime;
70 TMP_DECL;
71 TMP_MARK;
72
73 ASSERT (SIZ (n) > 0);
74 MPZ_TMP_INIT (nm, SIZ (n) + 1);
75 mpz_tdiv_q_2exp (nm, n, 1);
76
77 MPZ_TMP_INIT (x, SIZ (n) + 1);
78 MPZ_TMP_INIT (y, 2 * SIZ (n)); /* mpz_powm_ui needs excessive memory!!! */
79 MPZ_TMP_INIT (q, SIZ (n));
80
81 /* Find q and k, where q is odd and n = 1 + 2**k * q. */
82 k = mpn_scan1 (PTR (nm), 0);
83 mpz_tdiv_q_2exp (q, nm, k);
84 ++k;
85
86 /* BPSW test */
87 mpz_set_ui (x, 2);
88 is_prime = millerrabin (n, x, y, q, k) && mpz_stronglucas (n, x, y);
89
90 if (is_prime)
91 {
92 if (
93 #if GMP_BPSW_NOFALSEPOSITIVES_UPTO_64BITS
94 /* Consider numbers up to 2^64 that pass the BPSW test as primes. */
95 #if GMP_NUMB_BITS <= 64
96 SIZ (n) <= 64 / GMP_NUMB_BITS
97 #else
98 0
99 #endif
100 #if 64 % GMP_NUMB_BITS != 0
101 || SIZ (n) - 64 / GMP_NUMB_BITS == (PTR (n) [64 / GMP_NUMB_BITS] < CNST_LIMB(1) << 64 % GMP_NUMB_BITS)
102 #endif
103 #else
104 /* Consider numbers up to 35*2^46 that pass the BPSW test as primes.
105 This implementation was tested up to 2463*10^12 > 2^51+2^47+2^46 */
106 /* 2^5 < 35 = 0b100011 < 2^6 */
107 #define GMP_BPSW_LIMB_CONST CNST_LIMB(35)
108 #define GMP_BPSW_BITS_CONST (LOG2C(35) - 1)
109 #define GMP_BPSW_BITS_LIMIT (46 + GMP_BPSW_BITS_CONST)
110
111 #define GMP_BPSW_LIMBS_LIMIT (GMP_BPSW_BITS_LIMIT / GMP_NUMB_BITS)
112 #define GMP_BPSW_BITS_MOD (GMP_BPSW_BITS_LIMIT % GMP_NUMB_BITS)
113
114 #if GMP_NUMB_BITS <= GMP_BPSW_BITS_LIMIT
115 SIZ (n) <= GMP_BPSW_LIMBS_LIMIT
116 #else
117 0
118 #endif
119 #if GMP_BPSW_BITS_MOD >= GMP_BPSW_BITS_CONST
120 || SIZ (n) - GMP_BPSW_LIMBS_LIMIT == (PTR (n) [GMP_BPSW_LIMBS_LIMIT] < GMP_BPSW_LIMB_CONST << (GMP_BPSW_BITS_MOD - GMP_BPSW_BITS_CONST))
121 #else
122 #if GMP_BPSW_BITS_MOD != 0
123 || SIZ (n) - GMP_BPSW_LIMBS_LIMIT == (PTR (n) [GMP_BPSW_LIMBS_LIMIT] < GMP_BPSW_LIMB_CONST >> (GMP_BPSW_BITS_CONST - GMP_BPSW_BITS_MOD))
124 #else
125 #if GMP_NUMB_BITS > GMP_BPSW_BITS_CONST
126 || SIZ (nm) - GMP_BPSW_LIMBS_LIMIT + 1 == (PTR (nm) [GMP_BPSW_LIMBS_LIMIT - 1] < GMP_BPSW_LIMB_CONST << (GMP_NUMB_BITS - 1 - GMP_BPSW_BITS_CONST))
127 #endif
128 #endif
129 #endif
130
131 #undef GMP_BPSW_BITS_LIMIT
132 #undef GMP_BPSW_LIMB_CONST
133 #undef GMP_BPSW_BITS_CONST
134 #undef GMP_BPSW_LIMBS_LIMIT
135 #undef GMP_BPSW_BITS_MOD
136
137 #endif
138 )
139 is_prime = 2;
140 else
141 {
142 reps -= 24;
143 if (reps > 0)
144 {
145 gmp_randstate_t rstate;
146 /* (n-5)/2 */
147 mpz_sub_ui (nm, nm, 2L);
148 ASSERT (mpz_cmp_ui (nm, 1L) >= 0);
149
150 gmp_randinit_default (rstate);
151
152 do
153 {
154 /* 3 to (n-1)/2 inclusive, don't want 1, 0 or 2 */
155 mpz_urandomm (x, rstate, nm);
156 mpz_add_ui (x, x, 3L);
157
158 is_prime = millerrabin (n, x, y, q, k);
159 } while (--reps > 0 && is_prime);
160
161 gmp_randclear (rstate);
162 }
163 }
164 }
165 TMP_FREE;
166 return is_prime;
167 }
168
169 static int
170 mod_eq_m1 (mpz_srcptr x, mpz_srcptr m)
171 {
172 mp_size_t ms;
173 mp_srcptr mp, xp;
174
175 ms = SIZ (m);
176 if (SIZ (x) != ms)
177 return 0;
178 ASSERT (ms > 0);
179
180 mp = PTR (m);
181 xp = PTR (x);
182 ASSERT ((mp[0] - 1) == (mp[0] ^ 1)); /* n is odd */
183
184 if ((*xp ^ CNST_LIMB(1) ^ *mp) != CNST_LIMB(0)) /* xp[0] != mp[0] - 1 */
185 return 0;
186 else
187 {
188 int cmp;
189
190 --ms;
191 ++xp;
192 ++mp;
193
194 MPN_CMP (cmp, xp, mp, ms);
195
196 return cmp == 0;
197 }
198 }
199
200 static int
201 millerrabin (mpz_srcptr n, mpz_ptr x, mpz_ptr y,
202 mpz_srcptr q, mp_bitcnt_t k)
203 {
204 mpz_powm (y, x, q, n);
205
206 if (mpz_cmp_ui (y, 1L) == 0 || mod_eq_m1 (y, n))
207 return 1;
208
209 for (mp_bitcnt_t i = 1; i < k; ++i)
210 {
211 mpz_powm_ui (y, y, 2L, n);
212 if (mod_eq_m1 (y, n))
213 return 1;
214 }
215 return 0;
216 }