--- /dev/null
+/* 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 <stdio.h>
+#include<assert.h>
+#include <math.h>
+#include <stdlib.h>
+#include "giants.h"
+#include "tools.h"
+#include "ellproj.h"
+
+#define P_BREAK 32
+
+#ifdef _WIN32
+#include <string.h>
+#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; j<n; j++) {
+ p[j] = newgiant(MAX_DIGS);
+ }
+ return(p);
+}
+
+poly
+newpoly(int coeffs) {
+ poly pol;
+ pol = (poly) malloc(sizeof(polystruct));
+ pol->coe = (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; j<n; j++) gtog(a[j], b[j]);
+}
+
+void
+ptop(poly x, poly y)
+/* y := x. */
+{
+ y->deg = 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)||(kk<diffa)) continue;
+ gtog(b->coe[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");
+}
projects / apple / security.git / blobdiff