mrpower.c

大数运算库

/*
 *   MIRACL methods for modular exponentiation
 *   mrpower.c 
 *
 *   Copyright (c) 1988-1999 Shamus Software Ltd.
 */

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

void nres_powltr(_MIPD_ int x,big y,big w)
{ /* calculates w=x^y mod z using Left to Right Method   */
  /* uses only n^2 modular squarings, because x is small */
  /* Note: x is NOT an nresidue */
    int i,nb;
#ifndef MR_GENERIC_MT
    miracl *mr_mip=get_mip();
#endif

    if (mr_mip->ERNUM) return;
    copy(y,mr_mip->w1);

    MR_IN(86)

    zero(w);
    if (x==0) 
    {
        if (size(mr_mip->w1)==0) 
        { /* 0^0 = 1 */
            convert(_MIPP_ 1,w);
            nres(_MIPP_ w,w);
        }
        MR_OUT
        return;
    }

    convert(_MIPP_ 1,w);
    nres(_MIPP_ w,w);
    if (size(mr_mip->w1)==0) 
    {
        MR_OUT
        return;
    }
    if (size(mr_mip->w1)<0) mr_berror(_MIPP_ MR_ERR_NEG_POWER);
    if (mr_mip->ERNUM)
    {
        MR_OUT
        return;
    }

    if (mr_mip->base==mr_mip->base2)
    { 
        nb=logb2(_MIPP_ mr_mip->w1);
        convert(_MIPP_ x,w);
        nres(_MIPP_ w,w);
        if (nb>1) for (i=nb-2;i>=0;i--)
        { /* Left to Right binary method */

            if (mr_mip->user!=NULL) (*mr_mip->user)();

            nres_modmult(_MIPP_ w,w,w);
            if (mr_testbit(_MIPP_ mr_mip->w1,i))
            { /* this is quick */
                premult(_MIPP_ w,x,w);
                divide(_MIPP_ w,mr_mip->modulus,mr_mip->modulus);
            }
        }
    }    
    else
    {
        expint(_MIPP_ 2,logb2(_MIPP_ mr_mip->w1)-1,mr_mip->w2);
        while (size(mr_mip->w2)!=0)
        { /* Left to Right binary method */

            if (mr_mip->user!=NULL) (*mr_mip->user)();

            nres_modmult(_MIPP_ w,w,w);
            if (compare(mr_mip->w1,mr_mip->w2)>=0)
            {
                premult(_MIPP_ w,x,w);
                divide(_MIPP_ w,mr_mip->modulus,mr_mip->modulus);
                subtract(_MIPP_ mr_mip->w1,mr_mip->w2,mr_mip->w1);
            }
            subdiv(_MIPP_ mr_mip->w2,2,mr_mip->w2);
        }
    }
    MR_OUT
    return;
}

void nres_powmod(_MIPD_ big x,big y,big w)
{  /*  calculates w=x^y mod z, using m-residues       *
    *  See "Analysis of Sliding Window Techniques for *
    *  Exponentiation, C.K. Koc, Computers. Math. &   *
    *  Applic. Vol. 30 pp17-24 1995. Uses work-space  *
    *  variables for pre-computed table. */
    int i,j,k,t,nb,nbw,nzs,n;
    big table[16];
#ifndef MR_GENERIC_MT
    miracl *mr_mip=get_mip();
#endif
    if (mr_mip->ERNUM) return;

    copy(y,mr_mip->w1);
    copy(x,mr_mip->w3);

    MR_IN(84)
    zero(w);
    if (size(x)==0)
    {
       if (size(mr_mip->w1)==0)
       { /* 0^0 = 1 */
           convert(_MIPP_ 1,w);
           nres(_MIPP_ w,w);
       } 
       MR_OUT
       return;
    }

    convert(_MIPP_ 1,w);
    nres(_MIPP_ w,w);
    if (size(mr_mip->w1)==0) 
    {
        MR_OUT
        return;
    }

    if (size(mr_mip->w1)<0) mr_berror(_MIPP_ MR_ERR_NEG_POWER);

    if (mr_mip->ERNUM)
    {
        MR_OUT
        return;
    }
    if (mr_mip->base==mr_mip->base2)
    { /* build a table. Up to 5-bit sliding windows. Windows with
       * two adjacent 0 bits simply won't happen */
        table[0]=mr_mip->w3; table[1]=mr_mip->w4; table[2]=mr_mip->w5; table[3]=mr_mip->w16;
        table[4]=NULL;  table[5]=mr_mip->w6; table[6]=mr_mip->w17; table[7]=mr_mip->w8;
        table[8]=NULL;  table[9]=NULL;  table[10]=mr_mip->w9; table[11]=mr_mip->w10;
        table[12]=NULL; table[13]=mr_mip->w11; table[14]=mr_mip->w12; table[15]=mr_mip->w13;

        nres_modmult(_MIPP_ mr_mip->w3,mr_mip->w3,mr_mip->w2);  /* x^2 */
        n=15;
        j=0;
        do
        { /* pre-computations */
            t=1; k=j+1;
            while (table[k]==NULL) {k++; t++;}
            copy(table[j],table[k]);
            for (i=0;i<t;i++) nres_modmult(_MIPP_ table[k],mr_mip->w2,table[k]);
            j=k;
        } while (j<n);

        nb=logb2(_MIPP_ mr_mip->w1);
        copy(mr_mip->w3,w);
        if (nb>1) for (i=nb-2;i>=0;)
        { /* Left to Right method */

            if (mr_mip->user!=NULL) (*mr_mip->user)();

            n=mr_window(_MIPP_ mr_mip->w1,i,&nbw,&nzs); 
            for (j=0;j<nbw;j++)
                    nres_modmult(_MIPP_ w,w,w);
            if (n>0) nres_modmult(_MIPP_ w,table[n/2],w); 
            i-=nbw;
            if (nzs)
            {
                for (j=0;j<nzs;j++) nres_modmult(_MIPP_ w,w,w);
                i-=nzs;
            }
        }
    }
    else
    {
        copy(mr_mip->w3,mr_mip->w2);
        forever
        { /* "Russian peasant" Right-to-Left exponentiation */

            if (mr_mip->user!=NULL) (*mr_mip->user)();

            if (subdiv(_MIPP_ mr_mip->w1,2,mr_mip->w1)!=0)
                nres_modmult(_MIPP_ w,mr_mip->w2,w);
            if (mr_mip->ERNUM || size(mr_mip->w1)==0) break;
            nres_modmult(_MIPP_ mr_mip->w2,mr_mip->w2,mr_mip->w2);
        }
    }
    MR_OUT
}

