zzn12.cpp

密码算法程序包miracl contains all the source code for the latest version of MIRACL: a Multiprecision Int

/*
 *    MIRACL  C++ Implementation file zzn12.cpp
 *
 *    AUTHOR  : M. Scott
 *  
 *    PURPOSE : Implementation of class ZZn12  (Arithmetic over n^12)
 *
 * WARNING: This class has been cobbled together for a specific use with
 * the MIRACL library. It is not complete, and may not work in other 
 * applications
 *
 *    NOTE: - The irreducible polynomial is assumed to be of the form 
 *            x^6+i, where i is either 
 *      sqrt(-1) or 1+sqrt(-1) if p=3 mod 8
 *   or sqrt(-2) 1+sqrt(-2) if p=5 mod 8 or 7 mod 8
 *
 *    Copyright (c) 2006 Shamus Software Ltd.
 */

#include "zzn12.h"

using namespace std;

// Frobenius...

ZZn12& ZZn12::powq(const ZZn12* X)
{
    *this=conj(a[5])*X[5]+conj(a[4])*X[4]+conj(a[3])*X[3]+conj(a[2])*X[2]+conj(a[1])*X[1]+conj(a[0]);
    return *this;
}

int degree(const ZZn12& x)
{
    if (x.a[5]!=0) return 5;
    if (x.a[4]!=0) return 4;
    if (x.a[3]!=0) return 3;
    if (x.a[2]!=0) return 2;
    if (x.a[1]!=0) return 1;
    if (x.a[0]!=0) return 0;
    return -1;
}

void ZZn12::get(ZZn2* x)  
{for (int i=0;i<6;i++) x[i]=a[i]; } 

void ZZn12::get(ZZn2& x) 
{x=a[0]; }

