mrarth2.c

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

/*
 *   MIRACL arithmetic routines 2 - multiplying and dividing BIG NUMBERS.
 *   mrarth2.c
 *
 *   Copyright (c) 1988-2003 Shamus Software Ltd.
 */

#include "miracl.h"

/* If a number has more than this number of digits, then squaring is faster */

#define SQR_FASTER_THRESHOLD 5

mr_small normalise(_MIPD_ big x,big y)
{ /* normalise divisor */
    mr_small norm,r;
#ifdef MR_FP
    mr_small dres;
#endif
    int len;
#ifdef MR_OS_THREADS
    miracl *mr_mip=get_mip();
#endif

    MR_IN(4)

    if (x!=y) copy(x,y);
    len=(int)(y->len&MR_OBITS);
    if (mr_mip->base==0)
    {
#ifndef MR_NOFULLWIDTH
        if ((r=y->w[len-1]+1)==0) norm=1;
#ifdef MR_NOASM
        else norm=(mr_small)(((mr_large)1 << MIRACL)/r);
#else
        else norm=muldvm((mr_small)1,(mr_small)0,r,&r);
#endif
        if (norm!=1) mr_pmul(_MIPP_ y,norm,y);
#endif
    }
    else
    {
        norm=MR_DIV(mr_mip->base,(mr_small)(y->w[len-1]+1));   
        if (norm!=1) mr_pmul(_MIPP_ y,norm,y);
    }
    MR_OUT
    return norm;
}

