mrarth2.c

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

        }
        else for (i=0;i<xl;i++)
        { /* long multiplication */
       /* inline - substitutes for loop below */
#if INLINE_ASM == 1
            ASM cld
            ASM mov cx,yl
            ASM mov dx,i
            ASM shl dx,1
#ifdef MR_LMM
            ASM push ds
            ASM push es
            ASM les bx,DWORD PTR w0g
            ASM add bx,dx
            ASM lds si,DWORD PTR xg
            ASM add si,dx
            ASM mov di,[si]
            ASM lds si,DWORD PTR yg
#else
            ASM mov bx,w0g
            ASM add bx,dx
            ASM mov si,xg
            ASM add si,dx
            ASM mov di,[si]
            ASM mov si,yg
#endif
            ASM push bp
            ASM xor bp,bp
          tcl6:
            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 inc bx
            ASM inc bx
            ASM adc dx,0
            ASM mov bp,dx
            ASM loop tcl6

#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 cx,yl
            ASM mov dx,i
            ASM shl dx,2
#ifdef MR_LMM
            ASM push ds
            ASM push es
            ASM les bx,DWORD PTR w0g
            ASM add bx,dx
            ASM lds si,DWORD PTR xg
            ASM add si,dx
            ASM mov edi,[si]
            ASM lds si,DWORD PTR yg
#else           
            ASM mov bx,w0g
            ASM add bx,dx
            ASM mov si,xg
            ASM add si,dx
            ASM mov edi,[si]
            ASM mov si,yg
#endif
            ASM push ebp
            ASM xor ebp,ebp
          tcl6:
            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 tcl6

#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 ecx,yl
            ASM mov esi,i
            ASM shl esi,2
            ASM mov ebx,xg
            ASM add ebx,esi
            ASM mov edi,[ebx]
            ASM mov ebx,w0g
            ASM add ebx,esi
            ASM mov esi,yg
            ASM sub ebx,esi
            ASM sub ebx,4
            ASM push ebp
            ASM xor ebp,ebp
          tcl6:
            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 tcl6

            ASM mov [esi+ebx+4],ebp
            ASM pop ebp
#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),%%edi\n"
           "movl %3,%%ebx\n"
           "addl %%esi,%%ebx\n"
           "movl %4,%%esi\n"
           "subl %%esi,%%ebx\n"
           "subl $4,%%ebx\n"
           "pushl %%ebp\n"
           "xorl %%ebp,%%ebp\n"
         "tcl6:\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 tcl6\n"   

           "movl %%ebp,4(%%esi,%%ebx)\n"
           "popl %%ebp\n"
        :
        :"m"(yl),"m"(i),"m"(xg),"m"(w0g),"m"(yg)
        :"eax","edi","esi","ebx","ecx","edx","memory"
       );
#endif
#ifndef INLINE_ASM
            carry=0;
            for (j=0;j<yl;j++)
            { /* multiply each digit of y by x[i] */
#ifdef MR_NOASM 
                dble.d=(mr_large)x->w[i]*y->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],y->w[j],&carry,&w0->w[i+j]);
#endif
            }
            w0->w[yl+i]=carry;
#endif
        }