ZZn12& ZZn12::operator*=(const ZZn12& x)
{ 
    if (&x==this)
    {  // Karatsuba squaring
        if (a[5].iszero() && a[3].iszero() && a[1].iszero())
        {
            if (a[4].iszero() && a[2].iszero())
            {
                a[0]*=a[0];
                return *this;
            }

            ZZn2 r[6];

            r[0] = a[0]*a[0];
            r[2] = a[2]*a[2];
            r[4] = a[4]*a[4];

            r[1] = (a[0] + a[2]); r[1] *= r[1]; r[1] -= (r[0] + r[2]);
            r[3] = (a[2] + a[4]); r[3] *= r[3]; r[3] -= (r[2] + r[4]);
            r[5] = (a[0] + a[4]); r[5] *= r[5]; r[5] -= (r[0] + r[4]);
            r[2] += r[5];

            /* *** */

            r[1] -= txx(r[4]);
            r[0] -= txx(r[3]);
             
            a[0] = r[0]; a[1] = 0;
            a[2] = r[1]; a[3] = 0;
            a[4] = r[2]; a[5] = 0;
            return *this;
        }

        ZZn2 r[11],a02,a04,a13,a15,a23,a24,a35,wa,wb,wc;

        r[ 0] = a[0]*a[0];
        r[ 2] = a[1]*a[1];
        r[ 4] = a[2]*a[2];
        r[ 6] = a[3]*a[3];
        r[ 8] = a[4]*a[4];
        r[10] = a[5]*a[5];

        a02=a[0]+a[2];
        a04=a[0]+a[4];
        a13=a[1]+a[3];
        a15=a[1]+a[5];
        a23=a[2]+a[3];
        a24=a[2]+a[4];
        a35=a[3]+a[5];

        r[ 1] = (a[0] + a[1]); r[ 1] *= r[ 1]; wa=r[1]; r[ 1] -= (r[0]+r[2]);
        r[ 9] = (a[4] + a[5]); r[ 9] *= r[ 9]; wc=r[9]; r[ 9] -= (r[8]+r[10]);
 
        r[ 3] = a02+a13;  r[ 3] *= r[ 3];
        r[ 5] = a04+a15;  r[ 5] *= r[ 5];
        r[ 7] = a24+a35;  r[ 7] *= r[ 7];

        a02 *= a02; a02 -= (r[ 0] + r[ 4]);
        a04 *= a04; a04 -= (r[ 0] + r[ 8]);
        a13 *= a13; a13 -= (r[ 2] + r[ 6]);
        a23 *= a23; wb=a23; a23 -= (r[4]+r[6]);
        a15 *= a15; a15 -= (r[ 2] + r[10]);
        a24 *= a24; a24 -= (r[ 4] + r[ 8]);
        a35 *= a35; a35 -= (r[ 6] + r[10]);

        r[ 3] -=(wa + wb + a02 + a13);   
        r[ 5] -=(wa + wc + a04 + a15);
        r[ 7] -=(wb + wc + a24 + a35);

        r[ 8] += a35;
        r[ 6] += a15 + a24;
        r[ 5] += a23;
        r[ 4] += a04 + a13;
        r[ 2] += a02;


        r[4] -= txx(r[10]);
        r[3] -= txx(r[ 9]);
        r[2] -= txx(r[ 8]);
        r[1] -= txx(r[ 7]);
        r[0] -= txx(r[ 6]);
         
        a[0] = r[0]; a[1] = r[1];
        a[2] = r[2]; a[3] = r[3];
        a[4] = r[4]; a[5] = r[5];
        return *this;    
     }
     else
     { // Karatsuba multiplication

        if (a[5].iszero() && x.a[5].iszero()
         && a[3].iszero() && x.a[3].iszero()
         && a[1].iszero() && x.a[1].iszero())
        {
            if (a[4].iszero() && x.a[4].iszero()
             && a[2].iszero() && x.a[2].iszero())
            {
                a[0]*=x.a[0];
                return *this;
            }   

            ZZn2 r[6];

            r[0] = a[0]*x.a[0];
            r[2] = a[2]*x.a[2];
            r[4] = a[4]*x.a[4];

            r[1] = (a[0] + a[2])*(x.a[0] + x.a[2]) - (r[0] + r[2]);
            r[3] = (a[2] + a[4])*(x.a[2] + x.a[4]) - (r[2] + r[4]);
            r[5] = (a[0] + a[4])*(x.a[0] + x.a[4]) - (r[0] + r[4]);

            r[2] += r[5];


            r[1] -= txx(r[4]);
            r[0] -= txx(r[3]);
            
            a[0] = r[0]; a[1] = 0;
            a[2] = r[1]; a[3] = 0;
            a[4] = r[2]; a[5] = 0;
            return *this;
        }

        ZZn2 r[11],a02,a04,a13,a15,a23,a24,a35,wc,wa,wb,x02,x04,x13,x23,x15,x24,x35,za,zb,zc;

        r[ 0] = a[0]*x.a[0];
        r[ 2] = a[1]*x.a[1];
        r[ 4] = a[2]*x.a[2];
        r[ 6] = a[3]*x.a[3];
        r[ 8] = a[4]*x.a[4];
        r[10] = a[5]*x.a[5];

        a02=a[0]+a[2];
        a04=a[0]+a[4];
        a13=a[1]+a[3];
        a15=a[1]+a[5];
        a23=a[2]+a[3];
        a24=a[2]+a[4];
        a35=a[3]+a[5];

        x02=x.a[0] + x.a[2];
        x04=x.a[0] + x.a[4];
        x13=x.a[1] + x.a[3];
        x15=x.a[1] + x.a[5];
        x23=x.a[2] + x.a[3];
        x24=x.a[2] + x.a[4];
        x35=x.a[3] + x.a[5];

        za=(a02+a13)*(x02+x13);
        zb=(a04+a15)*(x04+x15);
        zc=(a24+a35)*(x24+x35);
 
        r[1] = (a[0] + a[1])*(x.a[0] + x.a[1]); wa=r[1]; r[1]-=(r[0]+r[2]);
        r[9] = (a[4] + a[5])*(x.a[4] + x.a[5]); wc=r[9]; r[9]-=(r[8]+r[10]);

        a02*=x02; a02 -= (r[ 0] + r[ 4]);
        a04*=x04; a04 -= (r[ 0] + r[ 8]);
        a13*=x13; a13 -= (r[ 2] + r[ 6]);
        a23*=x23; wb=a23; a23 -= (r[4]+r[6]);
        a15*=x15; a15 -= (r[ 2] + r[10]);
        a24*=x24; a24 -= (r[ 4] + r[ 8]);
        a35*=x35; a35 -= (r[ 6] + r[10]);

        r[ 3] = za - (wa + wb + a02 + a13);
        r[ 5] = zb - (wa + wc + a04 + a15);
        r[ 7] = zc - (wb + wc + a24 + a35);

        r[ 8] += a35;
        r[ 6] += a15 + a24;
        r[ 5] += a23;
        r[ 4] += a04 + a13;
        r[ 2] += a02;

        r[4] -= txx(r[10]);
        r[3] -= txx(r[ 9]);
        r[2] -= txx(r[ 8]);
        r[1] -= txx(r[ 7]);
        r[0] -= txx(r[ 6]);
        
        a[0] = r[0]; a[1] = r[1];
        a[2] = r[2]; a[3] = r[3];
        a[4] = r[4]; a[5] = r[5];
    }

    return *this;
}

