Mercurial > dropbear

comparison libtommath/bn_mp_prime_strong_lucas_selfridge.c @ 1739:13d834efc376 fuzz

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
merge from main
author Matt Johnston <matt@ucc.asn.au>
date Thu, 15 Oct 2020 19:55:15 +0800
parents 1051e4eea25a
children
comparison
equal deleted inserted replaced
1562:768ebf737aa0 1739:13d834efc376
1 #include "tommath_private.h"
2 #ifdef BN_MP_PRIME_STRONG_LUCAS_SELFRIDGE_C
3
4 /* LibTomMath, multiple-precision integer library -- Tom St Denis */
5 /* SPDX-License-Identifier: Unlicense */
6
7 /*
8 * See file bn_mp_prime_is_prime.c or the documentation in doc/bn.tex for the details
9 */
10 #ifndef LTM_USE_ONLY_MR
11
12 /*
13 * 8-bit is just too small. You can try the Frobenius test
14 * but that frobenius test can fail, too, for the same reason.
15 */
16 #ifndef MP_8BIT
17
18 /*
19 * multiply bigint a with int d and put the result in c
20 * Like mp_mul_d() but with a signed long as the small input
21 */
22 static mp_err s_mp_mul_si(const mp_int *a, int32_t d, mp_int *c)
23 {
24 mp_int t;
25 mp_err err;
26
27 if ((err = mp_init(&t)) != MP_OKAY) {
28 return err;
29 }
30
31 /*
32 * mp_digit might be smaller than a long, which excludes
33 * the use of mp_mul_d() here.
34 */
35 mp_set_i32(&t, d);
36 err = mp_mul(a, &t, c);
37 mp_clear(&t);
38 return err;
39 }
40 /*
41 Strong Lucas-Selfridge test.
42 returns MP_YES if it is a strong L-S prime, MP_NO if it is composite
43
44 Code ported from Thomas Ray Nicely's implementation of the BPSW test
45 at http://www.trnicely.net/misc/bpsw.html
46
47 Freeware copyright (C) 2016 Thomas R. Nicely <http://www.trnicely.net>.
48 Released into the public domain by the author, who disclaims any legal
49 liability arising from its use
50
51 The multi-line comments are made by Thomas R. Nicely and are copied verbatim.
52 Additional comments marked "CZ" (without the quotes) are by the code-portist.
53
54 (If that name sounds familiar, he is the guy who found the fdiv bug in the
55 Pentium (P5x, I think) Intel processor)
56 */
57 mp_err mp_prime_strong_lucas_selfridge(const mp_int *a, mp_bool *result)
58 {
59 /* CZ TODO: choose better variable names! */
60 mp_int Dz, gcd, Np1, Uz, Vz, U2mz, V2mz, Qmz, Q2mz, Qkdz, T1z, T2z, T3z, T4z, Q2kdz;
61 /* CZ TODO: Some of them need the full 32 bit, hence the (temporary) exclusion of MP_8BIT */
62 int32_t D, Ds, J, sign, P, Q, r, s, u, Nbits;
63 mp_err err;
64 mp_bool oddness;
65
66 *result = MP_NO;
67 /*
68 Find the first element D in the sequence {5, -7, 9, -11, 13, ...}
69 such that Jacobi(D,N) = -1 (Selfridge's algorithm). Theory
70 indicates that, if N is not a perfect square, D will "nearly
71 always" be "small." Just in case, an overflow trap for D is
72 included.
73 */
74
75 if ((err = mp_init_multi(&Dz, &gcd, &Np1, &Uz, &Vz, &U2mz, &V2mz, &Qmz, &Q2mz, &Qkdz, &T1z, &T2z, &T3z, &T4z, &Q2kdz,
76 NULL)) != MP_OKAY) {
77 return err;
78 }
79
80 D = 5;
81 sign = 1;
82
83 for (;;) {
84 Ds = sign * D;
85 sign = -sign;
86 mp_set_u32(&Dz, (uint32_t)D);
87 if ((err = mp_gcd(a, &Dz, &gcd)) != MP_OKAY) goto LBL_LS_ERR;
88
89 /* if 1 < GCD < N then N is composite with factor "D", and
90 Jacobi(D,N) is technically undefined (but often returned
91 as zero). */
92 if ((mp_cmp_d(&gcd, 1uL) == MP_GT) && (mp_cmp(&gcd, a) == MP_LT)) {
93 goto LBL_LS_ERR;
94 }
95 if (Ds < 0) {
96 Dz.sign = MP_NEG;
97 }
98 if ((err = mp_kronecker(&Dz, a, &J)) != MP_OKAY) goto LBL_LS_ERR;
99
100 if (J == -1) {
101 break;
102 }
103 D += 2;
104
105 if (D > (INT_MAX - 2)) {
106 err = MP_VAL;
107 goto LBL_LS_ERR;
108 }
109 }
110
111
112
113 P = 1; /* Selfridge's choice */
114 Q = (1 - Ds) / 4; /* Required so D = P*P - 4*Q */
115
116 /* NOTE: The conditions (a) N does not divide Q, and
117 (b) D is square-free or not a perfect square, are included by
118 some authors; e.g., "Prime numbers and computer methods for
119 factorization," Hans Riesel (2nd ed., 1994, Birkhauser, Boston),
120 p. 130. For this particular application of Lucas sequences,
121 these conditions were found to be immaterial. */
122
123 /* Now calculate N - Jacobi(D,N) = N + 1 (even), and calculate the
124 odd positive integer d and positive integer s for which
125 N + 1 = 2^s*d (similar to the step for N - 1 in Miller's test).
126 The strong Lucas-Selfridge test then returns N as a strong
127 Lucas probable prime (slprp) if any of the following
128 conditions is met: U_d=0, V_d=0, V_2d=0, V_4d=0, V_8d=0,
129 V_16d=0, ..., etc., ending with V_{2^(s-1)*d}=V_{(N+1)/2}=0
130 (all equalities mod N). Thus d is the highest index of U that
131 must be computed (since V_2m is independent of U), compared
132 to U_{N+1} for the standard Lucas-Selfridge test; and no
133 index of V beyond (N+1)/2 is required, just as in the
134 standard Lucas-Selfridge test. However, the quantity Q^d must
135 be computed for use (if necessary) in the latter stages of
136 the test. The result is that the strong Lucas-Selfridge test
137 has a running time only slightly greater (order of 10 %) than
138 that of the standard Lucas-Selfridge test, while producing
139 only (roughly) 30 % as many pseudoprimes (and every strong
140 Lucas pseudoprime is also a standard Lucas pseudoprime). Thus
141 the evidence indicates that the strong Lucas-Selfridge test is
142 more effective than the standard Lucas-Selfridge test, and a
143 Baillie-PSW test based on the strong Lucas-Selfridge test
144 should be more reliable. */
145
146 if ((err = mp_add_d(a, 1uL, &Np1)) != MP_OKAY) goto LBL_LS_ERR;
147 s = mp_cnt_lsb(&Np1);
148
149 /* CZ
150 * This should round towards zero because
151 * Thomas R. Nicely used GMP's mpz_tdiv_q_2exp()
152 * and mp_div_2d() is equivalent. Additionally:
153 * dividing an even number by two does not produce
154 * any leftovers.
155 */
156 if ((err = mp_div_2d(&Np1, s, &Dz, NULL)) != MP_OKAY) goto LBL_LS_ERR;
157 /* We must now compute U_d and V_d. Since d is odd, the accumulated
158 values U and V are initialized to U_1 and V_1 (if the target
159 index were even, U and V would be initialized instead to U_0=0
160 and V_0=2). The values of U_2m and V_2m are also initialized to
161 U_1 and V_1; the FOR loop calculates in succession U_2 and V_2,
162 U_4 and V_4, U_8 and V_8, etc. If the corresponding bits
163 (1, 2, 3, ...) of t are on (the zero bit having been accounted
164 for in the initialization of U and V), these values are then
165 combined with the previous totals for U and V, using the
166 composition formulas for addition of indices. */
167
168 mp_set(&Uz, 1uL); /* U=U_1 */
169 mp_set(&Vz, (mp_digit)P); /* V=V_1 */
170 mp_set(&U2mz, 1uL); /* U_1 */
171 mp_set(&V2mz, (mp_digit)P); /* V_1 */
172
173 mp_set_i32(&Qmz, Q);
174 if ((err = mp_mul_2(&Qmz, &Q2mz)) != MP_OKAY) goto LBL_LS_ERR;
175 /* Initializes calculation of Q^d */
176 mp_set_i32(&Qkdz, Q);
177
178 Nbits = mp_count_bits(&Dz);
179
180 for (u = 1; u < Nbits; u++) { /* zero bit off, already accounted for */
181 /* Formulas for doubling of indices (carried out mod N). Note that
182 * the indices denoted as "2m" are actually powers of 2, specifically
183 * 2^(ul-1) beginning each loop and 2^ul ending each loop.
184 *
185 * U_2m = U_m*V_m
186 * V_2m = V_m*V_m - 2*Q^m
187 */
188
189 if ((err = mp_mul(&U2mz, &V2mz, &U2mz)) != MP_OKAY) goto LBL_LS_ERR;
190 if ((err = mp_mod(&U2mz, a, &U2mz)) != MP_OKAY) goto LBL_LS_ERR;
191 if ((err = mp_sqr(&V2mz, &V2mz)) != MP_OKAY) goto LBL_LS_ERR;
192 if ((err = mp_sub(&V2mz, &Q2mz, &V2mz)) != MP_OKAY) goto LBL_LS_ERR;
193 if ((err = mp_mod(&V2mz, a, &V2mz)) != MP_OKAY) goto LBL_LS_ERR;
194
195 /* Must calculate powers of Q for use in V_2m, also for Q^d later */
196 if ((err = mp_sqr(&Qmz, &Qmz)) != MP_OKAY) goto LBL_LS_ERR;
197
198 /* prevents overflow */ /* CZ still necessary without a fixed prealloc'd mem.? */
199 if ((err = mp_mod(&Qmz, a, &Qmz)) != MP_OKAY) goto LBL_LS_ERR;
200 if ((err = mp_mul_2(&Qmz, &Q2mz)) != MP_OKAY) goto LBL_LS_ERR;
201
202 if (s_mp_get_bit(&Dz, (unsigned int)u) == MP_YES) {
203 /* Formulas for addition of indices (carried out mod N);
204 *
205 * U_(m+n) = (U_m*V_n + U_n*V_m)/2
206 * V_(m+n) = (V_m*V_n + D*U_m*U_n)/2
207 *
208 * Be careful with division by 2 (mod N)!
209 */
210 if ((err = mp_mul(&U2mz, &Vz, &T1z)) != MP_OKAY) goto LBL_LS_ERR;
211 if ((err = mp_mul(&Uz, &V2mz, &T2z)) != MP_OKAY) goto LBL_LS_ERR;
212 if ((err = mp_mul(&V2mz, &Vz, &T3z)) != MP_OKAY) goto LBL_LS_ERR;
213 if ((err = mp_mul(&U2mz, &Uz, &T4z)) != MP_OKAY) goto LBL_LS_ERR;
214 if ((err = s_mp_mul_si(&T4z, Ds, &T4z)) != MP_OKAY) goto LBL_LS_ERR;
215 if ((err = mp_add(&T1z, &T2z, &Uz)) != MP_OKAY) goto LBL_LS_ERR;
216 if (MP_IS_ODD(&Uz)) {
217 if ((err = mp_add(&Uz, a, &Uz)) != MP_OKAY) goto LBL_LS_ERR;
218 }
219 /* CZ
220 * This should round towards negative infinity because
221 * Thomas R. Nicely used GMP's mpz_fdiv_q_2exp().
222 * But mp_div_2() does not do so, it is truncating instead.
223 */
224 oddness = MP_IS_ODD(&Uz) ? MP_YES : MP_NO;
225 if ((err = mp_div_2(&Uz, &Uz)) != MP_OKAY) goto LBL_LS_ERR;
226 if ((Uz.sign == MP_NEG) && (oddness != MP_NO)) {
227 if ((err = mp_sub_d(&Uz, 1uL, &Uz)) != MP_OKAY) goto LBL_LS_ERR;
228 }
229 if ((err = mp_add(&T3z, &T4z, &Vz)) != MP_OKAY) goto LBL_LS_ERR;
230 if (MP_IS_ODD(&Vz)) {
231 if ((err = mp_add(&Vz, a, &Vz)) != MP_OKAY) goto LBL_LS_ERR;
232 }
233 oddness = MP_IS_ODD(&Vz) ? MP_YES : MP_NO;
234 if ((err = mp_div_2(&Vz, &Vz)) != MP_OKAY) goto LBL_LS_ERR;
235 if ((Vz.sign == MP_NEG) && (oddness != MP_NO)) {
236 if ((err = mp_sub_d(&Vz, 1uL, &Vz)) != MP_OKAY) goto LBL_LS_ERR;
237 }
238 if ((err = mp_mod(&Uz, a, &Uz)) != MP_OKAY) goto LBL_LS_ERR;
239 if ((err = mp_mod(&Vz, a, &Vz)) != MP_OKAY) goto LBL_LS_ERR;
240
241 /* Calculating Q^d for later use */
242 if ((err = mp_mul(&Qkdz, &Qmz, &Qkdz)) != MP_OKAY) goto LBL_LS_ERR;
243 if ((err = mp_mod(&Qkdz, a, &Qkdz)) != MP_OKAY) goto LBL_LS_ERR;
244 }
245 }
246
247 /* If U_d or V_d is congruent to 0 mod N, then N is a prime or a
248 strong Lucas pseudoprime. */
249 if (MP_IS_ZERO(&Uz) || MP_IS_ZERO(&Vz)) {
250 *result = MP_YES;
251 goto LBL_LS_ERR;
252 }
253
254 /* NOTE: Ribenboim ("The new book of prime number records," 3rd ed.,
255 1995/6) omits the condition V0 on p.142, but includes it on
256 p. 130. The condition is NECESSARY; otherwise the test will
257 return false negatives---e.g., the primes 29 and 2000029 will be
258 returned as composite. */
259
260 /* Otherwise, we must compute V_2d, V_4d, V_8d, ..., V_{2^(s-1)*d}
261 by repeated use of the formula V_2m = V_m*V_m - 2*Q^m. If any of
262 these are congruent to 0 mod N, then N is a prime or a strong
263 Lucas pseudoprime. */
264
265 /* Initialize 2*Q^(d*2^r) for V_2m */
266 if ((err = mp_mul_2(&Qkdz, &Q2kdz)) != MP_OKAY) goto LBL_LS_ERR;
267
268 for (r = 1; r < s; r++) {
269 if ((err = mp_sqr(&Vz, &Vz)) != MP_OKAY) goto LBL_LS_ERR;
270 if ((err = mp_sub(&Vz, &Q2kdz, &Vz)) != MP_OKAY) goto LBL_LS_ERR;
271 if ((err = mp_mod(&Vz, a, &Vz)) != MP_OKAY) goto LBL_LS_ERR;
272 if (MP_IS_ZERO(&Vz)) {
273 *result = MP_YES;
274 goto LBL_LS_ERR;
275 }
276 /* Calculate Q^{d*2^r} for next r (final iteration irrelevant). */
277 if (r < (s - 1)) {
278 if ((err = mp_sqr(&Qkdz, &Qkdz)) != MP_OKAY) goto LBL_LS_ERR;
279 if ((err = mp_mod(&Qkdz, a, &Qkdz)) != MP_OKAY) goto LBL_LS_ERR;
280 if ((err = mp_mul_2(&Qkdz, &Q2kdz)) != MP_OKAY) goto LBL_LS_ERR;
281 }
282 }
283 LBL_LS_ERR:
284 mp_clear_multi(&Q2kdz, &T4z, &T3z, &T2z, &T1z, &Qkdz, &Q2mz, &Qmz, &V2mz, &U2mz, &Vz, &Uz, &Np1, &gcd, &Dz, NULL);
285 return err;
286 }
287 #endif
288 #endif
289 #endif