www.pudn.com > miracl.zip > MODPOL.CPP
//
// Program to generate Modular Polynomials mod p, as required for fast
// implementations of the Schoof-Elkies-Atkins algorithm
// for counting points on an elliptic curve Y^2=X^3 + A.X + B mod p
//
// Implemented entirely from the description provided in:
// 1. "Distributed Computation of the number of points on an elliptic curve
// over a finite prime field", Buchmann, Mueller, & Shoup, SFB 124-TP D5
// Report 03/95, April 1995, Universitat des Saarlandes, and
// 2. "Counting the number of points on elliptic curves over finite fields
// of characteristic greater than three", Lehmann, Maurer, Mueller & Shoup,
// Proc. 1st Algorithmic Number Theory Symposium (ANTS), pp 60-70, 1994
//
// Both papers are available on the Web from Volker Mueller's home page
// www.informatik.th-darmstadt.de/TI/Mitarbeiter/vmueller.html
//
// Also strongly recommended is the book
//
// 3. "Elliptic Curves in Cryptography"
// by Blake, Seroussi and Smart, London Mathematical Society Lecture Note
// Series 265, Cambridge University Press. ISBN 0 521 65374 6
//
// The programs's output for each prime in the range is a bivariate polynomial
// in X and Y, which can optionally be stored to disk. Some informative output
// is generated just to convince you that it is still working, and to give an
// idea of progress.
//
// This program is a composite of the "mueller" and "process" applications.
// It generates the modular polynomials, pre-reduced wrt to a specified prime
// modulus. This may be the only feasible way to do it on a small computer
// system, for which the "mueller" application is too resource intensive.
//
// Although less memory intensive than "mueller", problems may still arise.
// See mueller.cpp for a description of the -s2, -s3 and -s6 flags
//
// .pol file format
// ,,<1st coef>,<1st power of X>,<1st power of Y>,<2nd coef>...
// Each polynomial ends wih zero powers of X and Y.
//
// For example
// modpol -d -f 2#512 0 500 -o test512.pol
//
// If appending to a file with the -a flag, make sure and use the same
// prime modulus as used to create the file originally - no check is made
//
// generates the test512.pol file directly, given the range 0 to 500 and
// using the first prime modulus it can find less than 2^512. This file can
// then be used directly with the "sea" application
//
// Copyright Shamus Software Ltd., 1999
//
#include
#include
#include
#include
#include "ps_zzn.h" // power series class
using namespace std;
extern int psN; // power series are modulo x^psN
BOOL fout;
BOOL append;
Miracl precision=20;
ofstream mueller;
// Code to parse formula in command line
// This code isn't mine, but its public domain
// Shamefully I forget the source
//
// NOTE: It may be necessary on some platforms to change the operators * and #
//
#define TIMES '*'
#define RAISE '#'
Big tt;
static char *ss;
void eval_power (Big& oldn,Big& n,char op)
{
if (op) n=pow(oldn,toint(n)); // power(oldn,size(n),n,n);
}
void eval_product (Big& oldn,Big& n,char op)
{
switch (op)
{
case TIMES:
n*=oldn;
break;
case '/':
n=oldn/n;
break;
case '%':
n=oldn%n;
}
}
void eval_sum (Big& oldn,Big& n,char op)
{
switch (op)
{
case '+':
n+=oldn;
break;
case '-':
n=oldn-n;
}
}
void eval (void)
{
Big oldn[3];
Big n;
int i;
char oldop[3];
char op;
char minus;
for (i=0;i<3;i++)
{
oldop[i]=0;
}
LOOP:
while (*ss==' ')
ss++;
if (*ss=='-') /* Unary minus */
{
ss++;
minus=1;
}
else
minus=0;
while (*ss==' ')
ss++;
if (*ss=='(' || *ss=='[' || *ss=='{') /* Number is subexpression */
{
ss++;
eval ();
n=tt;
}
else /* Number is decimal value */
{
for (i=0;ss[i]>='0' && ss[i]<='9';i++)
;
if (!i) /* No digits found */
{
cout << "Error - invalid number" << endl;
exit (20);
}
op=ss[i];
ss[i]=0;
n=atol(ss);
ss+=i;
*ss=op;
}
if (minus) n=-n;
do
op=*ss++;
while (op==' ');
if (op==0 || op==')' || op==']' || op=='}')
{
eval_power (oldn[2],n,oldop[2]);
eval_product (oldn[1],n,oldop[1]);
eval_sum (oldn[0],n,oldop[0]);
tt=n;
return;
}
else
{
if (op==RAISE)
{
eval_power (oldn[2],n,oldop[2]);
oldn[2]=n;
oldop[2]=RAISE;
}
else
{
if (op==TIMES || op=='/' || op=='%')
{
eval_power (oldn[2],n,oldop[2]);
oldop[2]=0;
eval_product (oldn[1],n,oldop[1]);
oldn[1]=n;
oldop[1]=op;
}
else
{
if (op=='+' || op=='-')
{
eval_power (oldn[2],n,oldop[2]);
oldop[2]=0;
eval_product (oldn[1],n,oldop[1]);
oldop[1]=0;
eval_sum (oldn[0],n,oldop[0]);
oldn[0]=n;
oldop[0]=op;
}
else /* Error - invalid operator */
{
cout << "Error - invalid operator" << endl;
exit (20);
}
}
}
}
goto LOOP;
}
//
// When summing the Zk^n 0<=k0 && (!first || !brackets)) cout << "+";
first=FALSE;
if (cf==1) cout << "Y";
else cout << cf << "*Y";
if (j!=1) cout << "^" << j;
}
cf=z.coeff(0);
if (fout) mueller << cf << "\n" << L+1-i << "\n" << 0 << endl;
if (brackets)
{
cout << "+" << cf << ")*X";
if (i!=L) cout << "^" << L+1-i ;
}
else
{
if (i==L+1)
{
cout << "+" << cf;
continue;
}
if (cf!=0)
{
if (cf==1) cout << "+X";
else cout << "+" << cf << "*X";
if (i!=L) cout << "^" << L+1-i ;
}
}
// all other coefficients should now be zero
if (z.coeff(L)!=0)
{ // check next coefficient is zero
cout << "\n\n Sanity Check Failed " << endl;
exit(0);
}
}
for (i=0;i<=L+1;i++) c[i].clear(); // reclaim space
for (i=0;i<=v;i++) jlt[i].clear();
cout << endl;
fft_reset();
}
int main(int argc,char **argv)
{
Big p;
miracl *mip=get_mip();
int i,j,s,lo,hi,sp,ip,skip;
int primes[200];
BOOL dir,gotP,gothi,gotlo;
argv++; argc--;
int Base;
if (argc<1)
{
cout << "Incorrect usage" << endl;
cout << "Program generates Modular Polynomials, for use by fast Schoof-Elkies-Atkins" << endl;
cout << "program for counting points on an elliptic curve" << endl;
cout << "modpol " << endl;
cout << "OR" << endl;
cout << "modpol " << endl;
cout << "where the numbers define a range. The program will find the" << endl;
cout << "Modular Polynomials for primes in this range wrt the specified modulus" << endl;
cout << "To input P in Hex, precede with -h" << endl;
cout << "To search downwards for a prime, use flag -d" << endl;
cout << "NOTE: Program is both memory and time intensive" << endl;
cout << "To skip \"difficult\" primes, use -s2, -s3 or -s6" << endl;
cout << "where -s2 skips most and -s6 skips least" << endl;
cout << "To output polynomials to a file use flag -o " << endl;
cout << "e.g. modpol -f 2#192-2#64-1 0 150 -o p192.pol" << endl;
cout << "Alternatively to append to a file use flag -a " << endl;
cout << "See source code file for details" << endl;
cout << "\nFreeware from Shamus Software, Dublin, Ireland" << endl;
cout << "Full C++ source code and MIRACL multiprecision library available" << endl;
cout << "http://indigo.ie/~mscott for details" << endl;
cout << "or email mscott@indigo.ie" << endl;
return 0;
}
if (argc<3)
{
cout << "Error in command line" << endl;
return 0;
}
ip=0;
skip=12;
fout=FALSE;
dir=gotP=gothi=gotlo=FALSE;
append=FALSE;
Base=10;
while (ipIOBASE=Base;
p=argv[ip++];
mip->IOBASE=10;
gotP=TRUE;
continue;
}
if (!gotlo)
{
lo=atoi(argv[ip++]);
gotlo=TRUE;
continue;
}
if (!gothi)
{
hi=atoi(argv[ip++]);
gothi=TRUE;
continue;
}
cout << "Error in command line" << endl;
return 0;
}
if (!gothi || !gotlo)
{
cout << "Error in command line" << endl;
return 0;
}
if (lo>hi || hi>1000)
{
cout << "Invalid range specified" << endl;
return 0;
}
gprime(1000); // get all primes < 1000
for (i=0;;i++)
{
sp=mip->PRIMES[i];
primes[i]=sp;
if (sp==0) break;
}
if (!prime(p))
{
int incr=0;
cout << "That number is not prime!" << endl;
if (dir)
{
cout << "Looking for next lower prime" << endl;
p-=1; incr++;
while (!prime(p)) { p-=1; incr++; }
cout << "Prime P = P-" << incr << endl;
}
else
{
cout << "Looking for next higher prime" << endl;
p+=1; incr++;
while (!prime(p)) { p+=1; incr++; }
cout << "Prime P = P+" << incr << endl;
}
cout << "Prime P = " << p << endl;
}
cout << "P mod 24 = " << p%24 << endl;
cout << "P is " << bits(p) << " bits long" << endl;
if (fout && !append) mueller << p << endl;
modulo(p); // Set prime modulus for ZZn type
for (j=0,i=1;;i++) // lets go
{
sp=primes[i];
if (sp==0) break;
if (sphi) break;
for (s=1;;s++)
if (s*(sp-1)%12==0) break;
if (s>=skip) continue;
j++;
cout << endl;
cout << "prime " << j << " = " << sp << " (s=" << s << ")" << endl;
mueller_pol(sp,s);
}
cout << endl;
if (j==0) cout << "No primes processed in the specified range" << endl;
if (j==1) cout << "One prime processed in the specified range" << endl;
if (j>1) cout << j << " primes processed in the specified range" << endl;
return 0;
}