void nres_powmodn(_MIPD_ int n,big *x,big *y,big w)
{
    int i,j,k,m,nb,ea;
    big *G;
#ifndef MR_GENERIC_MT
    miracl *mr_mip=get_mip();
#endif
    if (mr_mip->ERNUM) return;

    MR_IN(112)

    m=1<<n;
    G=(big *)mr_alloc(_MIPP_ m,sizeof(big));

/* 2^n - n - 1 modmults */
/* 4(n=3) 11(n=4) etc   */

    for (i=0,k=1;i<n;i++)
    {
        for (j=0; j < (1<<i) ;j++)
        {
            G[k]=mirvar(_MIPP_ 0);
            if (j==0) copy(x[i],G[k]);
            else      nres_modmult(_MIPP_ G[j],x[i],G[k]);
            k++;
        }
    }

    nb=0;
    for (j=0;j<n;j++) 
        if ((k=logb2(_MIPP_ y[j]))>nb) nb=k;

    convert(_MIPP_ 1,w);
    nres(_MIPP_ w,w);

    if (mr_mip->base==mr_mip->base2)
    {

        for (i=nb-1;i>=0;i--) 
        {
            if (mr_mip->user!=NULL) (*mr_mip->user)();
            ea=0;
            k=1;
            for (j=0;j<n;j++)
            {
                if (mr_testbit(_MIPP_ y[j],i)) ea+=k;
                k<<=1;
            }
            nres_modmult(_MIPP_ w,w,w);
            if (ea!=0) nres_modmult(_MIPP_ w,G[ea],w);
        }
    }
    else mr_berror(_MIPP_ MR_ERR_NOT_SUPPORTED);
    for (i=1;i<m;i++) mirkill(G[i]);
    mr_free(G);

    MR_OUT
}

void powmodn(_MIPD_ int n,big *x,big *y,big p,big w)
{/* w=x[0]^y[0].x[1]^y[1] .... x[n-1]^y[n-1] mod n */
    int j;
#ifndef MR_GENERIC_MT
    miracl *mr_mip=get_mip();
#endif
    if (mr_mip->ERNUM) return;

    MR_IN(113)

    prepare_monty(_MIPP_ p);
    for (j=0;j<n;j++) nres(_MIPP_ x[j],x[j]);
    nres_powmodn(_MIPP_ n,x,y,w);   
    for (j=0;j<n;j++) redc(_MIPP_ x[j],x[j]);

    redc(_MIPP_ w,w);

    MR_OUT
}

