//
// cl /O2 /GX bn.cpp zzn12.cpp zzn6a.cpp zzn2.cpp zzn.cpp ecn2.cpp ecn.cpp big.cpp ms32.lib
// Program to generate BN curves for use by ake12t.cpp and ake12.cpp
//

#include <iostream>
#include "big.h"
#include "ecn.h"
#include "ecn2.h"
#include "zzn12.h"

using namespace std;

Miracl precision=100;

int main()
{
    int ns;
    int sign;
    Big n,p,t,x,cube,y,b,eta,w,cf[4];
    miracl*mip=&precision;
    ECn P;
    ECn2 Q;
    ZZn2 x2,y2,xi;
    ZZn12 X,Y;
    ZZn6 A,B;
    
//    mip->IOBASE=16;

    x=pow((Big)2,62)+pow((Big)2,61)-1;  // x is this size to get right size for t, p and n
                                        // x is low hamming weight
    mip->IOBASE=16;
    sign=2;  // 1= positive, 2=negative for +/- x solutions
    ns=1;
    forever
    {
        forever
        {
            sign=3-sign;
            if (sign==1) x+=1;
          
            if (sign==1) p=36*pow(x,4)-36*pow(x,3)+24*x*x-6*x+1;
            else         p=36*pow(x,4)+36*pow(x,3)+24*x*x+6*x+1;

            if (p%6!=1) continue;   // check congruence conditions
            if (p%8!=7) continue;
            if (p%9!=4) continue;
            if (!prime(p)) continue;
            modulo(p);
            xi.set(1,1);
            if (pow(xi,(p-1)/2)==1) continue;
            if (pow(xi,(p-1)/3)==1) continue;  // make sure that x^6+c is irreducible

            t=6*x*x+1;
            n=p+1-t;
            if (prime(n)) break;
        }     
        
        cf[3]=1;
        cf[2]=6*x*x+1;
        cf[1]=36*x*x*x-18*x*x+12*x+1;
        cf[0]=36*x*x*x-30*x*x+18*x-2;

        b=0;
        forever
        {
            forever
            {
                b+=1;
                y=sqrt((b+1),p);
                if (y==0) continue;
                break;
            }
            ecurve(0,b,p,MR_AFFINE);

            P.set(1,y);
            if ((n*P).iszero()) break;
        }

        mip->TWIST=TRUE;
        while (!Q.set(randn2())) ;
        Q*=(p-1+t);
    
        Q.get(x2,y2);
        A.set((ZZn2)0,x2,(ZZn2)0);
        B.set((ZZn2)0,y2,(ZZn2)0);
        X.set(A); Y.set((ZZn6)0,B);

        if (!(n*Q).iszero() || Y*Y!=X*X*X +(ZZn12)b) continue;

        cout << "solution " << ns << endl;
        cout << "irreducible polynomial = X^6 + " << xi << endl;
        cout << "p=" << p << endl;
        cout << "p mod 72 = " << p%72 << endl;
        cout << "x mod 72 = " << x%72 << endl;
        if (sign==1) cout << "x= +" << x << endl;
        else         cout << "x= -" << x << endl;
        cout << "bits(p)= " << bits(p) << endl;
        cout << "n=" << n << endl;
        cout << "t=" << t << endl;
        cout << "coefficients of (p^4-p^2+1)/n to base p" << endl; 
        cout << "cf[0]= " << cf[0] << endl;
        cout << "cf[1]= " << cf[1] << endl;
        cout << "cf[2]=                  " << cf[2] << endl;
        cout << "cf[3]= " << cf[3] << endl;
        cout << "bits(t)= " << bits(t) << endl;
        cout << "ham(t-1) = " << ham(t-1) << " (small is better for Ate pairing)" << endl;
     //   cout << "ham(6*x)+ham(x)+1 = " << ham(6*x)+ham(x)+1 << endl;
     //   cout << "ham(6*x*x+1) = " << ham(6*x*x+1) << endl;
     //   if (ham(6*x)+ham(x)+1 > ham(6*x*x+1)) break;
        
        cout << "bits(n)= " << bits(n) << endl;
        cout << "ham(n-1) = " << ham(n-1) << " (small is better for Tate pairing)" << endl;

        cout << "E(Fp): y^2=x^3+" << b << endl;
        cout << "Point P= " << P << endl;
        cout << "Point Q= " << Q << endl << endl;
        ns++;

    }
    return 0;
}

