mrmonty.c

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

            ASM add [esi+ebx+4],edx
            ASM adc eax,0
            ASM mov delay_carry,eax

#endif
#if INLINE_ASM == 4
   ASM (
           "movl %0,%%ecx\n"
           "movl %1,%%esi\n"
           "shll $2,%%esi\n"
           "movl %2,%%ebx\n"
           "addl %%esi,%%ebx\n"
           "movl (%%ebx),%%eax\n"
           "mull %3\n"
           "movl %%eax,%%edi\n"
           "movl %4,%%esi\n"
           "subl %%esi,%%ebx\n"
           "subl $4,%%ebx\n"
           "pushl %%ebp\n"
           "xorl %%ebp,%%ebp\n"
        "m1:\n"
           "movl (%%esi),%%eax\n"
           "addl $4,%%esi\n" 
           "mull %%edi\n"
           "addl %%ebp,%%eax\n"
           "movl (%%esi,%%ebx),%%ebp\n"
           "adcl $0,%%edx\n"
           "addl %%eax,%%ebp\n" 
           "adcl $0,%%edx\n"
           "movl %%ebp,(%%esi,%%ebx)\n"
           "decl %%ecx\n"
           "movl %%edx,%%ebp\n"
           "jnz m1\n"   

           "popl %%ebp\n"
           "movl %5,%%eax\n"
           "addl %%eax,4(%%esi,%%ebx)\n"
           "movl $0,%%eax\n"
           "adcl $0,%%eax\n"
           "addl %%edx,4(%%esi,%%ebx)\n"
           "adcl $0,%%eax\n"
           "movl %%eax,%5\n"
       
        :
        :"m"(rn),"m"(i),"m"(w0g),"m"(ndash),"m"(mg),"m"(delay_carry)
        :"eax","edi","esi","ebx","ecx","edx","memory"
       );
#endif

#ifndef INLINE_ASM
/*        muldvd(w0->w[i],ndash,0,&m);    Note that after this time   */
        m=ndash*w0->w[i];
        carry=0;                       /* around the loop, w0[i]=0    */

        for (j=0;j<rn;j++)
        {
#ifdef MR_NOASM 
            dble.d=(mr_large)m*modulus->w[j]+carry+w0->w[i+j];
            w0->w[i+j]=dble.h[MR_BOT];
            carry=dble.h[MR_TOP];
#else
            muldvd2(m,modulus->w[j],&carry,&w0->w[i+j]);
#endif
        }
        w0->w[rn+i]+=delay_carry;
        if (w0->w[rn+i]<delay_carry) delay_carry=1;
        else delay_carry=0;
        w0->w[rn+i]+=carry;
        if (w0->w[rn+i]<carry) delay_carry=1; 
#endif
      }
#endif
    }
    else for (i=0;i<rn;i++) 
    {
#ifdef MR_FP_ROUNDING
        imuldiv(w0->w[i],ndash,0,mr_mip->base,mr_mip->inverse_base,&m);
#else
        muldiv(w0->w[i],ndash,0,mr_mip->base,&m);
#endif
        carry=0;
        for (j=0;j<rn;j++)
        {
#ifdef MR_NOASM 
          dbled=(mr_large)m*modulus->w[j]+carry+w0->w[i+j];
#ifdef MR_FP_ROUNDING
          carry=(mr_small)MR_LROUND(dbled*mr_mip->inverse_base);
#else
#ifndef MR_FP
          if (mr_mip->base==mr_mip->base2)
              carry=(mr_small)(dbled>>mr_mip->lg2b);
          else 
#endif  
              carry=(mr_small)MR_LROUND(dbled/mr_mip->base);
#endif
          w0->w[i+j]=(mr_small)(dbled-(mr_large)carry*mr_mip->base);  
#else
#ifdef MR_FP_ROUNDING
          carry=imuldiv(modulus->w[j],m,w0->w[i+j]+carry,mr_mip->base,mr_mip->inverse_base,&w0->w[i+j]);
#else
          carry=muldiv(modulus->w[j],m,w0->w[i+j]+carry,mr_mip->base,&w0->w[i+j]);
#endif
#endif
        }
        w0->w[rn+i]+=(delay_carry+carry);
        delay_carry=0;
        if (w0->w[rn+i]>=mr_mip->base)
        {
            w0->w[rn+i]-=mr_mip->base;
            delay_carry=1; 
        }
    }
    w0->w[rn2]=delay_carry;
    w0->len=rn2+1;
    mr_shift(_MIPP_ w0,(-rn),w0);
    mr_lzero(w0);
    
    if (compare(w0,modulus)>=0) mr_psub(_MIPP_ w0,modulus,w0);
    copy(w0,y);
    MR_OUT
}

