mrarth2.c

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

            ASM pop ds
#endif
            ASM pop edi
            ASM pop esi
#endif  
#if INLINE_ASM == 3
            ASM push esi
            ASM push edi
            ASM mov ebx,w0g
            ASM mov esi,k
            ASM shl esi,2
            ASM add ebx,esi
            ASM mov edx,[ebx+4]
            ASM mov eax,[ebx]
            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,[ebx-4]
            ASM jbe tcl12
          tcl11:
            ASM dec edi
            ASM add esi,ldy
            ASM jnc tcl10
          tcl12:
            ASM mov attemp,edi
            ASM pop edi
            ASM pop esi
#endif       
#if INLINE_ASM == 4
   ASM (
           "movl %1,%%ebx\n"
           "movl %2,%%esi\n"
           "shll $2,%%esi\n"
           "addl %%esi,%%ebx\n"
           "movl 4(%%ebx),%%edx\n"
           "movl (%%ebx),%%eax\n"
           "cmpl %3,%%edx\n"
           "jne tcl8\n"
           "movl $0xffffffff,%%edi\n"
           "movl %%eax,%%esi\n"
           "addl %3,%%esi\n"
           "jc tcl12\n"
           "jmp tcl10\n"
         "tcl8:\n"
           "divl %3\n"
           "movl %%eax,%%edi\n"
           "movl %%edx,%%esi\n"
         "tcl10:\n"
           "movl %4,%%eax\n"
           "mull %%edi\n"
           "cmpl %%esi,%%edx\n"
           "jb tcl12\n"
           "jne tcl11\n"
           "cmpl -4(%%ebx),%%eax\n"
           "jbe tcl12\n"
         "tcl11:\n"
           "decl %%edi\n"
           "addl %3,%%esi\n"
           "jnc tcl10\n"
         "tcl12:\n"
           "movl %%edi,%0\n"
        :"=m"(attemp)
        :"m"(w0g),"m"(k),"m"(ldy),"m"(sdy)
        :"eax","edi","esi","ebx","ecx","edx","memory"
       );
#endif
#ifndef INLINE_ASM
            carry=0;
            if (w0->w[k+1]==ldy) /* guess next quotient digit */
            {
                attemp=(mr_small)(-1);
                ra=ldy+w0->w[k];
                if (ra<ldy) carry=1;
            }
#ifdef MR_NOASM
            else
            {
                dble.h[MR_BOT]=w0->w[k];
                dble.h[MR_TOP]=w0->w[k+1];
                attemp=(mr_small)(dble.d/ldy);
                ra=(mr_small)(dble.d-(mr_large)attemp*ldy);
            }
#else
            else attemp=muldvm(w0->w[k+1],w0->w[k],ldy,&ra);
#endif
            while (carry==0)
            {
#ifdef MR_NOASM
                dble.d=(mr_large)attemp*sdy;
                r=dble.h[MR_BOT];
                tst=dble.h[MR_TOP];
#else
                tst=muldvd(sdy,attemp,(mr_small)0,&r);
#endif
                if (tst< ra || (tst==ra && r<=w0->w[k-1])) break;
                attemp--;  /* refine guess */
                ra+=ldy;
                if (ra<ldy) carry=1;
            }
