schoofs.c   [plain text]

/* 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 <stdio.h>
#include<assert.h>
#include <math.h>
#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; 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;
   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; 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;
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;
}