void nres_dotprod(_MIPD_ int n,big *x,big*y,big w)
{
    int i;
#ifdef MR_OS_THREADS
    miracl *mr_mip=get_mip();
#endif

    if (mr_mip->ERNUM) return;
    MR_IN(120)
    mr_mip->check=OFF;
    zero(mr_mip->w7);
    for (i=0;i<n;i++)
    {
        multiply(_MIPP_ x[i],y[i],mr_mip->w0);
        mr_padd(_MIPP_ mr_mip->w7,mr_mip->w0,mr_mip->w7);
    }
    mr_shift(_MIPP_ mr_mip->modulus,(int)mr_mip->modulus->len,mr_mip->w6);
        /* w6 = N.R */
    divide(_MIPP_ mr_mip->w7,mr_mip->w6,mr_mip->w6);
    redc(_MIPP_ mr_mip->w7,w);

    mr_mip->check=ON;
    MR_OUT
}

void nres_negate(_MIPD_ big x, big w)
{
#ifdef MR_OS_THREADS
    miracl *mr_mip=get_mip();
#endif
    if (mr_mip->ERNUM) return;

    MR_IN(92)
    if (size(x)==0) zero(w);
    else mr_psub(_MIPP_ mr_mip->modulus,x,w);    
    MR_OUT
}

void nres_modadd(_MIPD_ big x,big y,big w)
{ /* modular addition */
#ifdef MR_OS_THREADS
    miracl *mr_mip=get_mip();
#endif
#ifdef MR_COMBA
    if (mr_mip->ACTIVE)
    {
        comba_add(_MIPP_ x,y,w);
        return;
    }
    else
    {
#endif
        if (mr_mip->ERNUM) return;

        MR_IN(90)

        mr_padd(_MIPP_ x,y,mr_mip->w0);
        if (compare(mr_mip->w0,mr_mip->modulus)>=0) mr_psub(_MIPP_ mr_mip->w0,mr_mip->modulus,mr_mip->w0);
        copy(mr_mip->w0,w);

        MR_OUT
#ifdef MR_COMBA
    }
#endif
}

void nres_modsub(_MIPD_ big x,big y,big w)
{ /* modular subtraction */

#ifdef MR_OS_THREADS
    miracl *mr_mip=get_mip();
#endif
#ifdef MR_COMBA
    if (mr_mip->ACTIVE)
    {
        comba_sub(_MIPP_ x,y,w);
        return;
    }
    else
    {
#endif
        if (mr_mip->ERNUM) return;

        MR_IN(91)

        if (compare(x,y)>=0)
            mr_psub(_MIPP_ x,y,w);
        else
        {
            mr_psub(_MIPP_ y,x,w);
            mr_psub(_MIPP_ mr_mip->modulus,w,w);
        }
        MR_OUT
#ifdef MR_COMBA
    }
#endif

}

int nres_moddiv(_MIPD_ big x,big y,big w)
{ /* Modular division using n-residues w=x/y mod n */
    int gcd;
#ifdef MR_OS_THREADS
    miracl *mr_mip=get_mip();
#endif
    if (mr_mip->ERNUM) return 0;

    MR_IN(85)

    if (x==y)
    { /* Illegal parameter usage */
        mr_berror(_MIPP_ MR_ERR_BAD_PARAMETERS);
        MR_OUT
        
        return 0;
    }
    redc(_MIPP_ y,mr_mip->w6);
    gcd=xgcd(_MIPP_ mr_mip->w6,mr_mip->modulus,mr_mip->w6,mr_mip->w6,mr_mip->w6);
    
    if (gcd!=1) zero(w);
    else mad(_MIPP_ x,mr_mip->w6,x,mr_mip->modulus,mr_mip->modulus,w);

    MR_OUT
    return gcd;
}

void nres_premult(_MIPD_ big x,int k,big w)
{ /* multiply n-residue by small ordinary integer */
#ifdef MR_OS_THREADS
    miracl *mr_mip=get_mip();
#endif
    if (mr_mip->ERNUM) return;

    MR_IN(102)

    premult(_MIPP_ x,k,mr_mip->w0);
    divide(_MIPP_ mr_mip->w0,mr_mip->modulus,mr_mip->modulus);
    if (size(mr_mip->w0)<0) add(_MIPP_ mr_mip->w0,mr_mip->modulus,mr_mip->w0);
    copy(mr_mip->w0,w);

    MR_OUT
}