#endif    
            m=k-y0+1;
            if (attemp>0)
            { /* do partial subtraction */
                borrow=0;
    /*  inline - substitutes for loop below */
#if INLINE_ASM == 1
                ASM cld
                ASM mov cx,y0
                ASM mov si,m
                ASM shl si,1
                ASM mov di,attemp
#ifdef MR_LMM
                ASM push ds
                ASM push es
                ASM les bx,DWORD PTR w0g
                ASM add bx,si
                ASM sub bx,2
                ASM lds si,DWORD PTR yg
#else
                ASM mov bx,w0g
                ASM add bx,si
                ASM sub bx,2
                ASM mov si,yg
#endif
                ASM push bp
                ASM xor bp,bp

             tcl3:
                ASM lodsw
                ASM mul di
                ASM add ax,bp
                ASM adc dx,0
                ASM inc bx
                ASM inc bx
#ifdef MR_LMM
                ASM sub es:[bx],ax
#else
                ASM sub [bx],ax
#endif              
                ASM adc dx,0
                ASM mov bp,dx
                ASM loop tcl3

                ASM mov ax,bp
                ASM pop bp
#ifdef MR_LMM
                ASM pop es
                ASM pop ds
#endif
                ASM mov borrow,ax
#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
                ASM cld
                ASM mov cx,y0
                ASM mov si,m
                ASM shl si,2
                ASM mov edi,attemp
#ifdef MR_LMM
                ASM push ds
                ASM push es
                ASM les bx,DWORD PTR w0g
                ASM add bx,si
                ASM sub bx,4
                ASM lds si,DWORD PTR yg
#else
                ASM mov bx,w0g
                ASM add bx,si
                ASM sub bx,4
                ASM mov si,yg
#endif
                ASM push ebp
                ASM xor ebp,ebp

             tcl3:
                ASM lodsd
                ASM mul edi
                ASM add eax,ebp
                ASM adc edx,0
                ASM add bx,4
#ifdef MR_LMM
                ASM sub es:[bx],eax
#else
                ASM sub [bx],eax
#endif
                ASM adc edx,0
                ASM mov ebp,edx
                ASM loop tcl3

                ASM mov eax,ebp
                ASM pop ebp
#ifdef MR_LMM
                ASM pop es
                ASM pop ds
#endif
                ASM mov borrow,eax
                ASM pop edi
                ASM pop esi
#endif
#if INLINE_ASM == 3
                ASM push esi
                ASM push edi
                ASM mov ecx,y0
                ASM mov esi,m
                ASM shl esi,2
                ASM mov edi,attemp
                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

             tcl3:
                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 sub ebp,eax
                ASM adc edx,0
                ASM mov [esi+ebx],ebp
                ASM dec ecx
                ASM mov ebp,edx
                ASM jnz tcl3

                ASM mov eax,ebp
                ASM pop ebp
                ASM mov borrow,eax
                ASM pop edi
                ASM pop esi
#endif
#if INLINE_ASM == 4
   ASM (
           "movl %1,%%ecx\n"
           "movl %2,%%esi\n"
           "shll $2,%%esi\n"
           "movl %3,%%edi\n"
           "movl %4,%%ebx\n"
           "addl %%esi,%%ebx\n"
           "movl %5,%%esi\n"
           "subl %%esi,%%ebx\n"
           "subl $4,%%ebx\n"
           "pushl %%ebp\n"
           "xorl %%ebp,%%ebp\n"
         "tcl3:\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"
           "subl %%eax,%%ebp\n"
           "adcl $0,%%edx\n"
           "movl %%ebp,(%%esi,%%ebx)\n"
           "decl %%ecx\n"
           "movl %%edx,%%ebp\n"
           "jnz tcl3\n"
    
           "movl %%ebp,%%eax\n"
           "popl %%ebp\n"
           "movl %%eax,%0\n"
 
        :"=m"(borrow)
        :"m"(y0),"m"(m),"m"(attemp),"m"(w0g),"m"(yg)
        :"eax","edi","esi","ebx","ecx","edx","memory"
       );
#endif
#ifndef INLINE_ASM
                for (i=0;i<y0;i++)
                {
#ifdef MR_NOASM
                    dble.d=(mr_large)attemp*y->w[i]+borrow;
                    dig=dble.h[MR_BOT];
                    borrow=dble.h[MR_TOP];
#else
                  borrow=muldvd(attemp,y->w[i],borrow,&dig);
#endif
                  if (w0->w[m+i]<dig) borrow++;
                  w0->w[m+i]-=dig;
                }
#endif

                if (w0->w[k+1]<borrow)
                {  /* whoops! - over did it */
                    w0->w[k+1]=0;
                    carry=0;
                    for (i=0;i<y0;i++)
                    {  /* compensate for error ... */
                        psum=w0->w[m+i]+y->w[i]+carry;
                        if (psum>y->w[i]) carry=0;
                        if (psum<y->w[i]) carry=1;
                        w0->w[m+i]=psum;
                    }
                    attemp--;  /* ... and adjust guess */
                }
                else w0->w[k+1]-=borrow;
            }
            if (k==w00-1 && attemp==0) w00--;
            else if (y!=z) z->w[m]=attemp;
        }
