mrecgf2m.c

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

    }       
    if (n==2)
    {  /* in-line 2x2 */
        modmult2(_MIPP_ x[0],y[0],z[0]);
        modmult2(_MIPP_ x[1],y[1],z[2]);
        add2(x[0],x[1],t[0]);
        add2(y[0],y[1],t[1]);
        modmult2(_MIPP_ t[0],t[1],z[1]);
        add2(z[1],z[0],z[1]);
        add2(z[1],z[2],z[1]);
        zero(z[3]);
        return;
    }

    if (n==3)
    {
        modmult2(_MIPP_ x[0],y[0],z[0]); 
        modmult2(_MIPP_ x[1],y[1],z[2]); 
        modmult2(_MIPP_ x[2],y[2],z[4]); 
        add2(x[0],x[1],t[0]);
        add2(y[0],y[1],t[1]);
        modmult2(_MIPP_ t[0],t[1],z[1]); 
        add2(z[1],z[0],z[1]);  
        add2(z[1],z[2],z[1]);  
        add2(x[1],x[2],t[0]);
        add2(y[1],y[2],t[1]);
        modmult2(_MIPP_ t[0],t[1],z[3]); 
        add2(z[3],z[2],z[3]);
        add2(z[3],z[4],z[3]); 
        add2(x[0],x[2],t[0]);
        add2(y[0],y[2],t[1]);   
        modmult2(_MIPP_ t[0],t[1],t[0]); 
        add2(z[2],t[0],z[2]);
        add2(z[2],z[0],z[2]);
        add2(z[2],z[4],z[2]); 
        zero(z[5]);
        return;
    }

    if (n%2==0)
    {
        md=nd=n;
        md2=nd2=n/2;
    }
    else
    {
        nd=n+1;
        md=n-1;
        nd2=nd/2; md2=md/2;
    }

    for (m=0;m<nd2;m++)
    {
        copy(x[m],z[m]);
        copy(y[m],z[nd2+m]);
    }
    for (m=0;m<md2;m++)
    { 
        add2(z[m],x[nd2+m],z[m]);
        add2(z[nd2+m],y[nd2+m],z[nd2+m]);
    }

    karmul2_poly(_MIPP_ nd2,&t[nd],z,&z[nd2],t); 

    karmul2_poly(_MIPP_ nd2,&t[nd],x,y,z);  

    for (m=0;m<nd;m++) add2(t[m],z[m],t[m]);

    karmul2_poly(_MIPP_ md2,&t[nd],&x[nd2],&y[nd2],&z[nd]);

    for (m=0;m<md;m++) add2(t[m],z[nd+m],t[m]);
    for (m=0;m<nd;m++) add2(z[nd2+m],t[m],z[nd2+m]);
}

void karmul2_poly_upper(_MIPD_ int n,big *t,big *x,big *y,big *z)
{ /* n is large and even, and upper half of z is known already */
    int m,nd2,nd;
    nd2=n/2; nd=n;

    for (m=0;m<nd2;m++)
    { 
        add2(x[m],x[nd2+m],z[m]);
        add2(y[m],y[nd2+m],z[nd2+m]);
    }

    karmul2_poly(_MIPP_ nd2,&t[nd],z,&z[nd2],t); 

    karmul2_poly(_MIPP_ nd2,&t[nd],x,y,z);   /* only 2 karmuls needed! */

    for (m=0;m<nd;m++) add2(t[m],z[m],t[m]);

    for (m=0;m<nd2;m++) 
    {
        add2(z[nd+m],z[nd+nd2+m],z[nd+m]);
        add2(z[nd+m],t[nd2+m],z[nd+m]);
    }

    for (m=0;m<nd;m++) 
    {
        add2(t[m],z[nd+m],t[m]);
        add2(z[nd2+m],t[m],z[nd2+m]);
    }
}

#endif

/* Some in-line karatsuba down at the bottom... */

#ifndef MR_COMBA2

static void mr_bottom1(mr_small *x,mr_small *y,mr_small *z)
{
    z[1]=mr_mul2(x[0],y[0],&z[0]);
}

