Mercurial > dropbear

comparison libtommath/bn_mp_prime_strong_lucas_selfridge.c @ 1655:f52919ffd3b1

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