/* See "Fast Implementation of Elliptic Curve Arithmetic in GF(p^n)",
   Lim and Hwang */

ZZn12& ZZn12::invert(void)
{
    int degF,degG,degB,degC,i,j,n,d,g;
    ZZn2 BB[7],FF[7],CC[7],GG[7],alpha,beta,gamma;
    ZZn2 *B=BB,*C=CC,*F=FF,*G=GG,*T;

    ZZn12 t,keep=*this;

    /* *** */
    C[0]=1; F[0].set(1,1);  F[6]=1;  // F=x^6+i

    degF=6; degG=degree(*this); degC=0; degB=-1;

    if (degG==0)
    {
        a[0]=(ZZn)1/a[0];
        return *this;
    }

    for (i=0;i<6;i++)
    {
        G[i]=a[i];
        a[i]=0;
    }

    while (degF!=0)
    {
        if (degF<degG)
        {
            T=F; F=G; G=T; d=degF; degF=degG; degG=d;
            T=B; B=C; C=T; d=degB; degB=degC; degC=d;
        }
        j=degF-degG;
        alpha=G[degG]*G[degG];
        beta=F[degF]*G[degG];
        gamma=G[degG]*F[degF-1]-F[degF]*G[degG-1];

        for (i=0;i<=degF;i++)
        {
            F[i]*=alpha;
            if (i>=j-1) F[i]-=gamma*G[i-j+1];
            if (i>=j)   F[i]-=beta*G[i-j];
        }

        for (i=0;i<=degB || i<=degC+j;i++)
        {
            B[i]*=alpha;
            if (i>=j-1) B[i]-=gamma*C[i-j+1];
            if (i>=j)   B[i]-=beta*C[i-j];
        }        

        while (degF>=0 && F[degF]==0) degF--;

        if (degF==degG)
        {
            alpha=F[degF];
            for (i=0;i<=degF;i++)
            {
                F[i]*=G[degF];
                F[i]-=alpha*G[i];
            }
            for (i=0;i<=6-degF;i++) 
            {
                B[i]*=G[degF];
                B[i]-=alpha*C[i];
            } 
            while (degF>=0 && F[degF]==0) degF--;
        }

        degB=5; while (degB>=0 && B[degB]==0) degB--;
    
    }

    alpha=(ZZn2)1/F[0];
    for (i=0;i<=degB;i++)
        a[i]=alpha*B[i];

//    if (!(*this*keep).isunity()) cout << "inverse Failed" << endl;

    return *this;    
}

ZZn12& ZZn12::operator/=(const ZZn12& x)
{ // inversion 
 ZZn12 y=x;
 y.invert();
 *this *= y;
 return *this;
}

ZZn12 operator+(const ZZn12& x,const ZZn12& y) 
{ZZn12 w=x; for (int i=0;i<6;i++) w.a[i]+=y.a[i]; return w;} 

ZZn12 operator+(const ZZn12& x,const ZZn2& y) 
{ZZn12 w=x; w.a[0]+=y; return w; } //

ZZn12 operator-(const ZZn12& x,const ZZn12& y) 
{ZZn12 w=x; for (int i=0;i<6;i++) w.a[i]-=y.a[i]; return w; } 