static void mr_bottom2(mr_small *x,mr_small *y,mr_small *z)
{
    mr_small q0,r0,q1,r1,q2,r2;

    q0=mr_mul2(x[0],y[0],&r0);
    q1=mr_mul2(x[1],y[1],&r1);
    q2=mr_mul2((mr_small)(x[0]^x[1]),(mr_small)(y[0]^y[1]),&r2);

    z[0]=r0;
    z[1]=q0^r1^r0^r2;
    z[2]=q0^r1^q1^q2;
    z[3]=q1;
}

static void mr_bottom3(mr_small *x,mr_small *y,mr_small *z)
{ /* just 6 mr_muls... */
    mr_small q0,r0,q1,r1,q2,r2;
    mr_small a0,b0,a1,b1,a2,b2;

    q0=mr_mul2(x[0],y[0],&r0);
    q1=mr_mul2(x[1],y[1],&r1);
    q2=mr_mul2(x[2],y[2],&r2);

    a0=mr_mul2((mr_small)(x[0]^x[1]),(mr_small)(y[0]^y[1]),&b0);
    a1=mr_mul2((mr_small)(x[1]^x[2]),(mr_small)(y[1]^y[2]),&b1);
    a2=mr_mul2((mr_small)(x[0]^x[2]),(mr_small)(y[0]^y[2]),&b2);

    b0^=r0^r1;
    a0^=q0^q1;
    b1^=r1^r2;
    a1^=q1^q2;
    b2^=r0^r2;
    a2^=q0^q2;

    z[0]=r0;
    z[1]=q0^b0;
    z[2]=r1^a0^b2;
    z[3]=q1^b1^a2;
    z[4]=r2^a1;
    z[5]=q2;
}

static void mr_bottom4(mr_small *x,mr_small *y,mr_small *z)
{ /* unwound 4x4 karatsuba multiplication - only 9 muls */
    mr_small q0,r0,q1,r1,q2,r2,tx,ty;
    mr_small xx0,yy0,xx1,yy1,t0,t1,t2,t3;

    q0=mr_mul2(x[0],y[0],&r0);
    q1=mr_mul2(x[1],y[1],&r1);

    tx=x[0]^x[1];
    ty=y[0]^y[1];
    q2=mr_mul2(tx,ty,&r2);

    z[0]=r0;
    z[1]=q0^r1^r0^r2;
    z[2]=q0^r1^q1^q2;
    z[3]=q1;

    q0=mr_mul2(x[2],y[2],&r0);

    q1=mr_mul2(x[3],y[3],&r1);

    tx=x[2]^x[3];
    ty=y[2]^y[3];
    q2=mr_mul2(tx,ty,&r2);

    z[4]=r0;
    z[5]=q0^r1^r0^r2;
    z[6]=q0^r1^q1^q2;
    z[7]=q1;

    xx0=x[2]^x[0];
    yy0=y[2]^y[0];
    q0=mr_mul2(xx0,yy0,&r0);
   
    xx1=x[3]^x[1];
    yy1=y[3]^y[1];
    q1=mr_mul2(xx1,yy1,&r1);

    tx=xx0^xx1;
    ty=yy0^yy1;
    q2=mr_mul2(tx,ty,&r2);

    t0=z[0]^z[4]^r0;
    t1=z[1]^z[5]^q0^r1^r0^r2;
    t2=z[2]^z[6]^q0^r1^q1^q2;
    t3=z[3]^z[7]^q1; 

    z[2]^=t0;
    z[3]^=t1;
    z[4]^=t2;
    z[5]^=t3;
}