#endif
    }
    else
    {
        if (x==y && xl>SQR_FASTER_THRESHOLD)
        { /* squaring can be done nearly twice as fast */
            for (i=0;i<xl-1;i++)
            { /* long multiplication */
                carry=0;
                for (j=i+1;j<xl;j++)
                { /* Only do above the diagonal */
#ifdef MR_NOASM
                   dbled=(mr_large)x->w[i]*x->w[j]+w0->w[i+j]+carry;
  #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(x->w[i],x->w[j],w0->w[i+j]+carry,mr_mip->base,mr_mip->inverse_base,&w0->w[i+j]);
  #else
              carry=muldiv(x->w[i],x->w[j],w0->w[i+j]+carry,mr_mip->base,&w0->w[i+j]); 
  #endif
#endif
                }
                w0->w[xl+i]=carry;
            }
            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
                dbled=(mr_large)x->w[i]*x->w[i]+w0->w[ti]+carry;
#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[ti]=(mr_small)(dbled-(mr_large)carry*mr_mip->base);
#else

#ifdef MR_FP_ROUNDING
                carry=imuldiv(x->w[i],x->w[i],w0->w[ti]+carry,mr_mip->base,mr_mip->inverse_base,&w0->w[ti]);
#else
                carry=muldiv(x->w[i],x->w[i],w0->w[ti]+carry,mr_mip->base,&w0->w[ti]);
#endif

#endif
                w0->w[ti+1]+=carry;
                carry=0;
                if (w0->w[ti+1]>=mr_mip->base)
                {
                    carry=1;
                    w0->w[ti+1]-=mr_mip->base;
                }
            }
        }
        else for (i=0;i<xl;i++)
        { /* long multiplication */
            carry=0; 
            for (j=0;j<yl;j++)
            { /* multiply each digit of y by x[i] */
#ifdef MR_NOASM
                dbled=(mr_large)x->w[i]*y->w[j]+w0->w[i+j]+carry;

#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(x->w[i],y->w[j],w0->w[i+j]+carry,mr_mip->base,mr_mip->inverse_base,&w0->w[i+j]);
#else
                carry=muldiv(x->w[i],y->w[j],w0->w[i+j]+carry,mr_mip->base,&w0->w[i+j]);
#endif

#endif
            }
            w0->w[yl+i]=carry;
        }
    }
    w0->len=(sz|(xl+yl)); /* set length and sign of result */

    mr_lzero(w0);
    copy(w0,z);
    MR_OUT
}

