mrecgf2m.c

大数运算库

/*
 *   MIRACL routines for arithmetic over GF(2^m), and 
 *   implementation of Elliptic Curve Cryptography over GF(2^m) 
 *   mrecgf2m.c
 *
 *   Curve equation is Y^2 + XY = X^3 + A.X^2 + B
 *   where A is 0 or 1
 *
 *   For algorithms used, see IEEE P1363 Standard, Appendix A
 *   unless otherwise stated.
 * 
 *   The time-critical routines are the multiplication routine multiply2()
 *   and (for AFFINE co-ordinates), the modular inverse routine inverse2() 
 *   and the routines it calls.
 *
 *   No assembly language used.
 *
 *   Copyright (c) 2000-2001 Shamus Software Ltd.
 */

#include <stdio.h>
#include <miracl.h>

#ifndef MR_NOFULLWIDTH
                     /* This does not make sense using floating-point! */

#define M1 (MIRACL-1)
#define M2 (MIRACL-2)
#define M3 (MIRACL-3)
#define TOPBIT ((mr_small)1<<M1)
#define SECBIT ((mr_small)1<<M2)
#define THDBIT ((mr_small)1<<M3)
#define M8 (MIRACL-8)

#define KARATSUBA 2

static mr_small look[16]=
{0,(mr_small)1<<M8,(mr_small)4<<M8,(mr_small)5<<M8,(mr_small)16<<M8,
(mr_small)17<<M8,(mr_small)20<<M8,(mr_small)21<<M8,(mr_small)64<<M8,
(mr_small)65<<M8,(mr_small)68<<M8,(mr_small)69<<M8,(mr_small)80<<M8,
(mr_small)81<<M8,(mr_small)84<<M8,(mr_small)85<<M8};

/* This is extremely time-critical, and expensive */

/* wouldn't it be nice if instruction sets supported a 
   one cycle "carry-free" multiplication instruction ... */

