mrecgf2m.c

大数运算库

    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];
}

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

void multiply2(_MIPD_ big x,big y,big w)
{
    int i,j,xl,yl,ml;
    mr_small p,q;

#ifndef MR_GENERIC_MT
    miracl *mr_mip=get_mip();
#endif

    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);

    if (xl>=KARATSUBA && yl>=KARATSUBA)
    { /* Use Karatsuba */
        if (xl>yl) ml=xl;
        else       ml=yl;
     
        karmul2(ml,mr_mip->w7->w,x->w,y->w,w0->w);  /* w7 is work-space */

        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++)
    { /* slow old-fashioned O(n^2) way */
        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);
}

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];
        z->len=lm;
        for (i=lm;i<lz;i++)
            gz[i]=0;
        if (lx==ly) mr_lzero(z);
    }
}

static void remain2(_MIPD_ big y,big x)
{ /* generic "remainder" program. x%=y */
#ifndef MR_GENERIC_MT
    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;
}

static void gcd2(_MIPD_ big x,big y,big g)
{
#ifndef MR_GENERIC_MT
    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;

#ifndef MR_GENERIC_MT
    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;
    B=mr_mip->BB;
    C=mr_mip->CC;

    k1=1+M/MIRACL;       /* words from MSB to LSB */

    if (xl<=k1)
    {
      if (numbits(x)<=M) return;
    }

    rs1=M%MIRACL;
    ls1=MIRACL-rs1;

    if (M-A < MIRACL)
    { /* slow way */
        while (numbits(x)>=M+1)
        {
            copy(mr_mip->modulus,mr_mip->w7);
            shiftleftbits(mr_mip->w7,numbits(x)-M-1);
            add2(x,mr_mip->w7,x);    
        }
        return;
    }

    k2=1+(M-A)/MIRACL;   /* words from MSB to bit */
    rs2=(M-A)%MIRACL;
    ls2=MIRACL-rs2;

    if (B)
    { /* Pentanomial */
        k3=1+(M-B)/MIRACL;
        rs3=(M-B)%MIRACL;
        ls3=MIRACL-rs3;

        k4=1+(M-C)/MIRACL;
        rs4=(M-C)%MIRACL;
        ls4=MIRACL-rs4;
    }

    for (i=xl-1;i>=k1;i--)
    {
        if (rs1==0) gx[i-k1+1]^=gx[i];
        else
        {
            gx[i-k1+1]^=(gx[i]>>rs1);
            gx[i-k1]^=(gx[i]<<ls1);
        }
        if (rs2==0) gx[i-k2+1]^=gx[i];
        else
        {
            gx[i-k2+1]^=(gx[i]>>rs2);
            gx[i-k2]^=(gx[i]<<ls2);
        }
        if (B)
        {
            if (rs3==0) gx[i-k3+1]^=gx[i];
            else
            {
                gx[i-k3+1]^=(gx[i]>>rs3);
                gx[i-k3]^=(gx[i]<<ls3);
            }
            if (rs4==0) gx[i-k4+1]^=gx[i];
            else
            {
                gx[i-k4+1]^=(gx[i]>>rs4);
                gx[i-k4]^=(gx[i]<<ls4);
            }
        }
        gx[i]=0;
    }

    top=gx[k1-1]>>rs1;
    if (top!=0)
    {  
        gx[0]^=top;
        top<<=rs1;

        if (rs2==0) gx[k1-k2]^=top;
        else
        {
            gx[k1-k2]^=(top>>rs2);
            if (k1>k2) gx[k1-k2-1]^=(top<<ls2);
        }
        if (B)
        {
            if (rs3==0) gx[k1-k3]^=top;
            else
            {
                gx[k1-k3]^=(top>>rs3);
                if (k1>k3) gx[k1-k3-1]^=(top<<ls3);
            }
            if (rs4==0) gx[k1-k4]^=top;
            else
            {
                gx[k1-k4]^=(top>>rs4);
                if (k1>k4) gx[k1-k4-1]^=(top<<ls4);
            }
        }
       gx[k1-1]^=top;
    }

    if (gx[xl-1]==0) mr_lzero(x);
}

void incr2(big x,int n,big w)
{ /* increment x by small amount */
    if (x!=w) copy(x,w);
    if (n==0) return;
    if (w->len==0)
    {
        w->len=1;
        w->w[0]=n;
    }
    else
    {
        w->w[0]^=(mr_small)n;
        mr_lzero(w);
    }
}

void modsquare2(_MIPD_ big x,big w)
{ /* w=x*x mod f */
#ifndef MR_GENERIC_MT
    miracl *mr_mip=get_mip();
#endif
    square2(x,mr_mip->w0);
    reduce2(_MIPP_ mr_mip->w0,mr_mip->w0);
    copy(mr_mip->w0,w);
}

void modmult2(_MIPD_ big x,big y,big w)
{ /* w=x*y mod f */
#ifndef MR_GENERIC_MT
    miracl *mr_mip=get_mip();
#endif

    if (x==NULL || y==NULL)
    {
        zero(w);
        return;
    }

    if (x==y)
    {
        modsquare2(_MIPP_ x,w);
        return;
    }

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

    if (y->len==1)
    {
        if (y->w[0]==1)
        {
            copy(x,w);
            return;
        }
    } 
    
    multiply2(_MIPP_ x,y,mr_mip->w0);
    reduce2(_MIPP_ mr_mip->w0,mr_mip->w0);
    copy(mr_mip->w0,w);
}

void sqroot2(_MIPD_ big x,big y)
{ /* there are quicker ways, but its not time critical here... */
    int i;
#ifndef MR_GENERIC_MT
    miracl *mr_mip=get_mip();
#endif
    copy(x,y);
    for (i=1;i<mr_mip->M;i++)
        modsquare2(_MIPP_ y,y);

}

void power2(_MIPD_ big x,int m,big w)
{ /* w=x^m mod f. Could be optimised a lot, but not time critical for me */
#ifndef MR_GENERIC_MT
    miracl *mr_mip=get_mip();
#endif
    copy(x,mr_mip->w1);
    
    convert(_MIPP_ 1,w);
    forever
    {
        if (m%2!=0)
            modmult2(_MIPP_ w,mr_mip->w1,w);
        m/=2;
        if (m==0) break;
        modsquare2(_MIPP_ mr_mip->w1,mr_mip->w1);
    }
}

/* Schroeppel, Orman, O'Malley, Spatscheck    *
 * "Almost Inverse" algorithm, Crypto '95     *
 * More optimization here and in-lining would *
 * speed up AFFINE mode. I observe that       *
 * pentanomials would be more efficient if C  *
 * were greater                               */

BOOL inverse2(_MIPD_ big x,big w)
{
    mr_small lsw,*gw;
    int n,bits,step,k=0;
    int k1,k2,k3,k4,ls1,ls2,ls3,ls4,rs1,rs2,rs3,rs4;
    int M,A,B,C;
    big t;
#ifndef MR_GENERIC_MT
    miracl *mr_mip=get_mip();
#endif
    if (size(x)==0) return FALSE;