www.pudn.com > miracl.zip > MUELLER.CPP
//
// Program to generate Modular Polynomials, 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.
//
// .raw file format
// ,<1st coef>,<1st power of X>,<1st power of Y>,<2nd coef>...
// Each polynomial ends wih zero powers of X and Y.
//
// NOTE: This program is very hard on memory. In places "obvious" speed
// optimizations have not been applied, if they are memory intensive. But in
// any case 64Mb is really a minimal requirement to generate enough
// polynomials for serious work.
//
// The time/memory requirements for a particular prime L depend on the value
// s, defined as the smallest integer such that s.(L-1) is divisible by 12.
// The value of s can be 1, 2, 3 or 6. The bigger s, the worse the time/memory
// requirements, and the bigger the coefficients in the polynomial.
// The flags -s2, -s3, and -s6 cause the primes in the specified
// range to be skipped if, for example, s>=3.
//
// Four things can go wrong. An "Out of Space message means that you have run
// out of virtual memory. The "control panel" on your operating system should
// enable you to fix this. A "Number too big" error means that the value
// specified in the first parameter of the call to mirsys() has proven to be
// too small. This also means that you have pushed the program beyond the
// limits for which we have tested it (well done!). If your hard disk
// "trashes" continually, and processing slows to a crawl, you have run out of
// physical memory. Buy more memory for your computer, or a new computer.
// Another problem may be I/O buffer overflow. Use set_io_buffer_size() to
// specify a larger I/O buffer
//
// With 96 Mb of RAM, what works for me (so far, running over a long weekend)
// is
//
// mueller 0 180 -o mueller.raw
// mueller 180 300 -s6 -a mueller.raw
//
// and it should be possible to continue, for example, like
//
// mueller 300 400 -s3 -a mueller.raw
// mueller 400 500 -s2 -a mueller.raw
//
// This creates a file mueller.raw containing modular polynomials
// for 69 primes
//
// Of course different ranges can be covered simultaneously on multiple
// computers, if you have them.
//
// Alternatively, if these problems should prove insurmountable - see the
// "modpol" application
//
// This program would be a lot faster if the coefficients were calculated mod
// many 32-bit primes, and then combined via the CRT.
//
// Copyright Shamus Software Ltd., 1999
//
#include
#include
#include
#include
#include "ps_big.h" // power series class
using namespace std;
extern int psN; // power series are modulo x^psN
BOOL fout;
ofstream mueller;
//
// When summing the Zk^n 0<=k0 && (!first || !brackets)) cout << "+";
first=FALSE;
if (cf==1) cout << "Y";
if (cf==-1) cout << "-Y";
if (abs(cf)!=1) cout << cf << "*Y";
if (j!=1) cout << "^" << j;
}
cf=z.coeff(0);
if (fout) mueller << cf << "\n" << L+1-i << "\n" << 0 << endl;
if (cf>0) cout << "+";
if (brackets) cout << cf << ")*X";
else
{
if (i==L+1)
{
cout << cf;
continue;
}
if (cf==1) cout << "X";
if (cf==-1) cout << "-X";
if (abs(cf)!=1) 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)
{
miracl *mip;
int i,j,s,lo,hi,p,ip,skip;
int primes[200];
BOOL gothi,gotlo;
argv++; argc--;
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 << "mueller " << endl;
cout << "where the numbers define a range. The program will attempt to" << endl;
cout << "find the Modular Polynomials for all primes in this range" << endl;
cout << "To output polynomials to a file use flag -o " << endl;
cout << "e.g. mueller 0 10 -o mueller.raw" << endl;
cout << "Alternatively to append to a file use flag -a " << 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 << "See source code file for details" << endl;
cout << "Files generated from different ranges may be combined in an" << endl;
cout << "obvious way using a standard text editor" << 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<2)
{
cout << "Error in command line" << endl;
return 0;
}
ip=0;
skip=12;
fout=FALSE;
gothi=gotlo=FALSE;
while (iphi || hi>1000)
{
cout << "Invalid range specified" << endl;
return 0;
}
mip=mirsys(20,0);
gprime(1000); // get all primes < 1000
for (i=0;;i++)
{
p=mip->PRIMES[i];
primes[i]=p;
if (p==0) break;
}
mirexit();
for (j=0,i=1;;i++) // lets go
{
p=primes[i];
if (p==0) break;
if (phi) break;
for (s=1;;s++)
if (s*(p-1)%12==0) break;
if (s>=skip) continue;
// p*s/6 seems to be a fortuitous upper bound (?)
// on the size of integer coefficients generated.
// This MAY need to be increased if "number
// too big" errors accur.
mirsys(1+(p*s)/6,0);
set_io_buffer_size(4096);
j++;
cout << endl;
cout << "prime " << j << " = " << p << " (s=" << s << ")" << endl;
cout << 32*(1+(p*s)/6) << " bits reserved for each coefficient" << endl;
mueller_pol(p,s);
mirexit();
}
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;
}