mrecgf2m.c

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

/*
 *   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 <stdlib.h> 
#include "miracl.h"
#ifdef MR_STATIC
#include <string.h>
#endif

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

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

/*
#include <emmintrin.h>

static mr_small mr_mul2(mr_small a,mr_small b,mr_small *r)
{
    __m64 tt[4],xg,rg;
    mr_small q;

    tt[0]=_m_from_int(0);
    tt[1]=_m_from_int(a);
    tt[2]=_m_psllqi(tt[1],1);
    tt[3]=_m_pxor(tt[1],tt[2]);

    rg=tt[b&3]; 
    xg=tt[(b>>2)&3]; rg=_m_pxor(rg,_m_psllqi(xg,2));
    xg=tt[(b>>4)&3]; rg=_m_pxor(rg,_m_psllqi(xg,4));
    xg=tt[(b>>6)&3]; rg=_m_pxor(rg,_m_psllqi(xg,6));
    xg=tt[(b>>8)&3]; rg=_m_pxor(rg,_m_psllqi(xg,8));
    xg=tt[(b>>10)&3]; rg=_m_pxor(rg,_m_psllqi(xg,10));
    xg=tt[(b>>12)&3]; rg=_m_pxor(rg,_m_psllqi(xg,12));
    xg=tt[(b>>14)&3]; rg=_m_pxor(rg,_m_psllqi(xg,14));
    xg=tt[(b>>16)&3]; rg=_m_pxor(rg,_m_psllqi(xg,16));
    xg=tt[(b>>18)&3]; rg=_m_pxor(rg,_m_psllqi(xg,18));
    xg=tt[(b>>20)&3]; rg=_m_pxor(rg,_m_psllqi(xg,20));
    xg=tt[(b>>22)&3]; rg=_m_pxor(rg,_m_psllqi(xg,22));
    xg=tt[(b>>24)&3]; rg=_m_pxor(rg,_m_psllqi(xg,24));
    xg=tt[(b>>26)&3]; rg=_m_pxor(rg,_m_psllqi(xg,26));
    xg=tt[(b>>28)&3]; rg=_m_pxor(rg,_m_psllqi(xg,28));
    xg=tt[(b>>30)]; rg=_m_pxor(rg,_m_psllqi(xg,30));

    *r=_m_to_int(rg);
    q=_m_to_int(_m_psrlqi(rg,32));

    _m_empty();

    return q;
}

*/

/* wouldn't it be nice if instruction sets supported a 
   one cycle "carry-free" multiplication instruction ... 
   The SmartMips does - its called maddp                 */

#ifndef MR_COMBA2

static mr_small mr_mul2(mr_small a,mr_small b,mr_small *r)
{
    int k;
    mr_small kb,t[16];
    mr_small x,q,p;
    mr_utype tb0;
#if MIRACL > 32
    mr_utype tb1,tb2;
#endif
  
    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=(mr_utype)(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 UNWOUNDM
    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 UNWOUNDM
    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 UNWOUNDM
    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 UNWOUNDM

    q=p=(mr_small)0;
    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;
}

#endif

/*

Now inlined into square2(.) - see below

static mr_small mr_sqr2(mr_small a,mr_small *r)
{ 
    int i;
    mr_small t,q;
    static const 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};

	q=0; *r=0;
    
#if MIRACL == 8
#define UNWOUNDS
        t=look[a&0xF];
        a>>=4;
        *r>>=8;
        *r|=t;
        t=look[a&0xF];
        a>>=4;
        q>>=8;
        q|=t;
#endif

#if MIRACL == 16
#define UNWOUNDS
        t=look[a&0xF];
        a>>=4;
        *r>>=8;
        *r|=t;
        t=look[a&0xF];
        a>>=4;
        *r>>=8;
        *r|=t;
        t=look[a&0xF];
        a>>=4;
        q>>=8;
        q|=t;
        t=look[a&0xF];
        a>>=4;
        q>>=8;
        q|=t;
#endif

#if MIRACL == 32
#define UNWOUNDS

        t=look[a&0xF];
        a>>=4;
        *r>>=8;
        *r|=t;
        t=look[a&0xF];
        a>>=4;
        *r>>=8;
        *r|=t;
        t=look[a&0xF];
        a>>=4;
        *r>>=8;
        *r|=t;
        t=look[a&0xF];
        a>>=4;
        *r>>=8;
        *r|=t;
  
		t=look[a&0xF];
        a>>=4;
        q>>=8;
        q|=t;
        t=look[a&0xF];
        a>>=4;
        q>>=8;
        q|=t;
		t=look[a&0xF];
        a>>=4;
        q>>=8;
        q|=t;
        t=look[a&0xF];
        a>>=4;
        q>>=8;
        q|=t;

#endif

#ifndef UNWOUNDS

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

#endif

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

int degree2(big x)
{ /* returns -1 for x=0 */
    return (numbits(x)-1);
}

/*
static int zerobits(big x)
{ 
    int m,n,k;
    mr_small *gx,lsb,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;
        lsb=gx[m];
        while (!(bit&lsb))
        {
            n++; 
            bit<<=1;
        }
        break;
    }
    return (MIRACL*m+n);
}

static void shiftrightbits(big x,int m)
{
    int i,k=x->len;
    int w=m/MIRACL;
    int b=m%MIRACL;  
    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;
    }
    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,j,n,m;
    mr_small a,t,r,*gw;

    static const 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};

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

    for (i=n-1;i>=0;i--)
    {
        a=gw[i];

#if MIRACL == 8
#define UNWOUNDS
        r=0;
        t=look[a&0xF];
        a>>=4;
        r>>=8;
        r|=t;
        gw[i+i]=r; r=0;
        t=look[a&0xF];
        r>>=8;
        r|=t;
        gw[i+i+1]=r;
#endif

#if MIRACL == 16
#define UNWOUNDS
        r=0;
        t=look[a&0xF];
        a>>=4;
        r>>=8;
        r|=t;
        t=look[a&0xF];
        a>>=4;
        r>>=8;
        r|=t;
        gw[i+i]=r; r=0;
        t=look[a&0xF];
        a>>=4;
        r>>=8;
        r|=t;
        t=look[a&0xF];
        r>>=8;
        r|=t;
        gw[i+i+1]=r;
#endif

#if MIRACL == 32
#define UNWOUNDS
        gw[i+i]=(look[a&0xF]>>24)|(look[(a>>4)&0xF]>>16)|(look[(a>>8)&0xF]>>8)|look[(a>>12)&0xF]; 
        gw[i+i+1]=(look[(a>>16)&0xF]>>24)|(look[(a>>20)&0xF]>>16)|(look[(a>>24)&0xF]>>8)|look[(a>>28)];
#endif

#ifndef UNWOUNDS

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

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

#endif

    }
    if (gw[m-1]==0) 
    {
        w->len--;
        if (gw[m-2]==0)
            mr_lzero(w);
    }
}

/* Use karatsuba to multiply two polynomials with coefficients in GF(2^m) */

#ifndef MR_STATIC

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;