void karmul2(int n,mr_small *t,mr_small *x,mr_small *y,mr_small *z)
{ /* Karatsuba multiplication - note that n can be odd or even */
    int m,nd2,nd,md,md2;

    if (n<=4)
    {
        if (n==1)
        {
            mr_bottom1(x,y,z);
            return;
        }
        if (n==2)
        {   
            mr_bottom2(x,y,z);
            return;
        }
        if (n==3)
        {   
            mr_bottom3(x,y,z);
            return;
        }
        if (n==4)
        {   
            mr_bottom4(x,y,z);
            return;
        }
    }
    if (n%2==0)
    {
        md=nd=n;
        md2=nd2=n/2;
    }
    else
    {
        nd=n+1;
        md=n-1;
        nd2=nd/2; md2=md/2;
    }

    for (m=0;m<nd2;m++)
    {
        z[m]=x[m];
        z[nd2+m]=y[m];
    }
    for (m=0;m<md2;m++)
    {
        z[m]^=x[nd2+m];
        z[nd2+m]^=y[nd2+m];
    }

    karmul2(nd2,&t[nd],z,&z[nd2],t); 
    karmul2(nd2,&t[nd],x,y,z);  

    for (m=0;m<nd;m++) t[m]^=z[m];

    karmul2(md2,&t[nd],&x[nd2],&y[nd2],&z[nd]);

    for (m=0;m<md;m++) t[m]^=z[nd+m];
    for (m=0;m<nd;m++) z[nd2+m]^=t[m];
}

#endif

/* this is time-critical, so use karatsuba here, since addition is cheap *
 * and easy (no carries to worry about...)                               */

/* #define CLAIRE */

void multiply2(_MIPD_ big x,big y,big w)
{

#ifdef MR_OS_THREADS
    miracl *mr_mip=get_mip();
#endif

#ifdef MR_COMBA2
    comba_mult2(_MIPP_ x,y,w);
    return;
#else
    int i,j,xl,yl,ml;
#ifdef CLAIRE
    int d;
    mr_small hi,lo;
#endif
    mr_small p,q;

    big w0=mr_mip->w0;

    if (x==NULL || y==NULL)
    {
        zero(w);
        return;
    }
    if (x->len==0 || y->len==0)
    {
        zero(w);
        return;
    }

    xl=x->len;
    yl=y->len;
    zero(w0);

#ifdef CLAIRE

/* Comba method */

    w0->len=xl+yl;
    d=1+mr_mip->M/MIRACL;
    hi=lo=0;
    for (i=0;i<d;i++)
    {
        for (j=0;j<=i;j++)
        {
            q=mr_mul2(x->w[j],y->w[i-j],&p);
            hi^=q; lo^=p;
        }
        w0->w[i]=lo; lo=hi; hi=0;
    }
    for (i=d;i<2*d-1;i++)
    {
        for (j=i-d+1;j<d;j++)
        {
            q=mr_mul2(x->w[j],y->w[i-j],&p);
            hi^=q; lo^=p;
        }
        w0->w[i]=lo; lo=hi; hi=0;
    }
    w0->w[2*d-1]=lo;
    mr_lzero(w0);
    copy(w0,w);

#else

/* recommended method as mr_mul2 is so slow... */

    if (xl>=MR_KARATSUBA && yl>=MR_KARATSUBA)
    { 
        if (xl>yl) ml=xl;
        else       ml=yl;
     
        karmul2(ml,mr_mip->w7->w,x->w,y->w,w0->w);

        mr_mip->w7->len=w0->len=2*ml+1;
        mr_lzero(w0);
        mr_lzero(mr_mip->w7);    
        copy(w0,w);
        return;
    }

    w0->len=xl+yl;
    for (i=0;i<xl;i++)
    { 
        for (j=0;j<yl;j++)
        {
            q=mr_mul2(x->w[i],y->w[j],&p); 
            w0->w[i+j]^=p;
            w0->w[i+j+1]^=q;
        } 
    }
    mr_lzero(w0);
    copy(w0,w);
#endif
#endif
}

void add2(big x,big y,big z)
{ /* XOR x and y */
    int i,lx,ly,lz,lm;
    mr_small *gx,*gy,*gz;

    if (x==y)
    {
        zero(z);
        return;
    }
    if (y==NULL)
    {
        copy(x,z);
        return;
    }
    else if (x==NULL) 
    {
        copy(y,z);
        return;
    }

    if (x==z)
    {
        gy=y->w; gz=z->w;
        ly=y->len; lz=z->len;
        lm=lz; if (ly>lz) lm=ly;
        for (i=0;i<lm;i++)
            gz[i]^=gy[i];
        z->len=lm;
        if (gz[lm-1]==0) mr_lzero(z);
    }
    else
    {
        gx=x->w; gy=y->w; gz=z->w;
        lx=x->len; ly=y->len; lz=z->len;
        lm=lx; if (ly>lx) lm=ly;

        for (i=0;i<lm;i++)
            gz[i]=gx[i]^gy[i];
        for (i=lm;i<lz;i++)
            gz[i]=0;
        z->len=lm;
        if (gz[lm-1]==0) mr_lzero(z);
    }
}