void divide(_MIPD_ big x,big y,big z)
{  /*  divide two big numbers  z=x/y : x=x mod y  *
    *  returns quotient only if  divide(x,y,x)    *
    *  returns remainder only if divide(x,y,y)    */
    mr_small carry,attemp,ldy,sdy,ra,r,d,tst,psum;
#ifdef MR_FP
    mr_small dres;
#endif
    mr_unsign32 sx,sy,sz;
    mr_small borrow,dig,*w0g,*yg;
    int i,k,m,x0,y0,w00;
    big w0;

#ifdef MR_ITANIUM
    mr_small tm;
#endif
#ifdef MR_NOASM
    union doubleword dble;
    mr_large dbled;
    mr_large ldres;
#endif
    BOOL check;
#ifdef MR_OS_THREADS
    miracl *mr_mip=get_mip();
#endif
    if (mr_mip->ERNUM) return;
    w0=mr_mip->w0;

    MR_IN(6)

    if (x==y) mr_berror(_MIPP_ MR_ERR_BAD_PARAMETERS);
#ifdef MR_FLASH
    if (mr_notint(x) || mr_notint(y)) mr_berror(_MIPP_ MR_ERR_INT_OP);
#endif
    if (y->len==0) mr_berror(_MIPP_ MR_ERR_DIV_BY_ZERO);
    if (mr_mip->ERNUM)
    {
        MR_OUT
        return;
    }
    sx=(x->len&MR_MSBIT);   /*  extract signs ... */
    sy=(y->len&MR_MSBIT);
    sz=(sx^sy);
    x->len&=MR_OBITS;   /* ... and force operands to positive  */
    y->len&=MR_OBITS;
    x0=(int)x->len;
    y0=(int)y->len;
    copy(x,w0);
    w00=(int)w0->len;
    if (mr_mip->check && (w00-y0+1>mr_mip->nib))
    {
        mr_berror(_MIPP_ MR_ERR_OVERFLOW);
        MR_OUT
        return;
    }
    d=0;
    if (x0==y0)
    {
        if (x0==1) /* special case - x and y are both mr_smalls */
        { 
            d=MR_DIV(w0->w[0],y->w[0]);
            w0->w[0]=MR_REMAIN(w0->w[0],y->w[0]);
            mr_lzero(w0);
        }
        else if (MR_DIV(w0->w[x0-1],4)<y->w[x0-1])
        while (compare(w0,y)>=0)
        {  /* mr_small quotient - so do up to four subtracts instead */
            mr_psub(_MIPP_ w0,y,w0);
            d++;
        }
    }
    if (compare(w0,y)<0)
    {  /*  x less than y - so x becomes remainder */
        if (x!=z)  /* testing parameters */
        {
            copy(w0,x);
            if (x->len!=0) x->len|=sx;
        }
        if (y!=z)
        {
            zero(z);
            z->w[0]=d;
            if (d>0) z->len=(sz|1);
        }
        y->len|=sy;
        MR_OUT
        return;
    }

    if (y0==1)
    {  /* y is int - so use subdiv instead */
#ifdef MR_FP_ROUNDING
        r=mr_sdiv(_MIPP_ w0,y->w[0],mr_invert(y->w[0]),w0);
#else
        r=mr_sdiv(_MIPP_ w0,y->w[0],w0);
#endif
        if (y!=z)
        {
            copy(w0,z);
            z->len|=sz;
        }
        if (x!=z)
        {
            zero(x);
            x->w[0]=r;
            if (r>0) x->len=(sx|1);
        }
        y->len|=sy;
        MR_OUT
        return;
    }
    if (y!=z) zero(z);
    d=normalise(_MIPP_ y,y);
    check=mr_mip->check;
    mr_mip->check=OFF;

    if (mr_mip->base==0)
    {
#ifndef MR_NOFULLWIDTH
        if (d!=1) mr_pmul(_MIPP_ w0,d,w0);
        ldy=y->w[y0-1];
        sdy=y->w[y0-2];
        w0g=w0->w; yg=y->w;
        for (k=w00-1;k>=y0-1;k--)
        {  /* long division */

#if INLINE_ASM == 1
#ifdef MR_LMM
            ASM push ds
            ASM lds bx,DWORD PTR w0g
#else
            ASM mov bx,w0g
#endif
            ASM mov si,k
            ASM shl si,1
            ASM add bx,si
            ASM mov dx,[bx+2]
            ASM mov ax,[bx]
            ASM cmp dx,ldy
            ASM jne tcl8
            ASM mov di,0xffff
            ASM mov si,ax
            ASM add si,ldy
            ASM jc tcl12
            ASM jmp tcl10
          tcl8:
            ASM div WORD PTR ldy
            ASM mov di,ax
            ASM mov si,dx
          tcl10:
            ASM mov ax,sdy
            ASM mul di
            ASM cmp dx,si
            ASM jb tcl12
            ASM jne tcl11
            ASM cmp ax,[bx-2]
            ASM jbe tcl12
          tcl11:
            ASM dec di
            ASM add si,ldy
            ASM jnc tcl10
          tcl12:
            ASM mov attemp,di
#ifdef MR_LMM
            ASM pop ds
#endif
#endif
/* NOTE push and pop of esi/edi should not be necessary - Borland C bug *
 * These pushes are needed here even if register variables are disabled */
#if INLINE_ASM == 2
            ASM push esi
            ASM push edi
#ifdef MR_LMM
            ASM push ds
            ASM lds bx,DWORD PTR w0g
#else
            ASM mov bx,w0g
#endif
            ASM mov si,k
            ASM shl si,2
            ASM add bx,si
            ASM mov edx,[bx+4]
            ASM mov eax,[bx]
            ASM cmp edx,ldy
            ASM jne tcl8
            ASM mov edi,0xffffffff
            ASM mov esi,eax
            ASM add esi,ldy
            ASM jc tcl12
            ASM jmp tcl10
          tcl8:
            ASM div DWORD PTR ldy
            ASM mov edi,eax
            ASM mov esi,edx
          tcl10:
            ASM mov eax,sdy
            ASM mul edi
            ASM cmp edx,esi
            ASM jb tcl12
            ASM jne tcl11
            ASM cmp eax,[bx-4]
            ASM jbe tcl12
          tcl11:
            ASM dec edi
            ASM add esi,ldy
            ASM jnc tcl10
          tcl12:
            ASM mov attemp,edi
#ifdef MR_LMM