Line data Source code
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 0 : static mp_err s_mp_mul_si(const mp_int *a, int32_t d, mp_int *c)
23 : {
24 0 : mp_int t;
25 0 : mp_err err;
26 :
27 0 : if ((err = mp_init(&t)) != MP_OKAY) {
28 0 : 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 0 : mp_set_i32(&t, d);
36 0 : err = mp_mul(a, &t, c);
37 0 : mp_clear(&t);
38 0 : 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 0 : mp_err mp_prime_strong_lucas_selfridge(const mp_int *a, mp_bool *result)
58 : {
59 : /* CZ TODO: choose better variable names! */
60 0 : 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 0 : int32_t D, Ds, J, sign, P, Q, r, s, u, Nbits;
63 0 : mp_err err;
64 0 : mp_bool oddness;
65 :
66 0 : *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 0 : 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 0 : return err;
78 : }
79 :
80 0 : D = 5;
81 0 : sign = 1;
82 :
83 0 : for (;;) {
84 0 : Ds = sign * D;
85 0 : sign = -sign;
86 0 : mp_set_u32(&Dz, (uint32_t)D);
87 0 : 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 0 : if ((mp_cmp_d(&gcd, 1uL) == MP_GT) && (mp_cmp(&gcd, a) == MP_LT)) {
93 0 : goto LBL_LS_ERR;
94 : }
95 0 : if (Ds < 0) {
96 0 : Dz.sign = MP_NEG;
97 : }
98 0 : if ((err = mp_kronecker(&Dz, a, &J)) != MP_OKAY) goto LBL_LS_ERR;
99 :
100 0 : if (J == -1) {
101 0 : break;
102 : }
103 0 : D += 2;
104 :
105 0 : if (D > (INT_MAX - 2)) {
106 0 : err = MP_VAL;
107 0 : goto LBL_LS_ERR;
108 : }
109 : }
110 :
111 :
112 :
113 0 : P = 1; /* Selfridge's choice */
114 0 : 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 0 : if ((err = mp_add_d(a, 1uL, &Np1)) != MP_OKAY) goto LBL_LS_ERR;
147 0 : 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 0 : 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 0 : mp_set(&Uz, 1uL); /* U=U_1 */
169 0 : mp_set(&Vz, (mp_digit)P); /* V=V_1 */
170 0 : mp_set(&U2mz, 1uL); /* U_1 */
171 0 : mp_set(&V2mz, (mp_digit)P); /* V_1 */
172 :
173 0 : mp_set_i32(&Qmz, Q);
174 0 : if ((err = mp_mul_2(&Qmz, &Q2mz)) != MP_OKAY) goto LBL_LS_ERR;
175 : /* Initializes calculation of Q^d */
176 0 : mp_set_i32(&Qkdz, Q);
177 :
178 0 : Nbits = mp_count_bits(&Dz);
179 :
180 0 : 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 0 : if ((err = mp_mul(&U2mz, &V2mz, &U2mz)) != MP_OKAY) goto LBL_LS_ERR;
190 0 : if ((err = mp_mod(&U2mz, a, &U2mz)) != MP_OKAY) goto LBL_LS_ERR;
191 0 : if ((err = mp_sqr(&V2mz, &V2mz)) != MP_OKAY) goto LBL_LS_ERR;
192 0 : if ((err = mp_sub(&V2mz, &Q2mz, &V2mz)) != MP_OKAY) goto LBL_LS_ERR;
193 0 : 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 0 : 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 0 : if ((err = mp_mod(&Qmz, a, &Qmz)) != MP_OKAY) goto LBL_LS_ERR;
200 0 : if ((err = mp_mul_2(&Qmz, &Q2mz)) != MP_OKAY) goto LBL_LS_ERR;
201 :
202 0 : 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 0 : if ((err = mp_mul(&U2mz, &Vz, &T1z)) != MP_OKAY) goto LBL_LS_ERR;
211 0 : if ((err = mp_mul(&Uz, &V2mz, &T2z)) != MP_OKAY) goto LBL_LS_ERR;
212 0 : if ((err = mp_mul(&V2mz, &Vz, &T3z)) != MP_OKAY) goto LBL_LS_ERR;
213 0 : if ((err = mp_mul(&U2mz, &Uz, &T4z)) != MP_OKAY) goto LBL_LS_ERR;
214 0 : if ((err = s_mp_mul_si(&T4z, Ds, &T4z)) != MP_OKAY) goto LBL_LS_ERR;
215 0 : if ((err = mp_add(&T1z, &T2z, &Uz)) != MP_OKAY) goto LBL_LS_ERR;
216 0 : if (MP_IS_ODD(&Uz)) {
217 0 : 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 0 : oddness = MP_IS_ODD(&Uz) ? MP_YES : MP_NO;
225 0 : if ((err = mp_div_2(&Uz, &Uz)) != MP_OKAY) goto LBL_LS_ERR;
226 0 : if ((Uz.sign == MP_NEG) && (oddness != MP_NO)) {
227 0 : if ((err = mp_sub_d(&Uz, 1uL, &Uz)) != MP_OKAY) goto LBL_LS_ERR;
228 : }
229 0 : if ((err = mp_add(&T3z, &T4z, &Vz)) != MP_OKAY) goto LBL_LS_ERR;
230 0 : if (MP_IS_ODD(&Vz)) {
231 0 : if ((err = mp_add(&Vz, a, &Vz)) != MP_OKAY) goto LBL_LS_ERR;
232 : }
233 0 : oddness = MP_IS_ODD(&Vz) ? MP_YES : MP_NO;
234 0 : if ((err = mp_div_2(&Vz, &Vz)) != MP_OKAY) goto LBL_LS_ERR;
235 0 : if ((Vz.sign == MP_NEG) && (oddness != MP_NO)) {
236 0 : if ((err = mp_sub_d(&Vz, 1uL, &Vz)) != MP_OKAY) goto LBL_LS_ERR;
237 : }
238 0 : if ((err = mp_mod(&Uz, a, &Uz)) != MP_OKAY) goto LBL_LS_ERR;
239 0 : if ((err = mp_mod(&Vz, a, &Vz)) != MP_OKAY) goto LBL_LS_ERR;
240 :
241 : /* Calculating Q^d for later use */
242 0 : if ((err = mp_mul(&Qkdz, &Qmz, &Qkdz)) != MP_OKAY) goto LBL_LS_ERR;
243 0 : 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 0 : if (MP_IS_ZERO(&Uz) || MP_IS_ZERO(&Vz)) {
250 0 : *result = MP_YES;
251 0 : 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 0 : if ((err = mp_mul_2(&Qkdz, &Q2kdz)) != MP_OKAY) goto LBL_LS_ERR;
267 :
268 0 : for (r = 1; r < s; r++) {
269 0 : if ((err = mp_sqr(&Vz, &Vz)) != MP_OKAY) goto LBL_LS_ERR;
270 0 : if ((err = mp_sub(&Vz, &Q2kdz, &Vz)) != MP_OKAY) goto LBL_LS_ERR;
271 0 : if ((err = mp_mod(&Vz, a, &Vz)) != MP_OKAY) goto LBL_LS_ERR;
272 0 : if (MP_IS_ZERO(&Vz)) {
273 0 : *result = MP_YES;
274 0 : goto LBL_LS_ERR;
275 : }
276 : /* Calculate Q^{d*2^r} for next r (final iteration irrelevant). */
277 0 : if (r < (s - 1)) {
278 0 : if ((err = mp_sqr(&Qkdz, &Qkdz)) != MP_OKAY) goto LBL_LS_ERR;
279 0 : if ((err = mp_mod(&Qkdz, a, &Qkdz)) != MP_OKAY) goto LBL_LS_ERR;
280 0 : if ((err = mp_mul_2(&Qkdz, &Q2kdz)) != MP_OKAY) goto LBL_LS_ERR;
281 : }
282 : }
283 0 : LBL_LS_ERR:
284 0 : mp_clear_multi(&Q2kdz, &T4z, &T3z, &T2z, &T1z, &Qkdz, &Q2mz, &Qmz, &V2mz, &U2mz, &Vz, &Uz, &Np1, &gcd, &Dz, NULL);
285 0 : return err;
286 : }
287 : #endif
288 : #endif
289 : #endif
|
