X-Git-Url: https://git.saurik.com/apple/security.git/blobdiff_plain/80e2389990082500d76eb566d4946be3e786c3ef..d8f41ccd20de16f8ebe2ccc84d47bf1cb2b26bbb:/Security/libsecurity_cryptkit/lib/CurveParamDocs/schoof.c diff --git a/Security/libsecurity_cryptkit/lib/CurveParamDocs/schoof.c b/Security/libsecurity_cryptkit/lib/CurveParamDocs/schoof.c new file mode 100644 index 00000000..8c950d38 --- /dev/null +++ b/Security/libsecurity_cryptkit/lib/CurveParamDocs/schoof.c @@ -0,0 +1,1100 @@ +/* schoof.c + + Elliptic curve order calculator, for + + y^2 = x^3 + a x + b (mod p) + + Compile with: + + % cc -O schoof.c giants.c tools.c ellproj.c -lm -o schoof + + and run, entering p and the a,b parameters. + Eventually the curve order o is reported, together with + the twist order o'. The sum of these two orders is always + 2p + 2. A final check is performed to verify that a + random point (x,y,z) enjoys the annihilation + o * (x,y,z) = point-at-infinity. This is not absolutely + definitive, rather it is a necessary condition on the oder o + (i.e. it is a sanity check of sorts). + + * Change history: + + 18 Nov 99 REC fixed M. Scott bug (MAX_DIGS bumped by 2) + 5 Jul 99 REC Installed improved (base-4) power ladder + 5 Jul 99 REC Made use of binary segmentation square (per se) + 5 Jul 99 REC Improved memory usage + 2 May 99 REC Added initial primality check + 30 Apr 99 REC Added actual order-annihilation test + 29 Apr 99 REC Improvements due to A. Powell + 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 +#include "giants.h" +#include "tools.h" +#include "ellproj.h" + +#define P_BREAK 32 + +#ifdef _WIN32 +#include +#define bzero(D, n) memset(D, 0x00, n) +#define bcopy(S, D, n) memcpy(D, S, n) +#endif + +#define Q_MAX 256 /* Bits in largest primes handled. */ +#define L_CEIL 100 /* Bound on Schoof L values (not all needed in general). */ + + +typedef struct + { + int deg; + giant *coe; + } polystruct; +typedef polystruct *poly; + +extern int *pr; + +static int Q, L_MAX; +static int MAX_DIGS, MAX_COEFFS; + +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_CEIL+2], smonic; +static giant p, a, b; +static 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; + j = j * MAX_COEFFS; + globx = newgiant(j); + globy = newgiant(j); + s[0] = newpoly(MAX_COEFFS); + + for(j=1; j<=L_MAX+1; j++) { + s[j] = newpoly(j*(j+1)); + } + smonic = newpoly(MAX_COEFFS); + 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); +/* The next three need extra space for remaindering routines. */ + 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); + t3 = newgiant(MAX_DIGS); + t4 = newgiant(MAX_DIGS); + t5 = newgiant(MAX_DIGS); + cubic = newpoly(4); + 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; + 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]); + } + max = 0; + for(j=0; j <= y->deg; j++) { + if(bitlen(y->coe[j]) > max) max = bitlen(y->coe[j]); + } +} + +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 max_deg = 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) +{ int len, pos; + ptop(x, pbuff); + x->deg = 0; itog(1, x->coe[0]); + len = bitlen(n); + pos = 0; + while(1) { + if(bitval(n, pos++)) { + mulp(pbuff, x); + fullmod(x); + } + if(pos>=len) break; + mulp(pbuff, pbuff); + fullmod(pbuff); + } +} +*/ + +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]);} +} + +main(int argc, char **argv) { + int j, ct = 0, el, xmatch, ymatch; + int k, t; + int T[L_CEIL], P[L_CEIL], LL[L_CEIL]; + giant ss[L_CEIL]; + unsigned int ord, ordminus; + point_proj pt, pt2; + + p = newgiant(INFINITY); /* Also sets up internal giants' stacks. */ + j = ((Q_MAX+15)/16); + init_tools(2*j); + a = newgiant(j); + b = newgiant(j); + +entry: + printf("Give p > 3, a, b on separate lines:\n"); fflush(stdout); + gin(p); /* Field prime. */ + if((Q = bitlen(p)) > Q_MAX) { + fprintf(stderr, "p too large, larger than %d bits.\n", Q); + goto entry; + } + if(!prime_probable(p)) { + fprintf(stderr, "p is not but must be prime.\n", Q); + goto entry; + } + gin(a); gin(b); /* Curve parameters. */ + + t1 = newgiant(2*j); + t2 = newgiant(2*j); +/* Next, discriminant test for legitimacy of curve. */ + gtog(a, t1); squareg(t1); modg(p, t1); mulg(a, t1); modg(p, t1); + gshiftleft(2, t1); /* 4 a^3 mod p. */ + gtog(b, t2); squareg(t2); modg(p, t2); smulg(27, t2); + addg(t2, t1); modg(p, t1); + if(isZero(t1)) { + fprintf(stderr, "Discriminant FAILED\n"); + goto entry; + } + printf("Discriminant PASSED\n"); fflush(stdout); + +/* Next, find an efficient prime power array such that + Prod[powers] >= 16 p. */ + + /* Create minimal prime power array such that Prod[powers]^2 > 16p. */ + + gtog(p, t2); gshiftleft(4, t2); /* t2 := 16p. */ + + L_MAX = 3; + while(L_MAX <= L_CEIL-1) { + for(j=0; j <= L_MAX; j++) LL[j] = 0; + for(j=2; j <= L_MAX; j++) { + if(primeq(j)) { + LL[j] = 1; + k = j; + while(1) { + k *= j; + if(k <= L_MAX) { + LL[k] = 1; + LL[k/j] = 0; + } + else break; + } + } + } + itog(1, t1); + for(j=2; j <= L_MAX; j++) { + if(LL[j]) { smulg(j, t1); smulg(j, t1); } /* Building up t1^2. */ + } + if(gcompg(t1, t2) > 0) break; + ++L_MAX; + } + + printf("Initializing for prime powers:\n"); + for(j=2; j <= L_MAX; j++) { + if(LL[j]) printf("%d ", j); + } + printf("\n"); + fflush(stdout); + + + MAX_DIGS = (2 + (Q+15)/8); /* Size of (squared) coefficients. */ + MAX_COEFFS = ((L_MAX+1)*(L_MAX+2)); + + init_all(); + schain(L_MAX+1); + +for(L = 2; L <= L_MAX; L++) { + if(!LL[L]) continue; +printf("Resolving Schoof L = %d...\n", L); + P[ct] = L; /* Stuff another prime power. */ + init_recip(L); +// printf("s: "); polyout(s[L]); + gtog(p, aux2); + k = idivg(L, aux2); /* p (mod L). */ + +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; + 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; + } + 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++) { + free(s[j]); /* Just to clear memory. */ + smulg(P[j], t1); + } + + for(j=0; j < 2*ct; 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); gtog(t5, t3); /* Save order as t3. */ + printf("Twist order:\n"); + printf("o' = "); + addg(t2, t5); + addg(t2, t5); + gout(t5); +/* Next, verify the order. */ + printf("Verifying order o:...\n"); + init_ell_proj(MAX_DIGS); + pt = new_point_proj(MAX_DIGS); + pt2 = new_point_proj(MAX_DIGS); + itog(1,t2); + find_point_proj(pt, t2, a, b, p); + printf("A point on the curve y^2 = x^3 + a x + b (mod p) is:\n"); + printf("(x,y,z) = {\n"); gout(pt->x); printf(","); + gout(pt->y); printf(","); gout(pt->z); + printf("}\n"); + ell_mul_proj(pt, pt2, t3, a, p); + printf("A multiple is:\n"); + printf("o * (x,y,z) = {\n"); + gout(pt2->x); printf(",");gout(pt2->y); printf(",");gout(pt2->z); + printf("}\n"); + if(!isZero(pt2->z)) { + printf("TILT: multiple should be point-at-infinity.\n"); + exit(1); + } + printf("VERIFIED: multiple is indeed point-at-infinity.\n"); +}