/src/libgmp/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 2640*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, 2024 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 |
61 | | mod_eq_m1 (mpz_srcptr x, mpz_srcptr m) |
62 | 46.7k | { |
63 | 46.7k | mp_size_t ms; |
64 | 46.7k | mp_srcptr mp, xp; |
65 | | |
66 | 46.7k | ms = SIZ (m); |
67 | 46.7k | if (SIZ (x) != ms) |
68 | 6.06k | return 0; |
69 | 40.6k | ASSERT (ms > 0); |
70 | | |
71 | 40.6k | mp = PTR (m); |
72 | 40.6k | xp = PTR (x); |
73 | 40.6k | ASSERT ((mp[0] - 1) == (mp[0] ^ 1)); /* n is odd */ |
74 | | |
75 | 40.6k | if ((*xp ^ CNST_LIMB(1) ^ *mp) != CNST_LIMB(0)) /* xp[0] != mp[0] - 1 */ |
76 | 37.9k | return 0; |
77 | 2.67k | else |
78 | 2.67k | { |
79 | 2.67k | int cmp; |
80 | | |
81 | 2.67k | --ms; |
82 | 2.67k | ++xp; |
83 | 2.67k | ++mp; |
84 | | |
85 | 2.67k | MPN_CMP (cmp, xp, mp, ms); |
86 | | |
87 | 2.67k | return cmp == 0; |
88 | 2.67k | } |
89 | 40.6k | } |
90 | | |
91 | | static int |
92 | | millerrabin (mpz_srcptr n, mpz_ptr x, mpz_ptr y, |
93 | | mpz_srcptr q, mp_bitcnt_t k) |
94 | 24.4k | { |
95 | 24.4k | mpz_powm (y, x, q, n); |
96 | | |
97 | 24.4k | if (mpz_cmp_ui (y, 1L) == 0 || mod_eq_m1 (y, n)) |
98 | 2.41k | return 1; |
99 | | |
100 | 44.0k | for (mp_bitcnt_t i = 1; i < k; ++i) |
101 | 23.4k | { |
102 | 23.4k | mpz_powm_ui (y, y, 2L, n); |
103 | 23.4k | if (mod_eq_m1 (y, n)) |
104 | 1.50k | return 1; |
105 | 23.4k | } |
106 | 20.5k | return 0; |
107 | 22.0k | } |
108 | | |
109 | | int |
110 | | mpz_millerrabin (mpz_srcptr n, int reps) |
111 | 22.7k | { |
112 | 22.7k | mpz_t nm, x, y, q; |
113 | 22.7k | mp_bitcnt_t k; |
114 | 22.7k | int is_prime; |
115 | 22.7k | TMP_DECL; |
116 | 22.7k | TMP_MARK; |
117 | | |
118 | 22.7k | ASSERT (SIZ (n) > 0); |
119 | 22.7k | MPZ_TMP_INIT (nm, SIZ (n) + 1); |
120 | 22.7k | mpz_tdiv_q_2exp (nm, n, 1); |
121 | | |
122 | 22.7k | MPZ_TMP_INIT (x, SIZ (n) + 1); |
123 | 22.7k | MPZ_TMP_INIT (y, 2 * SIZ (n)); /* mpz_powm_ui needs excessive memory!!! */ |
124 | 22.7k | MPZ_TMP_INIT (q, SIZ (n)); |
125 | | |
126 | | /* Find q and k, where q is odd and n = 1 + 2**k * q. */ |
127 | 22.7k | k = mpn_scan1 (PTR (nm), 0); |
128 | 22.7k | mpz_tdiv_q_2exp (q, nm, k); |
129 | 22.7k | ++k; |
130 | | |
131 | | /* BPSW test */ |
132 | 22.7k | mpz_set_ui (x, 2); |
133 | 22.7k | is_prime = millerrabin (n, x, y, q, k) && mpz_stronglucas (n, x, y); |
134 | | |
135 | 22.7k | if (is_prime) |
136 | 2.13k | { |
137 | 2.13k | if ( |
138 | | #if GMP_BPSW_NOFALSEPOSITIVES_UPTO_64BITS |
139 | | /* Consider numbers up to 2^64 that pass the BPSW test as primes. */ |
140 | | #if GMP_NUMB_BITS <= 64 |
141 | | SIZ (n) <= 64 / GMP_NUMB_BITS |
142 | | #else |
143 | | 0 |
144 | | #endif |
145 | | #if 64 % GMP_NUMB_BITS != 0 |
146 | | || SIZ (n) - 64 / GMP_NUMB_BITS == (PTR (n) [64 / GMP_NUMB_BITS] < CNST_LIMB(1) << 64 % GMP_NUMB_BITS) |
147 | | #endif |
148 | | #else |
149 | | /* Consider numbers up to 75*2^45 that pass the BPSW test as primes. |
150 | | This implementation was tested up to 264*10^13 > 2^51+2^48+2^46+2^45 */ |
151 | | /* 2^6 < 75 = 0b1001011 < 2^7 */ |
152 | 2.13k | #define GMP_BPSW_LIMB_CONST CNST_LIMB(75) |
153 | 8.54k | #define GMP_BPSW_BITS_CONST (LOG2C(75) - 1) |
154 | 6.40k | #define GMP_BPSW_BITS_LIMIT (45 + GMP_BPSW_BITS_CONST) |
155 | | |
156 | 4.27k | #define GMP_BPSW_LIMBS_LIMIT (GMP_BPSW_BITS_LIMIT / GMP_NUMB_BITS) |
157 | 2.13k | #define GMP_BPSW_BITS_MOD (GMP_BPSW_BITS_LIMIT % GMP_NUMB_BITS) |
158 | | |
159 | | #if GMP_NUMB_BITS <= GMP_BPSW_BITS_LIMIT |
160 | | SIZ (n) <= GMP_BPSW_LIMBS_LIMIT |
161 | | #else |
162 | 2.13k | 0 |
163 | 2.13k | #endif |
164 | 2.13k | #if GMP_BPSW_BITS_MOD >= GMP_BPSW_BITS_CONST |
165 | 2.13k | || SIZ (n) - GMP_BPSW_LIMBS_LIMIT == (PTR (n) [GMP_BPSW_LIMBS_LIMIT] < GMP_BPSW_LIMB_CONST << (GMP_BPSW_BITS_MOD - GMP_BPSW_BITS_CONST)) |
166 | | #else |
167 | | #if GMP_BPSW_BITS_MOD != 0 |
168 | | || SIZ (n) - GMP_BPSW_LIMBS_LIMIT == (PTR (n) [GMP_BPSW_LIMBS_LIMIT] < GMP_BPSW_LIMB_CONST >> (GMP_BPSW_BITS_CONST - GMP_BPSW_BITS_MOD)) |
169 | | #else |
170 | | #if GMP_NUMB_BITS > GMP_BPSW_BITS_CONST |
171 | | || 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)) |
172 | | #endif |
173 | | #endif |
174 | | #endif |
175 | | |
176 | 2.13k | #undef GMP_BPSW_BITS_LIMIT |
177 | 2.13k | #undef GMP_BPSW_LIMB_CONST |
178 | 2.13k | #undef GMP_BPSW_BITS_CONST |
179 | 2.13k | #undef GMP_BPSW_LIMBS_LIMIT |
180 | 2.13k | #undef GMP_BPSW_BITS_MOD |
181 | | |
182 | 2.13k | #endif |
183 | 2.13k | ) |
184 | 259 | is_prime = 2; |
185 | 1.87k | else |
186 | 1.87k | { |
187 | 1.87k | reps -= 24; |
188 | 1.87k | if (reps > 0) |
189 | 1.78k | { |
190 | 1.78k | gmp_randstate_t rstate; |
191 | | /* (n-5)/2 */ |
192 | 1.78k | mpz_sub_ui (nm, nm, 2L); |
193 | 1.78k | ASSERT (mpz_cmp_ui (nm, 1L) >= 0); |
194 | | |
195 | 1.78k | gmp_randinit_default (rstate); |
196 | | |
197 | 1.78k | do |
198 | 1.78k | { |
199 | | /* 3 to (n-1)/2 inclusive, don't want 1, 0 or 2 */ |
200 | 1.78k | mpz_urandomm (x, rstate, nm); |
201 | 1.78k | mpz_add_ui (x, x, 3L); |
202 | | |
203 | 1.78k | is_prime = millerrabin (n, x, y, q, k); |
204 | 1.78k | } while (--reps > 0 && is_prime); |
205 | | |
206 | 1.78k | gmp_randclear (rstate); |
207 | 1.78k | } |
208 | 1.87k | } |
209 | 2.13k | } |
210 | 22.7k | TMP_FREE; |
211 | 22.7k | return is_prime; |
212 | 22.7k | } |