#endif
    }
    else
    {   /* have to do it the hard way */
        if (d!=1) mr_pmul(_MIPP_ w0,d,w0);
        ldy=y->w[y0-1];
        sdy=y->w[y0-2];

        for (k=w00-1;k>=y0-1;k--)
        {  /* long division */


            if (w0->w[k+1]==ldy) /* guess next quotient digit */
            {
                attemp=mr_mip->base-1;
                ra=ldy+w0->w[k];
            }
#ifdef MR_NOASM
            else 
            {
                dbled=(mr_large)w0->w[k+1]*mr_mip->base+w0->w[k];
                attemp=(mr_small)MR_LROUND(dbled/ldy);
                ra=(mr_small)(dbled-(mr_large)attemp*ldy);
            }
#else
            else attemp=muldiv(w0->w[k+1],mr_mip->base,w0->w[k],ldy,&ra);
#endif
            while (ra<mr_mip->base)
            {
#ifdef MR_NOASM
                dbled=(mr_large)sdy*attemp;
#ifdef MR_FP_ROUNDING
                tst=(mr_small)MR_LROUND(dbled*mr_mip->inverse_base);
#else
#ifndef MR_FP
                if (mr_mip->base==mr_mip->base2)
                    tst=(mr_small)(dbled>>mr_mip->lg2b);
                else 
#endif  
                    tst=(mr_small)MR_LROUND(dbled/mr_mip->base);
#endif
                r=(mr_small)(dbled-(mr_large)tst*mr_mip->base);
#else
#ifdef MR_FP_ROUNDING
                tst=imuldiv(sdy,attemp,(mr_small)0,mr_mip->base,mr_mip->inverse_base,&r); 
#else
                tst=muldiv(sdy,attemp,(mr_small)0,mr_mip->base,&r); 
#endif
#endif
                if (tst< ra || (tst==ra && r<=w0->w[k-1])) break;
                attemp--;  /* refine guess */
                ra+=ldy;
            }    
            m=k-y0+1;
            if (attemp>0)
            { /* do partial subtraction */
                borrow=0;
                for (i=0;i<y0;i++)
                {
#ifdef MR_NOASM
                  dbled=(mr_large)attemp*y->w[i]+borrow;
#ifdef MR_FP_ROUNDING
                  borrow=(mr_small)MR_LROUND(dbled*mr_mip->inverse_base);
#else
#ifndef MR_FP
                  if (mr_mip->base==mr_mip->base2)
                      borrow=(mr_small)(dbled>>mr_mip->lg2b);
                  else 
#endif  
                      borrow=(mr_small)MR_LROUND(dbled/mr_mip->base);
#endif
                  dig=(mr_small)(dbled-(mr_large)borrow*mr_mip->base);
#else
#ifdef MR_FP_ROUNDING
                  borrow=imuldiv(attemp,y->w[i],borrow,mr_mip->base,mr_mip->inverse_base,&dig);
#else
                  borrow=muldiv(attemp,y->w[i],borrow,mr_mip->base,&dig);
#endif
#endif
                  if (w0->w[m+i]<dig)
                  { /* set borrow */
                      borrow++;
                      w0->w[m+i]+=(mr_mip->base-dig);
                  }
                  else w0->w[m+i]-=dig;
                }
                if (w0->w[k+1]<borrow)
                {  /* whoops! - over did it */
                    w0->w[k+1]=0;
                    carry=0;
                    for (i=0;i<y0;i++)
                    {  /* compensate for error ... */
                        psum=w0->w[m+i]+y->w[i]+carry;
                        carry=0;
                        if (psum>=mr_mip->base)
                        {
                            carry=1;
                            psum-=mr_mip->base;
                        }
                        w0->w[m+i]=psum;
                    }
                    attemp--;  /* ... and adjust guess */
                }
                else
                    w0->w[k+1]-=borrow;
            }
            if (k==w00-1 && attemp==0) w00--;
            else if (y!=z) z->w[m]=attemp;
        }
    }
    if (y!=z) z->len=((w00-y0+1)|sz); /* set sign and length of result */

    w0->len=y0;

    mr_lzero(y);
    mr_lzero(z);

    if (x!=z)
    {
        mr_lzero(w0);
#ifdef MR_FP_ROUNDING
        if (d!=1) mr_sdiv(_MIPP_ w0,d,mr_invert(d),x);
#else
        if (d!=1) mr_sdiv(_MIPP_ w0,d,x);
#endif
        else copy(w0,x);
        if (x->len!=0) x->len|=sx;
    }
#ifdef MR_FP_ROUNDING
    if (d!=1) mr_sdiv(_MIPP_ y,d,mr_invert(d),y);
#else
    if (d!=1) mr_sdiv(_MIPP_ y,d,y);
#endif
    y->len|=sy;
    mr_mip->check=check;

    MR_OUT
}

BOOL divisible(_MIPD_ big x,big y)
{ /* returns y|x, that is TRUE if y divides x exactly */
#ifdef MR_OS_THREADS
    miracl *mr_mip=get_mip();
#endif
    if (mr_mip->ERNUM) return FALSE;

    MR_IN(87)

    copy (x,mr_mip->w0);
    divide(_MIPP_ mr_mip->w0,y,y);

    MR_OUT
    if (size(mr_mip->w0)==0) return TRUE;
    else                    return FALSE;
}     

void mad(_MIPD_ big x,big y,big z,big w,big q,big r)
{ /* Multiply, Add and Divide; q=(x*y+z)/w remainder r   *
   * returns remainder only if w=q, quotient only if q=r *
   * add done only if x, y and z are distinct.           */
#ifdef MR_OS_THREADS
    miracl *mr_mip=get_mip();
#endif
    if (mr_mip->ERNUM) return;

    MR_IN(24)

    mr_mip->check=OFF;           /* turn off some error checks */
    if (w==r)
    {
        mr_berror(_MIPP_ MR_ERR_BAD_PARAMETERS);
        MR_OUT
        return;
    }
    multiply(_MIPP_ x,y,mr_mip->w0);
    if (x!=z && y!=z) add(_MIPP_ mr_mip->w0,z,mr_mip->w0);

    divide(_MIPP_ mr_mip->w0,w,q);
    if (q!=r) copy(mr_mip->w0,r);
    mr_mip->check=ON;
    MR_OUT
}