1 #include <string.h>
2 #include <gmp.h>
3 #include "ptypes.h"
4
5 #include "bls75.h"
6 #include "primality.h"
7 #include "prime_iterator.h"
8 #include "pbrent63.h"
9 #include "squfof126.h"
10 #include "factor.h"
11 #include "simpqs.h"
12 #include "ecm.h"
13 #define _GMP_ECM_FACTOR(n, f, b1, ncurves) \
14 _GMP_ecm_factor_projective(n, f, b1, 0, ncurves)
15 #include "utility.h"
16
17 /*
18 * Lucas (1876): Given a completely factored n-1, if there exists an a s.t.
19 * a^(n-1) % n = 1
20 * a^((n-1/f) % n != 1 for ALL factors f of n-1
21 * then n is prime.
22 *
23 * PPBLS:, given n-1 = A*B, A > sqrt(n), if we can find an a s.t.
24 * a^A % n = 1
25 * gcd(a^(A/f)-1,n) = 1 for ALL factors f of A
26 * then n is prime.
27 *
28 * Generalized Pocklington: given n-1 = A*B, gcd(A,B)=1, A > sqrt(n), then if
29 * for each each factor f of A, there exists an a (1 < a < n-1) s.t.
30 * a^(n-1) % n = 1
31 * gcd(a^((n-1)/f)-1,n) = 1
32 * then n is prime.
33 *
34 * BLS T5: given n-1 = A*B, factored A, s=B/2A r=B mod (2A), and an a, then if:
35 * - A is even, B is odd, and AB=n-1 (all implied by n = odd and the above),
36 * - n < (A+1) * (2*A*A + (r-1) * A + 1)
37 * - for each each factor f of A, there exists an a (1 < a < n-1) s.t.
38 * - a^(n-1) % n = 1
39 * - gcd(a^((n-1)/f)-1,n) = 1 for ALL factors f of A
40 * then:
41 * if s = 0 or r*r - 8*s is not a perfect square
42 * n is prime
43 * else
44 * n is composite
45 *
46 * The generalized Pocklington test is also sometimes known as the
47 * Pocklington-Lehmer test. It's definitely an improvement over Lucas
48 * since we only have to find factors up to sqrt(n), _and_ we can choose
49 * a different 'a' value for each factor. This is corollary 1 from BLS75.
50 *
51 * BLS is the Brillhart-Lehmer-Selfridge 1975 theorem 5 (see link below).
52 * We can factor even less of n, and the test lets us kick out some
53 * composites early, without having to test n-3 different 'a' values.
54 *
55 * Once we've found the factors of n-1 (or enough of them), verification
56 * usually happens really fast. a=2 works for most, and few seem to require
57 * more than ~ log2(n). However all but BLS75 require testing all integers
58 * 1 < a < n-1 before answering in the negative, which is impractical.
59 *
60 * BLS75 theorem 7 is the final n-1 theorem and takes into account any
61 * knowledge that the remaining factor is not below a threshold B. Since
62 * we do initial trial division this helps. It is usually of only small
63 * benefit.
64 *
65 *
66 * AKS is not too hard to implement, but it's impractically slow.
67 *
68 * ECPP is very fast and definitely the best method for most numbers.
69 *
70 * BLS75: http://www.ams.org/journals/mcom/1975-29-130/S0025-5718-1975-0384673-1/S0025-5718-1975-0384673-1.pdf
71 *
72 */
73
74 /* Like all the primality functions:
75 * 2 = definitely prime, 1 = maybe prime, 0 = definitely composite
76 *
77 * You really should run is_prob_prime on n first, so we only have to run
78 * these tests on numbers that are very probably prime.
79 */
80
81
tfe(mpz_t f,mpz_t n,int effort)82 static int tfe(mpz_t f, mpz_t n, int effort)
83 {
84 int success = 0;
85 UV log2n = mpz_sizeinbase(n, 2);
86
87 if (mpz_cmp_ui(n,3) <= 0) {
88 mpz_set(f,n);
89 return 1;
90 }
91
92 if (effort == 0) {
93 if (!success && log2n <= 63) success = pbrent63(n, f, 1000000);
94 if (!success) success = (int)power_factor(n, f);
95 if (success) return success;
96 }
97 /* TODO: Use tinyqs effectively here, e.g. stage 3 for 50-90 bit */
98
99 switch (effort) {
100 case 0: success = _GMP_pminus1_factor(n, f, 400, 400);
101 break;
102 case 1: success = _GMP_pminus1_factor(n, f, 1000, 11000);
103 break;
104 case 2: { UV brent_rounds = (log2n <= 64) ? 100000 : 100000 / (log2n-63);
105 int final_B2 = 1000 * (150-(int)log2n);
106 if (log2n < 70) brent_rounds *= 3;
107 if (log2n < 80)
108 success = _GMP_ECM_FACTOR(n, f, 150, 5);
109 if (!success)
110 success = _GMP_pbrent_factor(n, f, 3, brent_rounds);
111 if (!success && final_B2 > 11000)
112 success = _GMP_pminus1_factor(n, f, 10000, final_B2);
113 } break;
114 case 3: success = _GMP_pminus1_factor(n, f, 20000, 200000);
115 if (!success && log2n <= 80) success = squfof126(n, f, 2000000);
116 break;
117 case 4: success = _GMP_ECM_FACTOR(n, f, 500, 30);
118 break;
119 case 5: success = _GMP_ECM_FACTOR(n, f, 2000, 20);
120 break;
121 case 6: if (log2n > 170) {
122 UV B1 = (log2n > 2500) ? 10000000 : 4000 * log2n;
123 success = _GMP_pminus1_factor(n, f, B1, 20*B1);
124 } break;
125 case 7: if (log2n > 210) { success = _GMP_ECM_FACTOR(n, f, 20000, 10); }
126 break;
127 case 8: if (log2n > 240) { success = _GMP_ECM_FACTOR(n, f, 40000, 10); }
128 break;
129 case 9: if (log2n > 240) { success = _GMP_ECM_FACTOR(n, f, 80000, 5); }
130 break;
131 case 10:if (log2n > 270) { success = _GMP_ECM_FACTOR(n, f,160000, 20); }
132 break;
133
134 /* QS for sizes 30-90 digits */
135 case 20:
136 case 21:{ UV log10n = mpz_sizeinbase(n, 10);
137 if (log10n >= 30 && log10n <= ((effort == 20) ? 54 : 90)) {
138 mpz_t farray[66];
139 int i, nfactors;
140 for (i = 0; i < 66; i++) mpz_init(farray[i]);
141 nfactors = _GMP_simpqs(n, farray);
142 /* TODO: Return all factors */
143 if (nfactors > 1) {
144 success = 1;
145 mpz_set(f, farray[nfactors-1]); /* Return largest */
146 }
147 for (i = 0; i < 66; i++) mpz_clear(farray[i]);
148 }
149 } break;
150
151 case 30:success = _GMP_pminus1_factor(n, f, 200000, 4000000);
152 break;
153
154 case 40:
155 case 41:
156 case 42:
157 case 43:
158 case 44:
159 case 45:
160 case 46:
161 case 47:
162 case 48:
163 case 49:
164 case 50:
165 case 51:
166 case 52:
167 case 53:
168 case 54:
169 case 55:
170 case 56:
171 case 57:{ UV B1 = UVCONST(1) << (13 + (effort-40));
172 success = _GMP_ECM_FACTOR(n, f, B1, 10);
173 } break;
174
175 default: break;
176 }
177
178 /* if (success) printf(" bls75 factored %lu-bit at effort %d\n", log2n, effort); */
179 return success;
180 }
181
182 /* F*R = n, F is factored part, R is remainder */
small_factor(mpz_t F,mpz_t R,UV B1)183 static void small_factor(mpz_t F, mpz_t R, UV B1)
184 {
185 PRIME_ITERATOR(iter);
186 UV tf;
187 for (tf = 2; tf < B1; tf = prime_iterator_next(&iter)) {
188 if (mpz_cmp_ui(R, tf*tf) < 0) break;
189 if (mpz_divisible_ui_p(R, tf)) {
190 do {
191 mpz_mul_ui(F, F, tf);
192 mpz_divexact_ui(R, R, tf);
193 } while (mpz_divisible_ui_p(R, tf));
194 }
195 }
196 prime_iterator_destroy(&iter);
197 }
198
199
200 typedef struct {
201 int cur;
202 int max;
203 mpz_t* stack;
204 } fstack_t;
205
206 #define FACTOR_STACK(name) fstack_t name = {0, 0, 0}
207
nstack(fstack_t * s)208 static int nstack(fstack_t* s) { return s->cur; }
push_fstack(fstack_t * s,mpz_t v)209 static void push_fstack(fstack_t* s, mpz_t v) {
210 if (s->stack == 0) New(0,s->stack, s->max = 10, mpz_t);
211 if (s->cur == s->max) Renew(s->stack, s->max += 10, mpz_t);
212 mpz_init_set(s->stack[(s->cur)++], v);
213 }
push_fstack_ui(fstack_t * s,unsigned long v)214 static void push_fstack_ui(fstack_t* s, unsigned long v) {
215 if (s->stack == 0) New(0,s->stack, s->max = 10, mpz_t);
216 if (s->cur == s->max) Renew(s->stack, s->max += 10, mpz_t);
217 mpz_init_set_ui(s->stack[(s->cur)++], v);
218 }
pop_fstack(mpz_t rv,fstack_t * s)219 static void pop_fstack(mpz_t rv, fstack_t* s) {
220 mpz_set(rv, s->stack[--(s->cur)]);
221 mpz_clear(s->stack[s->cur]);
222 }
clear_fstack(fstack_t * s)223 static void clear_fstack(fstack_t* s) {
224 while (s->cur > 0)
225 mpz_clear(s->stack[--(s->cur)]);
226 }
destroy_fstack(fstack_t * s)227 static void destroy_fstack(fstack_t* s) {
228 clear_fstack(s);
229 Safefree(s->stack);
230 s->stack = 0;
231 }
232
factor_out(mpz_t R,mpz_t F,mpz_t v)233 static void factor_out(mpz_t R, mpz_t F, mpz_t v) {
234 int ndiv = mpz_remove(R, R, v);
235 while (ndiv-- > 0)
236 mpz_mul(F, F, v);
237 }
factor_out_ui(mpz_t R,mpz_t F,unsigned long v)238 static void factor_out_ui(mpz_t R, mpz_t F, unsigned long v) {
239 while (mpz_divisible_ui_p(R, v)) {
240 mpz_mul_ui(F, F, v);
241 mpz_divexact_ui(R, R, v);
242 }
243 }
244
factor_test_ui(unsigned long f,mpz_t R,mpz_t F,fstack_t * s)245 static void factor_test_ui(unsigned long f, mpz_t R, mpz_t F, fstack_t* s) {
246 if (mpz_divisible_ui_p(R, f)) {
247 push_fstack_ui(s, f);
248 factor_out_ui(R, F, f);
249 }
250 }
251
252 typedef int (*bls_func_t)(mpz_t, int, char**);
253 typedef int (*limit_func_t)(mpz_t, mpz_t, mpz_t, UV, mpz_t,mpz_t,mpz_t,mpz_t);
254
handle_factor(mpz_t f,mpz_t R,mpz_t F,fstack_t * sf,fstack_t * sm,int effort,char ** prtext,int push_if_probable,bls_func_t func)255 static void handle_factor(mpz_t f, mpz_t R, mpz_t F,
256 fstack_t* sf, fstack_t* sm,
257 int effort, char** prtext,
258 int push_if_probable,
259 bls_func_t func) {
260 int pr = _GMP_BPSW(f);
261 if (pr == 1) { /* Try to prove */
262 if (effort > 1 || mpz_sizeinbase(f,2) < 200) {
263 pr = (*func)(f, effort, prtext);
264 }
265 }
266 if (pr == 2) {
267 push_fstack(sf, f);
268 factor_out(R, F, f);
269 } else if (pr == 0 || push_if_probable) {
270 push_fstack(sm, f);
271 }
272 }
handle_factor2(mpz_t f,mpz_t R,mpz_t F,fstack_t * sf,fstack_t * sp,fstack_t * sm,int effort,char ** prtext,bls_func_t func)273 static void handle_factor2(mpz_t f, mpz_t R, mpz_t F,
274 fstack_t* sf, fstack_t* sp, fstack_t* sm,
275 int effort, char** prtext,
276 bls_func_t func) {
277 int pr = _GMP_BPSW(f);
278 if (pr == 1) { /* Try to prove */
279 pr = (*func)(f, effort, prtext);
280 }
281 if (pr == 0) {
282 push_fstack(sm, f);
283 } else if (pr == 2) {
284 push_fstack(sf, f);
285 factor_out(R, F, f);
286 } else {
287 /* Save actual proof for later, but for now assume we can do it */
288 push_fstack(sp, f);
289 factor_out(R, F, f);
290 }
291 }
292
trim_factors(mpz_t F,mpz_t R,mpz_t n,mpz_t none,UV flim,fstack_t * fs,limit_func_t func,mpz_t t,mpz_t m,mpz_t r,mpz_t s)293 static void trim_factors(mpz_t F, mpz_t R, mpz_t n, mpz_t none, UV flim, fstack_t* fs, limit_func_t func, mpz_t t, mpz_t m, mpz_t r, mpz_t s) {
294 if (fs->cur > 1) {
295 int i;
296 mpz_set_ui(F, 1);
297 mpz_set(R, none);
298 for (i = 0; i < fs->cur; i++) {
299 if (i > 0 && func(n, F, R, flim, t, m, r, s))
300 break;
301 factor_out(R, F, fs->stack[i]);
302 }
303 /* Remove excess factors */
304 while (i < fs->cur)
305 pop_fstack(t, fs);
306 }
307 /* Verify Q[0] = 2 */
308 if (mpz_cmp_ui(fs->stack[0], 2) != 0)
309 croak("BLS75 internal error: 2 not at start of fstack");
310 /* r and s have been set by func */
311 }
312
313 /* Sort factors found from largest to smallest, but 2 must be at start. */
sort_and_trim_factors(int * fsp,mpz_t * fstack)314 static void sort_and_trim_factors(int* fsp, mpz_t* fstack)
315 {
316 int i, j;
317 for (i = 2; i < *fsp; i++)
318 for (j = i; j > 1 && mpz_cmp(fstack[j-1], fstack[j]) < 0; j--)
319 mpz_swap( fstack[j-1], fstack[j] );
320 for (i = 2; i < *fsp; i++) { /* Remove any duplicate factors */
321 if (mpz_cmp(fstack[i], fstack[i-1]) == 0) {
322 for (j = i+1; j < *fsp; j++)
323 mpz_set(fstack[j-1], fstack[j]);
324 *fsp -= 1;
325 }
326 }
327 }
fstack_sort_trim(fstack_t * s)328 static void fstack_sort_trim(fstack_t* s) {
329 sort_and_trim_factors(&(s->cur), s->stack);
330 }
331
332
333
334 /******************************************************************************/
335
336
bls_theorem5_limit(mpz_t n,mpz_t A,mpz_t B,UV dummy,mpz_t t,mpz_t y,mpz_t r,mpz_t s)337 static int bls_theorem5_limit(mpz_t n, mpz_t A, mpz_t B, UV dummy,
338 mpz_t t, mpz_t y, mpz_t r, mpz_t s)
339 {
340 mpz_mul(t, A, B);
341 mpz_add_ui(t, t, 1);
342 if (mpz_cmp(t, n) != 0) croak("BLS75 internal error: A*B != n-1\n");
343
344 mpz_mul_ui(t, A, 2);
345 mpz_tdiv_qr(s, r, B, t);
346
347 mpz_mul(y, t, A); /* y = 2*A*A */
348 mpz_sub_ui(t, r, 1); /* t = r-1 */
349 mpz_mul(t, t, A); /* t = A*(r-1) */
350 mpz_add(y, y, t); /* y = 2A^2 + A(r-1) */
351 mpz_add_ui(y, y, 1); /* y = 2A^2 + A(r-1) + 1 */
352 mpz_add_ui(t, A, 1); /* t = A+1 */
353 mpz_mul(y, y, t); /* y = (A+1)*(2A^2+(r-1)A+1) */
354
355 return (mpz_cmp(n, y) < 0) ? 1 : 0;
356 }
bls_theorem7_limit(mpz_t n,mpz_t F1,mpz_t R1,UV B1,mpz_t t,mpz_t y,mpz_t r,mpz_t s)357 static int bls_theorem7_limit(mpz_t n, mpz_t F1, mpz_t R1, UV B1,
358 mpz_t t, mpz_t y, mpz_t r, mpz_t s)
359 {
360 mpz_mul(t, F1, R1);
361 mpz_add_ui(t, t, 1);
362 if (mpz_cmp(t, n) != 0) croak("BLS75 internal error: F1*R1 != n-1\n");
363
364 mpz_mul_ui(t, F1, 2);
365 mpz_tdiv_qr(s, r, R1, t);
366
367 mpz_add(y, t, r);
368 mpz_sub_ui(y, y, B1);
369 mpz_mul(y, y, F1);
370 mpz_add_ui(y, y, 1); /* y = 2F1^2 + (r - B1)F1 + 1 */
371 mpz_mul_ui(t, F1, B1);
372 mpz_add_ui(t, t, 1);
373 mpz_mul(y, y, t); /* times (B1F1+1) */
374
375 return (mpz_cmp(n, y) < 0) ? 1 : 0;
376 }
bls_theorem17_limit(mpz_t n,mpz_t F2,mpz_t R2,UV dummy,mpz_t t,mpz_t y,mpz_t r,mpz_t s)377 static int bls_theorem17_limit(mpz_t n, mpz_t F2, mpz_t R2, UV dummy,
378 mpz_t t, mpz_t y, mpz_t r, mpz_t s)
379 {
380 mpz_mul(t, F2, R2);
381 mpz_sub_ui(t, t, 1);
382 if (mpz_cmp(t, n) != 0) croak("BLS75 internal error: F2*R2 != n+1\n");
383
384 mpz_mul_ui(t, F2, 2);
385 mpz_tdiv_qr(s, r, R2, t);
386 if (mpz_cmp(r, F2) >= 0) {
387 mpz_add_ui(s, s, 1);
388 mpz_sub(r, r, t);
389 }
390 /* Let m = 1 */
391 mpz_add_ui(y, t, 1);
392 mpz_abs(t, r);
393 mpz_sub(y, y, t);
394 mpz_mul(y, y, F2);
395 mpz_add_ui(y, y, 1); /* y = 2F2^2 + (m-r)F2 + 1 */
396 mpz_sub_ui(t, F2, 1);
397 mpz_mul(y, y, t); /* times (mF2-1) */
398
399 return (mpz_cmp(n, y) < 0) ? 1 : 0;
400 }
bls_theorem19_limit(mpz_t n,mpz_t F2,mpz_t R2,UV B2,mpz_t t,mpz_t y,mpz_t r,mpz_t s)401 static int bls_theorem19_limit(mpz_t n, mpz_t F2, mpz_t R2, UV B2,
402 mpz_t t, mpz_t y, mpz_t r, mpz_t s)
403 {
404 mpz_mul(t, F2, R2);
405 mpz_sub_ui(t, t, 1);
406 if (mpz_cmp(t, n) != 0) croak("BLS75 internal error: F2*R2 != n+1\n");
407
408 mpz_mul_ui(t, F2, 2);
409 mpz_tdiv_qr(s, r, R2, t);
410 if (mpz_cmp(r, F2) >= 0) {
411 mpz_add_ui(s, s, 1);
412 mpz_sub(r, r, t);
413 }
414
415 mpz_add_ui(y, t, B2);
416 mpz_abs(t, r);
417 mpz_sub(y, y, t);
418 mpz_mul(y, y, F2);
419 mpz_add_ui(y, y, 1); /* y = 2F2^2 + (B2 - r)F2 + 1 */
420 mpz_mul_ui(t, F2, B2);
421 mpz_sub_ui(t, t, 1);
422 mpz_mul(y, y, t); /* times (B2F2-1) */
423
424 return (mpz_cmp(n, y) < 0) ? 1 : 0;
425 }
426
427 /*
428 17: N < (m F2 - 1) ( 2 F2 F2 + m F2 - |r| F2 + 1 )
429 (III) test (l*F2+1) doesn't divide N for 1 .. m
430
431 19: N < (B2 F2 - 1) ( 2 F2 F2 + B2 F2 - |r| F2 + 1 )
432 (III) (IV) R2 factors > B2
433 */
434
bls_theorem20_limit(mpz_t n,mpz_t R1,mpz_t F1,mpz_t F2,UV B,UV m,mpz_t t,mpz_t g,mpz_t r,mpz_t s)435 static int bls_theorem20_limit(mpz_t n, mpz_t R1, mpz_t F1, mpz_t F2,
436 UV B, UV m,
437 mpz_t t, mpz_t g, mpz_t r, mpz_t s)
438 {
439 mpz_tdiv_q_2exp(t, F2, 1);
440 mpz_tdiv_qr(s, r, R1, t);
441
442 mpz_mul_ui(t, F1, B);
443 mpz_add_ui(g, t, 1);
444
445 mpz_mul_ui(t, F2, B);
446 mpz_sub_ui(t, t, 1);
447
448 if (mpz_cmp(t, g) > 0) mpz_set(g, t);
449
450 mpz_mul(t, F1, F2);
451 mpz_tdiv_q_2exp(t, t, 1);
452 mpz_mul_ui(t, t, B);
453 mpz_mul_ui(t, t, B);
454 mpz_add_ui(s, t, 1); /* s = B1*B2*F1*F2/2+1 */
455
456 mpz_mul(g, g, s);
457 if (mpz_cmp(n, g) < 0) {
458 mpz_set_ui(s, 0); /* Use s to signal whether we must test for m. */
459 return 1;
460 }
461
462 mpz_mul(t, F1, F2);
463 mpz_mul_ui(t, t, m);
464 mpz_tdiv_q_2exp(t, t, 1);
465 mpz_mul(r, r, F1); /* r *= F1; t += r, r /= F1; */
466 mpz_add(t, t, r);
467 mpz_divexact(r, r, F1);
468 mpz_add_ui(t, t, 1); /* t = m * F1 * F2/2 + r * F1 + 1 */
469 mpz_mul(g, s, t);
470 mpz_set_ui(s, 1);
471
472 return (mpz_cmp(n, g) < 0) ? 1 : 0;
473 }
474
475
476 /******************************************************************************/
477
478
479 /* (I) For each prime p_i dividing F1 [N-1 = F1R1] there exists an a_i
480 * such that N is a psp base a_i and gcd(A_i^{(N-1)/p_i}-1,N) = 1.
481 */
_verify_cond_I_p(mpz_t n,mpz_t pi,mpz_t ap,mpz_t t,int alimit,char * pspcache)482 static int _verify_cond_I_p(mpz_t n, mpz_t pi, mpz_t ap, mpz_t t, int alimit, char* pspcache)
483 {
484 int a, success = 0;
485 PRIME_ITERATOR(iter);
486
487 for (a = 2; !success && a <= alimit; a = prime_iterator_next(&iter)) {
488 int psp = -1;
489 mpz_set_ui(ap, a);
490
491 if (pspcache) psp = pspcache[a];
492 if (psp == -1) {
493 mpz_sub_ui(t, n, 1);
494 mpz_powm(t, ap, t, n);
495 psp = (mpz_cmp_ui(t, 1) == 0);
496 }
497 if (pspcache) pspcache[a] = psp;
498 if (!psp)
499 continue;
500
501 mpz_sub_ui(t, n, 1);
502 mpz_divexact(t, t, pi);
503 mpz_powm(t, ap, t, n);
504 mpz_sub_ui(t, t, 1);
505 mpz_gcd(t, t, n);
506 if (mpz_cmp_ui(t, 1) != 0)
507 continue;
508
509 success = 1; /* We found an a for this p */
510 }
511 prime_iterator_destroy(&iter);
512 return success;
513 }
514
515 /* (III) For each prime q_i dividing F2 [N+1 = F2R2] there exists a Lucas
516 * sequence U_k with discriminant D for which D/N = -1, N divides
517 * U_{N+1} and gcd(U_{(N+1)/q_i},N} = 1.
518 */
_verify_cond_III_q2(mpz_t n,mpz_t qi,IV p,IV q,mpz_t U,mpz_t V,mpz_t k,mpz_t t1,mpz_t t2)519 static int _verify_cond_III_q2(mpz_t n, mpz_t qi, IV p, IV q, mpz_t U, mpz_t V, mpz_t k, mpz_t t1, mpz_t t2)
520 {
521 mpz_add_ui(k, n, 1);
522 mpz_divexact(k, k, qi);
523 lucas_seq(U, V, n, p, q, k, t1, t2);
524 mpz_gcd(k, U, n);
525 return (mpz_cmp_ui(k, 1) == 0) ? 1 : 0;
526 }
527
528 #define MAXQV 50
529
_verify_cond_III_q(mpz_t n,mpz_t qi,IV * qv,int * pnumqv,IV * lp,IV * lq)530 static int _verify_cond_III_q(mpz_t n, mpz_t qi, IV* qv, int* pnumqv, IV* lp, IV* lq)
531 {
532 int i, numqv = *pnumqv, rval = 0;
533 IV d, p, q, startq;
534 mpz_t U, V, k, t1, t2;
535
536 mpz_init(U); mpz_init(V); mpz_init(k); mpz_init(t1); mpz_init(t2);
537
538 /* Try previous q values with (D|n)=-1 and U_(n+1) = 0 mod n */
539 for (i = 0; i < numqv; i++) {
540 q = qv[i];
541 p = (q % 2) ? 2 : 1;
542 if (_verify_cond_III_q2(n, qi, p, q, U, V, k, t1, t2)) {
543 rval = 2;
544 break;
545 }
546 }
547
548 if (rval == 0) {
549 /* Search for a q value */
550 startq = (numqv > 0) ? qv[numqv-1]+1 : 2;
551 for (q = startq; q < startq+1000; q++) {
552 if (mpz_cmp_ui(n, (unsigned long)q) <= 0) break;
553 p = (q % 2) ? 2 : 1;
554 d = p*p - 4*q;
555 mpz_set_si(t1, d);
556 if (mpz_jacobi(t1, n) != -1)
557 continue;
558 /* we have a d/p/q where d = -1. Check the first Lucas sequence. */
559 mpz_add_ui(k, n, 1);
560 lucas_seq(U, V, n, p, q, k, t1, t2);
561 if (mpz_sgn(U) != 0)
562 continue;
563 /* Passed first test, add to qv list */
564 if (numqv < MAXQV) {
565 qv[numqv] = q;
566 *pnumqv = ++numqv;
567 }
568 /* Verify second Lucas sequence */
569 if (_verify_cond_III_q2(n, qi, p, q, U, V, k, t1, t2)) {
570 rval = 2;
571 break;
572 }
573 }
574 }
575 mpz_clear(U); mpz_clear(V); mpz_clear(k); mpz_clear(t1); mpz_clear(t2);
576 if (lp) *lp = p;
577 if (lq) *lq = q;
578 return rval;
579 }
580 #if 0
581 /* N divides V_{(N+1)/2} and for q > 2 does not divide V{(N+1)/(2q)} */
582 static int _verify_theorem14_q(mpz_t n, mpz_t qi, IV* lastq, IV* lp, IV* lq)
583 {
584 int rval = 0;
585 IV d, p, q;
586 mpz_t U, V, k, t1, t2;
587
588 mpz_init(U); mpz_init(V); mpz_init(k); mpz_init(t1); mpz_init(t2);
589
590 if (lastq && *lastq > 0 && mpz_cmp_ui(qi,2) > 0) {
591 q = *lastq;
592 p = (q % 2) ? 2 : 1;
593 d = p*p - 4*q;
594 mpz_set_si(t1, d);
595 if (mpz_jacobi(t1, n) == -1) {
596 /* q passed the first tst. Do the second that depends on the factor. */
597 mpz_add_ui(k, n, 1);
598 mpz_tdiv_q_2exp(k, k, 1);
599 mpz_divexact(k, k, qi);
600 lucas_seq(U, V, n, p, q, k, t1, t2);
601 if (mpz_sgn(V) != 0) {
602 rval = 2;
603 }
604 }
605 }
606
607 if (!rval) {
608 for (q = (lastq && *lastq > 0) ? *lastq + 1 : 2; q < 1000; q++) {
609 p = (q % 2) ? 2 : 1;
610 d = p*p - 4*q;
611 mpz_set_si(t1, d);
612 if (mpz_jacobi(t1, n) != -1)
613 continue;
614
615 /* we have a d/p/q where d = -1. Check the Lucas sequence. */
616
617 mpz_add_ui(k, n, 1);
618 mpz_tdiv_q_2exp(k, k, 1);
619 lucas_seq(U, V, n, p, q, k, t1, t2);
620 if (mpz_sgn(V) == 0) {
621 if (mpz_cmp_ui(qi, 2) <= 0) {
622 rval = 2; break;
623 } else {
624 mpz_divexact(k, k, qi);
625 lucas_seq(U, V, n, p, q, k, t1, t2);
626 if (mpz_sgn(V) != 0) {
627 rval = 2; break;
628 }
629 }
630 }
631 }
632 }
633
634 mpz_clear(U); mpz_clear(V); mpz_clear(k); mpz_clear(t1); mpz_clear(t2);
635 if (lp) *lp = p;
636 if (lq) *lq = q;
637 if (lastq) *lastq = q;
638 return rval;
639 }
640 #endif
641
642
643
644 /******************************************************************************/
645
646
647
_GMP_primality_bls_nm1(mpz_t n,int effort,char ** prooftextptr)648 int _GMP_primality_bls_nm1(mpz_t n, int effort, char** prooftextptr)
649 {
650 mpz_t nm1, A, B, t, m, f, r, s;
651 FACTOR_STACK(fstack);
652 FACTOR_STACK(mstack);
653 int e, success = 1;
654 UV B1 = (mpz_sizeinbase(n,10) > 1000) ? 100000 : 10000;
655
656 /* We need to do this for BLS */
657 if (mpz_even_p(n)) return 0;
658
659 mpz_init(nm1);
660 mpz_sub_ui(nm1, n, 1);
661 mpz_init_set_ui(A, 1);
662 mpz_init_set(B, nm1);
663 mpz_init(m);
664 mpz_init(f);
665 mpz_init(t);
666 mpz_init(r);
667 mpz_init(s);
668
669 { /* Pull small factors out */
670 PRIME_ITERATOR(iter);
671 UV tf;
672 for (tf = 2; tf < B1; tf = prime_iterator_next(&iter)) {
673 if (mpz_cmp_ui(B, tf*tf) < 0) break;
674 factor_test_ui(tf, B, A, &fstack);
675 }
676 prime_iterator_destroy(&iter);
677 }
678
679 if (success && mpz_cmp_ui(B,1) > 0) {
680 mpz_set(f, B);
681 handle_factor(f, B, A, &fstack, &mstack, effort, prooftextptr, 1, &_GMP_primality_bls_nm1);
682 }
683
684 while (success) {
685
686 if ( (prooftextptr && bls_theorem5_limit(n, A, B, B1, t, m, r, s)) || (!prooftextptr && bls_theorem7_limit(n, A, B, B1, t, m, r, s)) )
687 break;
688
689 success = 0;
690 if (nstack(&mstack) == 0) /* If the stack is empty, we have failed. */
691 break;
692 pop_fstack(m, &mstack); /* pop a component off the stack */
693
694 for (e = 0; !success && e <= effort; e++)
695 success = tfe(f, m, e);
696
697 /* If we couldn't factor m and the stack is empty, we've failed. */
698 if ( (!success) && (nstack(&mstack) == 0) )
699 break;
700 /* Put the two factors f and m/f into the stacks, smallest first */
701 mpz_divexact(m, m, f);
702 if (mpz_cmp(m, f) < 0)
703 mpz_swap(m, f);
704 handle_factor(f, B, A, &fstack, &mstack, effort, prooftextptr, 0, &_GMP_primality_bls_nm1);
705 handle_factor(m, B, A, &fstack, &mstack, effort, prooftextptr, 0, &_GMP_primality_bls_nm1);
706 }
707
708 /* clear mstack since we don't care about it. Use to hold a values. */
709 clear_fstack(&mstack);
710
711 fstack_sort_trim(&fstack);
712
713 /* Shrink to smallest set and verify conditions. */
714 if (success > 0) {
715 if (prooftextptr)
716 trim_factors(A, B, n, nm1, B1, &fstack, &bls_theorem5_limit, t, m, r, s);
717 else
718 trim_factors(A, B, n, nm1, B1, &fstack, &bls_theorem7_limit, t, m, r, s);
719 /* Verify conditions */
720 success = 0;
721 if ( (prooftextptr && bls_theorem5_limit(n, A, B, B1, t, m, r, s)) || (!prooftextptr && bls_theorem7_limit(n, A, B, B1, t, m, r, s)) ) {
722 mpz_mul(t, r, r);
723 mpz_submul_ui(t, s, 8); /* t = r^2 - 8s */
724 /* N is prime if and only if s=0 OR t not a perfect square */
725 success = (mpz_sgn(s) == 0 || !mpz_perfect_square_p(t)) ? 1 : -1;
726 }
727 }
728
729 if (success > 0) {
730 int pcount, a;
731 int const alimit = (effort <= 2) ? 200 : 10000;
732 char afermat[10000+1];
733 mpz_t p, ap;
734
735 mpz_init(p);
736 mpz_init(ap);
737
738 /* Cache result that doesn't depend on factor */
739 for (a = 0; a <= alimit; a++) afermat[a] = -1;
740
741 for (pcount = 0; success && pcount < fstack.cur; pcount++) {
742 success = _verify_cond_I_p(n, fstack.stack[pcount], ap, t, alimit, afermat);
743 if (success)
744 push_fstack(&mstack, ap);
745 }
746 /* If we could not find 'a' values, then we should return 1 (maybe prime)
747 * since we did not perform an exhaustive search. It would be quite
748 * unusual to find a prime that didn't have an 'a' in the first 10,000
749 * primes, but it could happen. It's a "dubiously prime" :) */
750 if (!success && get_verbose_level() > 0)
751 printf("N-1 factored but failed to prove. Perhaps composite.\n");
752 mpz_clear(p);
753 mpz_clear(ap);
754 }
755
756 if (success > 0 && prooftextptr != 0) {
757 int i;
758 char *proofstr, *proofptr;
759 int curprooflen = (*prooftextptr == 0) ? 0 : strlen(*prooftextptr);
760 int fsp = nstack(&fstack);
761 int msp = nstack(&mstack);
762 int myprooflen = (5 + mpz_sizeinbase(n, 10)) * (2 + fsp + msp) + 200;
763
764 if (fsp != msp) croak("Different f and a counts\n");
765 New(0, proofstr, myprooflen + curprooflen + 1, char);
766 proofptr = proofstr;
767 proofptr += gmp_sprintf(proofptr, "Type BLS5\nN %Zd\n", n);
768 /* Q[0] is always 2 */
769 for (i = 1; i < fsp; i++)
770 proofptr += gmp_sprintf(proofptr, "Q[%d] %Zd\n", i, fstack.stack[i]);
771 /* A[i] only printed if not 2 */
772 for (i = 0; i < msp; i++)
773 if (mpz_cmp_ui(mstack.stack[i], 2) != 0)
774 proofptr += gmp_sprintf(proofptr, "A[%d] %Zd\n", i, mstack.stack[i]);
775 proofptr += gmp_sprintf(proofptr, "----\n");
776 /* Set or prepend */
777 if (*prooftextptr) {
778 proofptr += gmp_sprintf(proofptr, "\n");
779 strcat(proofptr, *prooftextptr);
780 Safefree(*prooftextptr);
781 }
782 *prooftextptr = proofstr;
783 }
784
785 destroy_fstack(&fstack);
786 destroy_fstack(&mstack);
787 mpz_clear(nm1);
788 mpz_clear(A);
789 mpz_clear(B);
790 mpz_clear(m);
791 mpz_clear(f);
792 mpz_clear(t);
793 mpz_clear(r);
794 mpz_clear(s);
795 if (success < 0) return 0;
796 if (success > 0) return 2;
797 return 1;
798 }
799
800
_GMP_primality_bls_np1(mpz_t n,int effort,char ** prooftextptr)801 int _GMP_primality_bls_np1(mpz_t n, int effort, char** prooftextptr)
802 {
803 mpz_t np1, F2, R2, t, m, f, r, s;
804 FACTOR_STACK(fstack);
805 FACTOR_STACK(mstack);
806 int e, success = 1;
807 UV B2 = (mpz_sizeinbase(n,10) > 1000) ? 100000 : 10000;
808
809 /* TODO: T19 doesn't seem right, and we're still
810 * passing some composites like 14299. */
811
812 /* We need to do this for BLS */
813 if (mpz_even_p(n)) return 0;
814
815 mpz_init(np1);
816 mpz_add_ui(np1, n, 1);
817 mpz_init_set_ui(F2, 1);
818 mpz_init_set(R2, np1);
819 mpz_init(m);
820 mpz_init(f);
821 mpz_init(t);
822 mpz_init(r);
823 mpz_init(s);
824
825 { /* Pull small factors out */
826 PRIME_ITERATOR(iter);
827 UV tf;
828 for (tf = 2; tf < B2; tf = prime_iterator_next(&iter)) {
829 if (mpz_cmp_ui(R2, tf*tf) < 0) break;
830 factor_test_ui(tf, R2, F2, &fstack);
831 }
832 prime_iterator_destroy(&iter);
833 }
834
835 /* printf("trial np1: %lu bits, %lu bits factored\n", mpz_sizeinbase(np1,2), mpz_sizeinbase(F2,2)); */
836
837 if (success && mpz_cmp_ui(R2,1) > 0) {
838 mpz_set(f, R2);
839 handle_factor(f, R2, F2, &fstack, &mstack, effort, prooftextptr, 1, &_GMP_primality_bls_np1);
840 }
841
842 while (success) {
843
844 if (bls_theorem17_limit(n, F2, R2, B2, t, m, r, s))
845 break;
846
847 success = 0;
848
849 if (nstack(&mstack) == 0) /* If the stack is empty, we have failed. */
850 break;
851 pop_fstack(m, &mstack); /* pop a component off the stack */
852
853 for (e = 0; !success && e <= effort; e++)
854 success = tfe(f, m, e);
855
856 /* If we couldn't factor m and the stack is empty, we've failed. */
857 if ( (!success) && (nstack(&mstack) == 0) )
858 break;
859 /* Put the two factors f and m/f into the stacks, smallest first */
860 mpz_divexact(m, m, f);
861 if (mpz_cmp(m, f) < 0)
862 mpz_swap(m, f);
863 handle_factor(f, R2, F2, &fstack, &mstack, effort, prooftextptr, 0, &_GMP_primality_bls_np1);
864 handle_factor(m, R2, F2, &fstack, &mstack, effort, prooftextptr, 0, &_GMP_primality_bls_np1);
865 }
866
867 /* clear mstack since we don't care about it. Use to hold a values. */
868 clear_fstack(&mstack);
869
870 /* printf("factor np1: %lu bits, %lu bits factored\n", mpz_sizeinbase(np1,2), mpz_sizeinbase(F2,2)); */
871
872 fstack_sort_trim(&fstack);
873
874 /* Shrink to smallest set and verify conditions. */
875 if (success > 0) {
876 trim_factors(F2, R2, n, np1, B2, &fstack, &bls_theorem17_limit, t, m, r, s);
877 /* Verify conditions */
878 success = 0;
879 if (bls_theorem17_limit(n, F2, R2, B2, t, m, r, s)) {
880 mpz_mul(t, r, r);
881 mpz_addmul_ui(t, s, 8); /* t = r^2 + 8s */
882 /* N is prime if and only if s=0 OR t not a perfect square */
883 success = (mpz_sgn(s) == 0 || !mpz_perfect_square_p(t)) ? 1 : -1;
884 }
885 }
886
887 if (success > 0) {
888 IV qv[MAXQV];
889 int pcount, numqv = 0;
890 for (pcount = 0; success && pcount < fstack.cur; pcount++) {
891 success = _verify_cond_III_q(n, fstack.stack[pcount], qv, &numqv, 0, 0);
892 }
893 /* If we meet theorem 17 limits, then no need to test R2 */
894 if (success && !bls_theorem17_limit(n, F2, R2, B2, t, m, r, s)) {
895 success = _verify_cond_III_q(n, R2, qv, &numqv, 0, 0);
896 }
897 if (!success && get_verbose_level() > 0)
898 printf("N+1 factored but failed to prove. Perhaps composite.\n");
899 }
900
901 /* TODO: Proof text */
902
903 destroy_fstack(&fstack);
904 destroy_fstack(&mstack);
905 mpz_clear(np1);
906 mpz_clear(F2);
907 mpz_clear(R2);
908 mpz_clear(m);
909 mpz_clear(f);
910 mpz_clear(t);
911 mpz_clear(r);
912 mpz_clear(s);
913 if (success < 0) return 0;
914 if (success > 0) return 2;
915 return 1;
916 }
917
918 /******************************************************************************/
919
920 #define PRINT_PCT 0
921
922 /* Will use one of:
923 * N-1 Corollary 1
924 * N-1 Theorem 5
925 * N-1 Theorem 7
926 * N+1 Corollary 8
927 * N+1 Theorem 17
928 * N+1 Theorem 19
929 * Comb Theorem 20
930 */
bls75_hybrid(mpz_t n,int effort,char ** prooftextptr)931 int bls75_hybrid(mpz_t n, int effort, char** prooftextptr)
932 {
933 mpz_t nm1, np1, F1, F2, R1, R2;
934 mpz_t r, s, t, u, f, c1, c2;
935 /* fstack: definite prime factors
936 * pstack: probable prime factors product of fstack and pstack = F
937 * mstack: composite remainders product of mstack = R
938 */
939 FACTOR_STACK(f1stack);
940 FACTOR_STACK(f2stack);
941 FACTOR_STACK(p1stack);
942 FACTOR_STACK(p2stack);
943 FACTOR_STACK(m1stack);
944 FACTOR_STACK(m2stack);
945 int pcount, e, success = 1;
946 int low_effort = (effort < 1) ? 0 : 1;
947 UV B1 = (effort < 2 && mpz_sizeinbase(n,2) < 160) ? 6000 :
948 (mpz_sizeinbase(n,2) < 1024) ? 20000 : 200000;
949 UV m = B1-1; /* m should be less than B1 */
950 #if PRINT_PCT
951 double trial_pct, prime_pct, fac_pct, fin_pct;
952 #endif
953
954 /* We need to do this for BLS */
955 if (mpz_even_p(n)) return 0;
956
957 mpz_init(nm1); mpz_sub_ui(nm1, n, 1);
958 mpz_init(np1); mpz_add_ui(np1, n, 1);
959
960 mpz_init_set_ui(F1, 1); mpz_init_set(R1, nm1);
961 mpz_init_set_ui(F2, 1); mpz_init_set(R2, np1);
962
963 mpz_init(r);
964 mpz_init(s);
965 mpz_init(u);
966 mpz_init(t);
967 mpz_init(f); mpz_init(c1); mpz_init(c2);
968
969 { /* Pull small factors out */
970 PRIME_ITERATOR(iter);
971 UV tf;
972 for (tf = 2; tf < B1; tf = prime_iterator_next(&iter)) {
973 /* Page 635 of BLS75 describes an optimization for divisibility
974 * testing. It seems slower than just doing two UI div tests. */
975 if (mpz_cmp_ui(R1, tf*tf) >= 0)
976 factor_test_ui(tf, R1, F1, &f1stack);
977 if (mpz_cmp_ui(R2, tf*tf) >= 0)
978 factor_test_ui(tf, R2, F2, &f2stack);
979 }
980 prime_iterator_destroy(&iter);
981 }
982
983 #if PRINT_PCT
984 trial_pct = (100.0 * (mpz_sizeinbase(F1,2) + mpz_sizeinbase(F2,2))) / (mpz_sizeinbase(nm1,2) + mpz_sizeinbase(np1,2));
985 printf("\n%6.2f .. ", trial_pct); fflush(stdout);
986 #endif
987
988 if ( bls_theorem7_limit(n, F1, R1, B1, t, u, r, s) ||
989 bls_theorem19_limit(n, F2, R2, B1, t, u, r, s) ||
990 bls_theorem20_limit(n, R1, F1, F2, B1, m, t, u, r, s) )
991 goto start_hybrid_proof;
992
993 if (mpz_cmp_ui(R1,1) > 0) {
994 mpz_set(f, R1);
995 handle_factor2(f, R1, F1, &f1stack, &p1stack, &m1stack, low_effort, prooftextptr, &bls75_hybrid);
996 }
997 if (mpz_cmp_ui(R2,1) > 0) {
998 mpz_set(f, R2);
999 handle_factor2(f, R2, F2, &f2stack, &p2stack, &m2stack, low_effort, prooftextptr, &bls75_hybrid);
1000 }
1001
1002 #if PRINT_PCT
1003 prime_pct = (100.0 * (mpz_sizeinbase(F1,2) + mpz_sizeinbase(F2,2))) / (mpz_sizeinbase(nm1,2) + mpz_sizeinbase(np1,2));
1004 printf("\n%6.2f .. ", prime_pct); fflush(stdout);
1005 #endif
1006
1007 while (1) {
1008 int d1, d2;
1009
1010 success = 1;
1011 if ( bls_theorem7_limit(n, F1, R1, B1, t, u, r, s) ||
1012 bls_theorem19_limit(n, F2, R2, B1, t, u, r, s) ||
1013 bls_theorem20_limit(n, R1, F1, F2, B1, m, t, u, r, s) )
1014 break;
1015
1016 success = 0;
1017 mpz_set_ui(c1, 0);
1018 mpz_set_ui(c2, 0);
1019 d1 = nstack(&m1stack) > 0;
1020 d2 = nstack(&m2stack) > 0;
1021 if (!d1 && !d2) break;
1022
1023 if (d1) pop_fstack(c1, &m1stack);
1024 if (d2) pop_fstack(c2, &m2stack);
1025 for (e = 0; !success && e <= effort; e++) {
1026 if (d1 && tfe(f, c1, e)) {
1027 if (d2) push_fstack(&m2stack, c2);
1028 mpz_set(u, c1);
1029 success = 1;
1030 } else if (d2 && tfe(f, c2, e)) {
1031 if (d1) push_fstack(&m1stack, c1);
1032 mpz_set(u, c2);
1033 success = 2;
1034 }
1035 }
1036
1037 /* No success for this set of composites. Move on. */
1038 if (!success) continue;
1039
1040 if (success == 1) {
1041 mpz_divexact(u, u, f);
1042 if (mpz_cmp(u, f) < 0)
1043 mpz_swap(u, f);
1044 handle_factor2(f, R1, F1, &f1stack, &p1stack, &m1stack, low_effort, prooftextptr, &bls75_hybrid);
1045 handle_factor2(u, R1, F1, &f1stack, &p1stack, &m1stack, low_effort, prooftextptr, &bls75_hybrid);
1046 } else if (success == 2) {
1047 mpz_divexact(u, u, f);
1048 if (mpz_cmp(u, f) < 0)
1049 mpz_swap(u, f);
1050 handle_factor2(f, R2, F2, &f2stack, &p2stack, &m2stack, low_effort, prooftextptr, &bls75_hybrid);
1051 handle_factor2(u, R2, F2, &f2stack, &p2stack, &m2stack, low_effort, prooftextptr, &bls75_hybrid);
1052 }
1053 #if PRINT_PCT
1054 fac_pct = (100.0 * (mpz_sizeinbase(F1,2) + mpz_sizeinbase(F2,2))) / (mpz_sizeinbase(nm1,2) + mpz_sizeinbase(np1,2));
1055 printf("%6.2f .. ", fac_pct); fflush(stdout);
1056 #endif
1057 }
1058
1059 start_hybrid_proof:
1060 mpz_mul(t, F1, R1); if (mpz_cmp(nm1, t) != 0) croak("Bad n-1 factor");
1061 mpz_mul(t, F2, R2); if (mpz_cmp(np1, t) != 0) croak("Bad n+1 factor");
1062
1063 /* We've done all the factoring we need or can. */
1064
1065 /* Finish proofs for p{1,2}stack as needed. */
1066 /* TODO: optimize for cases of both n-1 and n+1 working */
1067 if (nstack(&p1stack) > 0) {
1068 while (nstack(&p1stack) > 0) {
1069 int pr = 1;
1070 pop_fstack(f, &p1stack);
1071 if (effort > low_effort)
1072 pr = bls75_hybrid(f, effort, prooftextptr);
1073 if (pr == 0) croak("probable prime factor proved composite");
1074 else if (pr == 2) push_fstack(&f1stack, f); /* Proved, put on F stack */
1075 else factor_out(F1, R1, f); /* No proof. Move to R */
1076 }
1077 }
1078 if (nstack(&p2stack) > 0) {
1079 while (nstack(&p2stack) > 0) {
1080 int pr = 1;
1081 pop_fstack(f, &p2stack);
1082 if (effort > low_effort)
1083 pr = bls75_hybrid(f, effort, prooftextptr);
1084 if (pr == 0) croak("probable prime factor proved composite");
1085 else if (pr == 2) push_fstack(&f2stack, f); /* Proved, put on F stack */
1086 else factor_out(F2, R2, f); /* No proof. Move to R */
1087 }
1088 }
1089
1090 fstack_sort_trim(&f1stack);
1091 fstack_sort_trim(&f2stack);
1092
1093 #if PRINT_PCT
1094 fin_pct = (100.0 * (mpz_sizeinbase(F1,2) + mpz_sizeinbase(F2,2))) / (mpz_sizeinbase(nm1,2) + mpz_sizeinbase(np1,2));
1095 printf("%6.2f .. ", fin_pct); fflush(stdout);
1096 printf("\n"); fflush(stdout);
1097 #endif
1098
1099 /* Check the theorems we have available */
1100
1101 if (bls_theorem7_limit(n, F1, R1, B1, t, u, r, s)) {
1102 if (get_verbose_level() > 0) printf("BLS75 proof using N-1\n");
1103 trim_factors(F1, R1, n, nm1, B1, &f1stack, &bls_theorem7_limit, t, u, r, s);
1104 for (pcount = 0; success > 0 && pcount < f1stack.cur; pcount++)
1105 success = _verify_cond_I_p(n, f1stack.stack[pcount], u, t, 1000, 0);
1106 if (success > 0 && (mpz_mul(t, F1, F1), mpz_cmp(t,n) > 0))
1107 goto end_hybrid; /* Corollary 1, n-1 factored more than sqrt(n) */
1108 if (success > 0 && !bls_theorem5_limit(n, F1, R1, B1, t, u, r, s))
1109 success = _verify_cond_I_p(n, R1, u, t, 1000, 0);
1110 if (success > 0) {
1111 mpz_mul(t, r, r);
1112 mpz_submul_ui(t, s, 8); /* t = r^2 - 8s */
1113 /* N is prime if and only if s=0 OR t not a perfect square */
1114 success = (mpz_sgn(s) == 0 || !mpz_perfect_square_p(t)) ? 1 : -1;
1115 }
1116 goto end_hybrid; /* Theorem 5 or 7 */
1117 }
1118
1119
1120 if (bls_theorem19_limit(n, F2, R2, B1, t, u, r, s)) {
1121 IV qv[MAXQV];
1122 int numqv = 0;
1123 if (get_verbose_level() > 0) printf("BLS75 proof using N+1\n");
1124 trim_factors(F2, R2, n, np1, B1, &f2stack, &bls_theorem19_limit, t, u, r, s);
1125 for (pcount = 0; success > 0 && pcount < f2stack.cur; pcount++)
1126 success = _verify_cond_III_q(n, f2stack.stack[pcount], qv, &numqv, 0, 0);
1127 if (success > 0 && (mpz_mul(t,F2,F2), mpz_add_ui(t,t,1), mpz_cmp(t,n) > 0))
1128 goto end_hybrid; /* Corollary 8, n+1 factored more than sqrt(n)+1 */
1129 if (success > 0 && !bls_theorem17_limit(n, F2, R2, B1, t, u, r, s))
1130 success = _verify_cond_III_q(n, R2, qv, &numqv, 0, 0);
1131 if (success > 0) {
1132 mpz_mul(t, r, r);
1133 mpz_submul_ui(t, s, 8); /* t = r^2 - 8s */
1134 /* N is prime if and only if s=0 OR t not a perfect square */
1135 success = (mpz_sgn(s) == 0 || !mpz_perfect_square_p(t)) ? 1 : -1;
1136 }
1137 goto end_hybrid; /* Theorem 17 or 19 */
1138 }
1139
1140 if (get_verbose_level() > 0) printf("BLS75 proof using N-1 / N+1 (T20)\n");
1141
1142 /* Check N < B^3 F1*F2*F2/2 or N < B^3 F1*F1*F2/2 */
1143 success = bls_theorem20_limit(n, R1, F1, F2, B1, m, t, u, r, s);
1144
1145 /* Trim some factors from f2stack if possible */
1146 if (nstack(&f2stack) > 1) {
1147 int i;
1148 mpz_set_ui(F2, 1);
1149 mpz_set(R2, np1);
1150 for (i = 0; i < f2stack.cur; i++) {
1151 if (i > 0 && bls_theorem20_limit(n, R1, F1, F2, B1, m, t, u, r, s))
1152 break;
1153 factor_out(R2, F2, f2stack.stack[i]);
1154 }
1155 /* Remove excess factors */
1156 while (i < f2stack.cur)
1157 pop_fstack(t, &f2stack);
1158 /* Verify Q[0] = 2 */
1159 if (mpz_cmp_ui(f2stack.stack[0], 2) != 0)
1160 croak("BLS75 internal error: 2 not at start of fstack");
1161 }
1162
1163 /* Check lambda divisibility if needed */
1164 if (success > 0 && mpz_sgn(s)) {
1165 UV lambda;
1166 mpz_mul(t, F1, F2);
1167 mpz_tdiv_q_2exp(t, t, 1);
1168 mpz_mul(u, r, F1);
1169 mpz_add_ui(u, u, 1);
1170 for (lambda = 0; success > 0 && lambda < m; lambda++, mpz_add(u,u,t)) {
1171 if (lambda > 0 || mpz_sgn(r))
1172 if (mpz_divisible_p(n, u)) {
1173 /* Check that we found a non-trivial divisor */
1174 mpz_gcd(t, u, n);
1175 success = (mpz_cmp_ui(t,1) > 0 || mpz_cmp(t,n) < 0) ? -1 : 0;
1176 break;
1177 }
1178 }
1179 }
1180
1181 /* Verify (I) page 623 and (II) page 625 */
1182 if (success > 0) {
1183 for (pcount = 0; success > 0 && pcount < f1stack.cur; pcount++)
1184 success = _verify_cond_I_p(n, f1stack.stack[pcount], u, t, 1000, 0);
1185 if (success > 0)
1186 success = _verify_cond_I_p(n, R1, u, t, 1000, 0);
1187 }
1188
1189 /* Verify (III) page 631 and (IV) page 633 */
1190 if (success > 0) {
1191 IV qv[MAXQV];
1192 int numqv = 0;
1193 for (pcount = 0; success > 0 && pcount < f2stack.cur; pcount++)
1194 success = _verify_cond_III_q(n, f2stack.stack[pcount], qv, &numqv, 0, 0);
1195 if (success > 0)
1196 success = _verify_cond_III_q(n, R2, qv, &numqv, 0, 0);
1197 }
1198
1199 #if 0
1200 { double p1 = (100.0 * mpz_sizeinbase(F1,2) / mpz_sizeinbase(nm1,2));
1201 double p2 = (100.0 * mpz_sizeinbase(F2,2) / mpz_sizeinbase(np1,2));
1202 printf("%6.2f %6.2f\n", p1, p2); fflush(stdout); }
1203 //{ double pct = (100.0 * (mpz_sizeinbase(R1,2) + mpz_sizeinbase(R2,2))) / (mpz_sizeinbase(nm1,2) + mpz_sizeinbase(np1,2)); printf("%6.2f\n", 100.0-pct); fflush(stdout); }
1204 #endif
1205
1206 end_hybrid:
1207 destroy_fstack(&f1stack);
1208 destroy_fstack(&f2stack);
1209 destroy_fstack(&p1stack);
1210 destroy_fstack(&p2stack);
1211 destroy_fstack(&m1stack);
1212 destroy_fstack(&m2stack);
1213 mpz_clear(nm1); mpz_clear(np1);
1214 mpz_clear(F1); mpz_clear(F2);
1215 mpz_clear(R1); mpz_clear(R2);
1216 mpz_clear(r);
1217 mpz_clear(s);
1218 mpz_clear(u);
1219 mpz_clear(t);
1220 mpz_clear(f); mpz_clear(c1); mpz_clear(c2);
1221 if (success < 0) return 0;
1222 if (success > 0) return 2;
1223 return 1;
1224 }
1225
1226
1227 /* Given an n where we're factored n-1 down to p, check BLS theorem 3 */
_GMP_primality_bls_3(mpz_t n,mpz_t p,UV * reta)1228 int _GMP_primality_bls_3(mpz_t n, mpz_t p, UV* reta)
1229 {
1230 mpz_t nm1, m, t, t2;
1231 int rval = 0;
1232
1233 if (reta) *reta = 0;
1234 if (mpz_cmp_ui(n, 2) <= 0 || mpz_even_p(n) || mpz_even_p(p))
1235 return 0; /* n is <= 2, n is even, or p is even */
1236 if (!_GMP_is_prob_prime(p))
1237 return 0; /* p is not a probable prime */
1238
1239 mpz_init(nm1); mpz_init(m); mpz_init(t); mpz_init(t2);
1240 mpz_sub_ui(nm1, n, 1);
1241 mpz_divexact(m, nm1, p);
1242 mpz_mul(t, m, p);
1243 if (mpz_cmp(nm1, t) != 0)
1244 goto end_bls3; /* m*p != n+1 */
1245
1246 mpz_mul_ui(t, p, 2);
1247 mpz_add_ui(t, t, 1);
1248 mpz_sqrt(t2, n);
1249 if (mpz_cmp(t, t2) <= 0)
1250 goto end_bls3; /* 2p+1 <= sqrt(n) */
1251
1252 {
1253 /* N-1 = mp, p is an odd probable prime, and 2p+1 > sqrt(n).
1254 * Now find an 'a' where a^(n-1)/2 = -1 mod n, a^(m/2) != -1 mod n. */
1255 PRIME_ITERATOR(iter);
1256 UV const alimit = 1000;
1257 UV a;
1258 for (a = 2; a <= alimit; a = prime_iterator_next(&iter)) {
1259 /* should check kronecker(a,n) == -1 here */
1260 mpz_set_ui(t2, a);
1261 mpz_divexact_ui(t, m, 2);
1262 mpz_powm(t, t2, t, n); /* a^(m/2) mod n */
1263 if (mpz_cmp(t, nm1) == 0)
1264 continue;
1265 mpz_divexact_ui(t, nm1, 2);
1266 mpz_powm(t, t2, t, n); /* a^((n-1)/2) mod n */
1267 if (mpz_cmp(t, nm1) != 0)
1268 continue;
1269 rval = 2;
1270 if (reta) *reta = a;
1271 break;
1272 }
1273 prime_iterator_destroy(&iter);
1274 }
1275
1276 end_bls3:
1277 mpz_clear(nm1); mpz_clear(m); mpz_clear(t); mpz_clear(t2);
1278 return rval;
1279 }
1280
1281 /* Given an n where we're factored n+1 down to f, check BLS theorem 15 */
_GMP_primality_bls_15(mpz_t n,mpz_t f,IV * lp,IV * lq)1282 int _GMP_primality_bls_15(mpz_t n, mpz_t f, IV* lp, IV* lq)
1283 {
1284 mpz_t np1, m, t, t2;
1285 int rval = 0;
1286
1287 if (lp) *lp = 0;
1288 if (lq) *lq = 0;
1289 if (mpz_cmp_ui(n, 2) <= 0 || mpz_even_p(n) || mpz_even_p(f))
1290 return 0; /* n is <= 2, n is even, or f is even */
1291 if (!_GMP_is_prob_prime(f))
1292 return 0; /* f is not a probable prime */
1293
1294 mpz_init(np1); mpz_init(m); mpz_init(t); mpz_init(t2);
1295 mpz_add_ui(np1, n, 1);
1296 mpz_divexact(m, np1, f);
1297 mpz_mul(t, m, f);
1298 if (mpz_cmp(np1, t) != 0)
1299 goto end_bls15; /* m*f != n+1 */
1300
1301 mpz_mul_ui(t, f, 2);
1302 mpz_sub_ui(t, t, 1);
1303 mpz_sqrt(t2, n);
1304 if (mpz_cmp(t, t2) <= 0)
1305 goto end_bls15; /* 2f-1 <= sqrt(n) */
1306
1307 {
1308 /* N+1 = mf, f is an odd probable prime, and 2f-1 > sqrt(n).
1309 * Now find a Lucas sequence V_k with discriminant D s.t. D/N = -1
1310 * where N divides V_(N+1)/2 and N does not divide V_m/2. */
1311 IV d, p, q;
1312 mpz_t U, V, k;
1313
1314 mpz_init(U); mpz_init(V); mpz_init(k);
1315
1316 /* Primo gave me the idea of this p/q selection method */
1317 for (q = 2; q < 1000; q++) {
1318 p = (q % 2) ? 2 : 1;
1319 d = p*p - 4*q;
1320 mpz_set_si(t, d);
1321 if (mpz_jacobi(t, n) != -1)
1322 continue;
1323 /* we have a d/p/q where d = -1. Check the Lucas sequences. */
1324 mpz_divexact_ui(k, m, 2);
1325 lucas_seq(U, V, n, p, q, k, t, t2);
1326 if (mpz_sgn(V) != 0) {
1327 mpz_divexact_ui(k, np1, 2);
1328 lucas_seq(U, V, n, p, q, k, t, t2);
1329 if (mpz_sgn(V) == 0) {
1330 rval = 2;
1331 if (lp) *lp = p;
1332 if (lq) *lq = q;
1333 break;
1334 }
1335 }
1336 }
1337 mpz_clear(U); mpz_clear(V); mpz_clear(k);
1338 }
1339
1340 end_bls15:
1341 /* Somehow there is a tester getting 0 for LQ */
1342 if (rval && lq && *lq < 2) croak("Internal error in BLS15\n");
1343 mpz_clear(np1); mpz_clear(m); mpz_clear(t); mpz_clear(t2);
1344 return rval;
1345 }
1346
1347 /* Given an n, try using BLS75 theorem 15, N+1 = mq.
1348 * Note: this does _not_ prove n is prime! If it returns 1, then we have
1349 * found a q/D that satisfy theorem 15, but we leave proving q for the caller.
1350 */
_GMP_primality_bls_np1_split(mpz_t n,int effort,mpz_t q,IV * lp,IV * lq)1351 int _GMP_primality_bls_np1_split(mpz_t n, int effort, mpz_t q, IV* lp, IV* lq)
1352 {
1353 mpz_t np1, m, f, sqrtn, t;
1354 int e, success = 1;
1355 UV B1 = 2000;
1356
1357 /* We need to do this for BLS */
1358 if (mpz_even_p(n)) return 0;
1359
1360 mpz_init(np1); mpz_init(m); mpz_init(f); mpz_init(sqrtn); mpz_init(t);
1361 mpz_add_ui(np1, n, 1);
1362 mpz_set_ui(m, 1);
1363 mpz_set(q, np1);
1364 mpz_sqrt(sqrtn, n);
1365
1366 small_factor(m, q, B1);
1367
1368 while (success) {
1369 success = 0;
1370 mpz_mul_ui(t, q, 2);
1371 mpz_sub_ui(t, t, 1);
1372 if (mpz_cmp(t, sqrtn) <= 0)
1373 break;
1374 if (_GMP_is_prob_prime(q)) {
1375 success = 1;
1376 break;
1377 }
1378 for (e = 0; !success && e <= effort; e++)
1379 success = tfe(f, q, e);
1380 if (success) {
1381 mpz_divexact(q, q, f);
1382 if (mpz_cmp(q, f) < 0)
1383 mpz_swap(q, f);
1384 mpz_mul(m, m, f);
1385 }
1386 }
1387
1388 if (success)
1389 success = _GMP_primality_bls_15(n, q, lp, lq);
1390
1391 mpz_clear(np1);
1392 mpz_clear(m);
1393 mpz_clear(f);
1394 mpz_clear(sqrtn);
1395 mpz_clear(t);
1396 return success;
1397 }
1398
1399 /* Given an n, try using BLS75 theorem 3, N-1 = mp. */
_GMP_primality_bls_nm1_split(mpz_t n,int effort,mpz_t p,UV * reta)1400 int _GMP_primality_bls_nm1_split(mpz_t n, int effort, mpz_t p, UV *reta)
1401 {
1402 mpz_t nm1, m, f, sqrtn, t;
1403 int e, success = 1;
1404 UV B1 = 2000;
1405
1406 /* We need to do this for BLS */
1407 if (mpz_even_p(n)) return 0;
1408
1409 mpz_init(nm1); mpz_init(m); mpz_init(f); mpz_init(sqrtn); mpz_init(t);
1410 mpz_sub_ui(nm1, n, 1);
1411 mpz_set_ui(m, 1);
1412 mpz_set(p, nm1);
1413 mpz_sqrt(sqrtn, n);
1414
1415 small_factor(m, p, B1);
1416
1417 while (success) {
1418 success = 0;
1419 mpz_mul_ui(t, p, 2);
1420 mpz_add_ui(t, t, 1);
1421 if (mpz_cmp(t, sqrtn) <= 0)
1422 break;
1423 if (_GMP_is_prob_prime(p)) {
1424 success = 1;
1425 break;
1426 }
1427 for (e = 0; !success && e <= effort; e++)
1428 success = tfe(f, p, e);
1429 if (success) {
1430 mpz_divexact(p, p, f);
1431 if (mpz_cmp(p, f) < 0)
1432 mpz_swap(p, f);
1433 mpz_mul(m, m, f);
1434 }
1435 }
1436
1437 if (success)
1438 success = _GMP_primality_bls_3(n, p, reta);
1439
1440 mpz_clear(nm1);
1441 mpz_clear(m);
1442 mpz_clear(f);
1443 mpz_clear(sqrtn);
1444 mpz_clear(t);
1445 return success;
1446 }
1447