X-Git-Url: https://git.saurik.com/apple/security.git/blobdiff_plain/80e2389990082500d76eb566d4946be3e786c3ef..d8f41ccd20de16f8ebe2ccc84d47bf1cb2b26bbb:/Security/libsecurity_cryptkit/lib/CurveParamDocs/schoofs.c diff --git a/Security/libsecurity_cryptkit/lib/CurveParamDocs/schoofs.c b/Security/libsecurity_cryptkit/lib/CurveParamDocs/schoofs.c new file mode 100644 index 00000000..bf70b3ab --- /dev/null +++ b/Security/libsecurity_cryptkit/lib/CurveParamDocs/schoofs.c @@ -0,0 +1,1044 @@ +/* schoofs.c + + Elliptic curve order calculator, for + + y^2 = x^3 + a x + b (mod p) + + (NOTE: + This version has order sieving, triggering on the + initial b parameter and incrementing same. + Parameter details are described in schoof.c) + + Compile with: + + % cc -O schoofs.c giants.c tools.c -lm -o schoofs + + and run, entering the a,b parameters. + + * Change history: + + 20 Mar 01 (REC) Added binarysegsquare() and base-4 ladder + 20 Mar 01 (REC) Bumped MAX_DIGS as in schoof.c + 4 Feb 99 (REC) Sieving invoked. + 2 Feb 99 (REC) Added explicit CRT result + 12 Jan 99 (REC) Removed (hopefully) last of memory crashes + 20 Jan 98 (REC) Creation + + * c. 1998 Perfectly Scientific, Inc. + * All Rights Reserved. + * + * + *************************************************************/ + +#include +#include +#include +#include"giants.h" +#include "tools.h" + +#define P_BREAK 32 + + +#define Q 256 /* See schoof.c for explanation. */ +#define L_MAX 100 +#define MAX_DIGS (2 + (Q+15)/8) /* Size of (squared) coefficients. */ +#define MAX_COEFFS ((L_MAX+1)*(L_MAX+2)) + +typedef struct + { + int deg; + giant *coe; + } polystruct; +typedef polystruct *poly; + + +static giant *mcand, coe, tmp, err, aux, aux2, globx, globy, t1, t2, + t3, t4, t5; +static poly qscratch, rscratch, sscratch, tscratch, pbuff, pbuff2, + precip, cubic, powx, powy, kxn, kyn, kxd, kyd, + txn, txd, tyn, tyd, txn1, txd1, tyn1, tyd1, + nx, dx, ny, dy, mn, md, tmp1, tmp2, tmp3, tmp4; +static poly s[L_MAX+2], smonic; +static giant p, a, b, magcheck; +int L; + +void quickmodg(giant g, giant x) +{ int sgn = x->sign; + + if(sgn == 0) return; + if(sgn > 0) { + if(gcompg(x, g) >= 0) subg(g, x); + return; + } + addg(g,x); + return; +} + +int +log_2(int n) { + int c = 1, d = -1; + while(c <= n) { + c <<= 1; + ++d; + } + return(d); +} + +void +justifyp(poly x) { + int j; + for(j = x->deg; j >= 0; j--) { + if(!is0(x->coe[j])) break; + } + x->deg = (j>0)? j : 0; +} + +void +polyrem(poly x) { + int j; + for(j=0; j <= x->deg; j++) { + modg(p, x->coe[j]); + } + justifyp(x); +} + + +giant * +newa(int n) { + giant *p = (giant *)malloc(n*sizeof(giant)); + int j; + for(j=0; jcoe = (giant *)newa(coeffs); + return(pol); +} + +void +init_all() { + int j; + + j = (2*Q + log_2(MAX_COEFFS) + 32 + 15)/16; + globx = newgiant(MAX_COEFFS * j); + globy = newgiant(MAX_COEFFS * j); + + init_tools(MAX_DIGS); + powx = newpoly(MAX_COEFFS); + powy = newpoly(MAX_COEFFS); + kxn = newpoly(MAX_COEFFS); + kxd = newpoly(MAX_COEFFS); + kyn = newpoly(MAX_COEFFS); + kyd = newpoly(MAX_COEFFS); + txn = newpoly(MAX_COEFFS); + txd = newpoly(MAX_COEFFS); + tyn = newpoly(MAX_COEFFS); + tyd = newpoly(MAX_COEFFS); + txn1 = newpoly(MAX_COEFFS); + txd1 = newpoly(MAX_COEFFS); + tyn1 = newpoly(MAX_COEFFS); + tyd1 = newpoly(MAX_COEFFS); + nx = newpoly(MAX_COEFFS); + ny = newpoly(MAX_COEFFS); + dx = newpoly(MAX_COEFFS); + dy = newpoly(MAX_COEFFS); + mn = newpoly(MAX_COEFFS); + md = newpoly(MAX_COEFFS); + tmp1 = newpoly(MAX_COEFFS); + tmp2 = newpoly(MAX_COEFFS); + tmp3 = newpoly(MAX_COEFFS); + tmp4 = newpoly(MAX_COEFFS); + mcand = (giant *)newa(MAX_COEFFS); + + s[0] = newpoly(MAX_COEFFS); + + for(j=1; j<=L_MAX+1; j++) { + s[j] = newpoly(j*(j+1)); + } + smonic = newpoly(MAX_COEFFS); +/* The next three need extra space for remaindering routine. */ + qscratch = newpoly(2*MAX_COEFFS); + rscratch = newpoly(2*MAX_COEFFS); + pbuff = newpoly(2*MAX_COEFFS); + pbuff2 = newpoly(MAX_COEFFS); + sscratch = newpoly(MAX_COEFFS); + tscratch = newpoly(MAX_COEFFS); + tmp = newgiant(MAX_DIGS); + err = newgiant(MAX_DIGS); + coe = newgiant(MAX_DIGS); + aux = newgiant(MAX_DIGS); + aux2 = newgiant(MAX_DIGS); + t1 = newgiant(MAX_DIGS); + t2 = newgiant(MAX_DIGS); + t3 = newgiant(MAX_DIGS); + t4 = newgiant(MAX_DIGS); + t5 = newgiant(MAX_DIGS); + cubic = newpoly(4); + p = newgiant(MAX_DIGS); + a = newgiant(MAX_DIGS); + b = newgiant(MAX_DIGS); + magcheck = newgiant(MAX_DIGS); + precip = newpoly(MAX_COEFFS); +} + +void +atoa(giant *a, giant *b, int n) { + int j; + for(j=0; jdeg = x->deg; + atoa(x->coe, y->coe, y->deg+1); +} + +void +negp(poly y) +/* y := -y. */ +{ int j; + for(j=0; j <= y->deg; j++) { + negg(y->coe[j]); + quickmodg(p, y->coe[j]); + } +} + +int +iszer(giant a) { + + if(a->sign == 0) return(1); + return(0); + +} + +int +iszerop(poly x) { + if(x->deg == 0 && iszer(x->coe[0])) return 1; + return 0; +} + +int +is0(giant a) { + return(iszer(a)); +} + +int +is1(giant a) { + return(isone(a)); +} + + +void +addp(poly x, poly y) +/* y += x. */ +{ + int d = x->deg, j; + + if(y->deg > d) d = y->deg; + for(j = 0; j <= d; j++) { + if((j <= x->deg) && (j <= y->deg)) { + addg(x->coe[j], y->coe[j]); + quickmodg(p, y->coe[j]); + continue; + } + if((j <= x->deg) && (j > y->deg)) { + gtog(x->coe[j], y->coe[j]); + quickmodg(p, y->coe[j]); + continue; + } + } + y->deg = d; + justifyp(y); +} + +void +subp(poly x, poly y) +/* y -= x. */ +{ + int d = x->deg, j; + + if(y->deg > d) d = y->deg; + for(j = 0; j <= d; j++) { + if((j <= x->deg) && (j <= y->deg)) { + subg(x->coe[j], y->coe[j]); + quickmodg(p, y->coe[j]); + continue; + } + if((j <= x->deg) && (j > y->deg)) { + gtog(x->coe[j], y->coe[j]); + negg(y->coe[j]); + quickmodg(p, y->coe[j]); + continue; + } + } + y->deg = d; + justifyp(y); +} + + +void +grammarmulp(poly a, poly b) +/* b *= a. */ +{ + int dega = a->deg, degb = b->deg, deg = dega + degb; + register int d, kk, first, diffa; + + for(d=deg; d>=0; d--) { + diffa = d-dega; + itog(0, coe); + for(kk=0; kk<=d; kk++) { + if((kk>degb)||(kkcoe[kk], tmp); + mulg(a->coe[d-kk], tmp); + modg(p, tmp); + addg(tmp, coe); + quickmodg(p, coe); + } + gtog(coe, mcand[d]); + } + atoa(mcand, b->coe, deg+1); + b->deg = deg; + justifyp(b); +} + +void +atop(giant *a, poly x, int deg) +/* Copy array to poly, with monic option. */ +{ + int adeg = abs(deg); + x->deg = adeg; + atoa(a, x->coe, adeg); + if(deg < 0) { + itog(1, x->coe[adeg]); + } else { + gtog(a[adeg], x->coe[adeg]); + } +} + +void +just(giant g) { + while((g->n[g->sign-1] == 0) && (g->sign > 0)) --g->sign; +} + +void +unstuff_partial(giant g, poly y, int words){ + int j; + for(j=0; j < y->deg; j++) { + bcopy((g->n) + j*words, y->coe[j]->n, words*sizeof(short)); + y->coe[j]->sign = words; + just(y->coe[j]); + } +} + +void +stuff(poly x, giant g, int words) { + int deg = x->deg, j, coedigs; + + g->sign = words*(1 + deg); + for(j=0; j <= deg; j++) { + coedigs = (x->coe[j])->sign; + bcopy(x->coe[j]->n, (g->n) + j*words, coedigs*sizeof(short)); + bzero((g->n) + (j*words+coedigs), + sizeof(short)*(words-coedigs)); + } + just(g); +} + +int maxwords = 0; +void + +binarysegmul(poly x, poly y) { + int bits = bitlen(p), xwords, ywords, words; + + xwords = (2*bits + log_2(x->deg+1) + 32 + 15)/16; + ywords = (2*bits + log_2(y->deg+1) + 32 + 15)/16; + if(xwords > ywords) words = xwords; else words = ywords; +if(words > maxwords) { + maxwords = words; +// printf("m: %d\n", words); + fflush(stdout); +} + stuff(x, globx, words); + stuff(y, globy, words); + mulg(globx, globy); + gtog(y->coe[y->deg], globx); /* Save high coeff. */ + y->deg += x->deg; + gtog(globx, y->coe[y->deg]); /* Move high coeff. */ + unstuff_partial(globy, y, words); + mulg(x->coe[x->deg], y->coe[y->deg]); /* resolve high coeff. */ + justifyp(y); +} + +binarysegsquare(poly y) { + int bits = bitlen(p), words; + words = (2*bits + log_2(y->deg+1) + 32 + 15)/16; + stuff(y, globy, words); + squareg(globy); + gtog(y->coe[y->deg], globx); /* Save high coeff. */ + y->deg += y->deg; + gtog(globx, y->coe[y->deg]); /* Move high coeff. */ + unstuff_partial(globy, y, words); + mulg(y->coe[y->deg], y->coe[y->deg]); /* resolve high coeff. */ + justifyp(y); +} + + +void +assess(poly x, poly y){ + int max = 0, j; + for(j=0; j <= x->deg; j++) { + if(bitlen(x->coe[j]) > max) max = bitlen(x->coe[j]); + } +// printf("max: %d ", max); + max = 0; + for(j=0; j <= y->deg; j++) { + if(bitlen(y->coe[j]) > max) max = bitlen(y->coe[j]); + } +// printf("%d\n", max); + + +} + +int +pcompp(poly x, poly y) { + int j; + if(x->deg != y->deg) return 1; + for(j=0; j <= x->deg; j++) { + if(gcompg(x->coe[j], y->coe[j])) return 1; + } + return 0; +} + +int maxdeg = 0; + +void +mulp(poly x, poly y) +/* y *= x. */ +{ + int n, degx = x->deg, degy = y->deg; + +/* +if(degx > max_deg) { + max_deg = degx; printf("xdeg: %d\n", degx); +} + +if(degy > max_deg) { + max_deg = degy; printf("ydeg: %d\n", degy); +} +*/ + if((degx < P_BREAK) || (degy < P_BREAK)) { + grammarmulp(x,y); + justifyp(y); + return; + } + if(x==y) binarysegsquare(y); + else binarysegmul(x, y); +} + +void +revp(poly x) +/* Reverse the coefficients of x. */ +{ int j, deg = x->deg; + + for(j=0; j <= deg/2; j++) { + gtog(x->coe[deg-j], tmp); + gtog(x->coe[j], x->coe[deg-j]); + gtog(tmp, x->coe[j]); + } + justifyp(x); +} + +void +recipp(poly f, int deg) +/* f := 1/f, via newton method. */ +{ + int lim = deg + 1, prec = 1; + + sscratch->deg = 0; itog(1, sscratch->coe[0]); + itog(1, aux); + while(prec < lim) { + prec <<= 1; + if(prec > lim) prec = lim; + f->deg = prec-1; + ptop(sscratch, tscratch); + mulp(f, tscratch); + tscratch->deg = prec-1; + polyrem(tscratch); + subg(aux, tscratch->coe[0]); + quickmodg(p, tscratch->coe[0]); + mulp(sscratch, tscratch); + tscratch->deg = prec-1; + polyrem(tscratch); + subp(tscratch, sscratch); + sscratch->deg = prec-1; + polyrem(sscratch); + } + justifyp(sscratch); + ptop(sscratch, f); +} + +int +left_justifyp(poly x, int start) +/* Left-justify the polynomial, checking from position "start." */ +{ + int j, shift = 0; + + for(j = start; j <= x->deg; j++) { + if(!is0(x->coe[j])) { + shift = start; + break; + } + } + x->deg -= shift; + for(j=0; j<= x->deg; j++) { + gtog(x->coe[j+shift], x->coe[j]); + } + return(shift); +} + +void +remp(poly x, poly y, int mode) +/* y := x (mod y). + mode = 0 is normal operation, + = 1 jams a fixed reciprocal, + = 2 uses the fixed reciprocal. + */ +{ + int degx = x->deg, degy = y->deg, d, shift; + + if(degy == 0) { + y->deg = 0; + itog(0, y->coe[0]); + return; + } + d = degx - degy; + if(d < 0) { + ptop(x, y); + return; + } + revp(x); revp(y); + ptop(y, rscratch); + switch(mode) { + case 0: recipp(rscratch, d); + break; + case 1: recipp(rscratch, degy); /* degy -1. */ + ptop(rscratch, precip); + rscratch->deg = d; justifyp(rscratch); + break; + case 2: ptop(precip, rscratch); + rscratch->deg = d; justifyp(rscratch); + break; + } +/* Next, a limited-precision multiply. */ + if(d < degx) { x->deg = d; justifyp(x);} + mulp(x, rscratch); + rscratch->deg = d; + polyrem(rscratch); + justifyp(rscratch); + x->deg = degx; justifyp(x); + mulp(rscratch, y); + subp(x, y); + negp(y); polyrem(y); + shift = left_justifyp(y, d+1); + for(d = y->deg+1; d <= degx-shift; d++) itog(0, y->coe[d]); + y->deg = degx - shift; + revp(y); + revp(x); +} + +fullmod(poly x) { + polyrem(x); + ptop(smonic, s[0]); + remp(x, s[0], 2); + ptop(s[0], x); + polyrem(x); +} + +scalarmul(giant s, poly x) { + int j; + for(j=0; j <= x->deg; j++) { + mulg(s, x->coe[j]); + modg(p, x->coe[j]); + } +} + + +schain(int el) { + int j; + + s[0]->deg = 0; + itog(0, s[0]->coe[0]); + + s[1]->deg = 0; + itog(1, s[1]->coe[0]); + s[2]->deg = 0; + itog(2, s[2]->coe[0]); + + s[3]->deg = 4; + gtog(a, aux); mulg(a, aux); negg(aux); + gtog(aux, s[3]->coe[0]); + gtog(b, aux); smulg(12, aux); + gtog(aux, s[3]->coe[1]); + gtog(a, aux); smulg(6, aux); + gtog(aux, s[3]->coe[2]); + itog(0, s[3]->coe[3]); + itog(3, s[3]->coe[4]); + + s[4]->deg = 6; + gtog(a, aux); mulg(a, aux); mulg(a, aux); + gtog(b, tmp); mulg(b, tmp); smulg(8, tmp); addg(tmp, aux); + negg(aux); + gtog(aux, s[4]->coe[0]); + gtog(b, aux); mulg(a, aux); smulg(4, aux); negg(aux); + gtog(aux, s[4]->coe[1]); + gtog(a, aux); mulg(a, aux); smulg(5, aux); negg(aux); + gtog(aux, s[4]->coe[2]); + gtog(b, aux); smulg(20, aux); + gtog(aux, s[4]->coe[3]); + gtog(a, aux); smulg(5, aux); + gtog(aux, s[4]->coe[4]); + itog(0, s[4]->coe[5]); + itog(1, s[4]->coe[6]); + itog(4, aux); + scalarmul(aux, s[4]); + cubic->deg = 3; + itog(1, cubic->coe[3]); + itog(0, cubic->coe[2]); + gtog(a, cubic->coe[1]); + gtog(b, cubic->coe[0]); + for(j=5; j <= el; j++) { + if(j % 2 == 0) { + ptop(s[j/2 + 2], s[j]); mulp(s[j/2-1], s[j]); + polyrem(s[j]); mulp(s[j/2-1], s[j]); polyrem(s[j]); + ptop(s[j/2-2], s[0]); mulp(s[j/2+1], s[0]); polyrem(s[0]); + mulp(s[j/2+1], s[0]); polyrem(s[0]); + subp(s[0], s[j]); mulp(s[j/2], s[j]); polyrem(s[j]); + gtog(p, aux); iaddg(1, aux); gshiftright(1, aux); + scalarmul(aux, s[j]); + } else { + ptop(s[(j-1)/2+2], s[j]); + mulp(s[(j-1)/2], s[j]); polyrem(s[j]); + mulp(s[(j-1)/2], s[j]); polyrem(s[j]); + mulp(s[(j-1)/2], s[j]); polyrem(s[j]); + ptop(s[(j-1)/2-1], s[0]); + mulp(s[(j-1)/2 + 1], s[0]); polyrem(s[0]); + mulp(s[(j-1)/2 + 1], s[0]); polyrem(s[0]); + mulp(s[(j-1)/2 + 1], s[0]); polyrem(s[0]); + if(((j-1)/2) % 2 == 1) { + mulp(cubic, s[0]); polyrem(s[0]); + mulp(cubic, s[0]); polyrem(s[0]); + } else { + mulp(cubic, s[j]); polyrem(s[j]); + mulp(cubic, s[j]); polyrem(s[j]); + } +// polyout(s[1]); polyout(s[3]); polyout(s[0]); polyout(s[j]); + subp(s[0], s[j]); + polyrem(s[j]); + } + } +} + +init_recip(int el) { + int j; + ptop(s[el], smonic); + if(el == 2) { + mulp(cubic, smonic); polyrem(smonic); + } + gtog(smonic->coe[smonic->deg], aux); /* High coeff. */ + binvg(p, aux); + scalarmul(aux, smonic); /* s is now monic. */ + s[0]->deg = smonic->deg + 1; + for(j=0; j <= s[0]->deg; j++) itog(1, s[0]->coe[j]); + ptop(smonic, pbuff); + remp(s[0], pbuff, 1); /* Initialize reciprocal of s as precip. */ +} + +void powerpoly(poly x, giant n) +/* Base-4 window. */ +{ int pos, code; + ptop(x, pbuff); /* x. */ + ptop(pbuff, pbuff2); + mulmod(pbuff2, pbuff2); mulmod(pbuff, pbuff2); /* x^3. */ + pos = bitlen(n)-2; + while(pos >= 0) { + mulmod(x, x); + if(pos==0) { + if(bitval(n, pos) != 0) { + mulmod(pbuff, x); + } + break; + } + code = (bitval(n, pos) != 0) * 2 + (bitval(n, pos-1) != 0); + switch(code) { + case 0: mulmod(x,x); break; + case 1: mulmod(x,x); + mulmod(pbuff, x); + break; + case 2: mulmod(pbuff, x); + mulmod(x,x); break; + case 3: mulmod(x,x); mulmod(pbuff2, x); break; + } + pos -= 2; + } +} + +mulmod(poly x, poly y) { + mulp(x, y); fullmod(y); +} + +elldoublep(poly n1, poly d1, poly m1, poly c1, poly n0, poly d0, + poly m0, poly c0) { + + ptop(n1, mn); mulmod(n1, mn); + ptop(mn, pbuff); addp(mn, mn); addp(pbuff, mn); + fullmod(mn); + ptop(d1, pbuff); mulmod(d1, pbuff); + scalarmul(a, pbuff); addp(pbuff, mn); + fullmod(mn); + mulmod(c1, mn); + ptop(m1, md); addp(md, md); + mulmod(d1, md); mulmod(d1, md); mulmod(cubic, md); + + ptop(d1, n0); mulmod(mn, n0); mulmod(mn, n0); + mulmod(cubic, n0); + ptop(n1, pbuff); addp(pbuff, pbuff); fullmod(pbuff); + mulmod(md, pbuff); mulmod(md, pbuff); + subp(pbuff, n0); fullmod(n0); + ptop(md, d0); mulmod(md, d0); mulmod(d1, d0); + + ptop(mn, m0); mulmod(c1, m0); + ptop(d0, pbuff); mulmod(n1, pbuff); + ptop(n0, c0); mulmod(d1, c0); subp(c0, pbuff); + fullmod(pbuff); + mulmod(pbuff, m0); + ptop(m1, pbuff); mulmod(md, pbuff); + mulmod(d1, pbuff); mulmod(d0, pbuff); + subp(pbuff, m0); fullmod(m0); + + ptop(c1, c0); mulmod(md, c0); mulmod(d1, c0); mulmod(d0, c0); +} + +elladdp(poly n1, poly d1, poly m1, poly c1, poly n2, poly d2, poly m2, poly c2, poly n0, poly d0, poly m0, poly c0) { + ptop(m2, mn); mulmod(c1, mn); + ptop(m1, pbuff); mulmod(c2, pbuff); + subp(pbuff, mn); fullmod(mn); + mulmod(d1, mn); mulmod(d2, mn); + + ptop(n2, md); mulmod(d1, md); + ptop(n1, pbuff); mulmod(d2, pbuff); + subp(pbuff, md); fullmod(md); + mulmod(c1, md); mulmod(c2, md); + + ptop(cubic, n0); mulmod(mn, n0); mulmod(mn, n0); + mulmod(d1, n0); mulmod(d2, n0); + ptop(n1, pbuff); mulmod(d2, pbuff); + ptop(n2, d0); mulmod(d1, d0); + addp(d0, pbuff); mulmod(md, pbuff); mulmod(md, pbuff); + subp(pbuff, n0); fullmod(n0); + + ptop(md, d0); mulmod(md, d0); mulmod(d1, d0); mulmod(d2, d0); + + ptop(mn, m0); mulmod(c1, m0); + ptop(d0, pbuff); mulmod(n1, pbuff); + ptop(d1, c0); mulmod(n0, c0); + subp(c0, pbuff); fullmod(pbuff); + mulmod(pbuff, m0); + ptop(m1, pbuff); mulmod(md, pbuff); + mulmod(d0, pbuff); mulmod(d1, pbuff); + subp(pbuff, m0); fullmod(m0); + + ptop(md, c0); mulmod(c1, c0); mulmod(d0, c0); mulmod(d1, c0); + +} + +polyout(poly x) { + int j; + for(j=0; j <= x->deg; j++) {printf("%d: ",j); gout(x->coe[j]);} +} + +#define SIEVE_CUT 7 + +main(int argc, char **argv) { + int j, ct, el, xmatch, ymatch; + int k, k2, t, linit; + int T[L_MAX], P[L_MAX]; + giant ss[L_MAX]; + unsigned int ord, ordminus; + + printf("Initializing...\n"); fflush(stdout); + init_all(); + + for(j=0; j < L_MAX; j++) { /* Array to be used for CRT reconstruction. */ + ss[j] = NULL; + } + + printf("Give p, a, b on separate lines:\n"); fflush(stdout); + gin(p); /* Field prime. */ + if(bitlen(p) > Q) { + fprintf(stderr, "p too large, larger than %d bits.\n", Q); + exit(0); + } + + gin(a); gin(b); /* Curve parameters. */ + +newb: + printf("Checking discriminant for:\nb = "); + gout(b); + gtog(a, t1); squareg(t1); modg(p, t1); mulg(a, t1); modg(p, t1); + gshiftleft(2, t1); /* 4 a^3. */ + gtog(b, t2); squareg(t2); modg(p, t2); smulg(27, t2); + addg(t2, t1); modg(p, t1); + if(isZero(t1)) { + printf("Discriminant FAILED\n"); + iaddg(1, b); + goto newb; + } + printf("Discriminant PASSED\n"); + printf("Starting sieve process:\n"); + + + schain(SIEVE_CUT + 1); /* Do first piece of chain, for small-sieving. */ + linit = 0; + ct = 0; + itog(1, magcheck); +for(el = 2; el <= L_MAX; el += 1 + (el%2)) { + if(!primeq(el)) continue; + L = el; while(L*el <= 4) L *= el; +printf("Resolving Schoof L = %d...\n", L); + P[ct] = L; /* Stuff another prime power. */ + smulg(L, magcheck); + gtog(magcheck, tmp); squareg(tmp); gshiftright(2, tmp); + if(gcompg(tmp, p) > 0) break; /* Done...go to CRT phase. */ +if((L > SIEVE_CUT) && (linit == 0)) { + schain(L_MAX+1); + linit = 1; + } + init_recip(L); +// printf("s: "); polyout(s[L]); + gtog(p, aux2); + k = idivg(L, aux2); /* p (mod L). */ + gtog(p, aux2); + k2 = idivg(el, aux2); +// printf("power...\n"); + txd->deg = 0; itog(1, txd->coe[0]); + tyd->deg = 0; itog(1, tyd->coe[0]); + txn->deg = 1; itog(0, txn->coe[0]); itog(1, txn->coe[1]); + ptop(txn, kxn); + + gtog(p, aux2); + powerpoly(txn, aux2); /* x^p. */ +// printf("x^p done...\n"); + ptop(txn, powx); + powerpoly(powx, aux2); +// printf("x^p^2 done...\n"); + ptop(cubic, tyn); + gtog(p, aux2); itog(1, aux); subg(aux, aux2); + gshiftright(1, aux2); /* aux2 := (p-1)/2. */ + powerpoly(tyn, aux2); /* y^p. */ +// printf("y^p done...\n"); + ptop(tyn, powy); mulp(tyn, powy); fullmod(powy); + mulp(cubic, powy); fullmod(powy); + powerpoly(powy, aux2); + mulp(tyn, powy); fullmod(powy); +// printf("Powers done...\n"); + +// printf("pow" ); polyout(powx); polyout(powy); + ptop(txn, txn1); ptop(txd, txd1); /* Save t = 1 case. */ + ptop(tyn, tyn1); ptop(tyd, tyd1); +/* We now shall test + (powx, y powy) + k(kxn/kxd, y kyn/kyd) = t(txn/txd, y tyn/tyd) + */ + + if(k==1) { ptop(txd, kxd); ptop(txd, kyd); + ptop(txd, kyn); + } else { + ptop(s[k], kxd); mulp(s[k], kxd); fullmod(kxd); + if(k%2==0) { mulp(cubic, kxd); fullmod(kxd); } + mulp(kxd, kxn); fullmod(kxn); + ptop(s[k-1], pbuff); mulp(s[k+1], pbuff); fullmod(pbuff); + if(k%2==1) {mulp(cubic, pbuff); fullmod(pbuff); } + subp(pbuff, kxn); fullmod(kxn); + + ptop(s[k+2], kyn); mulp(s[k-1], kyn); fullmod(kyn); + mulp(s[k-1], kyn); fullmod(kyn); + if(k > 2) { + ptop(s[k-2], pbuff); mulp(s[k+1], pbuff); fullmod(pbuff); + mulp(s[k+1], pbuff); fullmod(pbuff); + subp(pbuff, kyn); fullmod(kyn); + } + ptop(s[k], kyd); mulp(s[k], kyd); fullmod(kyd); + mulp(s[k], kyd); fullmod(kyd); + if(k%2==0) { mulp(cubic, kyd); fullmod(kyd); + mulp(cubic, kyd); fullmod(kyd);} + itog(4, aux2); scalarmul(aux2, kyd); + } +//printf("kP: "); polyout(kxn); polyout(kxd); polyout(kyn); polyout(kyd); +/* Commence t = 0 check. */ +// printf("Checking t = %d ...\n", 0); +fflush(stdout); + + ptop(powx, pbuff); mulp(kxd, pbuff); + subp(kxn, pbuff); + fullmod(pbuff); + + xmatch = ymatch = 0; + if(iszerop(pbuff)) { + xmatch = 1; + /* Now check y coords. */ + if(L == 2) goto resolve; + ptop(powy, pbuff); mulp(kyd, pbuff); + addp(kyn, pbuff); + fullmod(pbuff); + if(iszerop(pbuff)) { + resolve: + printf("%d %d\n", L, 0); + T[ct++] = 0; + if((k2 + 1 - T[ct-1]) % el == 0) { + printf("TILT: %d\n", el); + el = 2; + iaddg(1, b); + goto newb; + } + continue; + } else ymatch = 1; + } +/* Combine pt1 and pt2. */ + if((xmatch == 1) && (ymatch == 1)) + elldoublep(powx, txd, powy, txd, nx, dx, ny, dy); + else + elladdp(powx, txd, powy, txd, kxn, kxd, kyn, kyd, nx, dx, ny, dy); + +/* Now {nx/dx, ny/dy} is (fixed) LHS. */ +// printf("add12: "); polyout(nx); polyout(dx); polyout(ny); polyout(dy); +/* Commence t > 0 check. */ + for(t=1; t <= L/2; t++) { +// printf("Checking t = %d ...\n", t); + if(t > 1) { /* Add (tx1, ty1) to (txn, tyn). */ + ptop(txn1, pbuff); mulmod(txd, pbuff); + ptop(txn, powx); mulmod(txd1, powx); + subp(powx, pbuff); fullmod(pbuff); + if(!iszerop(pbuff)) + elladdp(txn1, txd1, tyn1, tyd1, txn, txd, tyn, tyd, + tmp1, tmp2, tmp3, tmp4); + else elldoublep(txn, txd, tyn, tyd, + tmp1, tmp2, tmp3, tmp4); + ptop(tmp1, txn); ptop(tmp2, txd); + ptop(tmp3, tyn); ptop(tmp4, tyd); + } +// printf("tQ: "); polyout(txn); polyout(txd); polyout(tyn); polyout(tyd); + /* Next, check {nx/dx, ny/dy} =? {txn/txd, tyn/tyd}. */ + ptop(nx, pbuff); mulmod(txd, pbuff); + ptop(dx, powx); mulmod(txn, powx); + subp(powx, pbuff); fullmod(pbuff); + if(!iszerop(pbuff)) continue; + /* Next, check y. */ + // printf("y check!\n"); + ptop(ny, pbuff); mulmod(tyd, pbuff); + ptop(dy, powx); mulmod(tyn, powx); + subp(powx, pbuff); fullmod(pbuff); + if(iszerop(pbuff)) { + printf("%d %d\n", L, t); + T[ct++] = t; + } else { + printf("%d %d\n", L, L-t); + T[ct++] = L-t; + } + if((k2 + 1 - T[ct-1]) % el == 0) { + printf("TILT: %d\n", el); + el = 2; + iaddg(1, b); + goto newb; + } + + fflush(stdout); + break; + } +} + +/* Now, prime powers P[] and CRT residues T[] are intact. */ + printf("Prime powers L:\n"); + printf("{"); + for(j=0; j < ct-1; j++) { + printf("%d, ", P[j]); + } + printf("%d }\n", P[ct-1]); + + printf("Residues t (mod L):\n"); + printf("{"); + for(j=0; j < ct-1; j++) { + printf("%d, ", T[j]); + } + printf("%d }\n", T[ct-1]); + +/* Mathematica algorithm for order: +plis = {2^5, 3^3, 5^2, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47}; +tlis = {1, 26, 4, 2, 4, 11, 6, 5, 19, 22, 10, 16, 7, 22, 11}; +prod = Apply[Times, plis]; +prlis = prod/plis; +invlis = Table[PowerMod[prlis[[q]], -1, plis[[q]]],{q,1,Length[plis]}]; +p = 2^127 - 1; +t = Mod[tlis . (prlis * invlis), prod]; +ord = p + 1 - If[t^2 > 4p, t - prod, t] +*/ + + itog(1, t1); + for(j=0; j < ct; j++) { + smulg(P[j], t1); + } + + for(j=0; j < 2*ct; j++) { + if(!ss[j]) ss[j] = newgiant(MAX_DIGS); + } + + for(j=0; j < ct; j++) { + gtog(t1, ss[j]); + itog(P[j], t2); + divg(t2, ss[j]); + } + + for(j=0; j < ct; j++) { + gtog(ss[j], ss[j+ct]); + itog(P[j], t2); + invg(t2, ss[j+ct]); + } + + itog(0, t4); + for(j=0; j < ct; j++) { + itog(T[j], t5); + mulg(ss[j], t5); + mulg(ss[j+ct], t5); + addg(t5, t4); + } + modg(t1, t4); + gtog(p, t5); + iaddg(1, t5); + gtog(t4, t2); + squareg(t4); + gtog(p, t3); gshiftleft(2, t3); + if(gcompg(t4, t3) > 0) subg(t1, t2); + subg(t2, t5); + printf("Parameters:\n"); + printf("p = "); gout(p); + printf("a = "); gout(a); + printf("b = "); gout(b); + printf("Curve order:\n"); + printf("o = "); gout(t5); + printf("pprob: %d\n", prime_probable(t5)); + printf("Twist order:\n"); + printf("o' = "); + addg(t2, t5); + addg(t2, t5); + gout(t5); + printf("pprob: %d\n", prime_probable(t5)); + + iaddg(1,b); + goto newb; +}