1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
2 * gmpy_mpz_prp.c *
3 * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
4 * Python interface to the GMP or MPIR, MPFR, and MPC multiple precision *
5 * libraries. *
6 * *
7 * Copyright 2011 David Cleaver *
8 * *
9 * Copyright 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019, *
10 * 2020, 2021 Case Van Horsen *
11 * *
12 * The original file is available at: *
13 * <http://sourceforge.net/projects/mpzprp/files/> *
14 * *
15 * Modified by Case Van Horsen for inclusion into GMPY2. *
16 * *
17 * This file is part of GMPY2. *
18 * *
19 * GMPY2 is free software: you can redistribute it and/or modify it under *
20 * the terms of the GNU Lesser General Public License as published by the *
21 * Free Software Foundation, either version 3 of the License, or (at your *
22 * option) any later version. *
23 * *
24 * GMPY2 is distributed in the hope that it will be useful, but WITHOUT *
25 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or *
26 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public *
27 * License for more details. *
28 * *
29 * You should have received a copy of the GNU Lesser General Public *
30 * License along with GMPY2; if not, see <http://www.gnu.org/licenses/> *
31 * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
32
33 /* ******************************************************************
34 * mpz_prp: (also called a Fermat probable prime)
35 * A "probable prime" to the base a is a number n such that,
36 * (a,n)=1 and a^(n-1) = 1 mod n
37 * ******************************************************************/
38
39 PyDoc_STRVAR(doc_mpz_is_fermat_prp,
40 "is_fermat_prp(n,a) -> boolean\n\n"
41 "Return True if n is a Fermat probable prime to the base a.\n"
42 "Assuming:\n"
43 " gcd(n,a) == 1\n"
44 "Then a Fermat probable prime requires:\n"
45 " a**(n-1) == 1 (mod n)");
46
47 static PyObject *
GMPY_mpz_is_fermat_prp(PyObject * self,PyObject * args)48 GMPY_mpz_is_fermat_prp(PyObject *self, PyObject *args)
49 {
50 MPZ_Object *a = NULL, *n = NULL;
51 PyObject *result = NULL;
52 mpz_t res, nm1;
53
54 if (PyTuple_Size(args) != 2) {
55 TYPE_ERROR("is_fermat_prp() requires 2 integer arguments");
56 return NULL;
57 }
58
59 mpz_init(res);
60 mpz_init(nm1);
61
62 n = GMPy_MPZ_From_Integer(PyTuple_GET_ITEM(args, 0), NULL);
63 a = GMPy_MPZ_From_Integer(PyTuple_GET_ITEM(args, 1), NULL);
64 if (!a || !n) {
65 TYPE_ERROR("is_fermat_prp() requires 2 integer arguments");
66 goto cleanup;
67 }
68
69 /* Require a >= 2. */
70 if (mpz_cmp_ui(a->z, 2) < 0) {
71 VALUE_ERROR("is_fermat_prp() requires 'a' greater than or equal to 2");
72 goto cleanup;
73 }
74
75 /* Require n > 0. */
76 if (mpz_sgn(n->z) <= 0) {
77 VALUE_ERROR("is_fermat_prp() requires 'n' be greater than 0");
78 goto cleanup;
79 }
80
81 /* Check for n == 1 */
82 if (mpz_cmp_ui(n->z, 1) == 0) {
83 result = Py_False;
84 goto cleanup;
85 }
86
87 /* Handle n even. */
88 /* Should n even raise an exception? */
89 if (mpz_divisible_ui_p(n->z, 2)) {
90 if (mpz_cmp_ui(n->z, 2) == 0)
91 result = Py_True;
92 else
93 result = Py_False;
94 goto cleanup;
95 }
96
97 /* Check gcd(a,n) */
98 mpz_gcd(res, n->z, a->z);
99 if (mpz_cmp_ui(res, 1) > 0) {
100 VALUE_ERROR("is_fermat_prp() requires gcd(n,a) == 1");
101 goto cleanup;
102 }
103
104 mpz_set(nm1, n->z);
105 mpz_sub_ui(nm1, nm1, 1);
106 mpz_powm(res, a->z, nm1, n->z);
107
108 if (mpz_cmp_ui(res, 1) == 0)
109 result = Py_True;
110 else
111 result = Py_False;
112
113 cleanup:
114 Py_XINCREF(result);
115 mpz_clear(res);
116 mpz_clear(nm1);
117 Py_XDECREF((PyObject*)a);
118 Py_XDECREF((PyObject*)n);
119 return result;
120 }
121
122 /* *************************************************************************
123 * mpz_euler_prp: (also called a Solovay-Strassen probable prime)
124 * An "Euler probable prime" to the base a is an odd composite number n with,
125 * (a,n)=1 such that a^((n-1)/2)=(a/n) mod n [(a/n) is the Jacobi symbol]
126 * *************************************************************************/
127
128 PyDoc_STRVAR(doc_mpz_is_euler_prp,
129 "is_euler_prp(n,a) -> boolean\n\n"
130 "Return True if n is an Euler (also known as Solovay-Strassen)\n"
131 "probable prime to the base a.\n"
132 "Assuming:\n"
133 " gcd(n,a) == 1\n"
134 " n is odd\n"
135 "Then an Euler probable prime requires:\n"
136 " a**((n-1)/2) == 1 (mod n)");
137
138 static PyObject *
GMPY_mpz_is_euler_prp(PyObject * self,PyObject * args)139 GMPY_mpz_is_euler_prp(PyObject *self, PyObject *args)
140 {
141 MPZ_Object *a = NULL, *n = NULL;
142 PyObject *result = NULL;
143 mpz_t res, exp;
144 int ret;
145
146 if (PyTuple_Size(args) != 2) {
147 TYPE_ERROR("is_euler_prp() requires 2 integer arguments");
148 return NULL;
149 }
150
151 mpz_init(res);
152 mpz_init(exp);
153
154 n = GMPy_MPZ_From_Integer(PyTuple_GET_ITEM(args, 0), NULL);
155 a = GMPy_MPZ_From_Integer(PyTuple_GET_ITEM(args, 1), NULL);
156 if (!a || !n) {
157 TYPE_ERROR("is_euler_prp() requires 2 integer arguments");
158 goto cleanup;
159 }
160
161 /* Require a >= 2. */
162 if (mpz_cmp_ui(a->z, 2) < 0) {
163 VALUE_ERROR("is_euler_prp() requires 'a' greater than or equal to 2");
164 goto cleanup;
165 }
166
167 /* Require n > 0. */
168 if (mpz_sgn(n->z) <= 0) {
169 VALUE_ERROR("is_euler_prp() requires 'n' be greater than 0");
170 goto cleanup;
171 }
172
173 /* Check for n == 1 */
174 if (mpz_cmp_ui(n->z, 1) == 0) {
175 result = Py_False;
176 goto cleanup;
177 }
178
179 /* Handle n even. */
180 if (mpz_divisible_ui_p(n->z, 2)) {
181 if (mpz_cmp_ui(n->z, 2) == 0)
182 result = Py_True;
183 else
184 result = Py_False;
185 goto cleanup;
186 }
187
188 /* Check gcd(a,b) */
189 mpz_gcd(res, n->z, a->z);
190 if (mpz_cmp_ui(res, 1) > 0) {
191 VALUE_ERROR("is_euler_prp() requires gcd(n,a) == 1");
192 goto cleanup;
193 }
194
195 mpz_set(exp, n->z);
196 mpz_sub_ui(exp, exp, 1);
197 mpz_divexact_ui(exp, exp, 2);
198 mpz_powm(res, a->z, exp, n->z);
199
200 /* reuse exp to calculate jacobi(a,n) mod n */
201 ret = mpz_jacobi(a->z,n->z);
202 mpz_set(exp, n->z);
203 if (ret == -1)
204 mpz_sub_ui(exp, exp, 1);
205 else if (ret == 1)
206 mpz_add_ui(exp, exp, 1);
207 mpz_mod(exp, exp, n->z);
208
209 if (mpz_cmp(res, exp) == 0)
210 result = Py_True;
211 else
212 result = Py_False;
213
214 cleanup:
215 Py_XINCREF(result);
216 mpz_clear(res);
217 mpz_clear(exp);
218 Py_XDECREF((PyObject*)a);
219 Py_XDECREF((PyObject*)n);
220 return result;
221 }
222
223 /* *********************************************************************************************
224 * mpz_sprp: (also called a Miller-Rabin probable prime)
225 * A "strong probable prime" to the base a is an odd composite n = (2^r)*s+1 with s odd such that
226 * either a^s == 1 mod n, or a^((2^t)*s) == -1 mod n, for some integer t, with 0 <= t < r.
227 * *********************************************************************************************/
228
229 PyDoc_STRVAR(doc_mpz_is_strong_prp,
230 "is_strong_prp(n,a) -> boolean\n\n"
231 "Return True if n is an strong (also known as Miller-Rabin)\n"
232 "probable prime to the base a.\n"
233 "Assuming:\n"
234 " gcd(n,a) == 1\n"
235 " n is odd\n"
236 " n = s*(2**r) + 1, with s odd\n"
237 "Then a strong probable prime requires one of the following is true:\n"
238 " a**s == 1 (mod n)\n"
239 " or\n"
240 " a**(s*(2**t)) == -1 (mod n) for some t, 0 <= t < r.");
241
242 static PyObject *
GMPY_mpz_is_strong_prp(PyObject * self,PyObject * args)243 GMPY_mpz_is_strong_prp(PyObject *self, PyObject *args)
244 {
245 MPZ_Object *a = NULL, *n = NULL;
246 PyObject *result = NULL;
247 mpz_t s, nm1, mpz_test;
248 mp_bitcnt_t r = 0;
249
250 if (PyTuple_Size(args) != 2) {
251 TYPE_ERROR("is_strong_prp() requires 2 integer arguments");
252 return NULL;
253 }
254
255 mpz_init(s);
256 mpz_init(nm1);
257 mpz_init(mpz_test);
258
259 n = GMPy_MPZ_From_Integer(PyTuple_GET_ITEM(args, 0), NULL);
260 a = GMPy_MPZ_From_Integer(PyTuple_GET_ITEM(args, 1), NULL);
261 if (!a || !n) {
262 TYPE_ERROR("is_strong_prp() requires 2 integer arguments");
263 goto cleanup;
264 }
265
266 /* Require a >= 2. */
267 if (mpz_cmp_ui(a->z, 2) < 0) {
268 VALUE_ERROR("is_strong_prp() requires 'a' greater than or equal to 2");
269 goto cleanup;
270 }
271
272 /* Require n > 0. */
273 if (mpz_sgn(n->z) <= 0) {
274 VALUE_ERROR("is_strong_prp() requires 'n' be greater than 0");
275 goto cleanup;
276 }
277
278 /* Check for n == 1 */
279 if (mpz_cmp_ui(n->z, 1) == 0) {
280 result = Py_False;
281 goto cleanup;
282 }
283
284 /* Handle n even. */
285 if (mpz_divisible_ui_p(n->z, 2)) {
286 if (mpz_cmp_ui(n->z, 2) == 0)
287 result = Py_True;
288 else
289 result = Py_False;
290 goto cleanup;
291 }
292
293 /* Check gcd(a,b) */
294 mpz_gcd(s, n->z, a->z);
295 if (mpz_cmp_ui(s, 1) > 0) {
296 VALUE_ERROR("is_strong_prp() requires gcd(n,a) == 1");
297 goto cleanup;
298 }
299
300 mpz_set(nm1, n->z);
301 mpz_sub_ui(nm1, nm1, 1);
302
303 /* Find s and r satisfying: n-1=(2^r)*s, s odd */
304 r = mpz_scan1(nm1, 0);
305 mpz_fdiv_q_2exp(s, nm1, r);
306
307
308 /* Check a^((2^t)*s) mod n for 0 <= t < r */
309 mpz_powm(mpz_test, a->z, s, n->z);
310 if ((mpz_cmp_ui(mpz_test, 1) == 0) || (mpz_cmp(mpz_test, nm1) == 0)) {
311 result = Py_True;
312 goto cleanup;
313 }
314
315 while (--r) {
316 /* mpz_test = mpz_test^2%n */
317 mpz_mul(mpz_test, mpz_test, mpz_test);
318 mpz_mod(mpz_test, mpz_test, n->z);
319
320 if (mpz_cmp(mpz_test, nm1) == 0) {
321 result = Py_True;
322 goto cleanup;
323 }
324 }
325
326 result = Py_False;
327 cleanup:
328 Py_XINCREF(result);
329 mpz_clear(s);
330 mpz_clear(nm1);
331 mpz_clear(mpz_test);
332 Py_XDECREF((PyObject*)a);
333 Py_XDECREF((PyObject*)n);
334 return result;
335 }
336
337 /* *************************************************************************
338 * mpz_fibonacci_prp:
339 * A "Fibonacci probable prime" with parameters (P,Q), P > 0, Q=+/-1, is a
340 * composite n for which V_n == P mod n
341 * [V is the Lucas V sequence with parameters P,Q]
342 * *************************************************************************/
343
344 PyDoc_STRVAR(doc_mpz_is_fibonacci_prp,
345 "is_fibonacci_prp(n,p,q) -> boolean\n\n"
346 "Return True if n is an Fibonacci probable prime with parameters (p,q).\n"
347 "Assuming:\n"
348 " n is odd\n"
349 " p > 0, q = +/-1\n"
350 " p*p - 4*q != 0\n"
351 "Then a Fibonacci probable prime requires:\n"
352 " lucasv(p,q,n) == p (mod n).");
353
354 static PyObject *
GMPY_mpz_is_fibonacci_prp(PyObject * self,PyObject * args)355 GMPY_mpz_is_fibonacci_prp(PyObject *self, PyObject *args)
356 {
357 MPZ_Object *n = NULL, *p = NULL, *q = NULL;
358 PyObject *result = NULL;
359 mpz_t pmodn, zP;
360 /* used for calculating the Lucas V sequence */
361 mpz_t vl, vh, ql, qh, tmp;
362 mp_bitcnt_t s = 0, j = 0;
363
364 if (PyTuple_Size(args) != 3) {
365 TYPE_ERROR("is_fibonacci_prp() requires 3 integer arguments");
366 return NULL;
367 }
368
369 mpz_init(pmodn);
370 mpz_init(zP);
371 mpz_init(vl);
372 mpz_init(vh);
373 mpz_init(ql);
374 mpz_init(qh);
375 mpz_init(tmp);
376
377 n = GMPy_MPZ_From_Integer(PyTuple_GET_ITEM(args, 0), NULL);
378 p = GMPy_MPZ_From_Integer(PyTuple_GET_ITEM(args, 1), NULL);
379 q = GMPy_MPZ_From_Integer(PyTuple_GET_ITEM(args, 2), NULL);
380 if (!n || !p || !q) {
381 TYPE_ERROR("is_fibonacci_prp() requires 3 integer arguments");
382 goto cleanup;
383 }
384
385 /* Check if p*p - 4*q == 0. */
386
387 mpz_mul(tmp, p->z, p->z);
388 mpz_mul_ui(qh, q->z, 4);
389 mpz_sub(tmp, tmp, qh);
390 if (mpz_sgn(tmp) == 0) {
391 VALUE_ERROR("invalid values for p,q in is_fibonacci_prp()");
392 goto cleanup;
393 }
394
395 /* Verify q = +/-1 */
396
397 if ((mpz_cmp_si(q->z, 1) && mpz_cmp_si(q->z, -1)) || (mpz_sgn(p->z) <= 0)) {
398 VALUE_ERROR("invalid values for p,q in is_fibonacci_prp()");
399 goto cleanup;
400 }
401
402 /* Require n > 0. */
403 if (mpz_sgn(n->z) <= 0) {
404 VALUE_ERROR("is_fibonacci_prp() requires 'n' be greater than 0");
405 goto cleanup;
406 }
407
408 /* Check for n == 1 */
409 if (mpz_cmp_ui(n->z, 1) == 0) {
410 result = Py_False;
411 goto cleanup;
412 }
413
414 /* Handle n even. */
415 if (mpz_divisible_ui_p(n->z, 2)) {
416 if (mpz_cmp_ui(n->z, 2) == 0)
417 result = Py_True;
418 else
419 result = Py_False;
420 goto cleanup;
421 }
422
423 mpz_set(zP, p->z);
424 mpz_mod(pmodn, zP, n->z);
425
426 /* mpz_lucasvmod(res, p, q, n, n); */
427 mpz_set_si(vl, 2);
428 mpz_set(vh, p->z);
429 mpz_set_si(ql, 1);
430 mpz_set_si(qh, 1);
431 mpz_set_si(tmp,0);
432
433 s = mpz_scan1(n->z, 0);
434 for (j = mpz_sizeinbase(n->z,2)-1; j >= s+1; j--) {
435 /* ql = ql*qh (mod n) */
436 mpz_mul(ql, ql, qh);
437 mpz_mod(ql, ql, n->z);
438 if (mpz_tstbit(n->z,j) == 1) {
439 /* qh = ql*q */
440 mpz_mul(qh, ql, q->z);
441
442 /* vl = vh*vl - p*ql (mod n) */
443 mpz_mul(vl, vh, vl);
444 mpz_mul(tmp, ql, p->z);
445 mpz_sub(vl, vl, tmp);
446 mpz_mod(vl, vl, n->z);
447
448 /* vh = vh*vh - 2*qh (mod n) */
449 mpz_mul(vh, vh, vh);
450 mpz_mul_si(tmp, qh, 2);
451 mpz_sub(vh, vh, tmp);
452 mpz_mod(vh, vh, n->z);
453 }
454 else {
455 /* qh = ql */
456 mpz_set(qh, ql);
457
458 /* vh = vh*vl - p*ql (mod n) */
459 mpz_mul(vh, vh, vl);
460 mpz_mul(tmp, ql, p->z);
461 mpz_sub(vh, vh, tmp);
462 mpz_mod(vh, vh, n->z);
463
464 /* vl = vl*vl - 2*ql (mod n) */
465 mpz_mul(vl, vl, vl);
466 mpz_mul_si(tmp, ql, 2);
467 mpz_sub(vl, vl, tmp);
468 mpz_mod(vl, vl, n->z);
469 }
470 }
471 /* ql = ql*qh */
472 mpz_mul(ql, ql, qh);
473
474 /* qh = ql*q */
475 mpz_mul(qh, ql, q->z);
476
477 /* vl = vh*vl - p*ql */
478 mpz_mul(vl, vh, vl);
479 mpz_mul(tmp, ql, p->z);
480 mpz_sub(vl, vl, tmp);
481
482 /* ql = ql*qh */
483 mpz_mul(ql, ql, qh);
484
485 for (j = 1; j <= s; j++) {
486 /* vl = vl*vl - 2*ql (mod n) */
487 mpz_mul(vl, vl, vl);
488 mpz_mul_si(tmp, ql, 2);
489 mpz_sub(vl, vl, tmp);
490 mpz_mod(vl, vl, n->z);
491
492 /* ql = ql*ql (mod n) */
493 mpz_mul(ql, ql, ql);
494 mpz_mod(ql, ql, n->z);
495 }
496
497 /* vl contains our return value */
498 mpz_mod(vl, vl, n->z);
499
500 if (mpz_cmp(vl, pmodn) == 0)
501 result = Py_True;
502 else
503 result = Py_False;
504
505 cleanup:
506 Py_XINCREF(result);
507 mpz_clear(pmodn);
508 mpz_clear(zP);
509 mpz_clear(vl);
510 mpz_clear(vh);
511 mpz_clear(ql);
512 mpz_clear(qh);
513 mpz_clear(tmp);
514 Py_XDECREF((PyObject*)p);
515 Py_XDECREF((PyObject*)q);
516 Py_XDECREF((PyObject*)n);
517 return result;
518 }
519
520
521 /* *******************************************************************************
522 * mpz_lucas_prp:
523 * A "Lucas probable prime" with parameters (P,Q) is a composite n with D=P^2-4Q,
524 * (n,2QD)=1 such that U_(n-(D/n)) == 0 mod n [(D/n) is the Jacobi symbol]
525 * *******************************************************************************/
526
527 PyDoc_STRVAR(doc_mpz_is_lucas_prp,
528 "is_lucas_prp(n,p,q) -> boolean\n\n"
529 "Return True if n is a Lucas probable prime with parameters (p,q).\n"
530 "Assuming:\n"
531 " n is odd\n"
532 " D = p*p - 4*q, D != 0\n"
533 " gcd(n, 2*q*D) == 1\n"
534 "Then a Lucas probable prime requires:\n"
535 " lucasu(p,q,n - Jacobi(D,n)) == 0 (mod n)");
536
537 static PyObject *
GMPY_mpz_is_lucas_prp(PyObject * self,PyObject * args)538 GMPY_mpz_is_lucas_prp(PyObject *self, PyObject *args)
539 {
540 MPZ_Object *n = NULL, *p = NULL, *q = NULL;
541 PyObject *result = NULL;
542 mpz_t zD, res, index;
543 /* used for calculating the Lucas U sequence */
544 mpz_t uh, vl, vh, ql, qh, tmp;
545 mp_bitcnt_t s = 0, j = 0;
546 int ret;
547
548 if (PyTuple_Size(args) != 3) {
549 TYPE_ERROR("is_lucas_prp() requires 3 integer arguments");
550 return NULL;
551 }
552
553 mpz_init(zD);
554 mpz_init(res);
555 mpz_init(index);
556 mpz_init(uh);
557 mpz_init(vl);
558 mpz_init(vh);
559 mpz_init(ql);
560 mpz_init(qh);
561 mpz_init(tmp);
562
563 n = GMPy_MPZ_From_Integer(PyTuple_GET_ITEM(args, 0), NULL);
564 p = GMPy_MPZ_From_Integer(PyTuple_GET_ITEM(args, 1), NULL);
565 q = GMPy_MPZ_From_Integer(PyTuple_GET_ITEM(args, 2), NULL);
566 if (!n || !p || !q) {
567 TYPE_ERROR("is_lucas_prp() requires 3 integer arguments");
568 goto cleanup;
569 }
570
571 /* Check if p*p - 4*q == 0. */
572 mpz_mul(zD, p->z, p->z);
573 mpz_mul_ui(tmp, q->z, 4);
574 mpz_sub(zD, zD, tmp);
575 if (mpz_sgn(zD) == 0) {
576 VALUE_ERROR("invalid values for p,q in is_lucas_prp()");
577 goto cleanup;
578 }
579
580 /* Require n > 0. */
581 if (mpz_sgn(n->z) <= 0) {
582 VALUE_ERROR("is_lucas_prp() requires 'n' be greater than 0");
583 goto cleanup;
584 }
585
586 /* Check for n == 1 */
587 if (mpz_cmp_ui(n->z, 1) == 0) {
588 result = Py_False;
589 goto cleanup;
590 }
591
592 /* Handle n even. */
593 if (mpz_divisible_ui_p(n->z, 2)) {
594 if (mpz_cmp_ui(n->z, 2) == 0)
595 result = Py_True;
596 else
597 result = Py_False;
598 goto cleanup;
599 }
600
601 /* Check GCD */
602 mpz_mul(res, zD, q->z);
603 mpz_mul_ui(res, res, 2);
604 mpz_gcd(res, res, n->z);
605 if ((mpz_cmp(res, n->z) != 0) && (mpz_cmp_ui(res, 1) > 0)) {
606 VALUE_ERROR("is_lucas_prp() requires gcd(n,2*q*D) == 1");
607 goto cleanup;
608 }
609
610 /* index = n-(D/n), where (D/n) is the Jacobi symbol */
611 mpz_set(index, n->z);
612 ret = mpz_jacobi(zD, n->z);
613 if (ret == -1)
614 mpz_add_ui(index, index, 1);
615 else if (ret == 1)
616 mpz_sub_ui(index, index, 1);
617
618 /* mpz_lucasumod(res, p, q, index, n); */
619 mpz_set_si(uh, 1);
620 mpz_set_si(vl, 2);
621 mpz_set(vh, p->z);
622 mpz_set_si(ql, 1);
623 mpz_set_si(qh, 1);
624 mpz_set_si(tmp,0);
625
626 s = mpz_scan1(index, 0);
627 for (j = mpz_sizeinbase(index,2)-1; j >= s+1; j--) {
628 /* ql = ql*qh (mod n) */
629 mpz_mul(ql, ql, qh);
630 mpz_mod(ql, ql, n->z);
631 if (mpz_tstbit(index,j) == 1) {
632 /* qh = ql*q */
633 mpz_mul(qh, ql, q->z);
634
635 /* uh = uh*vh (mod n) */
636 mpz_mul(uh, uh, vh);
637 mpz_mod(uh, uh, n->z);
638
639 /* vl = vh*vl - p*ql (mod n) */
640 mpz_mul(vl, vh, vl);
641 mpz_mul(tmp, ql, p->z);
642 mpz_sub(vl, vl, tmp);
643 mpz_mod(vl, vl, n->z);
644
645 /* vh = vh*vh - 2*qh (mod n) */
646 mpz_mul(vh, vh, vh);
647 mpz_mul_si(tmp, qh, 2);
648 mpz_sub(vh, vh, tmp);
649 mpz_mod(vh, vh, n->z);
650 }
651 else {
652 /* qh = ql */
653 mpz_set(qh, ql);
654
655 /* uh = uh*vl - ql (mod n) */
656 mpz_mul(uh, uh, vl);
657 mpz_sub(uh, uh, ql);
658 mpz_mod(uh, uh, n->z);
659
660 /* vh = vh*vl - p*ql (mod n) */
661 mpz_mul(vh, vh, vl);
662 mpz_mul(tmp, ql, p->z);
663 mpz_sub(vh, vh, tmp);
664 mpz_mod(vh, vh, n->z);
665
666 /* vl = vl*vl - 2*ql (mod n) */
667 mpz_mul(vl, vl, vl);
668 mpz_mul_si(tmp, ql, 2);
669 mpz_sub(vl, vl, tmp);
670 mpz_mod(vl, vl, n->z);
671 }
672 }
673 /* ql = ql*qh */
674 mpz_mul(ql, ql, qh);
675
676 /* qh = ql*q */
677 mpz_mul(qh, ql, q->z);
678
679 /* uh = uh*vl - ql */
680 mpz_mul(uh, uh, vl);
681 mpz_sub(uh, uh, ql);
682
683 /* vl = vh*vl - p*ql */
684 mpz_mul(vl, vh, vl);
685 mpz_mul(tmp, ql, p->z);
686 mpz_sub(vl, vl, tmp);
687
688 /* ql = ql*qh */
689 mpz_mul(ql, ql, qh);
690
691 for (j = 1; j <= s; j++) {
692 /* uh = uh*vl (mod n) */
693 mpz_mul(uh, uh, vl);
694 mpz_mod(uh, uh, n->z);
695
696 /* vl = vl*vl - 2*ql (mod n) */
697 mpz_mul(vl, vl, vl);
698 mpz_mul_si(tmp, ql, 2);
699 mpz_sub(vl, vl, tmp);
700 mpz_mod(vl, vl, n->z);
701
702 /* ql = ql*ql (mod n) */
703 mpz_mul(ql, ql, ql);
704 mpz_mod(ql, ql, n->z);
705 }
706
707 /* uh contains our return value */
708 mpz_mod(res, uh, n->z);
709 if (mpz_cmp_ui(res, 0) == 0)
710 result = Py_True;
711 else
712 result = Py_False;
713
714 cleanup:
715 Py_XINCREF(result);
716 mpz_clear(zD);
717 mpz_clear(res);
718 mpz_clear(index);
719 mpz_clear(uh);
720 mpz_clear(vl);
721 mpz_clear(vh);
722 mpz_clear(ql);
723 mpz_clear(qh);
724 mpz_clear(tmp);
725 Py_XDECREF((PyObject*)p);
726 Py_XDECREF((PyObject*)q);
727 Py_XDECREF((PyObject*)n);
728 return result;
729 }
730
731 /* *********************************************************************************************
732 * mpz_stronglucas_prp:
733 * A "strong Lucas probable prime" with parameters (P,Q) is a composite n = (2^r)*s+(D/n), where
734 * s is odd, D=P^2-4Q, and (n,2QD)=1 such that either U_s == 0 mod n or V_((2^t)*s) == 0 mod n
735 * for some t, 0 <= t < r. [(D/n) is the Jacobi symbol]
736 * *********************************************************************************************/
737
738 PyDoc_STRVAR(doc_mpz_is_stronglucas_prp,
739 "is_strong_lucas_prp(n,p,q) -> boolean\n\n"
740 "Return True if n is a strong Lucas probable prime with parameters (p,q).\n"
741 "Assuming:\n"
742 " n is odd\n"
743 " D = p*p - 4*q, D != 0\n"
744 " gcd(n, 2*q*D) == 1\n"
745 " n = s*(2**r) + Jacobi(D,n), s odd\n"
746 "Then a strong Lucas probable prime requires:\n"
747 " lucasu(p,q,s) == 0 (mod n)\n"
748 " or\n"
749 " lucasv(p,q,s*(2**t)) == 0 (mod n) for some t, 0 <= t < r");
750
751 static PyObject *
GMPY_mpz_is_stronglucas_prp(PyObject * self,PyObject * args)752 GMPY_mpz_is_stronglucas_prp(PyObject *self, PyObject *args)
753 {
754 MPZ_Object *n = NULL, *p= NULL, *q = NULL;
755 PyObject *result = NULL;
756 mpz_t zD, s, nmj, res;
757 /* these are needed for the LucasU and LucasV part of this function */
758 mpz_t uh, vl, vh, ql, qh, tmp;
759 mp_bitcnt_t r = 0, j = 0;
760 int ret = 0;
761
762 if (PyTuple_Size(args) != 3) {
763 TYPE_ERROR("is_strong_lucas_prp() requires 3 integer arguments");
764 return NULL;
765 }
766
767 mpz_init(zD);
768 mpz_init(s);
769 mpz_init(nmj);
770 mpz_init(res);
771 mpz_init(uh);
772 mpz_init(vl);
773 mpz_init(vh);
774 mpz_init(ql);
775 mpz_init(qh);
776 mpz_init(tmp);
777
778 n = GMPy_MPZ_From_Integer(PyTuple_GET_ITEM(args, 0), NULL);
779 p = GMPy_MPZ_From_Integer(PyTuple_GET_ITEM(args, 1), NULL);
780 q = GMPy_MPZ_From_Integer(PyTuple_GET_ITEM(args, 2), NULL);
781 if (!n || !p || !q) {
782 TYPE_ERROR("is_strong_lucas_prp() requires 3 integer arguments");
783 goto cleanup;
784 }
785
786 /* Check if p*p - 4*q == 0. */
787 mpz_mul(zD, p->z, p->z);
788 mpz_mul_ui(tmp, q->z, 4);
789 mpz_sub(zD, zD, tmp);
790 if (mpz_sgn(zD) == 0) {
791 VALUE_ERROR("invalid values for p,q in is_strong_lucas_prp()");
792 goto cleanup;
793 }
794
795 /* Require n > 0. */
796 if (mpz_sgn(n->z) <= 0) {
797 VALUE_ERROR("is_strong_lucas_prp() requires 'n' be greater than 0");
798 goto cleanup;
799 }
800
801 /* Check for n == 1 */
802 if (mpz_cmp_ui(n->z, 1) == 0) {
803 result = Py_False;
804 goto cleanup;
805 }
806
807 /* Handle n even. */
808 if (mpz_divisible_ui_p(n->z, 2)) {
809 if (mpz_cmp_ui(n->z, 2) == 0)
810 result = Py_True;
811 else
812 result = Py_False;
813 goto cleanup;
814 }
815
816 /* Check GCD */
817 mpz_mul(res, zD, q->z);
818 mpz_mul_ui(res, res, 2);
819 mpz_gcd(res, res, n->z);
820 if ((mpz_cmp(res, n->z) != 0) && (mpz_cmp_ui(res, 1) > 0)) {
821 VALUE_ERROR("is_strong_lucas_prp() requires gcd(n,2*q*D) == 1");
822 goto cleanup;
823 }
824
825 /* nmj = n - (D/n), where (D/n) is the Jacobi symbol */
826 mpz_set(nmj, n->z);
827 ret = mpz_jacobi(zD, n->z);
828 if (ret == -1)
829 mpz_add_ui(nmj, nmj, 1);
830 else if (ret == 1)
831 mpz_sub_ui(nmj, nmj, 1);
832
833 r = mpz_scan1(nmj, 0);
834 mpz_fdiv_q_2exp(s, nmj, r);
835
836 /* make sure U_s == 0 mod n or V_((2^t)*s) == 0 mod n, for some t, 0 <= t < r */
837 mpz_set_si(uh, 1);
838 mpz_set_si(vl, 2);
839 mpz_set(vh, p->z);
840 mpz_set_si(ql, 1);
841 mpz_set_si(qh, 1);
842 mpz_set_si(tmp,0);
843
844 for (j = mpz_sizeinbase(s,2)-1; j >= 1; j--) {
845 /* ql = ql*qh (mod n) */
846 mpz_mul(ql, ql, qh);
847 mpz_mod(ql, ql, n->z);
848 if (mpz_tstbit(s,j) == 1) {
849 /* qh = ql*q */
850 mpz_mul(qh, ql, q->z);
851
852 /* uh = uh*vh (mod n) */
853 mpz_mul(uh, uh, vh);
854 mpz_mod(uh, uh, n->z);
855
856 /* vl = vh*vl - p*ql (mod n) */
857 mpz_mul(vl, vh, vl);
858 mpz_mul(tmp, ql, p->z);
859 mpz_sub(vl, vl, tmp);
860 mpz_mod(vl, vl, n->z);
861
862 /* vh = vh*vh - 2*qh (mod n) */
863 mpz_mul(vh, vh, vh);
864 mpz_mul_si(tmp, qh, 2);
865 mpz_sub(vh, vh, tmp);
866 mpz_mod(vh, vh, n->z);
867 }
868 else {
869 /* qh = ql */
870 mpz_set(qh, ql);
871
872 /* uh = uh*vl - ql (mod n) */
873 mpz_mul(uh, uh, vl);
874 mpz_sub(uh, uh, ql);
875 mpz_mod(uh, uh, n->z);
876
877 /* vh = vh*vl - p*ql (mod n) */
878 mpz_mul(vh, vh, vl);
879 mpz_mul(tmp, ql, p->z);
880 mpz_sub(vh, vh, tmp);
881 mpz_mod(vh, vh, n->z);
882
883 /* vl = vl*vl - 2*ql (mod n) */
884 mpz_mul(vl, vl, vl);
885 mpz_mul_si(tmp, ql, 2);
886 mpz_sub(vl, vl, tmp);
887 mpz_mod(vl, vl, n->z);
888 }
889 }
890 /* ql = ql*qh */
891 mpz_mul(ql, ql, qh);
892
893 /* qh = ql*q */
894 mpz_mul(qh, ql, q->z);
895
896 /* uh = uh*vl - ql */
897 mpz_mul(uh, uh, vl);
898 mpz_sub(uh, uh, ql);
899
900 /* vl = vh*vl - p*ql */
901 mpz_mul(vl, vh, vl);
902 mpz_mul(tmp, ql, p->z);
903 mpz_sub(vl, vl, tmp);
904
905 /* ql = ql*qh */
906 mpz_mul(ql, ql, qh);
907
908 mpz_mod(uh, uh, n->z);
909 mpz_mod(vl, vl, n->z);
910
911 /* uh contains LucasU_s and vl contains LucasV_s */
912 if ((mpz_cmp_ui(uh, 0) == 0) || (mpz_cmp_ui(vl, 0) == 0)) {
913 result = Py_True;
914 goto cleanup;
915 }
916
917 for (j = 1; j < r; j++) {
918 /* vl = vl*vl - 2*ql (mod n) */
919 mpz_mul(vl, vl, vl);
920 mpz_mul_si(tmp, ql, 2);
921 mpz_sub(vl, vl, tmp);
922 mpz_mod(vl, vl, n->z);
923
924 /* ql = ql*ql (mod n) */
925 mpz_mul(ql, ql, ql);
926 mpz_mod(ql, ql, n->z);
927
928 if (mpz_cmp_ui(vl, 0) == 0) {
929 result = Py_True;
930 goto cleanup;
931 }
932 }
933
934 result = Py_False;
935 cleanup:
936 Py_XINCREF(result);
937 mpz_clear(zD);
938 mpz_clear(s);
939 mpz_clear(nmj);
940 mpz_clear(res);
941 mpz_clear(uh);
942 mpz_clear(vl);
943 mpz_clear(vh);
944 mpz_clear(ql);
945 mpz_clear(qh);
946 mpz_clear(tmp);
947 Py_XDECREF((PyObject*)p);
948 Py_XDECREF((PyObject*)q);
949 Py_XDECREF((PyObject*)n);
950 return result;
951 }
952
953 /* *******************************************************************************************
954 * mpz_extrastronglucas_prp:
955 * Let U_n = LucasU(p,1), V_n = LucasV(p,1), and D=p^2-4.
956 * An "extra strong Lucas probable prime" to the base p is a composite n = (2^r)*s+(D/n), where
957 * s is odd and (n,2D)=1, such that either U_s == 0 mod n or V_s == +/-2 mod n, or
958 * V_((2^t)*s) == 0 mod n for some t with 0 <= t < r-1 [(D/n) is the Jacobi symbol]
959 * *******************************************************************************************/
960
961 PyDoc_STRVAR(doc_mpz_is_extrastronglucas_prp,
962 "is_extra_strong_lucas_prp(n,p) -> boolean\n\n"
963 "Return True if n is an extra strong Lucas probable prime with parameters\n"
964 "(p,1). Assuming:\n"
965 " n is odd\n"
966 " D = p*p - 4, D != 0\n"
967 " gcd(n, 2*D) == 1\n"
968 " n = s*(2**r) + Jacobi(D,n), s odd\n"
969 "Then an extra strong Lucas probable prime requires:\n"
970 " lucasu(p,1,s) == 0 (mod n)\n"
971 " or\n"
972 " lucasv(p,1,s) == +/-2 (mod n)\n"
973 " or\n"
974 " lucasv(p,1,s*(2**t)) == 0 (mod n) for some t, 0 <= t < r");
975
976 static PyObject *
GMPY_mpz_is_extrastronglucas_prp(PyObject * self,PyObject * args)977 GMPY_mpz_is_extrastronglucas_prp(PyObject *self, PyObject *args)
978 {
979 MPZ_Object *n = NULL, *p = NULL;
980 PyObject *result = NULL;
981 mpz_t zD, s, nmj, nm2, res;
982 /* these are needed for the LucasU and LucasV part of this function */
983 mpz_t uh, vl, vh, ql, qh, tmp;
984 mp_bitcnt_t r = 0, j = 0;
985 int ret = 0;
986
987 if (PyTuple_Size(args) != 2) {
988 TYPE_ERROR("is_extra_strong_lucas_prp() requires 2 integer arguments");
989 return NULL;
990 }
991
992 mpz_init(zD);
993 mpz_init(s);
994 mpz_init(nmj);
995 mpz_init(nm2);
996 mpz_init(res);
997 mpz_init(uh);
998 mpz_init(vl);
999 mpz_init(vh);
1000 mpz_init(ql);
1001 mpz_init(qh);
1002 mpz_init(tmp);
1003
1004 n = GMPy_MPZ_From_Integer(PyTuple_GET_ITEM(args, 0), NULL);
1005 p = GMPy_MPZ_From_Integer(PyTuple_GET_ITEM(args, 1), NULL);
1006 if (!n || !p) {
1007 TYPE_ERROR("is_extra_strong_lucas_prp() requires 2 integer arguments");
1008 goto cleanup;
1009 }
1010
1011 /* Check if p*p - 4 == 0. */
1012 mpz_mul(zD, p->z, p->z);
1013 mpz_sub_ui(zD, zD, 4);
1014 if (mpz_sgn(zD) == 0) {
1015 VALUE_ERROR("invalid value for p in is_extra_strong_lucas_prp()");
1016 goto cleanup;
1017 }
1018
1019 /* Require n > 0. */
1020 if (mpz_sgn(n->z) <= 0) {
1021 VALUE_ERROR("is_extra_strong_lucas_prp() requires 'n' be greater than 0");
1022 goto cleanup;
1023 }
1024
1025 /* Check for n == 1 */
1026 if (mpz_cmp_ui(n->z, 1) == 0) {
1027 result = Py_False;
1028 goto cleanup;
1029 }
1030
1031 /* Handle n even. */
1032 if (mpz_divisible_ui_p(n->z, 2)) {
1033 if (mpz_cmp_ui(n->z, 2) == 0)
1034 result = Py_True;
1035 else
1036 result = Py_False;
1037 goto cleanup;
1038 }
1039
1040 /* Check GCD */
1041 mpz_mul_ui(res, zD, 2);
1042 mpz_gcd(res, res, n->z);
1043 if ((mpz_cmp(res, n->z) != 0) && (mpz_cmp_ui(res, 1) > 0)) {
1044 VALUE_ERROR("is_extra_strong_lucas_prp() requires gcd(n,2*D) == 1");
1045 goto cleanup;
1046 }
1047
1048 /* nmj = n - (D/n), where (D/n) is the Jacobi symbol */
1049 mpz_set(nmj, n->z);
1050 ret = mpz_jacobi(zD, n->z);
1051 if (ret == -1)
1052 mpz_add_ui(nmj, nmj, 1);
1053 else if (ret == 1)
1054 mpz_sub_ui(nmj, nmj, 1);
1055
1056 r = mpz_scan1(nmj, 0);
1057 mpz_fdiv_q_2exp(s, nmj, r);
1058
1059 mpz_set(nm2, n->z);
1060 mpz_sub_ui(nm2, nm2, 2);
1061
1062 /* make sure that either U_s == 0 mod n or V_s == +/-2 mod n, or */
1063 /* V_((2^t)*s) == 0 mod n for some t with 0 <= t < r-1 */
1064 mpz_set_si(uh, 1);
1065 mpz_set_si(vl, 2);
1066 mpz_set(vh, p->z);
1067 mpz_set_si(ql, 1);
1068 mpz_set_si(qh, 1);
1069 mpz_set_si(tmp,0);
1070
1071 for (j = mpz_sizeinbase(s,2)-1; j >= 1; j--) {
1072 /* ql = ql*qh (mod n) */
1073 mpz_mul(ql, ql, qh);
1074 mpz_mod(ql, ql, n->z);
1075 if (mpz_tstbit(s,j) == 1) {
1076 /* qh = ql*q */
1077 mpz_set(qh, ql);
1078
1079 /* uh = uh*vh (mod n) */
1080 mpz_mul(uh, uh, vh);
1081 mpz_mod(uh, uh, n->z);
1082
1083 /* vl = vh*vl - p*ql (mod n) */
1084 mpz_mul(vl, vh, vl);
1085 mpz_mul(tmp, ql, p->z);
1086 mpz_sub(vl, vl, tmp);
1087 mpz_mod(vl, vl, n->z);
1088
1089 /* vh = vh*vh - 2*qh (mod n) */
1090 mpz_mul(vh, vh, vh);
1091 mpz_mul_si(tmp, qh, 2);
1092 mpz_sub(vh, vh, tmp);
1093 mpz_mod(vh, vh, n->z);
1094 }
1095 else {
1096 /* qh = ql */
1097 mpz_set(qh, ql);
1098
1099 /* uh = uh*vl - ql (mod n) */
1100 mpz_mul(uh, uh, vl);
1101 mpz_sub(uh, uh, ql);
1102 mpz_mod(uh, uh, n->z);
1103
1104 /* vh = vh*vl - p*ql (mod n) */
1105 mpz_mul(vh, vh, vl);
1106 mpz_mul(tmp, ql, p->z);
1107 mpz_sub(vh, vh, tmp);
1108 mpz_mod(vh, vh, n->z);
1109
1110 /* vl = vl*vl - 2*ql (mod n) */
1111 mpz_mul(vl, vl, vl);
1112 mpz_mul_si(tmp, ql, 2);
1113 mpz_sub(vl, vl, tmp);
1114 mpz_mod(vl, vl, n->z);
1115 }
1116 }
1117 /* ql = ql*qh */
1118 mpz_mul(ql, ql, qh);
1119
1120 /* qh = ql*q */
1121 mpz_set(qh, ql);
1122
1123 /* uh = uh*vl - ql */
1124 mpz_mul(uh, uh, vl);
1125 mpz_sub(uh, uh, ql);
1126
1127 /* vl = vh*vl - p*ql */
1128 mpz_mul(vl, vh, vl);
1129 mpz_mul(tmp, ql, p->z);
1130 mpz_sub(vl, vl, tmp);
1131
1132 /* ql = ql*qh */
1133 mpz_mul(ql, ql, qh);
1134
1135 mpz_mod(uh, uh, n->z);
1136 mpz_mod(vl, vl, n->z);
1137
1138 /* uh contains LucasU_s and vl contains LucasV_s */
1139 if ((mpz_cmp_ui(uh, 0) == 0) || (mpz_cmp_ui(vl, 0) == 0) ||
1140 (mpz_cmp(vl, nm2) == 0) || (mpz_cmp_si(vl, 2) == 0)) {
1141 result = Py_True;
1142 goto cleanup;
1143 }
1144
1145 for (j = 1; j < r-1; j++) {
1146 /* vl = vl*vl - 2*ql (mod n) */
1147 mpz_mul(vl, vl, vl);
1148 mpz_mul_si(tmp, ql, 2);
1149 mpz_sub(vl, vl, tmp);
1150 mpz_mod(vl, vl, n->z);
1151
1152 /* ql = ql*ql (mod n) */
1153 mpz_mul(ql, ql, ql);
1154 mpz_mod(ql, ql, n->z);
1155
1156 if (mpz_cmp_ui(vl, 0) == 0) {
1157 result = Py_True;
1158 goto cleanup;
1159 }
1160 }
1161
1162 result = Py_False;
1163 cleanup:
1164 Py_XINCREF(result);
1165 mpz_clear(zD);
1166 mpz_clear(s);
1167 mpz_clear(nmj);
1168 mpz_clear(nm2);
1169 mpz_clear(res);
1170 mpz_clear(uh);
1171 mpz_clear(vl);
1172 mpz_clear(vh);
1173 mpz_clear(ql);
1174 mpz_clear(qh);
1175 mpz_clear(tmp);
1176 Py_XDECREF((PyObject*)p);
1177 Py_XDECREF((PyObject*)n);
1178 return result;
1179 }
1180
1181 /* ***********************************************************************************************
1182 * mpz_selfridge_prp:
1183 * A "Lucas-Selfridge probable prime" n is a "Lucas probable prime" using Selfridge parameters of:
1184 * Find the first element D in the sequence {5, -7, 9, -11, 13, ...} such that Jacobi(D,n) = -1
1185 * Then use P=1 and Q=(1-D)/4 in the Lucas probable prime test.
1186 * Make sure n is not a perfect square, otherwise the search for D will only stop when D=n.
1187 * ***********************************************************************************************/
1188
1189 PyDoc_STRVAR(doc_mpz_is_selfridge_prp,
1190 "is_selfridge_prp(n) -> boolean\n\n"
1191 "Return True if n is a Lucas probable prime with Selfidge parameters\n"
1192 "(p,q). The Selfridge parameters are chosen by finding the first\n"
1193 "element D in the sequence {5, -7, 9, -11, 13, ...} such that\n"
1194 "Jacobi(D,n) == -1. Then let p=1 and q = (1-D)/4. Then perform\n"
1195 "a Lucas probable prime test.");
1196
1197 static PyObject *
GMPY_mpz_is_selfridge_prp(PyObject * self,PyObject * args)1198 GMPY_mpz_is_selfridge_prp(PyObject *self, PyObject *args)
1199 {
1200 MPZ_Object *n = NULL;
1201 PyObject *result = NULL, *temp = NULL;
1202 long d = 5, p = 1, q = 0, max_d = 1000000;
1203 int jacobi = 0;
1204 mpz_t zD;
1205
1206 if (PyTuple_Size(args) != 1) {
1207 TYPE_ERROR("is_selfridge_prp() requires 1 integer argument");
1208 return NULL;
1209 }
1210
1211 mpz_init(zD);
1212
1213 n = GMPy_MPZ_From_Integer(PyTuple_GET_ITEM(args, 0), NULL);
1214 if (!n) {
1215 TYPE_ERROR("is_selfridge_prp() requires 1 integer argument");
1216 goto cleanup;
1217 }
1218
1219 /* Require n > 0. */
1220 if (mpz_sgn(n->z) <= 0) {
1221 VALUE_ERROR("is_selfridge_prp() requires 'n' be greater than 0");
1222 goto cleanup;
1223 }
1224
1225 /* Check for n == 1 */
1226 if (mpz_cmp_ui(n->z, 1) == 0) {
1227 result = Py_False;
1228 goto cleanup;
1229 }
1230
1231 /* Handle n even. */
1232 if (mpz_divisible_ui_p(n->z, 2)) {
1233 if (mpz_cmp_ui(n->z, 2) == 0)
1234 result = Py_True;
1235 else
1236 result = Py_False;
1237 goto cleanup;
1238 }
1239
1240 mpz_set_ui(zD, d);
1241
1242 while (1) {
1243 jacobi = mpz_jacobi(zD, n->z);
1244
1245 /* if jacobi == 0, d is a factor of n, therefore n is composite... */
1246 /* if d == n, then either n is either prime or 9... */
1247 if (jacobi == 0) {
1248 if ((mpz_cmpabs(zD, n->z) == 0) && (mpz_cmp_ui(zD, 9) != 0)) {
1249 result = Py_True;
1250 goto cleanup;
1251 }
1252 else {
1253 result = Py_False;
1254 goto cleanup;
1255 }
1256 }
1257 if (jacobi == -1)
1258 break;
1259
1260 /* if we get to the 5th d, make sure we aren't dealing with a square... */
1261 if (d == 13) {
1262 if (mpz_perfect_square_p(n->z)) {
1263 result = Py_False;
1264 goto cleanup;
1265 }
1266 }
1267
1268 if (d < 0) {
1269 d *= -1;
1270 d += 2;
1271 }
1272 else {
1273 d += 2;
1274 d *= -1;
1275 }
1276
1277 /* make sure we don't search forever */
1278 if (d >= max_d) {
1279 VALUE_ERROR("appropriate value for D cannot be found in is_selfridge_prp()");
1280 goto cleanup;
1281 }
1282
1283 mpz_set_si(zD, d);
1284 }
1285
1286 q = (1-d)/4;
1287
1288 /* Since "O" is used, the refcount for n is incremented so deleting
1289 * temp will not delete n.
1290 */
1291 temp = Py_BuildValue("Oll", n, p, q);
1292 if (!temp)
1293 goto cleanup;
1294 result = GMPY_mpz_is_lucas_prp(NULL, temp);
1295 Py_DECREF(temp);
1296 goto return_result;
1297
1298 cleanup:
1299 Py_XINCREF(result);
1300 return_result:
1301 mpz_clear(zD);
1302 Py_DECREF((PyObject*)n);
1303 return result;
1304 }
1305
1306 /* *********************************************************************************************************
1307 * mpz_strongselfridge_prp:
1308 * A "strong Lucas-Selfridge probable prime" n is a "strong Lucas probable prime" using Selfridge parameters of:
1309 * Find the first element D in the sequence {5, -7, 9, -11, 13, ...} such that Jacobi(D,n) = -1
1310 * Then use P=1 and Q=(1-D)/4 in the strong Lucas probable prime test.
1311 * Make sure n is not a perfect square, otherwise the search for D will only stop when D=n.
1312 * **********************************************************************************************************/
1313
1314 PyDoc_STRVAR(doc_mpz_is_strongselfridge_prp,
1315 "is_strong_selfridge_prp(n) -> boolean\n\n"
1316 "Return True if n is a strong Lucas probable prime with Selfidge\n"
1317 "parameters (p,q). The Selfridge parameters are chosen by finding\n"
1318 "the first element D in the sequence {5, -7, 9, -11, 13, ...} such\n"
1319 "that Jacobi(D,n) == -1. Then let p=1 and q = (1-D)/4. Then perform\n"
1320 "a strong Lucas probable prime test.");
1321
1322 static PyObject *
GMPY_mpz_is_strongselfridge_prp(PyObject * self,PyObject * args)1323 GMPY_mpz_is_strongselfridge_prp(PyObject *self, PyObject *args)
1324 {
1325 MPZ_Object *n = NULL;
1326 PyObject *result = NULL, *temp = NULL;
1327 long d = 5, p = 1, q = 0, max_d = 1000000;
1328 int jacobi = 0;
1329 mpz_t zD;
1330
1331 if (PyTuple_Size(args) != 1) {
1332 TYPE_ERROR("is_strong_selfridge_prp() requires 1 integer argument");
1333 return NULL;
1334 }
1335
1336 mpz_init(zD);
1337
1338 n = GMPy_MPZ_From_Integer(PyTuple_GET_ITEM(args, 0), NULL);
1339 if (!n) {
1340 TYPE_ERROR("is_strong_selfridge_prp() requires 1 integer argument");
1341 goto cleanup;
1342 }
1343
1344 /* Require n > 0. */
1345 if (mpz_sgn(n->z) <= 0) {
1346 VALUE_ERROR("is_strong_selfridge_prp() requires 'n' be greater than 0");
1347 goto cleanup;
1348 }
1349
1350 /* Check for n == 1 */
1351 if (mpz_cmp_ui(n->z, 1) == 0) {
1352 result = Py_False;
1353 goto cleanup;
1354 }
1355
1356 /* Handle n even. */
1357 if (mpz_divisible_ui_p(n->z, 2)) {
1358 if (mpz_cmp_ui(n->z, 2) == 0)
1359 result = Py_True;
1360 else
1361 result = Py_False;
1362 goto cleanup;
1363 }
1364
1365
1366 mpz_set_ui(zD, d);
1367
1368 while (1) {
1369 jacobi = mpz_jacobi(zD, n->z);
1370
1371 /* if jacobi == 0, d is a factor of n, therefore n is composite... */
1372 /* if d == n, then either n is either prime or 9... */
1373 if (jacobi == 0) {
1374 if ((mpz_cmpabs(zD, n->z) == 0) && (mpz_cmp_ui(zD, 9) != 0)) {
1375 result = Py_True;
1376 goto cleanup;
1377 }
1378 else {
1379 result = Py_False;
1380 goto cleanup;
1381 }
1382 }
1383 if (jacobi == -1)
1384 break;
1385
1386 /* if we get to the 5th d, make sure we aren't dealing with a square... */
1387 if (d == 13) {
1388 if (mpz_perfect_square_p(n->z)) {
1389 result = Py_False;
1390 goto cleanup;
1391 }
1392 }
1393
1394 if (d < 0) {
1395 d *= -1;
1396 d += 2;
1397 }
1398 else {
1399 d += 2;
1400 d *= -1;
1401 }
1402
1403 /* make sure we don't search forever */
1404 if (d >= max_d) {
1405 VALUE_ERROR("appropriate value for D cannot be found in is_strong_selfridge_prp()");
1406 goto cleanup;
1407 }
1408
1409 mpz_set_si(zD, d);
1410 }
1411
1412 q = (1-d)/4;
1413
1414 /* Since "O" is used, the refcount for n is incremented so deleting
1415 * temp will not delete n.
1416 */
1417 temp = Py_BuildValue("Oll", n, p, q);
1418 if (!temp)
1419 goto cleanup;
1420 result = GMPY_mpz_is_stronglucas_prp(NULL, temp);
1421 Py_DECREF(temp);
1422 goto return_result;
1423
1424 cleanup:
1425 Py_XINCREF(result);
1426 return_result:
1427 mpz_clear(zD);
1428 Py_DECREF((PyObject*)n);
1429 return result;
1430 }
1431
1432 /* **********************************************************************************
1433 * mpz_bpsw_prp:
1434 * A "Baillie-Pomerance-Selfridge-Wagstaff probable prime" is a composite n such that
1435 * n is a strong probable prime to the base 2 and
1436 * n is a Lucas probable prime using the Selfridge parameters.
1437 * **********************************************************************************/
1438
1439 PyDoc_STRVAR(doc_mpz_is_bpsw_prp,
1440 "is_bpsw_prp(n) -> boolean\n\n"
1441 "Return True if n is a Baillie-Pomerance-Selfridge-Wagstaff probable \n"
1442 "prime. A BPSW probable prime passes the is_strong_prp() test with base\n"
1443 "2 and the is_selfridge_prp() test.\n");
1444
1445 static PyObject *
GMPY_mpz_is_bpsw_prp(PyObject * self,PyObject * args)1446 GMPY_mpz_is_bpsw_prp(PyObject *self, PyObject *args)
1447 {
1448 MPZ_Object *n = NULL;
1449 PyObject *result = NULL, *temp = NULL;
1450
1451 if (PyTuple_Size(args) != 1) {
1452 TYPE_ERROR("is_bpsw_prp() requires 1 integer argument");
1453 return NULL;
1454 }
1455
1456 n = GMPy_MPZ_From_Integer(PyTuple_GET_ITEM(args, 0), NULL);
1457 if (!n) {
1458 TYPE_ERROR("is_bpsw_prp() requires 1 integer argument");
1459 goto cleanup;
1460 }
1461
1462 /* Require n > 0. */
1463 if (mpz_sgn(n->z) <= 0) {
1464 VALUE_ERROR("is_bpsw_prp() requires 'n' be greater than 0");
1465 goto cleanup;
1466 }
1467
1468 /* Check for n == 1 */
1469 if (mpz_cmp_ui(n->z, 1) == 0) {
1470 result = Py_False;
1471 goto cleanup;
1472 }
1473
1474 /* Handle n even. */
1475 if (mpz_divisible_ui_p(n->z, 2)) {
1476 if (mpz_cmp_ui(n->z, 2) == 0)
1477 result = Py_True;
1478 else
1479 result = Py_False;
1480 goto cleanup;
1481 }
1482
1483 /* "O" is used to increment the reference to n so deleting temp won't
1484 * delete n.
1485 */
1486 temp = Py_BuildValue("Oi", n, 2);
1487 if (!temp)
1488 goto cleanup;
1489 result = GMPY_mpz_is_strong_prp(NULL, temp);
1490 Py_DECREF(temp);
1491 if (result == Py_False)
1492 goto return_result;
1493 /* Remember to ignore the preceding result */
1494 Py_DECREF(result);
1495
1496 temp = Py_BuildValue("(O)", n);
1497 if (!temp)
1498 goto cleanup;
1499 result = GMPY_mpz_is_selfridge_prp(NULL, temp);
1500 Py_DECREF(temp);
1501 goto return_result;
1502
1503 cleanup:
1504 Py_XINCREF(result);
1505 return_result:
1506 Py_DECREF((PyObject*)n);
1507 return result;
1508 }
1509
1510
1511 /* ****************************************************************************************
1512 * mpz_strongbpsw_prp:
1513 * A "strong Baillie-Pomerance-Selfridge-Wagstaff probable prime" is a composite n such that
1514 * n is a strong probable prime to the base 2 and
1515 * n is a strong Lucas probable prime using the Selfridge parameters.
1516 * ****************************************************************************************/
1517
1518 PyDoc_STRVAR(doc_mpz_is_strongbpsw_prp,
1519 "is_strong_bpsw_prp(n) -> boolean\n\n"
1520 "Return True if n is a strong Baillie-Pomerance-Selfridge-Wagstaff\n"
1521 "probable prime. A strong BPSW probable prime passes the is_strong_prp()\n"
1522 "test with base and the is_strong_selfridge_prp() test.\n");
1523
1524 static PyObject *
GMPY_mpz_is_strongbpsw_prp(PyObject * self,PyObject * args)1525 GMPY_mpz_is_strongbpsw_prp(PyObject *self, PyObject *args)
1526 {
1527 MPZ_Object *n = NULL;
1528 PyObject *result = NULL, *temp = NULL;
1529
1530 if (PyTuple_Size(args) != 1) {
1531 TYPE_ERROR("is_strong_bpsw_prp() requires 1 integer argument");
1532 return NULL;
1533 }
1534
1535 n = GMPy_MPZ_From_Integer(PyTuple_GET_ITEM(args, 0), NULL);
1536 if (!n) {
1537 TYPE_ERROR("is_strong_bpsw_prp() requires 1 integer argument");
1538 goto cleanup;
1539 }
1540
1541 /* Require n > 0. */
1542 if (mpz_sgn(n->z) <= 0) {
1543 VALUE_ERROR("is_strong_bpsw_prp() requires 'n' be greater than 0");
1544 goto cleanup;
1545 }
1546
1547 /* Check for n == 1 */
1548 if (mpz_cmp_ui(n->z, 1) == 0) {
1549 result = Py_False;
1550 goto cleanup;
1551 }
1552
1553 /* Handle n even. */
1554 if (mpz_divisible_ui_p(n->z, 2)) {
1555 if (mpz_cmp_ui(n->z, 2) == 0)
1556 result = Py_True;
1557 else
1558 result = Py_False;
1559 goto cleanup;
1560 }
1561
1562 /* "O" is used to increment the reference to n so deleting temp won't
1563 * delete n.
1564 */
1565 temp = Py_BuildValue("Oi", n, 2);
1566 if (!temp)
1567 goto cleanup;
1568 result = GMPY_mpz_is_strong_prp(NULL, temp);
1569 Py_DECREF(temp);
1570 if (result == Py_False)
1571 goto return_result;
1572 /* Remember to ignore the preceding result */
1573 Py_DECREF(result);
1574
1575 temp = Py_BuildValue("(O)", n);
1576 if (!temp)
1577 goto cleanup;
1578 result = GMPY_mpz_is_selfridge_prp(NULL, temp);
1579 Py_DECREF(temp);
1580 goto return_result;
1581
1582 cleanup:
1583 Py_XINCREF(result);
1584 return_result:
1585 Py_DECREF((PyObject*)n);
1586 return result;
1587 }
1588
1589