static void remain2(_MIPD_ big y,big x)
{ /* generic "remainder" program. x%=y */
#ifdef MR_OS_THREADS
    miracl *mr_mip=get_mip();
#endif
    int my=numbits(y);
    int mx=numbits(x);
    while (mx>=my)
    {
        copy(y,mr_mip->w7);
        shiftleftbits(mr_mip->w7,mx-my);
        add2(x,mr_mip->w7,x);    
        mx=numbits(x);
    }
    return;
}

void gcd2(_MIPD_ big x,big y,big g)
{
#ifdef MR_OS_THREADS
    miracl *mr_mip=get_mip();
#endif
    if (size(y)==0) 
    {
        copy(x,g);
        return;
    }
    copy(x,mr_mip->w1);
    copy(y,mr_mip->w2);
    forever
    {
        remain2(_MIPP_ mr_mip->w2,mr_mip->w1);
        if (size(mr_mip->w1)==0) break;
        copy(mr_mip->w1,mr_mip->w3);
        copy(mr_mip->w2,mr_mip->w1);
        copy(mr_mip->w3,mr_mip->w2);
    }
    copy(mr_mip->w2,g);
}


/* See "Elliptic Curves in Cryptography", Blake, Seroussi & Smart, 
   Cambridge University Press, 1999, page 20, for this fast reduction
   routine - algorithm II.9 */

void reduce2(_MIPD_ big y,big x)
{ /* reduction wrt the trinomial or pentanomial modulus        *
   * Note that this is linear O(n), and thus not time critical */
    int k1,k2,k3,k4,ls1,ls2,ls3,ls4,rs1,rs2,rs3,rs4,i;
    int M,A,B,C;
    int xl;
    mr_small top,*gx,w;

#ifdef MR_OS_THREADS
    miracl *mr_mip=get_mip();
#endif

    if (x!=y) copy(y,x);
    xl=x->len;
    gx=x->w;

    M=mr_mip->M;
    A=mr_mip->AA;
    if (A==0)
    {
        mr_berror(_MIPP_ MR_ERR_NO_BASIS);
        return;
    }
    B=mr_mip->BB;
    C=mr_mip->CC;

 
/* If optimizing agressively it makes sense to make this code specific to a particular field
   For example code like this can be optimized for the case 
   m=163. Note that the general purpose code involves lots of branches - these cause breaks in 
   the pipeline and they are slow. Further loop unrolling would be even faster...
*/

#if MIRACL == 32

    if (M==163 && A==7 && B==6 && C==3)
    {
        for (i=xl-1;i>=6;i--)
        {
            w=gx[i]; gx[i]=0;
            gx[i-5]^=((w>>3)^(w<<4)^(w<<3)^w);
            gx[i-6]^=(w<<29);              
            gx[i-4]^=((w>>28)^(w>>29));
        }
        top=gx[5]>>3;
        
        gx[0]^=top;
        top<<=3;
        gx[1]^=(top>>28);
        gx[0]^=(top<<4);
        gx[1]^=(top>>29);
        gx[0]^=(top<<3);
        gx[0]^=top;
        gx[5]^=top;
        
        x->len=6; 
        if (gx[5]==0) mr_lzero(x);
        
        return;
    }

    if (M==103 && A==9)
    {
        for (i=xl-1;i>=4;i--)
        {
            w=gx[i]; gx[i]=0;
            gx[i-3]^=((w>>7)^(w<<2));
            gx[i-4]^=(w<<25);
            gx[i-2]^=(w>>30);
        }
        top=gx[3]>>7;
        gx[0]^=top;
        top<<=7;
        gx[1]^=(top>>30);
        gx[0]^=(top<<2);
        gx[3]^=top;
        x->len=4;
        if (gx[3]==0) mr_lzero(x);

        return;
    }

    if (M==283 && A==119)
    {
        for (i=xl-1;i>=9;i--)