www.pudn.com > calc.rar > mpqsieve.c
/* mpqsieve.c */
/*
* Multiple Quadratic sieve method of factoring.
* Based on "The multiple polynomial quadratic sieve ", by R.D. Silverman,
* Mathematics of Computation, 48, 1987, 329-339.
*/
#include
#include
#include "integer.h"
#include "fun.h"
MPI *MPQS1(MPI *N)
/*
* The packed array version of MPQS(), done with a lot of help from Peter Adams.
*/
{
register unsigned int z;
unsigned int gap, h, i, j, t;
unsigned long a, a1, b, b1, l1, l2, r, e, k;
unsigned int TEST_ROWS;
unsigned int factored = 0, flag, found = 0;
unsigned int FBASE, TOLERANCE;
unsigned int CNUF;
unsigned long SRANGE, SSRANGE;
unsigned int *logn, *logQ;
unsigned long *r1, *r2, *soln, **M1, **M2, n;
unsigned int *exponents, *CTR;
MPI *N0, *X1, *X2, *X3, *X4, *D, *H0, *H1, *H2, *HH, *SAVE1, *Y;
MPI *B, *A, *Q, *Tmp, **H, **cofactor, *ROOT2N, *OVER2, *ROOTN;
MPI *P1, *P2, *T, *G, *srange, *xx, *hh;
int *INDEX, s;
long x, y, L, L1, L2, *base;
unsigned int m1cols, m2cols, m1bitcols, m2bitcols, index, ctr;
G = NULL;
n = LENGTHI(N);
printf("N = ");PRINTI(N);printf(" has %lu digits\n", n);
/* FBASE = no of odd primes in factor base. */
if (n <= 20)
{
FBASE = 150; SRANGE = 5000; TOLERANCE = 10;
}
else if (n > 20 && n <= 25)
{
FBASE = 150; SRANGE = 10000; TOLERANCE = 15;
}
else if (n > 25 && n <= 30)
{
FBASE = 250; SRANGE = 15000; TOLERANCE = 20;
}
else if (n > 30 && n <= 35)
{
FBASE = 400; SRANGE = 15000; TOLERANCE = 20;
}
else if (n > 35 && n <= 40)
{
FBASE = 500; SRANGE = 25000; TOLERANCE = 20;
}
else if (n > 40 && n <= 45)
{
FBASE = 900; SRANGE = 50000; TOLERANCE = 20;
}
else if (n > 45 && n <= 50)
{
FBASE = 1200; SRANGE = 100000; TOLERANCE = 22;
}
else /* n > 50 */
{
if (n > 61)
{
printf("Number of digits > 60.\n");
return(NULL);
}
FBASE = 3000; SRANGE = 180000; TOLERANCE = 22;
}
printf("FBASE = %u:", FBASE);
printf(" sieve_range = %lu:", SRANGE);
printf(" tolerance = %u:\n", TOLERANCE);
ROOTN = BIG_MTHROOT(N, 2);
TEST_ROWS = FBASE + 2;
INDEX = (int *)mmalloc((USL)((TEST_ROWS) * sizeof(int)));
INTSETI(INDEX, (USL)TEST_ROWS, -1);
CTR = (unsigned int *)ccalloc(TEST_ROWS, sizeof(unsigned int));
m1cols = (FBASE + 2);
m1bitcols = m1cols/O32 + 1;
M1 = idim2(TEST_ROWS, m1bitcols);
H = (MPI **)mmalloc((USL)((TEST_ROWS) * sizeof(MPI *)));
cofactor = (MPI **)mmalloc((USL)((TEST_ROWS) * sizeof(MPI *)));
base = (long *)mmalloc((USL)((FBASE + 2) * sizeof(long)));
base[0] = 2; base[FBASE + 1] = -1;
soln = (unsigned long *)mmalloc((USL)((FBASE + 1) * sizeof(unsigned long)));
/* soln[i] is a solution of x*x=n (mod base[i], 1 <= i <= FBASE. */
QERATOSTHENES(N, (USL *)base, soln, FBASE);
CNUF = BINARYB(N) / 2 + binary(SRANGE)- (TOLERANCE * binary((USL)(base[FBASE])) / 10);
r1 = (unsigned long *)mmalloc((USL)((FBASE + 1) * sizeof(unsigned long)));
r2 = (unsigned long *)mmalloc((USL)((FBASE + 1) * sizeof(unsigned long)));
/* we will sieve over -SRANGE <= x <= SRANGE. */
X1 = MULT_I(N, 2L);
ROOT2N = BIG_MTHROOT(X1, 2);
FREEMPI(X1);
SSRANGE = (USL)2 * SRANGE;
logQ = (unsigned int *)mmalloc((USL)(((unsigned int)SSRANGE + 1) * sizeof(unsigned int)));
logn = (unsigned int *)mmalloc((USL)((FBASE + 1) * sizeof(unsigned int)));
/* Table lookup: logn[i] = binary((USL)(base[i])), 1 <= i <= FBASE. */
k = N->V[0];
if (k % 8 == 1)
logn[0] = 3;
if (k % 8 == 5)
logn[0] = 2;
if (k % 4 == 3)
logn[0] = 1;
for (i = 1; i <= FBASE; i++)
logn[i] = binary((USL)(base[i]));
h = 0;
srange = CHANGE(SRANGE);
X3 = INT0(ROOT2N, srange); /* bug site */
FREEMPI(srange);
OVER2 = BIG_MTHROOT(X3, 2);
FREEMPI(X3);
/* calculating the polynomial Q(x) = A*x*x + 2*B*x + C */
Tmp = CHANGE((USL)(base[FBASE]));
if (RSV(OVER2, Tmp) <= 0)
{
T = OVER2;
OVER2 = ADD0_I(Tmp, (USL)1);
FREEMPI(T);
}
FREEMPI(Tmp);
while (factored < TEST_ROWS)
{
printf("changing polynomial\n");
if (h)
{
Tmp = OVER2;
hh = CHANGE(h);
OVER2 = ADD0I(OVER2, hh); /* bugsite */
FREEMPI(hh);
FREEMPI(Tmp);
}
D = PINCH_PRIME(OVER2, N, base, FBASE, &gap);
/*
D = PROB_PRIME(OVER2, N, &gap);
*/
h += gap + 2;
X1 = INT0_(D, (USL)4);
N0 = MOD0(N, D);
H0 = MPOWER(N0, X1, D);
FREEMPI(X1);
H1 = MULTM(N0, H0, D);
FREEMPI(N0);
A = MULTI(D, D);
X1 = ADD0_I(D, (USL)1);
X2 = INT0_(X1, (USL)2);
FREEMPI(X1);
Y = MULTM(X2, H0, D);
FREEMPI(H0);
FREEMPI(X2);
X1 = MULTI(H1, H1);
X2 = SUBI(N, X1);
FREEMPI(X1);
X3 = INT(X2, D);
FREEMPI(X2);
Tmp = MOD(X3, D);
FREEMPI(X3);
H2 = MULTM(Y, Tmp, D);
FREEMPI(Y);
FREEMPI(Tmp);
X4 = MULTI(H2, D);
FREEMPI(H2);
B = ADD0I(H1, X4);
FREEMPI(H1);
FREEMPI(X4);
SAVE1 = INVERSEM(D, N); /* 1/(D) (mod N) */
/*
printf("A = ");
PRINTI(A);
X1 = MULTI(B, B);
printf(", B = ");
PRINTI(B);
X2 = SUBI(X1, N);
C = INT(X2, A);
X3 = MOD(X2, A);
FREEMPI(X1);
FREEMPI(X2);
FREEMPI(X3);
printf(", C = ");
PRINTI(C);
printf("\n");
*/
/* C is freed later */
/*
printf("calculating the two roots r1[i], r2[i] (mod base[i]) of Q(x)= 0\n");
*/
r = (B->V[0]) % (USL)2 ? (USL)0 : (USL)1;
for (i = 1; i <= FBASE; i++)
{
/*
printf("mod base[%u], A, B, C = :%lu, %lu, %lu\n",i, MOD_(A,(USL)(base[i])), MOD_(B,(USL)(base[i])), MOD_(C,(USL)(base[i])));
*/
b = MOD0_(B, (USL)(base[i]));
b1 = MINUSm(b, (USL)(base[i]));
a = MOD0_(A, (USL)(base[i]));
a1 = INVERSEm(a, (USL)(base[i]));
r1[i] = MULTm(ADDm(b1, soln[i], (USL)(base[i])), a1, (USL)(base[i]));
r2[i] = MULTm(ADDm(b1, MINUSm(soln[i], (USL)(base[i])), (USL)(base[i])), a1, (USL)(base[i]));
/*
printf("a=%lu,a1 = %lu, b1 = %lu\n",a, a1, b1);
printf("soln[%u] = %lu\n", i, soln[i]);
printf("r1[%u] = %lu, r2[%u] = %lu ",i, r1[i], i,r2[i]);
printf("base[%u] = %lu\n", i, base[i]);
*/
}
/*
printf("finished calculating the two roots r1, r2 (mod base[i]) of Q(x)= 0\n");
*/
/* sieving over the range -SRANGE <= x <= SRANGE */
INTSETUI(logQ, (USL)SSRANGE + 1, 0);
L = SRANGE % (USL)2 == r ? -SRANGE : -SRANGE + 1;
for (y = L + SRANGE; y <= SSRANGE; y += 2)
logQ[y] += logn[0];
for (i = 1; i <= FBASE; i++)
{
l1 = (SRANGE + r1[i]) / (USL)(base[i]);
l2 = (SRANGE + r2[i]) / (USL)(base[i]);
L1 = r1[i] -(long)(l1 * (USL)(base[i]));
L2 = r2[i] - (long)(l2 * (USL)(base[i]));
/*
printf("base=%ld,r1=%lu,r2=%lu,l1=%lu,l2=%lu,L1=%ld,L2=%ld\n",base[i],r1[i],r2[i],l1,l2,L1,L2);
*/
for (y = L1 + SRANGE; y <= SSRANGE; y += base[i])
logQ[y] += logn[i];
for (y = L2 + SRANGE; y <= SSRANGE; y += base[i])
logQ[y] += logn[i];
}
/* Scanning the values logQ[x] to see if any exceed CNUF */
for (y = 0; y <= SSRANGE; y++)
{
if (logQ[y] >= CNUF)
{
x = y - SRANGE;
xx = CHANGEL(x);/* bugsite */
X1 = MULTI(A, xx);
FREEMPI(xx);
X2 = ADDI(X1, B);
X3 = MOD(X2, N);
H[factored] = MULTM(X3, SAVE1, N);
HH = MULTM(H[factored], H[factored], N);
if (RSV(X2, ROOTN) <= 0)
Q = SUBI(HH, N); /* Q->S = -1 */
else
Q = COPYI(HH); /* Q->S = 1 */
FREEMPI(HH);
FREEMPI(X1);
FREEMPI(X2);
FREEMPI(X3);
cofactor[factored] = FACTORQ1(Q, factored, (USL *)base, M1, FBASE, x, r1, r2, r, INDEX, CTR);
/* this factors Q over the factor base,
updates M1[factored][] and returns the cofactor. */
t = EQONEI(cofactor[factored]);
FREEMPI(cofactor[factored]);
if (t)
{
printf("Q[%u] = ",factored);PRINTI(Q);printf("\n");
cofactor[factored] = Q;
factored++;
}
else
{
FREEMPI(H[factored]);
FREEMPI(Q);
}
if (factored == TEST_ROWS)
break;
}
}
FREEMPI(SAVE1);
FREEMPI(A);
FREEMPI(B);
FREEMPI(D);
}
/*
printf("INDEX = ");
for (i = 0; i < TEST_ROWS; i++)
printf("%d,", INDEX[i]);
printf("\n");
printf("CTR = ");
for (i = 0; i < TEST_ROWS; i++)
printf("%u,", CTR[i]);
printf("\n");
printf("The packed array M1 is:\n");
printf("array M1 is:\n");
for (i = 0; i < TEST_ROWS; i++)
{
for (j = 0; j = 0; s--)
{
flag = 1;
for (j = 0; j < m1bitcols; j++)
if (M1[s][j])
{
flag = 0;
break;
}
if (flag == 0)
continue;
/* found a zero row, initialize exponent sums to zero */
INTSETUI(exponents, (USL)FBASE + 1, 0);
P1 = ONEI();
for (e = 0; e < TEST_ROWS; e++)
{
/*
f = (M2[s][e>>5] & ((USL)1<<(31-(e&31))))?1:0;
*/
if (get_element(M2,s,e))
{
T = MULTM(P1, H[e], N);
FREEMPI(P1);
P1 = T;
EXP_UPDATE(cofactor[e], (USL *)base, FBASE, exponents);
}
}
P2 = ONEI();
for (j = 0; j <= FBASE; j++)
{
if (exponents[j])
{
Y = CHANGE((USL)(base[j]));
T = MPOWER_M(Y,(unsigned long)(exponents[j] / 2), N);
FREEMPI(Y);
Tmp = P2;
P2 = MULTM(T, P2, N);
FREEMPI(T);
FREEMPI(Tmp);
}
}
T = SUBI(P1, P2);
G = GCD(T, N);
FREEMPI(T);
printf("P1 = "); PRINTI(P1);printf("\n");
printf("P2 = "); PRINTI(P2);printf("\n");
printf("GCD(P2 - P1, N) = "); PRINTI(G);printf("\n");
FREEMPI(P1);
FREEMPI(P2);
if (!EQONEI(G) && !EQUALI(G, N))
{
printf("FOUND A NON_TRIVIAL FACTOR: "); PRINTI(G);printf("\n");
found = 1;
break;
}
FREEMPI(G);
}
printf("N = "); PRINTI(N);printf("\n");
FREEMPI(ROOTN);
FREEMPI(ROOT2N);
FREEMPI(OVER2);
ffree((char *)base, (FBASE + 2) * sizeof(long));
ffree((char *)soln, (FBASE + 1) * sizeof(unsigned long));
ffree((char *)r1, (FBASE + 1) * sizeof(unsigned long));
ffree((char *)r2, (FBASE + 1) * sizeof(unsigned long));
ffree((char *)logn, (FBASE + 1) * sizeof(unsigned int));
ffree((char *)logQ, (SSRANGE + 1) * sizeof(unsigned int));
for (i = 0; i < factored; i++)
FREEMPI(H[i]);
ffree((char *)H, (TEST_ROWS) * sizeof(MPI *));
for (i = 0; i < TEST_ROWS; i++)
FREEMPI(cofactor[i]);
ffree((char *)cofactor, (TEST_ROWS) * sizeof(MPI *));
ffree((char *)INDEX, (TEST_ROWS) * sizeof(int));
ffree((char *)CTR, (TEST_ROWS) * sizeof(unsigned int));
ffree((char *)exponents, (FBASE + 1) * sizeof(unsigned int));
ifree2(M1, TEST_ROWS, m1bitcols);
ifree2(M2, TEST_ROWS, m2bitcols);
if (found)
return (G);
else
return (NULL);
}