Math-Prime-Util-GMP
view release on metacpan or search on metacpan
#include "squfof126.h"
#include "factor.h"
#include "simpqs.h"
#include "ecm.h"
#define _GMP_ECM_FACTOR(n, f, b1, ncurves) \
_GMP_ecm_factor_projective(n, f, b1, 0, ncurves)
#include "utility.h"
static const int maxuibits = (sizeof(unsigned long) > BITS_PER_WORD)
? BITS_PER_WORD : (8*sizeof(unsigned long));
/*
*
* Theorems where we are factoring n-1.
*
*
* Lucas (1876): Given a completely factored n-1, if there exists an a s.t.
* a^(n-1) % n = 1
* a^((n-1/f) % n != 1 for ALL factors f of n-1
* then n is prime.
*
* PPBLS:, given n-1 = A*B, A > sqrt(n), if we can find an a s.t.
* a^A % n = 1
* gcd(a^(A/f)-1,n) = 1 for ALL factors f of A
* then n is prime.
*
* Generalized Pocklington: given n-1 = A*B, gcd(A,B)=1, A > sqrt(n), then if
* for each each factor f of A, there exists an a (1 < a < n-1) s.t.
* a^(n-1) % n = 1
* gcd(a^((n-1)/f)-1,n) = 1
* then n is prime.
*
* BLS T5: given n-1 = A*B, factored A, s=B/2A r=B mod (2A), then if:
* - A is even, B is odd, and AB=n-1 (all implied by n = odd and the above),
* - n < (A+1) * (2*A*A + (r-1) * A + 1)
* - for each each factor f of A, there exists an a (1 < a < n-1) s.t.
* - a^(n-1) % n = 1
* - gcd(a^((n-1)/f)-1,n) = 1 for ALL factors f of A
* then:
* if s = 0 or r*r - 8*s is not a perfect square
* n is prime
* else
* n is composite
*
* BLS T7: given n-1 = A*B, factored A, a B1 where all factors of B are >= B1,
* s=B/2A r=B mod (2A), then if:
* - A is even, B is odd, and AB=n-1 (all implied by n = odd and the above),
* - n < (B1*A+1) * (2*A*A + (r-B1)*A + 1)
* - for each f in {B, factors of A}, there exists an a (1 < a < n-1) s.t.
* - a^(n-1) % n = 1
* - gcd(a^((n-1)/f)-1,n) = 1 for ALL factors f of A
* then:
* if s = 0 or r*r - 8*s is not a perfect square
* n is prime
* else
* n is composite
*
* The generalized Pocklington test is also sometimes known as the
* Pocklington-Lehmer test. It's definitely an improvement over Lucas
* since we only have to find factors up to sqrt(n), _and_ we can choose
* a different 'a' value for each factor. This is corollary 1 from BLS75.
*
* BLS T5 is the Brillhart-Lehmer-Selfridge 1975 theorem 5 (see link below).
* We can factor even less of n, and the test lets us kick out some
* composites early, without having to test n-3 different 'a' values.
*
* Once we've found the factors of n-1 (or enough of them), verification
* usually happens really fast. a=2 works for most, and few seem to require
* more than ~ log2(n). However all but BLS75 require testing all integers
* 1 < a < n-1 before answering in the negative, which is impractical.
*
* BLS75 theorem 7 is the final n-1 theorem and takes into account any
* knowledge that the remaining factor is not below a threshold B. Since
* we do initial trial division this helps. It is usually of only small
* benefit, but it costs nothing.
*
*
* AKS is not too hard to implement, but it's impractically slow.
*
* ECPP is very fast and definitely the best method for most numbers.
*
* APR-CL is very practical for numbers of a few hundred digits.
*
* BLS75: http://www.ams.org/journals/mcom/1975-29-130/S0025-5718-1975-0384673-1/S0025-5718-1975-0384673-1.pdf
*
*/
/*
* Theorems where we are factoring n+1.
*
* Corollary 8 is analagous to Corollary 1 (generalized Pocklington).
* Theorem 15 is analagous to Theorem 3. These are simple proofs for the case
* of n+1 (theorem 15) / n-1 (theorem 3) having a large prime factor. These
* are used by ECPP. T3 is slightly better than the original Proth theorem.
* Theorem 17 is analagous to Theorem 5.
* Theorem 19 is analagous to Theorem 7.
*
* Theorem 20 is a hybrid, combining theorems 7 and 19.
*/
/* Like all the primality functions:
* 2 = definitely prime, 1 = maybe prime, 0 = definitely composite
*
* You really should run is_prob_prime on n first, so we only have to run
* these tests on numbers that are very probably prime.
*/
static int tfe(mpz_t f, const mpz_t n, int effort)
{
int success = 0;
UV log2n = mpz_sizeinbase(n, 2);
if (mpz_cmp_ui(n,3) <= 0) {
mpz_set(f,n);
return 1;
}
/*
17: N < (m F2 - 1) ( 2 F2 F2 + m F2 - |r| F2 + 1 )
(III) test (l*F2+1) doesn't divide N for 1 .. m
19: N < (B2 F2 - 1) ( 2 F2 F2 + B2 F2 - |r| F2 + 1 )
(III) (IV) R2 factors > B2
*/
#if 0
static int 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)
{
mpz_mul(t, F2, R2);
mpz_sub_ui(t, t, 1);
if (mpz_cmp(t, n) != 0) croak("BLS75 internal error: F2*R2 != n+1\n");
mpz_mul_ui(t, F2, 2);
mpz_tdiv_qr(s, r, R2, t);
if (mpz_cmp(r, F2) >= 0) {
mpz_add_ui(s, s, 1);
mpz_sub(r, r, t);
}
/* Let m = 1 */
mpz_add_ui(y, t, 1);
mpz_abs(t, r);
mpz_sub(y, y, t);
mpz_mul(y, y, F2);
mpz_add_ui(y, y, 1); /* y = 2F2^2 + (m-r)F2 + 1 */
mpz_sub_ui(t, F2, 1);
mpz_mul(y, y, t); /* times (mF2-1) */
return (mpz_cmp(n, y) < 0) ? 1 : 0;
}
#endif
static int bls_theorem19_limit(const mpz_t n, mpz_t F2, mpz_t R2, UV B2,
mpz_t t, mpz_t y, mpz_t r, mpz_t s)
{
mpz_mul(t, F2, R2);
mpz_sub_ui(t, t, 1);
if (mpz_cmp(t, n) != 0) croak("BLS75 internal error: F2*R2 != n+1\n");
mpz_mul_ui(t, F2, 2);
mpz_tdiv_qr(s, r, R2, t);
if (mpz_cmp(r, F2) >= 0) {
mpz_add_ui(s, s, 1);
mpz_sub(r, r, t);
}
mpz_add_ui(y, t, B2);
mpz_abs(t, r);
mpz_sub(y, y, t);
mpz_mul(y, y, F2);
mpz_add_ui(y, y, 1); /* y = 2F2^2 + (B2 - r)F2 + 1 */
mpz_mul_ui(t, F2, B2);
mpz_sub_ui(t, t, 1);
mpz_mul(y, y, t); /* times (B2F2-1) */
return (mpz_cmp(n, y) < 0) ? 1 : 0;
}
static int bls_corollary11_limit(const mpz_t n, mpz_t R1, mpz_t F1, mpz_t F2, UV B,
mpz_t t, mpz_t g, mpz_t r, mpz_t s)
{
if (mpz_cmp(F1,F2) >= 0) {
mpz_tdiv_q_2exp(t, F2, 1);
mpz_mul(t, t, F1);
mpz_mul(t, t, F1);
} else {
mpz_tdiv_q_2exp(t, F1, 1);
mpz_mul(t, t, F2);
mpz_mul(t, t, F2);
}
mpz_ui_pow_ui(g, B, 3);
mpz_mul(t, t, g);
return (mpz_cmp(n, t) < 0);
}
static int bls_theorem20_limit(const 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)
{
int m_used = 0;
if (bls_corollary11_limit(n,R1,F1,F2,B,t,g,r,s)) {
mpz_set_ui(s,0); /* No test for m needed */
return 1;
}
mpz_mul_ui(t, F1, B);
mpz_add_ui(g, t, 1);
mpz_mul_ui(t, F2, B);
mpz_sub_ui(t, t, 1);
if (mpz_cmp(t, g) > 0) mpz_set(g, t);
mpz_tdiv_q_2exp(t, F2, 1);
mpz_tdiv_qr(s, r, R1, t);
mpz_mul(t, F1, F2);
mpz_mul_ui(t, t, m);
mpz_tdiv_q_2exp(t, t, 1);
mpz_mul(r, r, F1); /* r *= F1; t += r */
mpz_add(t, t, r);
mpz_add_ui(t, t, 1); /* t = m * F1 * F2/2 + r * F1 + 1 */
if (mpz_cmp(t, g) > 0) {
m_used = 1;
mpz_set(g, t);
}
mpz_mul(t, F1, F2);
mpz_tdiv_q_2exp(t, t, 1);
mpz_mul_ui(t, t, B);
mpz_mul_ui(t, t, B);
mpz_add_ui(s, t, 1); /* s = B1*B2*F1*F2/2+1 */
mpz_mul(g, g, s);
mpz_set_ui(s, m_used); /* Use s to signal whether we must test for m. */
return (mpz_cmp(n, g) < 0);
}
/******************************************************************************/
/* (I) For each prime p_i dividing F1 [N-1 = F1R1] there exists an a_i
* such that N is a psp base a_i and gcd(a_i^{(N-1)/p_i}-1,N) = 1.
*/
static int _verify_cond_I_p(const mpz_t n, mpz_t pi, mpz_t ap, mpz_t t, int alimit, signed char* pspcache)
{
int a, success = 0;
PRIME_ITERATOR(iter);
for (a = 2; !success && a <= alimit; a = prime_iterator_next(&iter)) {
int psp = pspcache ? pspcache[a] : -1;
mpz_set_ui(ap, a);
if (psp == -1) {
mpz_sub_ui(t, n, 1);
mpz_powm(t, ap, t, n);
psp = (mpz_cmp_ui(t, 1) == 0);
}
if (!psp)
return -1; /* We failed a Fermat test, n is composite. */
/* TODO: optimize for cases of both n-1 and n+1 working */
if (nstack(&p1stack) > 0) {
while (nstack(&p1stack) > 0) {
int pr = 1;
pop_fstack(f, &p1stack);
if (effort > low_effort)
pr = BLS_primality(f, effort, prooftextptr);
if (pr == 0) croak("probable prime factor proved composite");
else if (pr == 2) push_fstack(&f1stack, f); /* Proved, put on F stack */
else factor_out(F1, R1, f); /* No proof. Move to R */
}
}
if (nstack(&p2stack) > 0) {
while (nstack(&p2stack) > 0) {
int pr = 1;
pop_fstack(f, &p2stack);
if (effort > low_effort)
pr = BLS_primality(f, effort, prooftextptr);
if (pr == 0) croak("probable prime factor proved composite");
else if (pr == 2) push_fstack(&f2stack, f); /* Proved, put on F stack */
else factor_out(F2, R2, f); /* No proof. Move to R */
}
}
if (ev) gmp_printf("start tidy: N %Zd F2 %Zd R2 %Zd B1 %lu\n", n, F2, R2, B1);
fstack_tidy(&f1stack, B1);
fstack_tidy(&f2stack, B1);
if (ev) gmp_printf("end tidy: N %Zd F2 %Zd R2 %Zd B1 %lu\n", n, F2, R2, B1);
#if PRINT_PCT
fin_pct = (100.0 * (mpz_sizeinbase(F1,2) + mpz_sizeinbase(F2,2))) / (mpz_sizeinbase(nm1,2) + mpz_sizeinbase(np1,2));
printf("%6.2f .. ", fin_pct); fflush(stdout);
printf("\n"); fflush(stdout);
#endif
if (ev) gmp_printf("check: N %Zd F2 %Zd R2 %Zd B1 %lu\n", n, F2, R2, B1);
/* Check the theorems we have available */
/* If we can do a standard n-1 proof, do that. */
if (bls_theorem7_limit(n, F1, R1, B1, t, u, r, s)) {
if (get_verbose_level() > 0) printf("BLS75 proof using N-1\n");
if (ev) gmp_printf("N %Zd F1 %Zd R1 %Zd B1 %lu\n", n, F1, R1, B1);
trim_factors(F1, R1, n, nm1, B1, &f1stack, &bls_theorem7_limit, t, u, r, s);
if (ev) gmp_printf("N %Zd F1 %Zd R1 %Zd B1 %lu\n", n, F1, R1, B1);
for (pcount = 0; success > 0 && pcount < f1stack.cur; pcount++)
success = _verify_cond_I_p(n, f1stack.stack[pcount], u, t, 1000, 0);
if (success > 0 && (mpz_mul(t, F1, F1), mpz_cmp(t,n) > 0))
goto end_hybrid; /* Corollary 1, n-1 factored more than sqrt(n) */
if (success > 0 && !bls_theorem5_limit(n, F1, R1, B1, t, u, r, s))
success = _verify_cond_I_p(n, R1, u, t, 1000, 0);
if (success > 0) {
mpz_mul(t, r, r);
mpz_submul_ui(t, s, 8); /* t = r^2 - 8s */
/* N is prime if and only if s=0 OR t not a perfect square */
success = (mpz_sgn(s) == 0 || !mpz_perfect_square_p(t)) ? 1 : -1;
}
goto end_hybrid; /* Theorem 5 or 7 */
}
/* Rather than looking at theorem 19 for an N+1 proof, we'll just go to
* theorem 20 (or corollary 11). T20 is faster. */
/* Check N < B^3 F1*F2*F2/2 or N < B^3 F1*F1*F2/2 */
success = bls_theorem20_limit(n, R1, F1, F2, B1, m, t, u, r, s);
if (get_verbose_level() > 0) printf("BLS75 proof using N-1 / N+1 (T20)\n");
/* Trim some factors from f2stack if possible */
if (nstack(&f2stack) > 1) {
int i;
mpz_set_ui(F2, 1);
mpz_set(R2, np1);
for (i = 0; i < f2stack.cur; i++) {
if (i > 0 && bls_theorem20_limit(n, R1, F1, F2, B1, m, t, u, r, s))
break;
factor_out(R2, F2, f2stack.stack[i]);
}
/* Remove excess factors */
while (i < f2stack.cur)
pop_fstack(t, &f2stack);
/* Verify Q[0] = 2 */
if (mpz_cmp_ui(f2stack.stack[0], 2) != 0)
croak("BLS75 internal error: 2 not at start of fstack");
}
/* Check lambda divisibility if needed */
if (success > 0 && mpz_sgn(s)) {
UV lambda;
mpz_mul(t, F1, F2);
mpz_tdiv_q_2exp(t, t, 1);
mpz_mul(u, r, F1);
mpz_add_ui(u, u, 1);
for (lambda = 0; success > 0 && lambda < m; lambda++, mpz_add(u,u,t)) {
if (lambda > 0 || mpz_sgn(r))
if (mpz_divisible_p(n, u)) {
/* Check that we found a non-trivial divisor */
mpz_gcd(t, u, n);
success = (mpz_cmp_ui(t,1) > 0 && mpz_cmp(t,n) < 0) ? -1 : 0;
break;
}
}
}
/* Verify (I) page 623 and (II) page 625 */
if (success > 0) {
for (pcount = 0; success > 0 && pcount < f1stack.cur; pcount++)
success = _verify_cond_I_p(n, f1stack.stack[pcount], u, t, 1000, 0);
if (success > 0)
success = _verify_cond_I_p(n, R1, u, t, 1000, 0);
}
/* Verify (III) page 631 and (IV) page 633 */
if (success > 0) {
success = _test_III_IV(n, np1, R2, &f2stack, t, u, r, s);
}
#if 0
{ double p1 = (100.0 * mpz_sizeinbase(F1,2) / mpz_sizeinbase(nm1,2));
double p2 = (100.0 * mpz_sizeinbase(F2,2) / mpz_sizeinbase(np1,2));
printf("%6.2f %6.2f\n", p1, p2); fflush(stdout); }
//{ 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); }
( run in 1.224 second using v1.01-cache-2.11-cpan-df04353d9ac )