/src/gmp/mpz/millerrabin.c
Line | Count | Source (jump to first uncovered line) |
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 | 0 | { |
67 | 0 | mpz_t nm, x, y, q; |
68 | 0 | mp_bitcnt_t k; |
69 | 0 | int is_prime; |
70 | 0 | TMP_DECL; |
71 | 0 | TMP_MARK; |
72 | |
|
73 | 0 | ASSERT (SIZ (n) > 0); |
74 | 0 | MPZ_TMP_INIT (nm, SIZ (n) + 1); |
75 | 0 | mpz_tdiv_q_2exp (nm, n, 1); |
76 | |
|
77 | 0 | MPZ_TMP_INIT (x, SIZ (n) + 1); |
78 | 0 | MPZ_TMP_INIT (y, 2 * SIZ (n)); /* mpz_powm_ui needs excessive memory!!! */ |
79 | 0 | MPZ_TMP_INIT (q, SIZ (n)); |
80 | | |
81 | | /* Find q and k, where q is odd and n = 1 + 2**k * q. */ |
82 | 0 | k = mpn_scan1 (PTR (nm), 0); |
83 | 0 | mpz_tdiv_q_2exp (q, nm, k); |
84 | 0 | ++k; |
85 | | |
86 | | /* BPSW test */ |
87 | 0 | mpz_set_ui (x, 2); |
88 | 0 | is_prime = millerrabin (n, x, y, q, k) && mpz_stronglucas (n, x, y); |
89 | |
|
90 | 0 | if (is_prime) |
91 | 0 | { |
92 | 0 | 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 | 0 | #define GMP_BPSW_LIMB_CONST CNST_LIMB(35) |
108 | 0 | #define GMP_BPSW_BITS_CONST (LOG2C(35) - 1) |
109 | 0 | #define GMP_BPSW_BITS_LIMIT (46 + GMP_BPSW_BITS_CONST) |
110 | |
|
111 | 0 | #define GMP_BPSW_LIMBS_LIMIT (GMP_BPSW_BITS_LIMIT / GMP_NUMB_BITS) |
112 | 0 | #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 | 0 |
118 | 0 | #endif |
119 | 0 | #if GMP_BPSW_BITS_MOD >= GMP_BPSW_BITS_CONST |
120 | 0 | || 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 | 0 | #undef GMP_BPSW_BITS_LIMIT |
132 | 0 | #undef GMP_BPSW_LIMB_CONST |
133 | 0 | #undef GMP_BPSW_BITS_CONST |
134 | 0 | #undef GMP_BPSW_LIMBS_LIMIT |
135 | 0 | #undef GMP_BPSW_BITS_MOD |
136 | |
|
137 | 0 | #endif |
138 | 0 | ) |
139 | 0 | is_prime = 2; |
140 | 0 | else |
141 | 0 | { |
142 | 0 | reps -= 24; |
143 | 0 | if (reps > 0) |
144 | 0 | { |
145 | 0 | gmp_randstate_t rstate; |
146 | | /* (n-5)/2 */ |
147 | 0 | mpz_sub_ui (nm, nm, 2L); |
148 | 0 | ASSERT (mpz_cmp_ui (nm, 1L) >= 0); |
149 | |
|
150 | 0 | gmp_randinit_default (rstate); |
151 | |
|
152 | 0 | do |
153 | 0 | { |
154 | | /* 3 to (n-1)/2 inclusive, don't want 1, 0 or 2 */ |
155 | 0 | mpz_urandomm (x, rstate, nm); |
156 | 0 | mpz_add_ui (x, x, 3L); |
157 | |
|
158 | 0 | is_prime = millerrabin (n, x, y, q, k); |
159 | 0 | } while (--reps > 0 && is_prime); |
160 | |
|
161 | 0 | gmp_randclear (rstate); |
162 | 0 | } |
163 | 0 | } |
164 | 0 | } |
165 | 0 | TMP_FREE; |
166 | 0 | return is_prime; |
167 | 0 | } |
168 | | |
169 | | static int |
170 | | mod_eq_m1 (mpz_srcptr x, mpz_srcptr m) |
171 | 0 | { |
172 | 0 | mp_size_t ms; |
173 | 0 | mp_srcptr mp, xp; |
174 | |
|
175 | 0 | ms = SIZ (m); |
176 | 0 | if (SIZ (x) != ms) |
177 | 0 | return 0; |
178 | 0 | ASSERT (ms > 0); |
179 | |
|
180 | 0 | mp = PTR (m); |
181 | 0 | xp = PTR (x); |
182 | 0 | ASSERT ((mp[0] - 1) == (mp[0] ^ 1)); /* n is odd */ |
183 | |
|
184 | 0 | if ((*xp ^ CNST_LIMB(1) ^ *mp) != CNST_LIMB(0)) /* xp[0] != mp[0] - 1 */ |
185 | 0 | return 0; |
186 | 0 | else |
187 | 0 | { |
188 | 0 | int cmp; |
189 | |
|
190 | 0 | --ms; |
191 | 0 | ++xp; |
192 | 0 | ++mp; |
193 | |
|
194 | 0 | MPN_CMP (cmp, xp, mp, ms); |
195 | |
|
196 | 0 | return cmp == 0; |
197 | 0 | } |
198 | 0 | } |
199 | | |
200 | | static int |
201 | | millerrabin (mpz_srcptr n, mpz_ptr x, mpz_ptr y, |
202 | | mpz_srcptr q, mp_bitcnt_t k) |
203 | 0 | { |
204 | 0 | mpz_powm (y, x, q, n); |
205 | |
|
206 | 0 | if (mpz_cmp_ui (y, 1L) == 0 || mod_eq_m1 (y, n)) |
207 | 0 | return 1; |
208 | | |
209 | 0 | for (mp_bitcnt_t i = 1; i < k; ++i) |
210 | 0 | { |
211 | 0 | mpz_powm_ui (y, y, 2L, n); |
212 | 0 | if (mod_eq_m1 (y, n)) |
213 | 0 | return 1; |
214 | 0 | } |
215 | 0 | return 0; |
216 | 0 | } |