static mr_small mr_mul2(mr_small a,mr_small b,mr_small *r)
{
    int k;
    mr_small kb,x,q,p,t[16];
    mr_utype tb0,tb1,tb2;
    q=p=(mr_small)0;

    if (b<a) {x=a; a=b; b=x;}   /* make sure a is smaller */

    if (a<16)
    { /* a is small */
        t[0]=0; t[1]=b;
        p=t[a&1];
        a>>=1;
        x=t[a&1];
        p^=x<<1;
        q^=x>>M1;
        a>>=1;
        x=t[a&1];
        p^=x<<2;
        q^=x>>M2;
        a>>=1;
        x=t[a&1];
        p^=x<<3;
        q^=x>>M3;
        *r=p;
        return q;
    }

    kb=b;

#if MIRACL <= 32
    t[0]=0;               /* small look up table */
    t[3]=t[2]=a<<1;       /* it can overflow.... */
    t[1]=t[2]>>1;
    t[3]^=t[1];

    tb0=(a&TOPBIT);       /* remember top bit    */
    tb0>>=M1;             /* all ones if top bit is one */
#else
    t[0]=0;               /* larger look-up table */
    t[8]=a<<3;
    t[4]=t[8]>>1;
    t[2]=t[4]>>1;
    t[1]=t[2]>>1;
    t[3]=t[5]=t[7]=t[9]=t[11]=t[13]=t[15]=t[1];
    t[3]^=t[2];
    t[5]^=t[4];
    t[9]^=t[8]; 
    t[6]=t[3]<<1;
    t[7]^=t[6];
    t[10]=t[5]<<1;
    t[11]^=t[10];
    t[12]=t[6]<<1;
    t[13]^=t[12];
    t[14]=t[7]<<1;
    t[15]^=t[14];

    tb0=(a&TOPBIT);       /* remember top bits  */
    tb0>>=M1;             /* all bits one, if this bit is set in a */
    tb1=(a&SECBIT)<<1;      
    tb1>>=M1;
    tb2=(a&THDBIT)<<2;
    tb2>>=M1;
#endif

#if MIRACL == 16
#define UNWOUND
    p=q=t[b&3];                       q>>=2;
    x=t[(b>>2)&3];  q^=x; p^=(x<<2);  q>>=2;   
    x=t[(b>>4)&3];  q^=x; p^=(x<<4);  q>>=2;   
    x=t[(b>>6)&3];  q^=x; p^=(x<<6);  q>>=2;
    x=t[(b>>8)&3];  q^=x; p^=(x<<8);  q>>=2;
    x=t[(b>>10)&3]; q^=x; p^=(x<<10); q>>=2;
    x=t[(b>>12)&3]; q^=x; p^=(x<<12); q>>=2;
    x=t[(b>>14)];   q^=x; p^=(x<<14); q>>=2;
#endif

#if MIRACL == 32
#define UNWOUND
    p=q=t[b&3];                       q>>=2;
    x=t[(b>>2)&3];  q^=x; p^=(x<<2);  q>>=2;   /* 8 ASM 80386 instructions */
    x=t[(b>>4)&3];  q^=x; p^=(x<<4);  q>>=2;   /* but only 4 ARM instructions! */
    x=t[(b>>6)&3];  q^=x; p^=(x<<6);  q>>=2;
    x=t[(b>>8)&3];  q^=x; p^=(x<<8);  q>>=2;
    x=t[(b>>10)&3]; q^=x; p^=(x<<10); q>>=2;
    x=t[(b>>12)&3]; q^=x; p^=(x<<12); q>>=2;
    x=t[(b>>14)&3]; q^=x; p^=(x<<14); q>>=2;
    x=t[(b>>16)&3]; q^=x; p^=(x<<16); q>>=2;
    x=t[(b>>18)&3]; q^=x; p^=(x<<18); q>>=2;
    x=t[(b>>20)&3]; q^=x; p^=(x<<20); q>>=2;
    x=t[(b>>22)&3]; q^=x; p^=(x<<22); q>>=2;
    x=t[(b>>24)&3]; q^=x; p^=(x<<24); q>>=2;
    x=t[(b>>26)&3]; q^=x; p^=(x<<26); q>>=2;
    x=t[(b>>28)&3]; q^=x; p^=(x<<28); q>>=2;
    x=t[(b>>30)];   q^=x; p^=(x<<30); q>>=2;
#endif

#if MIRACL == 64
#define UNWOUND
    p=q=t[b&0xf];                       q>>=4;
    x=t[(b>>4)&0xf];  q^=x; p^=(x<<4);  q>>=4;
    x=t[(b>>8)&0xf];  q^=x; p^=(x<<8);  q>>=4;
    x=t[(b>>12)&0xf]; q^=x; p^=(x<<12); q>>=4;
    x=t[(b>>16)&0xf]; q^=x; p^=(x<<16); q>>=4;
    x=t[(b>>20)&0xf]; q^=x; p^=(x<<20); q>>=4;
    x=t[(b>>24)&0xf]; q^=x; p^=(x<<24); q>>=4;
    x=t[(b>>28)&0xf]; q^=x; p^=(x<<28); q>>=4;
    x=t[(b>>32)&0xf]; q^=x; p^=(x<<32); q>>=4;
    x=t[(b>>36)&0xf]; q^=x; p^=(x<<36); q>>=4;
    x=t[(b>>40)&0xf]; q^=x; p^=(x<<40); q>>=4;
    x=t[(b>>44)&0xf]; q^=x; p^=(x<<44); q>>=4;
    x=t[(b>>48)&0xf]; q^=x; p^=(x<<48); q>>=4;
    x=t[(b>>52)&0xf]; q^=x; p^=(x<<52); q>>=4;
    x=t[(b>>56)&0xf]; q^=x; p^=(x<<56); q>>=4;
    x=t[(b>>60)];     q^=x; p^=(x<<60); q>>=4;

#endif

#ifndef UNWOUND
    for (k=0;k<MIRACL;k+=8)
    {                 
        q^=(t[b&3]);   
        b>>=2; 
        p>>=2; 
        p|=q<<M2;
        q>>=2;

        q^=(t[b&3]);
        b>>=2;
        p>>=2;
        p|=q<<M2;
        q>>=2;

        q^=(t[b&3]);
        b>>=2;
        p>>=2;
        p|=q<<M2;
        q>>=2;

        q^=(t[b&3]);
        b>>=2;
        p>>=2;
        p|=q<<M2;
        q>>=2;
    }
#endif

#if MIRACL <= 32
    p^=(tb0&(kb<<M1));       /* compensate for top bit */
    q^=(tb0&(kb>>1));        /* don't break pipeline.. */
#else
    p^=(tb0&(kb<<M1));
    q^=(tb0&(kb>>1));
    p^=(tb1&(kb<<M2));
    q^=(tb1&(kb>>2));
    p^=(tb2&(kb<<M3));
    q^=(tb2&(kb>>3));
#endif

    *r=p;
    return q;
}