void multiply(_MIPD_ big x,big y,big z)
{  /*  multiply two big numbers: z=x.y  */
    int i,xl,yl,j,ti;
    mr_small carry,*xg,*yg,*w0g;

#ifdef MR_ITANIUM
    mr_small tm;
#endif
    mr_unsign32 sz;
    big w0;
#ifdef MR_NOASM
    union doubleword dble;
    mr_large dbled;
    mr_large ldres;
#endif
#ifdef MR_OS_THREADS
    miracl *mr_mip=get_mip();
#endif
    if (mr_mip->ERNUM) return;
    if (y->len==0 || x->len==0) 
    {
        zero(z);
        return;
    }
    w0=mr_mip->w0;    /* local pointer */

    MR_IN(5)

#ifdef MR_FLASH
    if (mr_notint(x) || mr_notint(y))
    {
        mr_berror(_MIPP_ MR_ERR_INT_OP);
        MR_OUT
        return;
    }
#endif
    sz=((x->len&MR_MSBIT)^(y->len&MR_MSBIT));
    xl=(int)(x->len&MR_OBITS);
    yl=(int)(y->len&MR_OBITS);
    zero(w0);
    if (mr_mip->check && xl+yl>mr_mip->nib)
    {
        mr_berror(_MIPP_ MR_ERR_OVERFLOW);
        MR_OUT
        return;
    }
    if (mr_mip->base==0)
    {
#ifndef MR_NOFULLWIDTH
        xg=x->w; yg=y->w; w0g=w0->w; 
        if (x==y && xl>SQR_FASTER_THRESHOLD)    
                             /* extra hassle make it not    */
                             /* worth it for small numbers */
        { /* fast squaring */
            for (i=0;i<xl-1;i++)
            {  /* long multiplication */
           /* inline - substitutes for loop below */
#if INLINE_ASM == 1
                ASM cld
                ASM mov dx,i
                ASM mov cx,xl
                ASM sub cx,dx
                ASM dec cx
                ASM shl dx,1
#ifdef MR_LMM
                ASM push ds
                ASM push es
                ASM les bx,DWORD PTR w0g
                ASM lds si,DWORD PTR xg
                ASM add si,dx
                ASM mov di,[si]
#else
                ASM mov bx,w0g
                ASM mov si,xg
                ASM add si,dx
                ASM mov di,[si]
#endif
                ASM add bx,dx
                ASM add bx,dx 
                ASM add bx,2
                ASM add si,2
                ASM push bp
                ASM xor bp,bp
              tcl4:
                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 tcl4

#ifdef MR_LMM
                ASM mov es:[bx],bp
                ASM pop bp
                ASM pop es
                ASM pop ds
#else
                ASM mov [bx],bp
                ASM pop bp
#endif
#endif
#if INLINE_ASM == 2
                ASM cld
                ASM mov dx,i
                ASM mov cx,xl
                ASM sub cx,dx
                ASM dec cx
                ASM shl dx,2
#ifdef MR_LMM
                ASM push ds
                ASM push es
                ASM les bx,DWORD PTR w0g
                ASM lds si,DWORD PTR xg
                ASM add si,dx
                ASM mov edi,[si]
#else
                ASM mov bx,w0g
                ASM mov si,xg
                ASM add si,dx
                ASM mov edi,[si]
#endif
                ASM add bx,dx
                ASM add bx,dx
                ASM add bx,4
                ASM add si,4   
                ASM push ebp
                ASM xor ebp,ebp
              tcl4:
                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 tcl4

#ifdef MR_LMM
                ASM mov es:[bx],ebp
                ASM pop ebp
                ASM pop es
                ASM pop ds
#else
                ASM mov [bx],ebp
                ASM pop ebp
#endif
#endif
#if INLINE_ASM == 3
                ASM mov esi,i
                ASM mov ecx,xl
                ASM sub ecx,esi
                ASM dec ecx
                ASM shl esi,2
                ASM mov edx, xg
                ASM mov ebx,edx
                ASM add ebx,esi
                ASM mov edi,[ebx]
                ASM mov ebx,w0g
                ASM add ebx,esi
                ASM add esi,edx
                ASM sub ebx,edx
                ASM add esi,4  
                ASM sub ebx,4
                ASM push ebp
                ASM xor ebp,ebp
              tcl4:
                ASM mov eax,[esi]   /* optimized for Pentium */
                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 tcl4

                ASM mov [esi+ebx+4],ebp
                ASM pop ebp
#endif
#if INLINE_ASM == 4
   ASM (
           "movl %0,%%esi\n"
           "movl %1,%%ecx\n"
           "subl %%esi,%%ecx\n"
           "decl %%ecx\n"
           "shll $2,%%esi\n"
           "movl %2,%%edx\n"
           "movl %%edx,%%ebx\n"
           "addl %%esi,%%ebx\n"
           "movl (%%ebx),%%edi\n"
           "movl %3,%%ebx\n"
           "addl %%esi,%%ebx\n"
           "addl %%edx,%%esi\n"
           "subl %%edx,%%ebx\n"
           "addl $4,%%esi\n"
           "subl $4,%%ebx\n"
           "pushl %%ebp\n"
           "xorl %%ebp,%%ebp\n"
         "tcl4:\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 tcl4\n"

           "movl %%ebp,4(%%esi,%%ebx)\n"
           "popl %%ebp\n"
           
        :
        :"m"(i),"m"(xl),"m"(xg),"m"(w0g)
        :"eax","edi","esi","ebx","ecx","edx","memory"
       );
#endif
#ifndef INLINE_ASM
                carry=0;
                for (j=i+1;j<xl;j++)
                { /* Only do above the diagonal */
#ifdef MR_NOASM
                    dble.d=(mr_large)x->w[i]*x->w[j]+carry+w0->w[i+j];
                    w0->w[i+j]=dble.h[MR_BOT];
                    carry=dble.h[MR_TOP];
#else
                    muldvd2(x->w[i],x->w[j],&carry,&w0->w[i+j]);
#endif
                }
                w0->w[xl+i]=carry;
#endif
            }
#if INLINE_ASM == 1
            ASM mov cx,xl
            ASM shl cx,1
#ifdef MR_LMM
            ASM push ds
            ASM push es
            ASM les bx,DWORD PTR w0g  
#else
            ASM mov bx,w0g
#endif
          tcl5:
#ifdef MR_LMM
            ASM rcl WORD PTR es:[bx],1
#else
            ASM rcl WORD PTR [bx],1
#endif
            ASM inc bx
            ASM inc bx
            ASM loop tcl5

            ASM cld
            ASM mov cx,xl
#ifdef MR_LMM
            ASM les di,DWORD PTR w0g
            ASM lds si,DWORD PTR xg
#else
            ASM mov di,w0g
            ASM mov si,xg
#endif
       
            ASM xor bx,bx
          tcl7:
            ASM lodsw
            ASM mul ax
            ASM add ax,bx
            ASM adc dx,0
#ifdef MR_LMM
            ASM add es:[di],ax
#else
            ASM add [di],ax
#endif
            ASM adc dx,0
            ASM xor bx,bx
            ASM inc di
            ASM inc di
#ifdef MR_LMM
            ASM add es:[di],dx
#else
            ASM add [di],dx
#endif
            ASM adc bx,0
            ASM inc di
            ASM inc di
            ASM loop tcl7
#ifdef MR_LMM
            ASM pop es
            ASM pop ds
#endif
#endif
#if INLINE_ASM == 2
            ASM mov cx,xl
            ASM shl cx,1
#ifdef MR_LMM
            ASM push ds
            ASM push es
            ASM les bx,DWORD PTR w0g
#else
            ASM mov bx,w0g
#endif
          tcl5:
#ifdef MR_LMM
            ASM rcl DWORD PTR es:[bx],1
#else
            ASM rcl DWORD PTR [bx],1
#endif
            ASM inc bx
            ASM inc bx
            ASM inc bx
            ASM inc bx
            ASM loop tcl5

            ASM cld
            ASM mov cx,xl
#ifdef MR_LMM
            ASM les di,DWORD PTR w0g
            ASM lds si,DWORD PTR xg
#else
            ASM mov di,w0g
            ASM mov si,xg
#endif
            ASM xor ebx,ebx
          tcl7:
            ASM lodsd
            ASM mul eax
            ASM add eax,ebx
            ASM adc edx,0
#ifdef MR_LMM
            ASM add es:[di],eax
#else
            ASM add [di],eax
#endif
            ASM adc edx,0
            ASM xor ebx,ebx
            ASM add di,4
#ifdef MR_LMM
            ASM add es:[di],edx
#else
            ASM add [di],edx
#endif
            ASM adc ebx,0
            ASM add di,4
            ASM loop tcl7
#ifdef MR_LMM
            ASM pop es
            ASM pop ds
#endif
#endif
#if INLINE_ASM == 3
            ASM mov ecx,xl
            ASM shl ecx,1
            ASM mov edi,w0g
          tcl5:
            ASM rcl DWORD PTR [edi],1
            ASM inc edi
            ASM inc edi
            ASM inc edi
            ASM inc edi
            ASM loop tcl5

            ASM mov ecx,xl
            ASM mov esi,xg
            ASM mov edi,w0g
            ASM xor ebx,ebx
          tcl7:
            ASM mov eax,[esi]
            ASM add esi,4
            ASM mul eax
            ASM add eax,ebx
            ASM adc edx,0
            ASM add [edi],eax
            ASM adc edx,0
            ASM xor ebx,ebx
            ASM add edi,4
            ASM add [edi],edx
            ASM adc ebx,0
            ASM add edi,4
            ASM dec ecx
            ASM jnz tcl7
#endif
#if INLINE_ASM == 4
   ASM (
           "movl %0,%%ecx\n"
           "shll $1,%%ecx\n"
           "movl %1,%%edi\n"
         "tcl5:\n"
           "rcll $1,(%%edi)\n"
           "incl %%edi\n"
           "incl %%edi\n"
           "incl %%edi\n"
           "incl %%edi\n"
           "loop tcl5\n"

           "movl %0,%%ecx\n"
           "movl %2,%%esi\n"
           "movl %1,%%edi\n"
           "xorl %%ebx,%%ebx\n"
         "tcl7:\n"
           "movl (%%esi),%%eax\n"
           "addl $4,%%esi\n"
           "mull %%eax\n"
           "addl %%ebx,%%eax\n"
           "adcl $0,%%edx\n"
           "addl %%eax,(%%edi)\n"
           "adcl $0,%%edx\n"
           "xorl %%ebx,%%ebx\n"
           "addl $4,%%edi\n"
           "addl %%edx,(%%edi)\n"
           "adcl $0,%%ebx\n"
           "addl $4,%%edi\n"
           "decl %%ecx\n"
           "jnz tcl7\n"                       
        :
        :"m"(xl),"m"(w0g),"m"(xg)
        :"eax","edi","esi","ebx","ecx","edx","memory"
       );
#endif
#ifndef INLINE_ASM
            w0->len=xl+xl-1;
            mr_padd(_MIPP_ w0,w0,w0);     /* double it */
            carry=0;
            for (i=0;i<xl;i++)
            { /* add in squared elements */
                ti=i+i;
#ifdef MR_NOASM               
                dble.d=(mr_large)x->w[i]*x->w[i]+carry+w0->w[ti];
                w0->w[ti]=dble.h[MR_BOT];
                carry=dble.h[MR_TOP];
#else
                muldvd2(x->w[i],x->w[i],&carry,&w0->w[ti]);
#endif
                w0->w[ti+1]+=carry;
                if (w0->w[ti+1]<carry) carry=1;
                else                   carry=0;
            }
#endif