void nres_powmod2(_MIPD_ big x,big y,big a,big b,big w)
{ /* finds w = x^y.a^b mod n. Fast for some cryptosystems */ 
    int i,j,nb,nb2,nbw,nzs,n;
    big table[16];
#ifndef MR_GENERIC_MT
    miracl *mr_mip=get_mip();
#endif
    if (mr_mip->ERNUM) return;
    copy(y,mr_mip->w1);
    copy(x,mr_mip->w2);
    copy(b,mr_mip->w3);
    copy(a,mr_mip->w4);
    zero(w);
    if (size(mr_mip->w2)==0 || size(mr_mip->w4)==0) return;

    MR_IN(99)

    convert(_MIPP_ 1,w);
    nres(_MIPP_ w,w);
    if (size(mr_mip->w1)==0 && size(mr_mip->w3)==0) 
    {
        MR_OUT
        return;
    }
    if (size(mr_mip->w1)<0 || size(mr_mip->w3)<0) mr_berror(_MIPP_ MR_ERR_NEG_POWER);
    if (mr_mip->ERNUM)
    {
        MR_OUT
        return;
    }
     
    if (mr_mip->base==mr_mip->base2)
    { /* uses 2-bit sliding window. This is 25% faster! */

        nres_modmult(_MIPP_ mr_mip->w2,mr_mip->w4,mr_mip->w5);
        nres_modmult(_MIPP_ mr_mip->w2,mr_mip->w2,mr_mip->w12);
        nres_modmult(_MIPP_ mr_mip->w4,mr_mip->w4,mr_mip->w13);
        nres_modmult(_MIPP_ mr_mip->w4,mr_mip->w13,mr_mip->w16);
        nres_modmult(_MIPP_ mr_mip->w2,mr_mip->w13,mr_mip->w6);
        nres_modmult(_MIPP_ mr_mip->w6,mr_mip->w4,mr_mip->w17);
        nres_modmult(_MIPP_ mr_mip->w4,mr_mip->w12,mr_mip->w8);
        nres_modmult(_MIPP_ mr_mip->w2,mr_mip->w12,mr_mip->w9);
        nres_modmult(_MIPP_ mr_mip->w4,mr_mip->w9,mr_mip->w10);
        nres_modmult(_MIPP_ mr_mip->w16,mr_mip->w12,mr_mip->w11);
        nres_modmult(_MIPP_ mr_mip->w9,mr_mip->w13,mr_mip->w12);
        nres_modmult(_MIPP_ mr_mip->w10,mr_mip->w13,mr_mip->w13);

        table[0]=NULL; table[1]=mr_mip->w4;  table[2]=mr_mip->w2;   table[3]=mr_mip->w5;
        table[4]=NULL; table[5]=mr_mip->w16; table[6]=mr_mip->w6;   table[7]=mr_mip->w17;
        table[8]=NULL; table[9]=mr_mip->w8;  table[10]=mr_mip->w9;  table[11]=mr_mip->w10;
        table[12]=NULL;table[13]=mr_mip->w11;table[14]=mr_mip->w12; table[15]=mr_mip->w13;
        nb=logb2(_MIPP_ mr_mip->w1);
        if ((nb2=logb2(_MIPP_ mr_mip->w3))>nb) nb=nb2;

        for (i=nb-1;i>=0;)
        { 
            if (mr_mip->user!=NULL) (*mr_mip->user)();

            n=mr_window2(_MIPP_ mr_mip->w1,mr_mip->w3,i,&nbw,&nzs);
            for (j=0;j<nbw;j++)
                nres_modmult(_MIPP_ w,w,w);
            if (n>0) nres_modmult(_MIPP_ w,table[n],w);
            i-=nbw;
            if (nzs) 
            {
                nres_modmult(_MIPP_ w,w,w);
                i--;
            } 
        }

    }
    else
    {
        nres_modmult(_MIPP_ mr_mip->w2,mr_mip->w4,mr_mip->w5);

        if (compare(mr_mip->w1,mr_mip->w3)>=0)
             expint(_MIPP_ 2,logb2(_MIPP_ mr_mip->w1)-1,mr_mip->w6);
        else expint(_MIPP_ 2,logb2(_MIPP_ mr_mip->w3)-1,mr_mip->w6);
        while (size(mr_mip->w6)!=0)
        {
            if (mr_mip->user!=NULL) (*mr_mip->user)();

            nres_modmult(_MIPP_ w,w,w);
            if (compare(mr_mip->w1,mr_mip->w6)>=0)
            {
                if (compare(mr_mip->w3,mr_mip->w6)>=0)
                {
                     nres_modmult(_MIPP_ w,mr_mip->w5,w);
                     subtract(_MIPP_ mr_mip->w3,mr_mip->w6,mr_mip->w3);
                }

                else nres_modmult(_MIPP_ w,mr_mip->w2,w);
                subtract(_MIPP_ mr_mip->w1,mr_mip->w6,mr_mip->w1);
            }
            else
            {
                if (compare(mr_mip->w3,mr_mip->w6)>=0)
                {
                     nres_modmult(_MIPP_ w,mr_mip->w4,w);
                     subtract(_MIPP_ mr_mip->w3,mr_mip->w6,mr_mip->w3);
                }
            }