static mr_small mr_sqr2(mr_small a,mr_small *r)
{ /* squaring is linear O(n), and not time-critical */
    int i;
    mr_small t,q;
    q=0; *r=0;
    
  /* simply insert a 0 between bits - look-up table is a bit faster */

    for (i=0;i<MIRACL/2;i+=4)
    {
        t=look[a&0xF];
        a>>=4;
        *r>>=8;
        *r|=t;
    }

    for (i=0;i<MIRACL/2;i+=4)
    {
        t=look[a&0xF];
        a>>=4;
        q>>=8;
        q|=t;
    }

    return q;
}

static int numbits(big x)
{ /* return degree of x */
    mr_small *gx=x->w,bit=TOPBIT;
    int m,k=x->len;
    if (k==0) return 0;
    m=k*MIRACL;
    while (!(gx[k-1]&bit))
    {
        m--;
        bit>>=1;
    }
    return m;
}

static int zerobits(big x)
{ /* return number of zero bits at the end of x */
    int m,n,k;
    mr_small *gx,bit=1;
    k=x->len; 
    if (k==0) return (-1);
    gx=x->w;
    for (m=0;m<k;m++)
    {
        if (gx[m]==0) continue;
        n=0;
        while (!(bit&gx[m]))
        {
            n++; 
            bit<<=1;
        }
        break;
    }
    return (MIRACL*m+n);
}

static void shiftrightbits(big x,int m)
{
    int i,k=x->len;
    int w=m/MIRACL;  /* words */
    int b=m%MIRACL;  /* bits  */
    mr_small *gx=x->w;
    if (k==0 || m==0) return;
    if (w>0)
    {
        for (i=0;i<k-w;i++)
            gx[i]=gx[i+w];
        for (i=k-w;i<k;i++) gx[i]=0;
        x->len-=w;
    }
/* time critical */
    if (b!=0) 
    {
        for (i=0;i<k-w-1;i++)
            gx[i]=(gx[i]>>b)|(gx[i+1]<<(MIRACL-b));   
        gx[k-w-1]>>=b;
        if (gx[k-w-1]==0) x->len--;
    }
}

static void shiftleftbits(big x,int m)
{
    int i,k=x->len;
    mr_small j; 
    int w=m/MIRACL;  /* words */
    int b=m%MIRACL;  /* bits  */
    mr_small *gx=x->w;
    if (k==0 || m==0) return;
    if (w>0)
    {
        for (i=k+w-1;i>=w;i--)
            gx[i]=gx[i-w];
        for (i=w-1;i>=0;i--) gx[i]=0;
        x->len+=w;
    }
/* time critical */
    if (b!=0) 
    {
        j=gx[k+w-1]>>(MIRACL-b);
        if (j!=0)
        {
            x->len++;
            gx[k+w]=j;
        }
        for (i=k+w-1;i>w;i--)
        {
            gx[i]=(gx[i]<<b)|(gx[i-1]>>(MIRACL-b));
        }
        gx[w]<<=b;
    }
}

static void square2(big x,big w)
{ /* w=x*x where x can be NULL so be careful */
    int i,n,m;
    mr_small a,b,*gw;

    if (x!=w) copy(x,w);
    n=w->len;
    m=n+n;
    w->len=m;
    gw=w->w; 

    for (i=n-1;i>=0;i--)
    {
        a=mr_sqr2(gw[i],&b);
        gw[i+i]=b;
        gw[i+i+1]=a;
    }
    if (gw[m-1]==0) mr_lzero(w);
}


/* Use katatsuba to multiply two polynomial with coefficients in GF(2^m) */

void karmul2_poly(_MIPD_ int n,big *t,big *x,big *y,big *z)
{
    int m,nd2,nd,md,md2;                          
    if (n==1) 
    { /* finished */
        modmult2(_MIPP_ *x,*y,*z);
        zero(z[1]);
        return;
    }       
    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%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]);
    }
}

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

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