ZZn12 operator-(const ZZn12& x,const ZZn2& y) 
{ZZn12 w=x; w.a[0]-=y; return w; } //

ZZn12 operator-(const ZZn12& x) 
{ZZn12 w; for (int i=0;i<6;i++) w.a[i]=-x.a[i]; return w; } 

ZZn12 operator*(const ZZn12& x,const ZZn12& y)
{ZZn12 w=x;  w*=y;  return w;}

ZZn12 operator*(const ZZn12& x,const ZZn2& y)
{ZZn12 w=x; for (int i=0;i<6;i++) w.a[i]*=y; return w;} //

ZZn12 operator*(const ZZn2& y,const ZZn12& x)
{ZZn12 w=x; for (int i=0;i<6;i++) w.a[i]*=y; return w;} //

ZZn12 operator*(const ZZn12& x,int y)
{ZZn12 w=x; for (int i=0;i<6;i++) w.a[i]*=y; return w;}
                                         
ZZn12 operator*(int y,const ZZn12& x)
{ZZn12 w=x; for (int i=0;i<6;i++) w.a[i]*=y; return w;}

ZZn12 operator/(const ZZn12& x,const ZZn12& y)
{ZZn12 w=x; w/=y; return w;}

ZZn12 operator/(const ZZn12& x,const ZZn2& y)
{ZZn12 w=x; ZZn2 j=(ZZn)1/y; for (int i=0;i<6;i++) w.a[i]*=j; return w;} //

ZZn12 randn12(void)
{ZZn12 w; for (int i=0;i<6;i++) w.a[i]=randn2(); return w;}

#ifndef MR_NO_STANDARD_IO

ostream& operator<<(ostream& s,ZZn12& b)
{
    int i;
    ZZn2 x[6];
    b.get(x);
    s << "[";
    for (i=0;i<=4;i++)
        s << x[i] << ",";
    s << x[5] << "]";
    return s;    
}

#endif

// Left to right method - with windows

ZZn12 pow(const ZZn12* t,const ZZn12& x,const Big& k)
{
    int i,j,nb,n,nbw,nzs;
    ZZn12 u=x;

    if (k==0) return (ZZn12)1;
    if (k==1) return u;

    nb=bits(k);
    if (nb>1) for (i=nb-2;i>=0;)
    {
        n=window(k,i,&nbw,&nzs);
        for (j=0;j<nbw;j++) u*=u;
        if (n>0) u*=t[n/2];
        i-=nbw;
        if (nzs)
        {
            for (j=0;j<nzs;j++) u*=u;
            i-=nzs;
        }
    }
    return u;
}

void precompute(const ZZn12& x,ZZn12* t)
{
    int i;
    ZZn12 u2,u=x;
    u2=(u*u);
    t[0]=u;
   
    for (i=1;i<16;i++)
        t[i]=u2*t[i-1];

}

ZZn12 pow(const ZZn12& x,const Big& k)
{
    ZZn12 u,t[16];

    if (k==0) return (ZZn12)1;
    u=x;
    if (k==1) return u;
//
// Prepare table for windowing
//
    precompute(u,t);
    return pow(t,u,k);
}

// standard MIRACL multi-exponentiation

ZZn12 pow(int n,const ZZn12* x,const Big* b)
{
    int k,j,i,m,nb,ea;
    ZZn12 *G;
    ZZn12 r;
    m=1<<n;
    G=new ZZn12[m];

 // precomputation
    
    for (i=0,k=1;i<n;i++)
    {
        for (j=0; j < (1<<i) ;j++)
        {
            if (j==0)   G[k]=x[i];
            else        G[k]=G[j]*x[i];      
            k++;
        }
    }

    nb=0;
    for (j=0;j<n;j++) 
        if ((k=bits(b[j]))>nb) nb=k;

    r=1;
    for (i=nb-1;i>=0;i--) 
    {
        ea=0;
        k=1;
        for (j=0;j<n;j++)
        {
            if (bit(b[j],i)) ea+=k;
            k<<=1;
        }
        r*=r;
        if (ea!=0) r*=G[ea];
    }
    delete [] G;
    return r;
}