From 65b2c999c088ce5ae28990797b4a80b72d6df6a2 Mon Sep 17 00:00:00 2001 From: David Cleaver Date: Sat, 31 Jan 2015 15:23:14 +0000 Subject: [PATCH] Remove unused code from aprcl and clean-up usage of aprcl in ecm git-svn-id: svn://scm.gforge.inria.fr/svnroot/ecm/trunk@2571 404564d9-a503-0410-82bf-e18ce2cf3989 --- aprtcle/aprcl.c | 52 +- aprtcle/mpz_aprcl.c | 1183 +------------------------------------------ aprtcle/mpz_aprcl.h | 109 +--- auxi.c | 29 -- candi.c | 18 +- ecm-ecm.h | 20 + 6 files changed, 83 insertions(+), 1328 deletions(-) diff --git a/aprtcle/aprcl.c b/aprtcle/aprcl.c index 9efad4ce..9cab38cf 100644 --- a/aprtcle/aprcl.c +++ b/aprtcle/aprcl.c @@ -1,4 +1,4 @@ -/* Copyright 2012,2013 David Cleaver +/* Copyright 2012-2015 David Cleaver * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by @@ -83,7 +83,7 @@ int main(int argc, char* argv[]) printf(" Where is a number to test for primality\n"); printf(" %s -inp \n", argv[0]); printf(" Where contains one or more numbers to test for primality\n"); - printf(" Note: This program returns the APR-CL and BPSW status of the input number(s).\n"); + printf(" Note: This program returns the APR-CL and MPZ PRP status of the input number(s).\n"); printf(" Note: This program will also print out timing information for each check.\n"); printf(" Note: All input numbers are assumed to be base-10 numbers.\n"); return 0; @@ -124,9 +124,9 @@ int main(int argc, char* argv[]) printf("===============================================================================\n"); dig = b10_digits(test); gmp_printf("Testing: %Zd (%d digits)\n", test, dig); fflush(stdout); - printf("Running the BPSW PRP test...\n"); fflush(stdout); + printf("Running the MPZ PRP test with 5 iterations...\n"); fflush(stdout); t0 = ttime(0); - ret_b = mpz_bpsw_prp(test); + ret_b = mpz_probab_prime_p(test, 5); t1 = ttime(t0); printf("Running the APRCL prime test...\n"); fflush(stdout); @@ -134,13 +134,13 @@ int main(int argc, char* argv[]) ret_a = mpz_aprtcle(test, 1); t2 = ttime(t0); - printf("\n BPSW took %.4f seconds\n", t1); - if (ret_b == PRP_COMPOSITE) - printf(" BPSW says the number is COMPOSITE\n"); - else if (ret_b == PRP_PRP) - printf(" BPSW says the number is PRP\n"); - else if (ret_b == PRP_PRIME) - printf(" BPSW says the number is PRIME\n"); + printf("\n MPZ PRP took %.4f seconds\n", t1); + if (ret_b == APRTCLE_COMPOSITE) + printf(" MPZ PRP says the number is COMPOSITE\n"); + else if (ret_b == APRTCLE_PRP) + printf(" MPZ PRP says the number is PRP\n"); + else if (ret_b == APRTCLE_PRIME) + printf(" MPZ PRP says the number is PRIME\n"); printf("APRCL took %.4f seconds\n", t2); if (ret_a == APRTCLE_COMPOSITE) @@ -150,11 +150,11 @@ int main(int argc, char* argv[]) else if (ret_a == APRTCLE_PRIME) printf("APRCL says the number is PRIME\n"); - if ((ret_b == PRP_COMPOSITE && ret_a != APRTCLE_COMPOSITE) || - (ret_b != PRP_COMPOSITE && ret_a == APRTCLE_COMPOSITE)) + if ((ret_b == APRTCLE_COMPOSITE && ret_a != APRTCLE_COMPOSITE) || + (ret_b != APRTCLE_COMPOSITE && ret_a == APRTCLE_COMPOSITE)) { printf(" *** ATTENTION *** ATTENTION *** ATTENTION *** ATTENTION ***\n"); - printf("BPSW and APRCL do not agree on the status of this number!!!\n"); + printf("MPZ PRP and APRCL do not agree on the status of this number!!!\n"); printf("Please report this to http://www.mersenneforum.org/showthread.php?t=18353\n"); gmp_printf("N = %Zd\n", test); } @@ -172,9 +172,9 @@ int main(int argc, char* argv[]) dig = b10_digits(test); gmp_printf("Testing: %Zd (%d digits)\n", test, dig); fflush(stdout); - printf("Running the BPSW PRP test...\n"); fflush(stdout); + printf("Running the MPZ PRP test with 5 iterations...\n"); fflush(stdout); t0 = ttime(0); - ret_b = mpz_bpsw_prp(test); + ret_b = mpz_probab_prime_p(test, 5); t1 = ttime(t0); printf("Running the APRCL prime test...\n"); fflush(stdout); @@ -182,13 +182,13 @@ int main(int argc, char* argv[]) ret_a = mpz_aprtcle(test, 1); t2 = ttime(t0); - printf("\n BPSW took %.4f seconds\n", t1); - if (ret_b == PRP_COMPOSITE) - printf(" BPSW says the number is COMPOSITE\n"); - else if (ret_b == PRP_PRP) - printf(" BPSW says the number is PRP\n"); - else if (ret_b == PRP_PRIME) - printf(" BPSW says the number is PRIME\n"); + printf("\n MPZ PRP took %.4f seconds\n", t1); + if (ret_b == APRTCLE_COMPOSITE) + printf(" MPZ PRP says the number is COMPOSITE\n"); + else if (ret_b == APRTCLE_PRP) + printf(" MPZ PRP says the number is PRP\n"); + else if (ret_b == APRTCLE_PRIME) + printf(" MPZ PRP says the number is PRIME\n"); printf("APRCL took %.4f seconds\n", t2); if (ret_a == APRTCLE_COMPOSITE) @@ -198,11 +198,11 @@ int main(int argc, char* argv[]) else if (ret_a == APRTCLE_PRIME) printf("APRCL says the number is PRIME\n"); - if ((ret_b == PRP_COMPOSITE && ret_a != APRTCLE_COMPOSITE) || - (ret_b != PRP_COMPOSITE && ret_a == APRTCLE_COMPOSITE)) + if ((ret_b == APRTCLE_COMPOSITE && ret_a != APRTCLE_COMPOSITE) || + (ret_b != APRTCLE_COMPOSITE && ret_a == APRTCLE_COMPOSITE)) { printf(" *** ATTENTION *** ATTENTION *** ATTENTION *** ATTENTION ***\n"); - printf("BPSW and APRCL do not agree on the status of this number!!!\n"); + printf("MPZ PRP and APRCL do not agree on the status of this number!!!\n"); printf("Please report this to http://www.mersenneforum.org/showthread.php?t=18353\n"); gmp_printf("N = %Zd\n", test); } diff --git a/aprtcle/mpz_aprcl.c b/aprtcle/mpz_aprcl.c index 4f7e976d..65e70411 100644 --- a/aprtcle/mpz_aprcl.c +++ b/aprtcle/mpz_aprcl.c @@ -1,4 +1,4 @@ -/* Copyright 2011,2012,2013 David Cleaver +/* Copyright 2011-2015 David Cleaver * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by @@ -18,8 +18,6 @@ * v1.0 Posted to SourceForge on 2013/07/04 * * v1.1 Posted to SourceForge on 2013/12/27 - * [The following fix was recommended by Dana Jacobsen and verified by Jon Grantham] - * - Bug fix: Removed unnecessary vl==0 check in mpz_extrastronglucas_prp * [The following improvements/fixes were recommended by Laurent Desnogues in 2013/08] * - Speed improvement 1: Removed extraneous NormalizeJS calls in ARPCL * - Speed improvement 2: Removed/consolidated calls to mpz_mod in APRCL @@ -27,10 +25,6 @@ * - Bug fix: Final test in APRCL routine is now correct */ -/* - * The PRP functions presented here are based on the paper: - * Grantham, Jon. Frobenius Pseudoprimes. Math. Comp. 70 (2001), 873-891. - */ #include #include @@ -43,1148 +37,6 @@ typedef long long s64_t; typedef unsigned long long u64_t; #endif -/* ****************************************************************** - * mpz_prp: (also called a Fermat pseudoprime) - * A "pseudoprime" to the base a is a composite number n such that, - * (a,n)=1 and a^(n-1) = 1 mod n - * ******************************************************************/ -int mpz_prp(mpz_t n, mpz_t a) -{ - mpz_t res; - mpz_t nm1; - - if (mpz_cmp_ui(a, 2) < 0) - return PRP_ERROR; - - if (mpz_cmp_ui(n, 2) < 0) - return PRP_COMPOSITE; - - if (mpz_divisible_ui_p(n, 2)) - { - if (mpz_cmp_ui(n, 2) == 0) - return PRP_PRIME; - else - return PRP_COMPOSITE; - } - - mpz_init_set_ui(res, 0); - mpz_gcd(res, n, a); - - if (mpz_cmp_ui(res, 1) > 0) - { - mpz_clear(res); - return PRP_COMPOSITE; - } - - mpz_init_set(nm1, n); - mpz_sub_ui(nm1, nm1, 1); - mpz_powm(res, a, nm1, n); - - if (mpz_cmp_ui(res, 1) == 0) - { - mpz_clear(res); - mpz_clear(nm1); - return PRP_PRP; - } - else - { - mpz_clear(res); - mpz_clear(nm1); - return PRP_COMPOSITE; - } - -}/* method mpz_prp */ - - -/* ************************************************************************* - * mpz_euler_prp: (also called a Solovay-Strassen pseudoprime) - * An "Euler pseudoprime" to the base a is an odd composite number n with, - * (a,n)=1 such that a^((n-1)/2)=(a/n) mod n [(a/n) is the Jacobi symbol] - * *************************************************************************/ -int mpz_euler_prp(mpz_t n, mpz_t a) -{ - mpz_t res; - mpz_t exp; - int ret = 0; - - if (mpz_cmp_ui(a, 2) < 0) - return PRP_ERROR; - - if (mpz_cmp_ui(n, 2) < 0) - return PRP_COMPOSITE; - - if (mpz_divisible_ui_p(n, 2)) - { - if (mpz_cmp_ui(n, 2) == 0) - return PRP_PRIME; - else - return PRP_COMPOSITE; - } - - mpz_init_set_ui(res, 0); - mpz_gcd(res, n, a); - - if (mpz_cmp_ui(res, 1) > 0) - { - mpz_clear(res); - return PRP_COMPOSITE; - } - - mpz_init_set(exp, n); - - mpz_sub_ui(exp, exp, 1); /* exp = n-1 */ - mpz_divexact_ui(exp, exp, 2); /* exp = (n-1)/2 */ - mpz_powm(res, a, exp, n); - - /* reuse exp to calculate jacobi(a,n) mod n */ - ret = mpz_jacobi(a,n); - mpz_set(exp, n); - if (ret == -1) - mpz_sub_ui(exp, exp, 1); - else if (ret == 1) - mpz_add_ui(exp, exp, 1); - mpz_mod(exp, exp, n); - - if (mpz_cmp(res, exp) == 0) - { - mpz_clear(res); - mpz_clear(exp); - return PRP_PRP; - } - else - { - mpz_clear(res); - mpz_clear(exp); - return PRP_COMPOSITE; - } - -}/* method mpz_euler_prp */ - - -/* ********************************************************************************************* - * mpz_sprp: (also called a Miller-Rabin pseudoprime) - * A "strong pseudoprime" to the base a is an odd composite n = (2^r)*s+1 with s odd such that - * either a^s == 1 mod n, or a^((2^t)*s) == -1 mod n, for some integer t, with 0 <= t < r. - * *********************************************************************************************/ -int mpz_sprp(mpz_t n, mpz_t a) -{ - mpz_t s; - mpz_t nm1; - mpz_t mpz_test; - unsigned long int r = 0; - - if (mpz_cmp_ui(a, 2) < 0) - return PRP_ERROR; - - if (mpz_cmp_ui(n, 2) < 0) - return PRP_COMPOSITE; - - if (mpz_divisible_ui_p(n, 2)) - { - if (mpz_cmp_ui(n, 2) == 0) - return PRP_PRIME; - else - return PRP_COMPOSITE; - } - - mpz_init_set_ui(mpz_test, 0); - mpz_init_set_ui(s, 0); - mpz_init_set(nm1, n); - mpz_sub_ui(nm1, nm1, 1); - - /***********************************************/ - /* Find s and r satisfying: n-1=(2^r)*s, s odd */ - r = mpz_scan1(nm1, 0); - mpz_fdiv_q_2exp(s, nm1, r); - - - /******************************************/ - /* Check a^((2^t)*s) mod n for 0 <= t < r */ - mpz_powm(mpz_test, a, s, n); - if ( (mpz_cmp_ui(mpz_test, 1) == 0) || (mpz_cmp(mpz_test, nm1) == 0) ) - { - mpz_clear(s); - mpz_clear(nm1); - mpz_clear(mpz_test); - return PRP_PRP; - } - - while ( --r ) - { - /* mpz_test = mpz_test^2%n */ - mpz_mul(mpz_test, mpz_test, mpz_test); - mpz_mod(mpz_test, mpz_test, n); - - if (mpz_cmp(mpz_test, nm1) == 0) - { - mpz_clear(s); - mpz_clear(nm1); - mpz_clear(mpz_test); - return PRP_PRP; - } - } - - mpz_clear(s); - mpz_clear(nm1); - mpz_clear(mpz_test); - return PRP_COMPOSITE; - -}/* method mpz_sprp */ - - -/* ************************************************************************* - * mpz_fibonacci_prp: - * A "Fibonacci pseudoprime" with parameters (P,Q), P > 0, Q=+/-1, is a - * composite n for which V_n == P mod n - * [V is the Lucas V sequence with parameters P,Q] - * *************************************************************************/ -int mpz_fibonacci_prp(mpz_t n, long int p, long int q) -{ - mpz_t pmodn, zP; - mpz_t vl, vh, ql, qh, tmp; /* used for calculating the Lucas V sequence */ - int s = 0, j = 0; - - if (p*p-4*q == 0) - return PRP_ERROR; - - if (((q != 1) && (q != -1)) || (p <= 0)) - return PRP_ERROR; - - if (mpz_cmp_ui(n, 2) < 0) - return PRP_COMPOSITE; - - if (mpz_divisible_ui_p(n, 2)) - { - if (mpz_cmp_ui(n, 2) == 0) - return PRP_PRIME; - else - return PRP_COMPOSITE; - } - - mpz_init_set_ui(zP, p); - mpz_init(pmodn); - mpz_mod(pmodn, zP, n); - - /* mpz_lucasvmod(res, p, q, n, n); */ - mpz_init_set_si(vl, 2); - mpz_init_set_si(vh, p); - mpz_init_set_si(ql, 1); - mpz_init_set_si(qh, 1); - mpz_init_set_si(tmp,0); - - s = mpz_scan1(n, 0); - for (j = mpz_sizeinbase(n,2)-1; j >= s+1; j--) - { - /* ql = ql*qh (mod n) */ - mpz_mul(ql, ql, qh); - mpz_mod(ql, ql, n); - if (mpz_tstbit(n,j) == 1) - { - /* qh = ql*q */ - mpz_mul_si(qh, ql, q); - - /* vl = vh*vl - p*ql (mod n) */ - mpz_mul(vl, vh, vl); - mpz_mul_si(tmp, ql, p); - mpz_sub(vl, vl, tmp); - mpz_mod(vl, vl, n); - - /* vh = vh*vh - 2*qh (mod n) */ - mpz_mul(vh, vh, vh); - mpz_mul_si(tmp, qh, 2); - mpz_sub(vh, vh, tmp); - mpz_mod(vh, vh, n); - } - else - { - /* qh = ql */ - mpz_set(qh, ql); - - /* vh = vh*vl - p*ql (mod n) */ - mpz_mul(vh, vh, vl); - mpz_mul_si(tmp, ql, p); - mpz_sub(vh, vh, tmp); - mpz_mod(vh, vh, n); - - /* vl = vl*vl - 2*ql (mod n) */ - mpz_mul(vl, vl, vl); - mpz_mul_si(tmp, ql, 2); - mpz_sub(vl, vl, tmp); - mpz_mod(vl, vl, n); - } - } - /* ql = ql*qh */ - mpz_mul(ql, ql, qh); - - /* qh = ql*q */ - mpz_mul_si(qh, ql, q); - - /* vl = vh*vl - p*ql */ - mpz_mul(vl, vh, vl); - mpz_mul_si(tmp, ql, p); - mpz_sub(vl, vl, tmp); - - /* ql = ql*qh */ - mpz_mul(ql, ql, qh); - - for (j = 1; j <= s; j++) - { - /* vl = vl*vl - 2*ql (mod n) */ - mpz_mul(vl, vl, vl); - mpz_mul_si(tmp, ql, 2); - mpz_sub(vl, vl, tmp); - mpz_mod(vl, vl, n); - - /* ql = ql*ql (mod n) */ - mpz_mul(ql, ql, ql); - mpz_mod(ql, ql, n); - } - - mpz_mod(vl, vl, n); /* vl contains our return value */ - - if (mpz_cmp(vl, pmodn) == 0) - { - mpz_clear(zP); - mpz_clear(pmodn); - mpz_clear(vl); - mpz_clear(vh); - mpz_clear(ql); - mpz_clear(qh); - mpz_clear(tmp); - return PRP_PRP; - } - mpz_clear(zP); - mpz_clear(pmodn); - mpz_clear(vl); - mpz_clear(vh); - mpz_clear(ql); - mpz_clear(qh); - mpz_clear(tmp); - return PRP_COMPOSITE; - -}/* method mpz_fibonacci_prp */ - - -/* ******************************************************************************* - * mpz_lucas_prp: - * A "Lucas pseudoprime" with parameters (P,Q) is a composite n with D=P^2-4Q, - * (n,2QD)=1 such that U_(n-(D/n)) == 0 mod n [(D/n) is the Jacobi symbol] - * *******************************************************************************/ -int mpz_lucas_prp(mpz_t n, long int p, long int q) -{ - mpz_t zD; - mpz_t res; - mpz_t index; - mpz_t uh, vl, vh, ql, qh, tmp; /* used for calculating the Lucas U sequence */ - int s = 0, j = 0; - int ret = 0; - long int d = p*p - 4*q; - - if (d == 0) /* Does not produce a proper Lucas sequence */ - return PRP_ERROR; - - if (mpz_cmp_ui(n, 2) < 0) - return PRP_COMPOSITE; - - if (mpz_divisible_ui_p(n, 2)) - { - if (mpz_cmp_ui(n, 2) == 0) - return PRP_PRIME; - else - return PRP_COMPOSITE; - } - - mpz_init(index); - mpz_init_set_si(zD, d); - mpz_init(res); - - mpz_mul_si(res, zD, q); - mpz_mul_ui(res, res, 2); - mpz_gcd(res, res, n); - if ((mpz_cmp(res, n) != 0) && (mpz_cmp_ui(res, 1) > 0)) - { - mpz_clear(zD); - mpz_clear(res); - mpz_clear(index); - return PRP_COMPOSITE; - } - - /* index = n-(D/n), where (D/n) is the Jacobi symbol */ - mpz_set(index, n); - ret = mpz_jacobi(zD, n); - if (ret == -1) - mpz_add_ui(index, index, 1); - else if (ret == 1) - mpz_sub_ui(index, index, 1); - - /* mpz_lucasumod(res, p, q, index, n); */ - mpz_init_set_si(uh, 1); - mpz_init_set_si(vl, 2); - mpz_init_set_si(vh, p); - mpz_init_set_si(ql, 1); - mpz_init_set_si(qh, 1); - mpz_init_set_si(tmp,0); - - s = mpz_scan1(index, 0); - for (j = mpz_sizeinbase(index,2)-1; j >= s+1; j--) - { - /* ql = ql*qh (mod n) */ - mpz_mul(ql, ql, qh); - mpz_mod(ql, ql, n); - if (mpz_tstbit(index,j) == 1) - { - /* qh = ql*q */ - mpz_mul_si(qh, ql, q); - - /* uh = uh*vh (mod n) */ - mpz_mul(uh, uh, vh); - mpz_mod(uh, uh, n); - - /* vl = vh*vl - p*ql (mod n) */ - mpz_mul(vl, vh, vl); - mpz_mul_si(tmp, ql, p); - mpz_sub(vl, vl, tmp); - mpz_mod(vl, vl, n); - - /* vh = vh*vh - 2*qh (mod n) */ - mpz_mul(vh, vh, vh); - mpz_mul_si(tmp, qh, 2); - mpz_sub(vh, vh, tmp); - mpz_mod(vh, vh, n); - } - else - { - /* qh = ql */ - mpz_set(qh, ql); - - /* uh = uh*vl - ql (mod n) */ - mpz_mul(uh, uh, vl); - mpz_sub(uh, uh, ql); - mpz_mod(uh, uh, n); - - /* vh = vh*vl - p*ql (mod n) */ - mpz_mul(vh, vh, vl); - mpz_mul_si(tmp, ql, p); - mpz_sub(vh, vh, tmp); - mpz_mod(vh, vh, n); - - /* vl = vl*vl - 2*ql (mod n) */ - mpz_mul(vl, vl, vl); - mpz_mul_si(tmp, ql, 2); - mpz_sub(vl, vl, tmp); - mpz_mod(vl, vl, n); - } - } - /* ql = ql*qh */ - mpz_mul(ql, ql, qh); - - /* qh = ql*q */ - mpz_mul_si(qh, ql, q); - - /* uh = uh*vl - ql */ - mpz_mul(uh, uh, vl); - mpz_sub(uh, uh, ql); - - /* vl = vh*vl - p*ql */ - mpz_mul(vl, vh, vl); - mpz_mul_si(tmp, ql, p); - mpz_sub(vl, vl, tmp); - - /* ql = ql*qh */ - mpz_mul(ql, ql, qh); - - for (j = 1; j <= s; j++) - { - /* uh = uh*vl (mod n) */ - mpz_mul(uh, uh, vl); - mpz_mod(uh, uh, n); - - /* vl = vl*vl - 2*ql (mod n) */ - mpz_mul(vl, vl, vl); - mpz_mul_si(tmp, ql, 2); - mpz_sub(vl, vl, tmp); - mpz_mod(vl, vl, n); - - /* ql = ql*ql (mod n) */ - mpz_mul(ql, ql, ql); - mpz_mod(ql, ql, n); - } - - mpz_mod(res, uh, n); /* uh contains our return value */ - - mpz_clear(zD); - mpz_clear(index); - mpz_clear(uh); - mpz_clear(vl); - mpz_clear(vh); - mpz_clear(ql); - mpz_clear(qh); - mpz_clear(tmp); - - if (mpz_cmp_ui(res, 0) == 0) - { - mpz_clear(res); - return PRP_PRP; - } - else - { - mpz_clear(res); - return PRP_COMPOSITE; - } - -}/* method mpz_lucas_prp */ - - -/* ********************************************************************************************* - * mpz_stronglucas_prp: - * A "strong Lucas pseudoprime" with parameters (P,Q) is a composite n = (2^r)*s+(D/n), where - * 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 - * for some t, 0 <= t < r. [(D/n) is the Jacobi symbol] - * *********************************************************************************************/ -int mpz_stronglucas_prp(mpz_t n, long int p, long int q) -{ - mpz_t zD; - mpz_t s; - mpz_t nmj; /* n minus jacobi(D/n) */ - mpz_t res; - mpz_t uh, vl, vh, ql, qh, tmp; /* these are needed for the LucasU and LucasV part of this function */ - long int d = p*p - 4*q; - unsigned long int r = 0; - int ret = 0; - unsigned long int j = 0; - - if (d == 0) /* Does not produce a proper Lucas sequence */ - return PRP_ERROR; - - if (mpz_cmp_ui(n, 2) < 0) - return PRP_COMPOSITE; - - if (mpz_divisible_ui_p(n, 2)) - { - if (mpz_cmp_ui(n, 2) == 0) - return PRP_PRIME; - else - return PRP_COMPOSITE; - } - - mpz_init_set_si(zD, d); - mpz_init(res); - - mpz_mul_si(res, zD, q); - mpz_mul_ui(res, res, 2); - mpz_gcd(res, res, n); - if ((mpz_cmp(res, n) != 0) && (mpz_cmp_ui(res, 1) > 0)) - { - mpz_clear(zD); - mpz_clear(res); - return PRP_COMPOSITE; - } - - mpz_init(s); - mpz_init(nmj); - - /* nmj = n - (D/n), where (D/n) is the Jacobi symbol */ - mpz_set(nmj, n); - ret = mpz_jacobi(zD, n); - if (ret == -1) - mpz_add_ui(nmj, nmj, 1); - else if (ret == 1) - mpz_sub_ui(nmj, nmj, 1); - - r = mpz_scan1(nmj, 0); - mpz_fdiv_q_2exp(s, nmj, r); - - /* make sure U_s == 0 mod n or V_((2^t)*s) == 0 mod n, for some t, 0 <= t < r */ - mpz_init_set_si(uh, 1); - mpz_init_set_si(vl, 2); - mpz_init_set_si(vh, p); - mpz_init_set_si(ql, 1); - mpz_init_set_si(qh, 1); - mpz_init_set_si(tmp,0); - - for (j = mpz_sizeinbase(s,2)-1; j >= 1; j--) - { - /* ql = ql*qh (mod n) */ - mpz_mul(ql, ql, qh); - mpz_mod(ql, ql, n); - if (mpz_tstbit(s,j) == 1) - { - /* qh = ql*q */ - mpz_mul_si(qh, ql, q); - - /* uh = uh*vh (mod n) */ - mpz_mul(uh, uh, vh); - mpz_mod(uh, uh, n); - - /* vl = vh*vl - p*ql (mod n) */ - mpz_mul(vl, vh, vl); - mpz_mul_si(tmp, ql, p); - mpz_sub(vl, vl, tmp); - mpz_mod(vl, vl, n); - - /* vh = vh*vh - 2*qh (mod n) */ - mpz_mul(vh, vh, vh); - mpz_mul_si(tmp, qh, 2); - mpz_sub(vh, vh, tmp); - mpz_mod(vh, vh, n); - } - else - { - /* qh = ql */ - mpz_set(qh, ql); - - /* uh = uh*vl - ql (mod n) */ - mpz_mul(uh, uh, vl); - mpz_sub(uh, uh, ql); - mpz_mod(uh, uh, n); - - /* vh = vh*vl - p*ql (mod n) */ - mpz_mul(vh, vh, vl); - mpz_mul_si(tmp, ql, p); - mpz_sub(vh, vh, tmp); - mpz_mod(vh, vh, n); - - /* vl = vl*vl - 2*ql (mod n) */ - mpz_mul(vl, vl, vl); - mpz_mul_si(tmp, ql, 2); - mpz_sub(vl, vl, tmp); - mpz_mod(vl, vl, n); - } - } - /* ql = ql*qh */ - mpz_mul(ql, ql, qh); - - /* qh = ql*q */ - mpz_mul_si(qh, ql, q); - - /* uh = uh*vl - ql */ - mpz_mul(uh, uh, vl); - mpz_sub(uh, uh, ql); - - /* vl = vh*vl - p*ql */ - mpz_mul(vl, vh, vl); - mpz_mul_si(tmp, ql, p); - mpz_sub(vl, vl, tmp); - - /* ql = ql*qh */ - mpz_mul(ql, ql, qh); - - mpz_mod(uh, uh, n); - mpz_mod(vl, vl, n); - - /* uh contains LucasU_s and vl contains LucasV_s */ - if ((mpz_cmp_ui(uh, 0) == 0) || (mpz_cmp_ui(vl, 0) == 0)) - { - mpz_clear(zD); - mpz_clear(s); - mpz_clear(nmj); - mpz_clear(res); - mpz_clear(uh); - mpz_clear(vl); - mpz_clear(vh); - mpz_clear(ql); - mpz_clear(qh); - mpz_clear(tmp); - return PRP_PRP; - } - - for (j = 1; j < r; j++) - { - /* vl = vl*vl - 2*ql (mod n) */ - mpz_mul(vl, vl, vl); - mpz_mul_si(tmp, ql, 2); - mpz_sub(vl, vl, tmp); - mpz_mod(vl, vl, n); - - /* ql = ql*ql (mod n) */ - mpz_mul(ql, ql, ql); - mpz_mod(ql, ql, n); - - if (mpz_cmp_ui(vl, 0) == 0) - { - mpz_clear(zD); - mpz_clear(s); - mpz_clear(nmj); - mpz_clear(res); - mpz_clear(uh); - mpz_clear(vl); - mpz_clear(vh); - mpz_clear(ql); - mpz_clear(qh); - mpz_clear(tmp); - return PRP_PRP; - } - } - - mpz_clear(zD); - mpz_clear(s); - mpz_clear(nmj); - mpz_clear(res); - mpz_clear(uh); - mpz_clear(vl); - mpz_clear(vh); - mpz_clear(ql); - mpz_clear(qh); - mpz_clear(tmp); - return PRP_COMPOSITE; - -}/* method mpz_stronglucas_prp */ - - -/* ******************************************************************************************* - * mpz_extrastronglucas_prp: - * Let U_n = LucasU(p,1), V_n = LucasV(p,1), and D=p^2-4. - * An "extra strong Lucas pseudoprime" to the base p is a composite n = (2^r)*s+(D/n), where - * s is odd and (n,2D)=1, such that either U_s == 0 mod n and V_s == +/-2 mod n, or - * V_((2^t)*s) == 0 mod n for some t with 0 <= t < r-1 [(D/n) is the Jacobi symbol] - * *******************************************************************************************/ -int mpz_extrastronglucas_prp(mpz_t n, long int p) -{ - mpz_t zD; - mpz_t s; - mpz_t nmj; /* n minus jacobi(D/n) */ - mpz_t res; - mpz_t uh, vl, vh, ql, qh, tmp; /* these are needed for the LucasU and LucasV part of this function */ - long int d = p*p - 4; - long int q = 1; - unsigned long int r = 0; - int ret = 0; - unsigned long int j = 0; - - if (d == 0) /* Does not produce a proper Lucas sequence */ - return PRP_ERROR; - - if (mpz_cmp_ui(n, 2) < 0) - return PRP_COMPOSITE; - - if (mpz_divisible_ui_p(n, 2)) - { - if (mpz_cmp_ui(n, 2) == 0) - return PRP_PRIME; - else - return PRP_COMPOSITE; - } - - mpz_init_set_si(zD, d); - mpz_init(res); - - mpz_mul_ui(res, zD, 2); - mpz_gcd(res, res, n); - if ((mpz_cmp(res, n) != 0) && (mpz_cmp_ui(res, 1) > 0)) - { - mpz_clear(zD); - mpz_clear(res); - return PRP_COMPOSITE; - } - - mpz_init(s); - mpz_init(nmj); - - /* nmj = n - (D/n), where (D/n) is the Jacobi symbol */ - mpz_set(nmj, n); - ret = mpz_jacobi(zD, n); - if (ret == -1) - mpz_add_ui(nmj, nmj, 1); - else if (ret == 1) - mpz_sub_ui(nmj, nmj, 1); - - r = mpz_scan1(nmj, 0); - mpz_fdiv_q_2exp(s, nmj, r); - - /* make sure that either (U_s == 0 mod n and V_s == +/-2 mod n), or */ - /* V_((2^t)*s) == 0 mod n for some t with 0 <= t < r-1 */ - mpz_init_set_si(uh, 1); - mpz_init_set_si(vl, 2); - mpz_init_set_si(vh, p); - mpz_init_set_si(ql, 1); - mpz_init_set_si(qh, 1); - mpz_init_set_si(tmp,0); - - for (j = mpz_sizeinbase(s,2)-1; j >= 1; j--) - { - /* ql = ql*qh (mod n) */ - mpz_mul(ql, ql, qh); - mpz_mod(ql, ql, n); - if (mpz_tstbit(s,j) == 1) - { - /* qh = ql*q */ - mpz_mul_si(qh, ql, q); - - /* uh = uh*vh (mod n) */ - mpz_mul(uh, uh, vh); - mpz_mod(uh, uh, n); - - /* vl = vh*vl - p*ql (mod n) */ - mpz_mul(vl, vh, vl); - mpz_mul_si(tmp, ql, p); - mpz_sub(vl, vl, tmp); - mpz_mod(vl, vl, n); - - /* vh = vh*vh - 2*qh (mod n) */ - mpz_mul(vh, vh, vh); - mpz_mul_si(tmp, qh, 2); - mpz_sub(vh, vh, tmp); - mpz_mod(vh, vh, n); - } - else - { - /* qh = ql */ - mpz_set(qh, ql); - - /* uh = uh*vl - ql (mod n) */ - mpz_mul(uh, uh, vl); - mpz_sub(uh, uh, ql); - mpz_mod(uh, uh, n); - - /* vh = vh*vl - p*ql (mod n) */ - mpz_mul(vh, vh, vl); - mpz_mul_si(tmp, ql, p); - mpz_sub(vh, vh, tmp); - mpz_mod(vh, vh, n); - - /* vl = vl*vl - 2*ql (mod n) */ - mpz_mul(vl, vl, vl); - mpz_mul_si(tmp, ql, 2); - mpz_sub(vl, vl, tmp); - mpz_mod(vl, vl, n); - } - } - /* ql = ql*qh */ - mpz_mul(ql, ql, qh); - - /* qh = ql*q */ - mpz_mul_si(qh, ql, q); - - /* uh = uh*vl - ql */ - mpz_mul(uh, uh, vl); - mpz_sub(uh, uh, ql); - - /* vl = vh*vl - p*ql */ - mpz_mul(vl, vh, vl); - mpz_mul_si(tmp, ql, p); - mpz_sub(vl, vl, tmp); - - /* ql = ql*qh */ - mpz_mul(ql, ql, qh); - - mpz_mod(uh, uh, n); - mpz_mod(vl, vl, n); - - /* tmp = n-2, for the following comparison */ - mpz_sub_ui(tmp, n, 2); - - /* uh contains LucasU_s and vl contains LucasV_s */ - if ((mpz_cmp_ui(uh, 0) == 0) && ((mpz_cmp(vl, tmp) == 0) || (mpz_cmp_si(vl, 2) == 0))) - { - mpz_clear(zD); - mpz_clear(s); - mpz_clear(nmj); - mpz_clear(res); - mpz_clear(uh); - mpz_clear(vl); - mpz_clear(vh); - mpz_clear(ql); - mpz_clear(qh); - mpz_clear(tmp); - return PRP_PRP; - } - - for (j = 1; j < r-1; j++) - { - /* vl = vl*vl - 2*ql (mod n) */ - mpz_mul(vl, vl, vl); - mpz_mul_si(tmp, ql, 2); - mpz_sub(vl, vl, tmp); - mpz_mod(vl, vl, n); - - /* ql = ql*ql (mod n) */ - mpz_mul(ql, ql, ql); - mpz_mod(ql, ql, n); - - if (mpz_cmp_ui(vl, 0) == 0) - { - mpz_clear(zD); - mpz_clear(s); - mpz_clear(nmj); - mpz_clear(res); - mpz_clear(uh); - mpz_clear(vl); - mpz_clear(vh); - mpz_clear(ql); - mpz_clear(qh); - mpz_clear(tmp); - return PRP_PRP; - } - } - - mpz_clear(zD); - mpz_clear(s); - mpz_clear(nmj); - mpz_clear(res); - mpz_clear(uh); - mpz_clear(vl); - mpz_clear(vh); - mpz_clear(ql); - mpz_clear(qh); - mpz_clear(tmp); - return PRP_COMPOSITE; - -}/* method mpz_extrastronglucas_prp */ - - -/* *********************************************************************************************** - * mpz_selfridge_prp: - * A "Lucas-Selfridge pseudoprime" n is a "Lucas pseudoprime" using Selfridge parameters of: - * Find the first element D in the sequence {5, -7, 9, -11, 13, ...} such that Jacobi(D,n) = -1 - * Then use P=1 and Q=(1-D)/4 in the Lucas pseudoprime test. - * Make sure n is not a perfect square, otherwise the search for D will only stop when D=n. - * ***********************************************************************************************/ -int mpz_selfridge_prp(mpz_t n) -{ - long int d = 5, p = 1, q = 0; - int max_d = 1000000; - int jacobi = 0; - mpz_t zD; - - if (mpz_cmp_ui(n, 2) < 0) - return PRP_COMPOSITE; - - if (mpz_divisible_ui_p(n, 2)) - { - if (mpz_cmp_ui(n, 2) == 0) - return PRP_PRIME; - else - return PRP_COMPOSITE; - } - - mpz_init_set_ui(zD, d); - - while (1) - { - jacobi = mpz_jacobi(zD, n); - - /* if jacobi == 0, d is a factor of n, therefore n is composite... */ - /* if d == n, then either n is either prime or 9... */ - if (jacobi == 0) - { - if ((mpz_cmpabs(zD, n) == 0) && (mpz_cmp_ui(zD, 9) != 0)) - { - mpz_clear(zD); - return PRP_PRIME; - } - else - { - mpz_clear(zD); - return PRP_COMPOSITE; - } - } - if (jacobi == -1) - break; - - /* if we get to the 5th d, make sure we aren't dealing with a square... */ - if (d == 13) - { - if (mpz_perfect_square_p(n)) - { - mpz_clear(zD); - return PRP_COMPOSITE; - } - } - - if (d < 0) - { - d *= -1; - d += 2; - } - else - { - d += 2; - d *= -1; - } - - /* make sure we don't search forever */ - if (d >= max_d) - { - mpz_clear(zD); - return PRP_ERROR; - } - - mpz_set_si(zD, d); - } - mpz_clear(zD); - - q = (1-d)/4; - - return mpz_lucas_prp(n, p, q); - -}/* method mpz_selfridge_prp */ - - -/* ********************************************************************************************************* - * mpz_strongselfridge_prp: - * A "strong Lucas-Selfridge pseudoprime" n is a "strong Lucas pseudoprime" using Selfridge parameters of: - * Find the first element D in the sequence {5, -7, 9, -11, 13, ...} such that Jacobi(D,n) = -1 - * Then use P=1 and Q=(1-D)/4 in the strong Lucase pseudoprime test. - * Make sure n is not a perfect square, otherwise the search for D will only stop when D=n. - * **********************************************************************************************************/ -int mpz_strongselfridge_prp(mpz_t n) -{ - long int d = 5, p = 1, q = 0; - int max_d = 1000000; - int jacobi = 0; - mpz_t zD; - - if (mpz_cmp_ui(n, 2) < 0) - return PRP_COMPOSITE; - - if (mpz_divisible_ui_p(n, 2)) - { - if (mpz_cmp_ui(n, 2) == 0) - return PRP_PRIME; - else - return PRP_COMPOSITE; - } - - mpz_init_set_ui(zD, d); - - while (1) - { - jacobi = mpz_jacobi(zD, n); - - /* if jacobi == 0, d is a factor of n, therefore n is composite... */ - /* if d == n, then either n is either prime or 9... */ - if (jacobi == 0) - { - if ((mpz_cmpabs(zD, n) == 0) && (mpz_cmp_ui(zD, 9) != 0)) - { - mpz_clear(zD); - return PRP_PRIME; - } - else - { - mpz_clear(zD); - return PRP_COMPOSITE; - } - } - if (jacobi == -1) - break; - - /* if we get to the 5th d, make sure we aren't dealing with a square... */ - if (d == 13) - { - if (mpz_perfect_square_p(n)) - { - mpz_clear(zD); - return PRP_COMPOSITE; - } - } - - if (d < 0) - { - d *= -1; - d += 2; - } - else - { - d += 2; - d *= -1; - } - - /* make sure we don't search forever */ - if (d >= max_d) - { - mpz_clear(zD); - return PRP_ERROR; - } - - mpz_set_si(zD, d); - } - mpz_clear(zD); - - q = (1-d)/4; - - return mpz_stronglucas_prp(n, p, q); - -}/* method mpz_strongselfridge_prp */ - - -/* ********************************************************************************** - * mpz_bpsw_prp: - * A "Baillie-Pomerance-Selfridge-Wagstaff pseudoprime" is a composite n such that - * n is a strong pseudoprime to the base 2 and - * n is a Lucas pseudoprime using the Selfridge parameters. - * **********************************************************************************/ -int mpz_bpsw_prp(mpz_t n) -{ - int ret = 0; - mpz_t two; - - mpz_init_set_ui(two, 2); - - ret = mpz_sprp(n, two); - mpz_clear(two); - - /* with a base of 2, mpz_sprp, won't return PRP_ERROR */ - /* so, only check for PRP_COMPOSITE or PRP_PRIME here */ - if ((ret == PRP_COMPOSITE) || (ret == PRP_PRIME)) - return ret; - - return mpz_selfridge_prp(n); - -}/* method mpz_bpsw_prp */ - - -/* **************************************************************************************** - * mpz_strongbpsw_prp: - * A "strong Baillie-Pomerance-Selfridge-Wagstaff pseudoprime" is a composite n such that - * n is a strong pseudoprime to the base 2 and - * n is a strong Lucas pseudoprime using the Selfridge parameters. - * ****************************************************************************************/ -int mpz_strongbpsw_prp(mpz_t n) -{ - int ret = 0; - mpz_t two; - - mpz_init_set_ui(two, 2); - - ret = mpz_sprp(n, two); - mpz_clear(two); - - /* with a base of 2, mpz_sprp won't return PRP_ERROR */ - /* so, only check for PRP_COMPOSITE or PRP_PRIME here */ - if ((ret == PRP_COMPOSITE) || (ret == PRP_PRIME)) - return ret; - - return mpz_strongselfridge_prp(n); - -}/* method mpz_strongbpsw_prp */ - - -/* ************************ - * mpz_frobenius_prp - * ************************/ - -/* ************************ - * mpz_strongfrobenius_prp - * ************************/ - -/* ************************ - * mpz_elliptic_prp - * ************************/ - - -/********************************************/ -/********************************************/ -/* End of probable primality (prp) routines */ -/********************************************/ -/********************************************/ - - -/***********************************/ -/***********************************/ -/* Start of prime proving routines */ -/***********************************/ -/***********************************/ /* ********************************************************************************** * APR-CL (also known as APRT-CLE) is a prime proving algorithm developed by: @@ -1644,7 +496,6 @@ int mpz_aprtcle(mpz_t N, int verbose) int IV, InvX, LEVELnow, NP, PK, PL, PM, SW, VK, TestedQs, TestingQs; int QQ, T1, T3, U1, U3, V1, V3; int break_this = 0; - int zz = 0; /* make sure the input is >= 2 and odd */ @@ -1674,8 +525,8 @@ int mpz_aprtcle(mpz_t N, int verbose) if (NumberLength > 7000) { if (verbose >= APRTCLE_VERBOSE2) - printf(" Info: Number too large, returning BPSW(N)\n"); - return mpz_bpsw_prp(N); + printf(" Info: Number too large, returning mpz_probab_prime_p(N)\n"); + return mpz_probab_prime_p(N, 1); // mpz_probab_prime_p mpz_bpsw_prp } allocate_vars(); @@ -1722,7 +573,7 @@ int mpz_aprtcle(mpz_t N, int verbose) if (i == LEVELmax) { /* too big */ free_vars(); - return mpz_bpsw_prp(N); + return mpz_probab_prime_p(N, 1); } LEVELnow = i; TestingQs = j; @@ -1734,11 +585,19 @@ int mpz_aprtcle(mpz_t N, int verbose) { for (i = 0; i < NP; i++) { - P = aiP[i+zz]; + P = aiP[i]; while (T%P != 0) { - zz++; - P = aiP[i+zz]; + i++; + if (i >= NP) break; + P = aiP[i]; + } + if (i >= NP) + { /* too big */ + free_vars(); + if (verbose >= APRTCLE_VERBOSE2) + printf(" ERROR: Not enough primes in P to fully process T, returning mpz_probab_prime_p(N)\n"); + return mpz_probab_prime_p(N, 1); } SW = TestedQs = 0; @@ -2270,9 +1129,9 @@ int mpz_aprtcle(mpz_t N, int verbose) if (LEVELnow == LEVELmax) { if (verbose >= APRTCLE_VERBOSE2) - {printf("APR-CL cannot tell: lvlnow == lvlmax, returning BPSW(N)\n"); fflush(stdout);} + {printf("APR-CL cannot tell: lvlnow == lvlmax, returning mpz_probab_prime_p(N)\n"); fflush(stdout);} free_vars(); - return mpz_bpsw_prp(N); /* Cannot tell */ + return mpz_probab_prime_p(N, 1); /* Cannot tell */ } T = aiT[LEVELnow]; NP = aiNP[LEVELnow]; @@ -2298,9 +1157,9 @@ int mpz_aprtcle(mpz_t N, int verbose) } } /* end for J */ if (verbose >= APRTCLE_VERBOSE2) - {printf("Failed: APR-CL error, returning BPSW(N)\n"); fflush(stdout);} + {printf("Failed: APR-CL error, returning mpz_probab_prime_p(N)\n"); fflush(stdout);} free_vars(); - return mpz_bpsw_prp(N); /* Program error */ + return mpz_probab_prime_p(N, 1); /* Program error */ } /* end if */ break; } /* end for (;;) */ @@ -2335,9 +1194,9 @@ int mpz_aprtcle(mpz_t N, int verbose) } /* End for U */ /* This should never be reached. */ if (verbose >= APRTCLE_VERBOSE2) - {printf("Failed: APR-CL error with final test, returning BPSW(N)\n"); fflush(stdout);} + {printf("Failed: APR-CL error with final test, returning mpz_probab_prime_p(N)\n"); fflush(stdout);} free_vars(); - return mpz_bpsw_prp(N); /* Program error */ + return mpz_probab_prime_p(N, 1); /* Program error */ } } diff --git a/aprtcle/mpz_aprcl.h b/aprtcle/mpz_aprcl.h index 34425a4e..c0839058 100644 --- a/aprtcle/mpz_aprcl.h +++ b/aprtcle/mpz_aprcl.h @@ -1,4 +1,4 @@ -/* Copyright 2011,2012,2013 David Cleaver +/* Copyright 2011-2015 David Cleaver * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by @@ -25,113 +25,6 @@ typedef unsigned long long u64_t; #include "jacobi_sum.h" -/* - * The PRP functions presented here are based on the paper: - * Grantham, Jon. Frobenius Pseudoprimes. Math. Comp. 70 (2001), 873-891. - */ - -/*************************************************************/ -/*************************************************************/ -/* These are the definitions for the probable prime routines */ -/*************************************************************/ -/*************************************************************/ -#define PRP_ERROR -1 -#define PRP_COMPOSITE 0 -#define PRP_PRP 1 -#define PRP_PRIME 2 - -/* ****************************************************************** - * mpz_prp: (also called a Fermat pseudoprime) - * A "pseudoprime" to the base a is a composite number n such that, - * (a,n)=1 and a^(n-1) = 1 mod n - * ******************************************************************/ -int mpz_prp(mpz_t n, mpz_t a); - -/* ************************************************************************* - * mpz_euler_prp: (also called a Solovay-Strassen pseudoprime) - * An "Euler pseudoprime" to the base a is an odd composite number n with, - * (a,n)=1 such that a^((n-1)/2)=(a/n) mod n [(a/n) is the Jacobi symbol] - * *************************************************************************/ -int mpz_euler_prp(mpz_t n, mpz_t a); - -/* ********************************************************************************************* - * mpz_sprp: (also called a Miller-Rabin pseudoprime) - * A "strong pseudoprime" to the base a is an odd composite n = (2^r)*s+1 with s odd such that - * either a^s == 1 mod n, or a^((2^t)*s) == -1 mod n, for some integer t, with 0 <= t < r. - * *********************************************************************************************/ -int mpz_sprp(mpz_t n, mpz_t a); - -/* ************************************************************************* - * mpz_fibonacci_prp: - * A "Fibonacci pseudoprime" with parameters (P,Q), P > 0, Q=+/-1, is a - * composite n for which V_n == P mod n - * [V is the Lucas V sequence with parameters P,Q] - * *************************************************************************/ -int mpz_fibonacci_prp(mpz_t n, long int p, long int q); - -/* ******************************************************************************* - * mpz_lucas_prp: - * A "Lucas pseudoprime" with parameters (P,Q) is a composite n with D=P^2-4Q, - * (n,2QD)=1 such that U_(n-(D/n)) == 0 mod n [(D/n) is the Jacobi symbol] - * *******************************************************************************/ -int mpz_lucas_prp(mpz_t n, long int p, long int q); - -/* ********************************************************************************************* - * mpz_stronglucas_prp: - * A "strong Lucas pseudoprime" with parameters (P,Q) is a composite n = (2^r)*s+(D/n), where - * 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 - * for some t, 0 <= t < r. [(D/n) is the Jacobi symbol] - * *********************************************************************************************/ -int mpz_stronglucas_prp(mpz_t n, long int p, long int q); - -/* ******************************************************************************************* - * mpz_extrastronglucas_prp: - * Let U_n = LucasU(p,1), V_n = LucasV(p,1), and D=p^2-4. - * An "extra strong Lucas pseudoprime" to the base p is a composite n = (2^r)*s+(D/n), where - * s is odd and (n,2D)=1, such that either U_s == 0 mod n and V_s == +/-2 mod n, or - * V_((2^t)*s) == 0 mod n for some t with 0 <= t < r-1 [(D/n) is the Jacobi symbol] - * *******************************************************************************************/ -int mpz_extrastronglucas_prp(mpz_t n, long int p); - -/* *********************************************************************************************** - * mpz_selfridge_prp: - * A "Lucas-Selfridge pseudoprime" n is a "Lucas pseudoprime" using Selfridge parameters of: - * Find the first element D in the sequence {5, -7, 9, -11, 13, ...} such that Jacobi(D,n) = -1 - * Then use P=1 and Q=(1-D)/4 in the Lucas pseudoprime test. - * Make sure n is not a perfect square, otherwise the search for D will only stop when D=n. - * ***********************************************************************************************/ -int mpz_selfridge_prp(mpz_t n); - -/* ********************************************************************************************************* - * mpz_strongselfridge_prp: - * A "strong Lucas-Selfridge pseudoprime" n is a "strong Lucas pseudoprime" using Selfridge parameters of: - * Find the first element D in the sequence {5, -7, 9, -11, 13, ...} such that Jacobi(D,n) = -1 - * Then use P=1 and Q=(1-D)/4 in the strong Lucase pseudoprime test. - * Make sure n is not a perfect square, otherwise the search for D will only stop when D=n. - * **********************************************************************************************************/ -int mpz_strongselfridge_prp(mpz_t n); - -/* ********************************************************************************** - * mpz_bpsw_prp: - * A "Baillie-Pomerance-Selfridge-Wagstaff pseudoprime" is a composite n such that - * n is a strong pseudoprime to the base 2 and - * n is a Lucas pseudoprime using the Selfridge parameters. - * **********************************************************************************/ -int mpz_bpsw_prp(mpz_t n); - -/* **************************************************************************************** - * mpz_strongbpsw_prp: - * A "strong Baillie-Pomerance-Selfridge-Wagstaff pseudoprime" is a composite n such that - * n is a strong pseudoprime to the base 2 and - * n is a strong Lucas pseudoprime using the Selfridge parameters. - * ****************************************************************************************/ -int mpz_strongbpsw_prp(mpz_t n); - - - - - - /*******************************************************/ /*******************************************************/ /* These are the definitions for the APRT-CLE routines */ diff --git a/auxi.c b/auxi.c index 742081be..fe98df71 100644 --- a/auxi.c +++ b/auxi.c @@ -24,18 +24,6 @@ along with this program; see the file COPYING. If not, see #include "ecm-impl.h" #include "ecm-ecm.h" -#ifdef HAVE_APRCL -#include "aprtcle/mpz_aprcl.h" -#define ECM_FAC_PRIME APRTCLE_PRIME -#define ECM_FAC_PRP APRTCLE_PRP -#else -#define mpz_aprtcle(x,y) mpz_probab_prime_p(x,PROBAB_PRIME_TESTS) -#define ECM_FAC_PRIME 2 /* mpz_probab_prime_p returns 2 when a number is - definitely prime */ -#define ECM_FAC_PRP 1 /* mpz_probab_prime_p returns 1 when a number is - a probable prime */ -#endif - #ifdef HAVE_GWNUM /* For GWNUM_VERSION */ #include "gwnum.h" @@ -231,23 +219,6 @@ process_newfactor (mpz_t g, int result, mpcandi_t *n, int method, if (verbose >= 1) { -#ifdef HAVE_APRCL -#define APRCL_CUTOFF 400 /* for more than APRCL_CUTOFF digits, print - progress */ -#define APRCL_CUTOFF2 1000 /* for more than APRCL_CUTOFF2 digits, perform - a pseudo-primality test */ - if (n->ndigits > APRCL_CUTOFF2) - cofactor_is_prime = n->isPrp; /* from mpcandi_t_addfoundfactor */ - else if (n->ndigits > APRCL_CUTOFF) - { - printf ("Proving primality of %u-digit cofactor may take a while...\n", n->ndigits); - /* the second argument here tells aprtcle to print primality proving progress info */ - cofactor_is_prime = mpz_aprtcle (n->n, 1); - printf ("\n"); - } - else - cofactor_is_prime = mpz_aprtcle (n->n, 0); -#endif cofactor_is_prime = n->isPrp; /* from mpcandi_t_addfoundfactor */ if (cofactor_is_prime == ECM_FAC_PRIME) diff --git a/candi.c b/candi.c index efad0a04..4c357e07 100644 --- a/candi.c +++ b/candi.c @@ -108,12 +108,19 @@ mpcandi_t_add_candidate (mpcandi_t *n, mpz_t c, const char *cpExpr, strcpy (n->cpExpr, cpExpr); } mpz_set (n->n, c); + n->ndigits = nb_digits (c); if (primetest) - n->isPrp = mpz_probab_prime_p (c, PROBAB_PRIME_TESTS); + { + if (n->ndigits < APRCL_CUTOFF) + n->isPrp = mpz_aprtcle (c, 0); + if (n->ndigits < APRCL_CUTOFF2) + n->isPrp = mpz_aprtcle (c, 1); + else + n->isPrp = mpz_probab_prime_p (c, PROBAB_PRIME_TESTS); + } else n->isPrp = 0; /* there is a candidate there now, and the user did not tell us to prp it, so assume it is composite */ - n->ndigits = nb_digits (c); #if defined (CANDI_DEBUG) Candi_Validate("Post mpcandi_t_add_candidate", n); @@ -147,7 +154,12 @@ mpcandi_t_addfoundfactor (mpcandi_t *n, mpz_t f, int displaywarning) /* remove f from n->n */ mpz_divexact (n->n, n->n, f); n->ndigits = nb_digits (n->n); - n->isPrp = mpz_probab_prime_p (n->n, PROBAB_PRIME_TESTS); + if (n->ndigits < APRCL_CUTOFF) + n->isPrp = mpz_aprtcle (n->n, 0); + if (n->ndigits < APRCL_CUTOFF2) + n->isPrp = mpz_aprtcle (n->n, 1); + else + n->isPrp = mpz_probab_prime_p (n->n, PROBAB_PRIME_TESTS); if (n->cpExpr != NULL) { /* If there is an expression, then lets preserve it */ diff --git a/ecm-ecm.h b/ecm-ecm.h index 190ee672..0278c1af 100644 --- a/ecm-ecm.h +++ b/ecm-ecm.h @@ -175,4 +175,24 @@ long PeakMemusage (void); #define ABS(x) ((x) >= 0 ? (x) : -(x)) +#ifdef HAVE_APRCL +#include "aprtcle/mpz_aprcl.h" +#define ECM_FAC_PRIME APRTCLE_PRIME +#define ECM_FAC_PRP APRTCLE_PRP +#else +#define mpz_aprtcle(x,y) mpz_probab_prime_p(x,PROBAB_PRIME_TESTS) +#define ECM_FAC_PRIME 2 /* mpz_probab_prime_p and aprcl returns 2 when + a number is definitely prime */ +#define ECM_FAC_PRP 1 /* mpz_probab_prime_p and aprcl returns 1 when + a number is a probable prime */ +#endif + +/* Cutoff value: Show APRCL progress only if n has more digits than cutoff */ +#define APRCL_CUTOFF 400 /* for more than APRCL_CUTOFF digits, print + progress */ +/* Cutoff value: Use APRCL if n has fewer digits than cutoff2 + Use mpz_probab_prime_p if n has more digits than cutoff2 */ +#define APRCL_CUTOFF2 1000 /* for more than APRCL_CUTOFF2 digits, perform + a pseudo-primality test */ + #endif /* _ECM_ECM_H */