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 *
6 * LibTomMath is a library that provides multiple-precision
7 * integer arithmetic as well as number theoretic functionality.
8 *
9 * The library was designed directly after the MPI library by
10 * Michael Fromberger but has been written from scratch with
11 * additional optimizations in place.
12 *
13 * SPDX-License-Identifier: Unlicense
14 */
15
16 /*
17 * See file bn_mp_prime_is_prime.c or the documentation in doc/bn.tex for the details
18 */
19 #ifndef LTM_USE_FIPS_ONLY
20
21 /*
22 * 8-bit is just too small. You can try the Frobenius test
23 * but that frobenius test can fail, too, for the same reason.
24 */
25 #ifndef MP_8BIT
26
27 /*
28 * multiply bigint a with int d and put the result in c
29 * Like mp_mul_d() but with a signed long as the small input
30 */
s_mp_mul_si(const mp_int * a,long d,mp_int * c)31 static int s_mp_mul_si(const mp_int *a, long d, mp_int *c)
32 {
33 mp_int t;
34 int err, neg = 0;
35
36 if ((err = mp_init(&t)) != MP_OKAY) {
37 return err;
38 }
39 if (d < 0) {
40 neg = 1;
41 d = -d;
42 }
43
44 /*
45 * mp_digit might be smaller than a long, which excludes
46 * the use of mp_mul_d() here.
47 */
48 if ((err = mp_set_long(&t, (unsigned long) d)) != MP_OKAY) {
49 goto LBL_MPMULSI_ERR;
50 }
51 if ((err = mp_mul(a, &t, c)) != MP_OKAY) {
52 goto LBL_MPMULSI_ERR;
53 }
54 if (neg == 1) {
55 c->sign = (a->sign == MP_NEG) ? MP_ZPOS: MP_NEG;
56 }
57 LBL_MPMULSI_ERR:
58 mp_clear(&t);
59 return err;
60 }
61
62
63
64 /*
65 Strong Lucas-Selfridge test.
66 returns MP_YES if it is a strong L-S prime, MP_NO if it is composite
67
68 Code ported from Thomas Ray Nicely's implementation of the BPSW test
69 at http://www.trnicely.net/misc/bpsw.html
70
71 Freeware copyright (C) 2016 Thomas R. Nicely <http://www.trnicely.net>.
72 Released into the public domain by the author, who disclaims any legal
73 liability arising from its use
74
75 The multi-line comments are made by Thomas R. Nicely and are copied verbatim.
76 Additional comments marked "CZ" (without the quotes) are by the code-portist.
77
78 (If that name sounds familiar, he is the guy who found the fdiv bug in the
79 Pentium (P5x, I think) Intel processor)
80 */
mp_prime_strong_lucas_selfridge(const mp_int * a,int * result)81 int mp_prime_strong_lucas_selfridge(const mp_int *a, int *result)
82 {
83 /* CZ TODO: choose better variable names! */
84 mp_int Dz, gcd, Np1, Uz, Vz, U2mz, V2mz, Qmz, Q2mz, Qkdz, T1z, T2z, T3z, T4z, Q2kdz;
85 /* CZ TODO: Some of them need the full 32 bit, hence the (temporary) exclusion of MP_8BIT */
86 int D, Ds, J, sign, P, Q, r, s, u, Nbits;
87 int e;
88 int isset;
89
90 *result = MP_NO;
91
92 /*
93 Find the first element D in the sequence {5, -7, 9, -11, 13, ...}
94 such that Jacobi(D,N) = -1 (Selfridge's algorithm). Theory
95 indicates that, if N is not a perfect square, D will "nearly
96 always" be "small." Just in case, an overflow trap for D is
97 included.
98 */
99
100 if ((e = mp_init_multi(&Dz, &gcd, &Np1, &Uz, &Vz, &U2mz, &V2mz, &Qmz, &Q2mz, &Qkdz, &T1z, &T2z, &T3z, &T4z, &Q2kdz,
101 NULL)) != MP_OKAY) {
102 return e;
103 }
104
105 D = 5;
106 sign = 1;
107
108 for (;;) {
109 Ds = sign * D;
110 sign = -sign;
111 if ((e = mp_set_long(&Dz, (unsigned long)D)) != MP_OKAY) {
112 goto LBL_LS_ERR;
113 }
114 if ((e = mp_gcd(a, &Dz, &gcd)) != MP_OKAY) {
115 goto LBL_LS_ERR;
116 }
117 /* if 1 < GCD < N then N is composite with factor "D", and
118 Jacobi(D,N) is technically undefined (but often returned
119 as zero). */
120 if ((mp_cmp_d(&gcd, 1uL) == MP_GT) && (mp_cmp(&gcd, a) == MP_LT)) {
121 goto LBL_LS_ERR;
122 }
123 if (Ds < 0) {
124 Dz.sign = MP_NEG;
125 }
126 if ((e = mp_kronecker(&Dz, a, &J)) != MP_OKAY) {
127 goto LBL_LS_ERR;
128 }
129
130 if (J == -1) {
131 break;
132 }
133 D += 2;
134
135 if (D > (INT_MAX - 2)) {
136 e = MP_VAL;
137 goto LBL_LS_ERR;
138 }
139 }
140
141 P = 1; /* Selfridge's choice */
142 Q = (1 - Ds) / 4; /* Required so D = P*P - 4*Q */
143
144 /* NOTE: The conditions (a) N does not divide Q, and
145 (b) D is square-free or not a perfect square, are included by
146 some authors; e.g., "Prime numbers and computer methods for
147 factorization," Hans Riesel (2nd ed., 1994, Birkhauser, Boston),
148 p. 130. For this particular application of Lucas sequences,
149 these conditions were found to be immaterial. */
150
151 /* Now calculate N - Jacobi(D,N) = N + 1 (even), and calculate the
152 odd positive integer d and positive integer s for which
153 N + 1 = 2^s*d (similar to the step for N - 1 in Miller's test).
154 The strong Lucas-Selfridge test then returns N as a strong
155 Lucas probable prime (slprp) if any of the following
156 conditions is met: U_d=0, V_d=0, V_2d=0, V_4d=0, V_8d=0,
157 V_16d=0, ..., etc., ending with V_{2^(s-1)*d}=V_{(N+1)/2}=0
158 (all equalities mod N). Thus d is the highest index of U that
159 must be computed (since V_2m is independent of U), compared
160 to U_{N+1} for the standard Lucas-Selfridge test; and no
161 index of V beyond (N+1)/2 is required, just as in the
162 standard Lucas-Selfridge test. However, the quantity Q^d must
163 be computed for use (if necessary) in the latter stages of
164 the test. The result is that the strong Lucas-Selfridge test
165 has a running time only slightly greater (order of 10 %) than
166 that of the standard Lucas-Selfridge test, while producing
167 only (roughly) 30 % as many pseudoprimes (and every strong
168 Lucas pseudoprime is also a standard Lucas pseudoprime). Thus
169 the evidence indicates that the strong Lucas-Selfridge test is
170 more effective than the standard Lucas-Selfridge test, and a
171 Baillie-PSW test based on the strong Lucas-Selfridge test
172 should be more reliable. */
173
174 if ((e = mp_add_d(a, 1uL, &Np1)) != MP_OKAY) {
175 goto LBL_LS_ERR;
176 }
177 s = mp_cnt_lsb(&Np1);
178
179 /* CZ
180 * This should round towards zero because
181 * Thomas R. Nicely used GMP's mpz_tdiv_q_2exp()
182 * and mp_div_2d() is equivalent. Additionally:
183 * dividing an even number by two does not produce
184 * any leftovers.
185 */
186 if ((e = mp_div_2d(&Np1, s, &Dz, NULL)) != MP_OKAY) {
187 goto LBL_LS_ERR;
188 }
189 /* We must now compute U_d and V_d. Since d is odd, the accumulated
190 values U and V are initialized to U_1 and V_1 (if the target
191 index were even, U and V would be initialized instead to U_0=0
192 and V_0=2). The values of U_2m and V_2m are also initialized to
193 U_1 and V_1; the FOR loop calculates in succession U_2 and V_2,
194 U_4 and V_4, U_8 and V_8, etc. If the corresponding bits
195 (1, 2, 3, ...) of t are on (the zero bit having been accounted
196 for in the initialization of U and V), these values are then
197 combined with the previous totals for U and V, using the
198 composition formulas for addition of indices. */
199
200 mp_set(&Uz, 1uL); /* U=U_1 */
201 mp_set(&Vz, (mp_digit)P); /* V=V_1 */
202 mp_set(&U2mz, 1uL); /* U_1 */
203 mp_set(&V2mz, (mp_digit)P); /* V_1 */
204
205 if (Q < 0) {
206 Q = -Q;
207 if ((e = mp_set_long(&Qmz, (unsigned long)Q)) != MP_OKAY) {
208 goto LBL_LS_ERR;
209 }
210 if ((e = mp_mul_2(&Qmz, &Q2mz)) != MP_OKAY) {
211 goto LBL_LS_ERR;
212 }
213 /* Initializes calculation of Q^d */
214 if ((e = mp_set_long(&Qkdz, (unsigned long)Q)) != MP_OKAY) {
215 goto LBL_LS_ERR;
216 }
217 Qmz.sign = MP_NEG;
218 Q2mz.sign = MP_NEG;
219 Qkdz.sign = MP_NEG;
220 Q = -Q;
221 } else {
222 if ((e = mp_set_long(&Qmz, (unsigned long)Q)) != MP_OKAY) {
223 goto LBL_LS_ERR;
224 }
225 if ((e = mp_mul_2(&Qmz, &Q2mz)) != MP_OKAY) {
226 goto LBL_LS_ERR;
227 }
228 /* Initializes calculation of Q^d */
229 if ((e = mp_set_long(&Qkdz, (unsigned long)Q)) != MP_OKAY) {
230 goto LBL_LS_ERR;
231 }
232 }
233
234 Nbits = mp_count_bits(&Dz);
235 for (u = 1; u < Nbits; u++) { /* zero bit off, already accounted for */
236 /* Formulas for doubling of indices (carried out mod N). Note that
237 * the indices denoted as "2m" are actually powers of 2, specifically
238 * 2^(ul-1) beginning each loop and 2^ul ending each loop.
239 *
240 * U_2m = U_m*V_m
241 * V_2m = V_m*V_m - 2*Q^m
242 */
243
244 if ((e = mp_mul(&U2mz, &V2mz, &U2mz)) != MP_OKAY) {
245 goto LBL_LS_ERR;
246 }
247 if ((e = mp_mod(&U2mz, a, &U2mz)) != MP_OKAY) {
248 goto LBL_LS_ERR;
249 }
250 if ((e = mp_sqr(&V2mz, &V2mz)) != MP_OKAY) {
251 goto LBL_LS_ERR;
252 }
253 if ((e = mp_sub(&V2mz, &Q2mz, &V2mz)) != MP_OKAY) {
254 goto LBL_LS_ERR;
255 }
256 if ((e = mp_mod(&V2mz, a, &V2mz)) != MP_OKAY) {
257 goto LBL_LS_ERR;
258 }
259 /* Must calculate powers of Q for use in V_2m, also for Q^d later */
260 if ((e = mp_sqr(&Qmz, &Qmz)) != MP_OKAY) {
261 goto LBL_LS_ERR;
262 }
263 /* prevents overflow */ /* CZ still necessary without a fixed prealloc'd mem.? */
264 if ((e = mp_mod(&Qmz, a, &Qmz)) != MP_OKAY) {
265 goto LBL_LS_ERR;
266 }
267 if ((e = mp_mul_2(&Qmz, &Q2mz)) != MP_OKAY) {
268 goto LBL_LS_ERR;
269 }
270
271 if ((isset = mp_get_bit(&Dz, u)) == MP_VAL) {
272 e = isset;
273 goto LBL_LS_ERR;
274 }
275 if (isset == MP_YES) {
276 /* Formulas for addition of indices (carried out mod N);
277 *
278 * U_(m+n) = (U_m*V_n + U_n*V_m)/2
279 * V_(m+n) = (V_m*V_n + D*U_m*U_n)/2
280 *
281 * Be careful with division by 2 (mod N)!
282 */
283
284 if ((e = mp_mul(&U2mz, &Vz, &T1z)) != MP_OKAY) {
285 goto LBL_LS_ERR;
286 }
287 if ((e = mp_mul(&Uz, &V2mz, &T2z)) != MP_OKAY) {
288 goto LBL_LS_ERR;
289 }
290 if ((e = mp_mul(&V2mz, &Vz, &T3z)) != MP_OKAY) {
291 goto LBL_LS_ERR;
292 }
293 if ((e = mp_mul(&U2mz, &Uz, &T4z)) != MP_OKAY) {
294 goto LBL_LS_ERR;
295 }
296 if ((e = s_mp_mul_si(&T4z, (long)Ds, &T4z)) != MP_OKAY) {
297 goto LBL_LS_ERR;
298 }
299 if ((e = mp_add(&T1z, &T2z, &Uz)) != MP_OKAY) {
300 goto LBL_LS_ERR;
301 }
302 if (mp_isodd(&Uz) != MP_NO) {
303 if ((e = mp_add(&Uz, a, &Uz)) != MP_OKAY) {
304 goto LBL_LS_ERR;
305 }
306 }
307 /* CZ
308 * This should round towards negative infinity because
309 * Thomas R. Nicely used GMP's mpz_fdiv_q_2exp().
310 * But mp_div_2() does not do so, it is truncating instead.
311 */
312 if ((e = mp_div_2(&Uz, &Uz)) != MP_OKAY) {
313 goto LBL_LS_ERR;
314 }
315 if ((Uz.sign == MP_NEG) && (mp_isodd(&Uz) != MP_NO)) {
316 if ((e = mp_sub_d(&Uz, 1uL, &Uz)) != MP_OKAY) {
317 goto LBL_LS_ERR;
318 }
319 }
320 if ((e = mp_add(&T3z, &T4z, &Vz)) != MP_OKAY) {
321 goto LBL_LS_ERR;
322 }
323 if (mp_isodd(&Vz) != MP_NO) {
324 if ((e = mp_add(&Vz, a, &Vz)) != MP_OKAY) {
325 goto LBL_LS_ERR;
326 }
327 }
328 if ((e = mp_div_2(&Vz, &Vz)) != MP_OKAY) {
329 goto LBL_LS_ERR;
330 }
331 if ((Vz.sign == MP_NEG) && (mp_isodd(&Vz) != MP_NO)) {
332 if ((e = mp_sub_d(&Vz, 1uL, &Vz)) != MP_OKAY) {
333 goto LBL_LS_ERR;
334 }
335 }
336 if ((e = mp_mod(&Uz, a, &Uz)) != MP_OKAY) {
337 goto LBL_LS_ERR;
338 }
339 if ((e = mp_mod(&Vz, a, &Vz)) != MP_OKAY) {
340 goto LBL_LS_ERR;
341 }
342 /* Calculating Q^d for later use */
343 if ((e = mp_mul(&Qkdz, &Qmz, &Qkdz)) != MP_OKAY) {
344 goto LBL_LS_ERR;
345 }
346 if ((e = mp_mod(&Qkdz, a, &Qkdz)) != MP_OKAY) {
347 goto LBL_LS_ERR;
348 }
349 }
350 }
351
352 /* If U_d or V_d is congruent to 0 mod N, then N is a prime or a
353 strong Lucas pseudoprime. */
354 if ((mp_iszero(&Uz) != MP_NO) || (mp_iszero(&Vz) != MP_NO)) {
355 *result = MP_YES;
356 goto LBL_LS_ERR;
357 }
358
359 /* NOTE: Ribenboim ("The new book of prime number records," 3rd ed.,
360 1995/6) omits the condition V0 on p.142, but includes it on
361 p. 130. The condition is NECESSARY; otherwise the test will
362 return false negatives---e.g., the primes 29 and 2000029 will be
363 returned as composite. */
364
365 /* Otherwise, we must compute V_2d, V_4d, V_8d, ..., V_{2^(s-1)*d}
366 by repeated use of the formula V_2m = V_m*V_m - 2*Q^m. If any of
367 these are congruent to 0 mod N, then N is a prime or a strong
368 Lucas pseudoprime. */
369
370 /* Initialize 2*Q^(d*2^r) for V_2m */
371 if ((e = mp_mul_2(&Qkdz, &Q2kdz)) != MP_OKAY) {
372 goto LBL_LS_ERR;
373 }
374
375 for (r = 1; r < s; r++) {
376 if ((e = mp_sqr(&Vz, &Vz)) != MP_OKAY) {
377 goto LBL_LS_ERR;
378 }
379 if ((e = mp_sub(&Vz, &Q2kdz, &Vz)) != MP_OKAY) {
380 goto LBL_LS_ERR;
381 }
382 if ((e = mp_mod(&Vz, a, &Vz)) != MP_OKAY) {
383 goto LBL_LS_ERR;
384 }
385 if (mp_iszero(&Vz) != MP_NO) {
386 *result = MP_YES;
387 goto LBL_LS_ERR;
388 }
389 /* Calculate Q^{d*2^r} for next r (final iteration irrelevant). */
390 if (r < (s - 1)) {
391 if ((e = mp_sqr(&Qkdz, &Qkdz)) != MP_OKAY) {
392 goto LBL_LS_ERR;
393 }
394 if ((e = mp_mod(&Qkdz, a, &Qkdz)) != MP_OKAY) {
395 goto LBL_LS_ERR;
396 }
397 if ((e = mp_mul_2(&Qkdz, &Q2kdz)) != MP_OKAY) {
398 goto LBL_LS_ERR;
399 }
400 }
401 }
402 LBL_LS_ERR:
403 mp_clear_multi(&Q2kdz, &T4z, &T3z, &T2z, &T1z, &Qkdz, &Q2mz, &Qmz, &V2mz, &U2mz, &Vz, &Uz, &Np1, &gcd, &Dz, NULL);
404 return e;
405 }
406 #endif
407 #endif
408 #endif
409
410 /* ref: $Format:%D$ */
411 /* git commit: $Format:%H$ */
412 /* commit time: $Format:%ai$ */
413