void nres_modmult(_MIPD_ big x,big y,big w)
{ /* Modular multiplication using n-residues w=x*y mod n */

#ifdef MR_OS_THREADS
    miracl *mr_mip=get_mip();
#endif
#ifdef MR_COMBA
    if (mr_mip->ACTIVE)
    {
        if (x==y) comba_square(_MIPP_ x,mr_mip->w0);
        else      comba_mult(_MIPP_ x,y,mr_mip->w0);
        comba_redc(_MIPP_ mr_mip->w0,w);
    }
    else
    {
#endif
#ifdef MR_KCM
    if (mr_mip->ACTIVE)
    {
        if (x==y) kcm_sqr(_MIPP_ x,mr_mip->w0);
        else      kcm_mul(_MIPP_ x,y,mr_mip->w0);
        kcm_redc(_MIPP_ mr_mip->w0,w);
    }
    else
    { 
#endif
#ifdef MR_PENTIUM
    if (mr_mip->ACTIVE)
    {
        if (x==y) fastmodsquare(_MIPP_ x,w);
        else      fastmodmult(_MIPP_ x,y,w);
    }
    else
    { 
#endif
        if (mr_mip->ERNUM) return;

        MR_IN(83)

        mr_mip->check=OFF;
        multiply(_MIPP_ x,y,mr_mip->w0);
        redc(_MIPP_ mr_mip->w0,w);
        mr_mip->check=ON;
        MR_OUT
#ifdef MR_COMBA
}
#endif
#ifdef MR_KCM
}
#endif
#ifdef MR_PENTIUM
}
#endif

}

/* Montgomery's trick for finding multiple   *
 * simultaneous modular inverses             *
 * Based on the observation that             *
 *           1/x = yz*(1/xyz)                *
 *           1/y = xz*(1/xyz)                *
 *           1/z = xy*(1/xyz)                *
 * Why are all of Peter Montgomery's clever  *
 * algorithms always described as "tricks" ??*/


BOOL nres_double_inverse(_MIPD_ big x,big y,big w,big z)
{ /* find y=1/x mod n and z=1/w mod n */
  /* 1/x = w/xw, and 1/w = x/xw       */
#ifdef MR_OS_THREADS
    miracl *mr_mip=get_mip();
#endif
    MR_IN(145)

    nres_modmult(_MIPP_ x,w,mr_mip->w6);  /* xw */

    if (size(mr_mip->w6)==0)
    {
        mr_berror(_MIPP_ MR_ERR_DIV_BY_ZERO);
        MR_OUT
        return FALSE;
    }
    redc(_MIPP_ mr_mip->w6,mr_mip->w6);
    redc(_MIPP_ mr_mip->w6,mr_mip->w6);
    xgcd(_MIPP_ mr_mip->w6,mr_mip->modulus,mr_mip->w6,mr_mip->w6,mr_mip->w6);

    nres_modmult(_MIPP_ w,mr_mip->w6,mr_mip->w5);
    nres_modmult(_MIPP_ x,mr_mip->w6,z);
    copy(mr_mip->w5,y);

    MR_OUT
    return TRUE;
}

BOOL nres_multi_inverse(_MIPD_ int m,big *x,big *w)
{ /* find w[i]=1/x[i] mod n, for i=0 to m-1 *
   * x and w MUST be distinct               */
    int i;
#ifdef MR_OS_THREADS
    miracl *mr_mip=get_mip();
#endif
    if (m==0) return TRUE;
    if (m<0) return FALSE;
    MR_IN(118)

    if (x==w)
    {
        mr_berror(_MIPP_ MR_ERR_BAD_PARAMETERS);
        MR_OUT
        return FALSE;
    }

    if (m==1)
    {
        convert(_MIPP_ 1,w[0]);
        nres(_MIPP_ w[0],w[0]);
        nres_moddiv(_MIPP_ w[0],x[0],w[0]);
        MR_OUT
        return TRUE;
    }

    convert(_MIPP_ 1,w[0]);
    copy(x[0],w[1]);
    for (i=2;i<m;i++)
        nres_modmult(_MIPP_ w[i-1],x[i-1],w[i]); 

    nres_modmult(_MIPP_ w[m-1],x[m-1],mr_mip->w6);  /* y=x[0]*x[1]*x[2]....x[m-1] */
    if (size(mr_mip->w6)==0)
    {
        mr_berror(_MIPP_ MR_ERR_DIV_BY_ZERO);
        MR_OUT
        return FALSE;
    }

    redc(_MIPP_ mr_mip->w6,mr_mip->w6);
    redc(_MIPP_ mr_mip->w6,mr_mip->w6);
    xgcd(_MIPP_ mr_mip->w6,mr_mip->modulus,mr_mip->w6,mr_mip->w6,mr_mip->w6);

/* Now y=1/y */

    copy(x[m-1],mr_mip->w5);
    nres_modmult(_MIPP_ w[m-1],mr_mip->w6,w[m-1]);

    for (i=m-2;;i--)
    {
        if (i==0)
        {
            nres_modmult(_MIPP_ mr_mip->w5,mr_mip->w6,w[0]);
            break;
        }
        nres_modmult(_MIPP_ w[i],mr_mip->w5,w[i]);
        nres_modmult(_MIPP_ w[i],mr_mip->w6,w[i]);
        nres_modmult(_MIPP_ mr_mip->w5,x[i],mr_mip->w5);
    }

    MR_OUT 
    return TRUE;   
}