mrmonty.c

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

/*
 *   MIRACL Montgomery's method for modular arithmetic without division.
 *   mrmonty.c 
 *
 *   Programs to implement Montgomery's method
 *   See "Modular Multiplication Without Trial Division", Math. Comp. 
 *   Vol 44, Number 170, April 1985, Pages 519-521
 *   NOTE - there is an important correction to this paper mentioned as a
 *   footnote in  "Speeding the Pollard and Elliptic Curve Methods", 
 *   Math. Comput., Vol. 48, January 1987, 243-264
 *
 *   The advantage of this approach is that no division required in order
 *   to compute a modular reduction - useful if division is slow
 *   e.g. on a SPARC processor, or a DSP.
 *
 *   NOTE:- nres_modmult() *might* use mr_mip->w7
 *   
 *   The disadvantage is that numbers must first be converted to an internal
 *   "n-residue" form.
 *
 *   Copyright (c) 1988-2001 Shamus Software Ltd.
 */

#include <stdlib.h> 
#include "miracl.h"

#define mr_abs(x)  ((x)<0? (-(x)) : (x))

void kill_monty(_MIPDO_ )
{
#ifdef MR_OS_THREADS
    miracl *mr_mip=get_mip();
#endif
    zero(mr_mip->modulus);
#ifdef MR_KCM
    zero(mr_mip->big_ndash);
#endif
}

mr_small prepare_monty(_MIPD_ big n)
{ /* prepare Montgomery modulus */ 
#ifdef MR_KCM
    int nl;
#endif
#ifdef MR_PENTIUM
    mr_small ndash;
    mr_small base;
    mr_small magic=13835058055282163712.0;   
    int control=0x1FFF;
#endif
#ifdef MR_OS_THREADS
    miracl *mr_mip=get_mip();
#endif
    if (mr_mip->ERNUM) return (mr_small)0;
/* Is it set-up already? */
    if (size(mr_mip->modulus)!=0)
        if (compare(n,mr_mip->modulus)==0) return mr_mip->ndash;

    MR_IN(80)

    if (size(n)<=2) 
    {
        mr_berror(_MIPP_ MR_ERR_BAD_MODULUS);
        MR_OUT
        return (mr_small)0;
    }

    zero(mr_mip->w6);
    zero(mr_mip->w16);

/* set a small negative QNR (on the assumption that n is prime!) */

    mr_mip->pmod8=remain(_MIPP_ n,8);
    switch (mr_mip->pmod8)
    {
    case 0:
    case 1:
    case 2:
    case 4:
    case 6:
        mr_mip->qnr=0;  /* none defined */
        break;
    case 3:
        mr_mip->qnr=-1;
        break;
    case 5:
        mr_mip->qnr=-2;
        break;
    case 7:
        mr_mip->qnr=-2;
        break;
    }

#ifdef MR_PENTIUM

mr_mip->ACTIVE=FALSE;
if (mr_mip->base!=0)
    if (MR_PENTIUM==n->len) mr_mip->ACTIVE=TRUE;
    if (MR_PENTIUM<0)
    {
        if (n->len<=(-MR_PENTIUM)) mr_mip->ACTIVE=TRUE;
        if (logb2(_MIPP_ n)%mr_mip->lg2b==0) mr_mip->ACTIVE=FALSE;
    }
#endif

    mr_mip->MONTY=ON;
#ifdef MR_COMBA
    mr_mip->ACTIVE=FALSE;
    if (MR_COMBA==n->len && mr_mip->base==0) 
    {
        mr_mip->ACTIVE=TRUE;
#ifdef MR_SPECIAL
        mr_mip->MONTY=OFF;      /* "special" modulus reduction */
#endif                          /* implemented in mrcomba.c    */
    }

#endif

    if (!mr_mip->MONTY)
    { /* Montgomery arithmetic is turned off */
        copy(n,mr_mip->modulus);
        mr_mip->ndash=0;
        MR_OUT
        return (mr_small)0;
    }

#ifdef MR_KCM
  
/* test for base==0 & n->len=MR_KCM.2^x */

    mr_mip->ACTIVE=FALSE;
    if (mr_mip->base==0)
    {
        nl=(int)n->len;
        while (nl>=MR_KCM)
        {
            if (nl==MR_KCM)
            {
                mr_mip->ACTIVE=TRUE;
                break;
            }
            if (nl%2!=0) break;
            nl/=2;
        }
    }  
    if (mr_mip->ACTIVE)
    {
        mr_mip->w6->len=n->len+1;
        mr_mip->w6->w[n->len]=1;
        if (xgcd(_MIPP_ n,mr_mip->w6,mr_mip->w14,mr_mip->w14,mr_mip->w14)!=1)
        { /* problems */
            mr_berror(_MIPP_ MR_ERR_BAD_MODULUS);
            MR_OUT
            return (mr_small)0;
        }
    }
    else
    {
#endif
        mr_mip->w6->len=2;
        mr_mip->w6->w[0]=0;
        mr_mip->w6->w[1]=1;    /* w6 = base */
        mr_mip->w16->len=1;
        mr_mip->w16->w[0]=n->w[0];  /* w16 = n mod base */
        if (xgcd(_MIPP_ mr_mip->w16,mr_mip->w6,mr_mip->w14,mr_mip->w14,mr_mip->w14)!=1)
        { /* problems */
            mr_berror(_MIPP_ MR_ERR_BAD_MODULUS);
            MR_OUT
            return (mr_small)0;
        }
#ifdef MR_KCM
    }
    copy(mr_mip->w14,mr_mip->big_ndash);
#endif
    mr_mip->ndash=mr_mip->base-mr_mip->w14->w[0]; /* = N' mod b */
    copy(n,mr_mip->modulus);
#ifdef MR_PENTIUM
/* prime the FP stack */
    if (mr_mip->ACTIVE)
    {
        ndash=mr_mip->ndash;
        base=mr_mip->base;
        magic *=base;
        ASM
        {
            finit
            fldcw WORD PTR control
            fld QWORD PTR ndash
            fld1
            fld QWORD PTR base
            fdiv
            fld QWORD PTR magic
        }
    }
#endif
    MR_OUT

    return mr_mip->ndash;
}

