modpol.cpp

miracl-大数运算库,大家使用有什么问题请多多提意见

// Newton's Identities - Calculate coefficients from Power Sums
//
// from a Web page somewhere and 3. page 54
//
    for (i=1;i<=L+1;i++)
    {
        cout << setw(3) << i << flush;
        c[i]=0;
        for (j=1;j<=i;j++)
              c[i]+=(ps[j]*c[i-j]);   // cheap, but lots of them
        c[i]=(-c[i])/i;
        cout << "\b\b\b" << flush;
    }

    cout << setw(3) << L+1 << endl;

    for (i=0;i<=L+1;i++) ps[i].clear(); // reclaim space
//
// Get powers of j(Lt)^i, i=1 to v
// These will be needed to determine the exponent of Y in each
// coefficient of the Modular Polynomial
//
    jlt[0]=1;
    jlt[1]=klein;

    for (i=2;i<=v;i++)
        jlt[i]=jlt[i-1]*klein;    // cheap
//
// Find Modular Polynomial, format it, and output
//
// 2. page 5, middle "Hl(X) = ..."
//
    cout << "\nG" << L << "(X,Y) = X^" << L+1 ;

    if (fout) 
    {
        mueller << L << endl;
        mueller << 1 << "\n" << L+1 << "\n" << 0 << endl;
    }

    for (i=1;i<=L+1;i++)
    {
        ZZn cf;
        BOOL brackets,first;
        first=TRUE;
        brackets=FALSE;
        z=c[i];            

// idea is to reduce this to an integer
// by subtracting j(Lt)^k as necessary
// The power of k required is then
// the coefficient of Y^k in G(X,Y)
// The first coefficient of c[i] tells us which j(Lt)^k to try

        if (z.first()!=0) 
        {
            brackets=TRUE;
            cout << "+(" ;
        }
   // coefficient may be a polynomial in Y
        while (z.first()!=0)    
        {
            int j=(-z.first()/L);    // index into jlt
            cf=z.coeff(z.first());   // get coefficient to be cancelled
            if (fout) mueller << cf << "\n" << L+1-i << "\n" << j << endl;
            if (cf==0) break;
            z-=(jlt[j]*cf);
            if (cf>0 && (!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 <prime modulus P> <low number> <high number>" << endl;
        cout << "OR" << endl;
        cout << "modpol <formula for P>   <low number> <high number>" << 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 <filename>" << 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 <filename>" << 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 [e-mail]" << 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 (ip<argc)
    {
        if (!gotP && strcmp(argv[ip],"-f")==0)
        {
            ip++;
            if (!gotP && ip<argc)
            {

                ss=argv[ip++];
                tt=0;
                eval();
                p=tt;
                gotP=TRUE;
                continue;
            }
            else
            {
                cout << "Error in command line" << endl;
                return 0;
            }
        }

        if (strcmp(argv[ip],"-d")==0)
        {
            ip++;
            dir=TRUE;
            continue;
        }
        if (skip==12 && strcmp(argv[ip],"-s2")==0)
        {
            ip++;
            skip=2;
            continue;
        }
        if (skip==12 && strcmp(argv[ip],"-s3")==0)
        {
            ip++;
            skip=3;
            continue;
        }
        if (skip==12 && strcmp(argv[ip],"-s6")==0)
        {
            ip++;
            skip=6;
            continue;
        }

        if (!fout && strcmp(argv[ip],"-o")==0)
        {
            ip++;
            if (ip<argc)
            {
                fout=TRUE;
                append=FALSE;
                mueller.open(argv[ip++]);
                continue;
            }
            else
            {
                cout << "Error in command line" << endl;
                return 0;
            }
        }
        if (!fout && strcmp(argv[ip],"-a")==0)
        {
            ip++;
            if (ip<argc)
            {
                fout=TRUE;
                append=TRUE;
                mueller.open(argv[ip++],ios::app);
                continue;
            }
            else
            {
                cout << "Error in command line" << endl;
                return 0;
            }
        }
        if (strcmp(argv[ip],"-h")==0)
        {
            ip++;
            Base=16;
            continue;
        }
        if (!gotP)
        {
            mip->IOBASE=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 (sp<lo) continue;
        if (sp>hi) 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;   
}