void nres(_MIPD_ big x,big y)
{ /* convert x to n-residue format */
#ifdef MR_OS_THREADS
    miracl *mr_mip=get_mip();
#endif
    if (mr_mip->ERNUM) return;

    MR_IN(81)

    if (size(mr_mip->modulus)==0)
    {
        mr_berror(_MIPP_ MR_ERR_NO_MODULUS);
        MR_OUT
        return;
    }
    copy(x,y);
    divide(_MIPP_ y,mr_mip->modulus,mr_mip->modulus);
    if (size(y)<0) add(_MIPP_ y,mr_mip->modulus,y);
    if (!mr_mip->MONTY) 
    {
        MR_OUT
        return;
    }
    mr_mip->check=OFF;

    mr_shift(_MIPP_ y,(int)mr_mip->modulus->len,mr_mip->w0);
    divide(_MIPP_ mr_mip->w0,mr_mip->modulus,mr_mip->modulus);
    mr_mip->check=ON;
    copy(mr_mip->w0,y);

    MR_OUT
}

void redc(_MIPD_ big x,big y)
{ /* Montgomery's REDC function p. 520 */
  /* also used to convert n-residues back to normal form */
    mr_small carry,delay_carry,m,ndash,*w0g,*mg;

#ifdef MR_ITANIUM
    mr_small tm;
#endif
    int i,j,rn,rn2;
    big w0,modulus;
#ifdef MR_NOASM
    union doubleword dble;
    mr_large dbled,ldres;
#endif
#ifdef MR_OS_THREADS
    miracl *mr_mip=get_mip();
#endif
    if (mr_mip->ERNUM) return;

    MR_IN(82)

    w0=mr_mip->w0;        /* get these into local variables (for inline assembly) */
    modulus=mr_mip->modulus;
    ndash=mr_mip->ndash;

    copy(x,w0);
    if (!mr_mip->MONTY)
    {
        divide(_MIPP_ w0,modulus,modulus);
        copy(w0,y);
        MR_OUT
        return;
    }
    delay_carry=0;
    rn=(int)modulus->len;
    rn2=rn+rn;
    if (mr_mip->base==0) 
    {
#ifndef MR_NOFULLWIDTH
      mg=modulus->w;
      w0g=w0->w;
      for (i=0;i<rn;i++)
      {
       /* inline - substitutes for loop below */
#if INLINE_ASM == 1
            ASM cld
            ASM mov cx,rn
            ASM mov si,i
            ASM shl si,1
#ifdef MR_LMM
            ASM push ds
            ASM push es
            ASM les bx,DWORD PTR w0g
            ASM add bx,si
            ASM mov ax,es:[bx]
            ASM mul WORD PTR ndash
            ASM mov di,ax
            ASM lds si,DWORD PTR mg
#else
            ASM mov bx,w0g
            ASM add bx,si
            ASM mov ax,[bx]
            ASM mul WORD PTR ndash
            ASM mov di,ax
            ASM mov si,mg
#endif
            ASM push bp
            ASM xor bp,bp
          m1:
            ASM lodsw
            ASM mul di
            ASM add ax,bp
            ASM adc dx,0
#ifdef MR_LMM
            ASM add es:[bx],ax
#else
            ASM add [bx],ax
#endif
            ASM adc dx,0
            ASM inc bx
            ASM inc bx
            ASM mov bp,dx
            ASM loop m1

            ASM pop bp
            ASM mov ax,delay_carry     
#ifdef MR_LMM
            ASM add es:[bx],ax
            ASM mov ax,0
            ASM adc ax,0
            ASM add es:[bx],dx
            ASM pop es
            ASM pop ds
#else
            ASM add [bx],ax
            ASM mov ax,0
            ASM adc ax,0
            ASM add [bx],dx
#endif
            ASM adc ax,0
            ASM mov delay_carry,ax
#endif
#if INLINE_ASM == 2
            ASM cld
            ASM mov cx,rn
            ASM mov si,i
            ASM shl si,2
#ifdef MR_LMM
            ASM push ds
            ASM push es
            ASM les bx,DWORD PTR w0g
            ASM add bx,si
            ASM mov eax,es:[bx]
            ASM mul DWORD PTR ndash
            ASM mov edi,eax
            ASM lds si,DWORD PTR mg
#else
            ASM mov bx,w0g
            ASM add bx,si
            ASM mov eax,[bx]
            ASM mul DWORD PTR ndash
            ASM mov edi,eax
            ASM mov si,mg
#endif
            ASM push ebp
            ASM xor ebp,ebp
          m1:
            ASM lodsd
            ASM mul edi
            ASM add eax,ebp
            ASM adc edx,0
#ifdef MR_LMM
            ASM add es:[bx],eax
#else
            ASM add [bx],eax
#endif
            ASM adc edx,0
            ASM add bx,4
            ASM mov ebp,edx
            ASM loop m1

            ASM pop ebp
            ASM mov eax,delay_carry    
#ifdef MR_LMM
            ASM add es:[bx],eax
            ASM mov eax,0
            ASM adc eax,0
            ASM add es:[bx],edx
            ASM pop es
            ASM pop ds
#else 
            ASM add [bx],eax
            ASM mov eax,0
            ASM adc eax,0
            ASM add [bx],edx
#endif
            ASM adc eax,0
            ASM mov delay_carry,eax

#endif
#if INLINE_ASM == 3
            ASM mov ecx,rn
            ASM mov esi,i
            ASM shl esi,2
            ASM mov ebx,w0g
            ASM add ebx,esi
            ASM mov eax,[ebx]
            ASM mul DWORD PTR ndash
            ASM mov edi,eax
            ASM mov esi,mg
            ASM sub ebx,esi
            ASM sub ebx,4
            ASM push ebp
            ASM xor ebp,ebp
          m1:
            ASM mov eax,[esi]
            ASM add esi,4
            ASM mul edi
            ASM add eax,ebp
            ASM mov ebp,[esi+ebx]
            ASM adc edx,0
            ASM add ebp,eax
            ASM adc edx,0
            ASM mov [esi+ebx],ebp
            ASM dec ecx
            ASM mov ebp,edx
            ASM jnz m1

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