mrmuldv.any

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

        swap   D1            ;get remainder
        move.l rp(A6),A0     ;get address for remainder
        move   D1,(A0)       ;store remainder
    }
}

unsigned int muldvd(a,b,c,rp)
unsigned int a,b,c,*rp;
{
    asm
    {
    ;
    ; Version of muldvd for Apple MAC
    ;
        move   a(A6),D1      ;get a
        mulu   b(a6),D1      ;multiply by b
        clr.l  D0
        move   c(A6),D0      ;get c
        add.l  D0,D1         ;add in c
        move.l D1,D0         
        swap   D0            ;return quotient in D0
        move.l rp(A6),A0     ;get address for remainder
        move   D1,(A0)       ;store remainder
    }
}


**********************************************************************


#
# 68020+ versions for Next, and for new 32-bit Macs
# Parameters come off the stack
#

.globl _muldiv,_muldvm,_muldvd

_muldiv:
    movel sp@(4),d0
    mulul sp@(8),d1:d0
    addl  sp@(12),d0
    negxl d1          # tricky stuff!
    negl  d1
    divul sp@(16),d1:d0
    movel sp@(20),a0
    movel d1,a0@
    rts

_muldvm:
    movel sp@(4),d1
    movel sp@(8),d0
    divul sp@(12),d1:d0
    movel sp@(16),a0
    movel d1,a0@
    rts

_muldvd:
    movel sp@(4),d1
    mulul sp@(8),d0:d1
    addl  sp@(12),d1
    negxl d0
    negl  d0
    movel sp@(16),a0
    movel d1,a0@
    rts



*************************************************************************


unsigned int muldiv(a,b,c,m,rp)
unsigned int a,b,c,m,*rp;
{
    asm
    {
    ;
    ; 32016 processor version for BBC Master Scientific
    ; with 32-bit int, by Dudley Long, Rutherford-Appleton Labs.
    ; No muldvm() or muldvd()
    ;
        movd   a,0           ;move a to R0
        meid   b,0           ;multiply  by b, result extended
        addd   c,0           ;add c to extended number in R0 & R1
        addcd  #0,1
        deid   m,0           ;divide by m
        movd   0,0(rp)       ;remainder to *rp
        movd   1,0           ;quotient returned in R0
    }
}


*******************************************************************


; 
; MOTE! This code is obsolete. Newer ARMs support a 32x32 UMULL instruction
; The ARM compiler supports a long long type, so a C only version may be 
; faster
;
; Acorn ARM Risc version (32-bit) for Archimedes micro
; Wingpass Macro Assembler
; Use with mirdef.h32
;
.INCLUDE "A.REGNAMES"

.AREA C$$code, .CODE, .READONLY

muldiv::
         MOV     ip, sp             ;standard linkage
         STMFD   sp!, {v1-v4}

         CMPS    a2,#0x80000000     ;check for b=MAXBASE
         MOVEQ   v3,a1,LSL #31      ;this idea is quicker because
         MOVEQ   v4,a1,LSR #1       ;of ARM barrel shifting capability
         BEQ     addin
         MOV     v1,a1,LSR #16      ;do it the hard way
         MOV     v2,a2,LSR #16
         BIC     a1,a1,v1,LSL #16
         BIC     a2,a2,v2,LSL #16
         MUL     v3,a1,a2           ;form partial products of a*b
         MUL     v4,v1,v2
         SUB     v1,v1,a1
         SUB     v2,a2,v2
         MLA     v1,v2,v1,v3        ;look - only 3 MULs!
         ADD     v1,v1,v4
         ADDS    v3,v3,v1,LSL #16
         ADC     v4,v4,v1,LSR #16
addin:
         ADDS    v3,v3,a3           ;now add in c
         ADCCS   v4,v4,#0

         CMPS    a4,#0x80000000     ;check for m=MAXBASE
         MOVEQ   a1,v3,LSR #31
         ADDEQ   a1,a1,v4,LSL #1
         BICEQ   v4,v3,#0x80000000
         BEQ     leave
         MOV     a1,#0              ;do long division by m

divlp:

.REPEAT 32                          ;2xfaster than a loop!
         MOVS    v3,v3,ASL #1       ;get next bit into carry
         ADC     v4,v4,v4           ;accumulate remainder
         CMPS    v4,a4
         SUBCS   v4,v4,a4
         ADC     a1,a1,a1           ;accumulate quotient
.ENDREPEAT

leave:
         LDR     v3,[ip]
         STR     v4,[v3]       ;store remainder
         LDMFD   sp!, {v1-v4}
         MOVS pc,lr

muldvm::
         STMFD   sp!, {v1-v2}

         MOV     v2,a1              ;'multiply' by 2^32
         MOV     v1,a2              ;add in c

         MOV     a1,#0              ;do long division by m

.REPEAT 32                          ;2xfaster than a loop!
         MOVS    v1,v1,ASL #1       ;get next bit into carry
         ADCS    v2,v2,v2           ;accumulate remainder
         CMPCCS  v2,a3
         SUBCS   v2,v2,a3
         ADC     a1,a1,a1           ;accumulate quotient
.ENDREPEAT

         STR     v2,[a4]       ;store remainder
         LDMFD   sp!, {v1-v2}
         MOVS    pc,lr


muldvd::
         STMFD   sp!, {v1-v2}

         MOV     ip,a1,LSR #16      ;do it the hard way
         MOV     v2,a2,LSR #16
         BIC     a1,a1,ip,LSL #16
         BIC     a2,a2,v2,LSL #16
         MUL     v1,a1,a2           ;form partial products of a*b
         MUL     a2,ip,a2
         MUL     a1,v2,a1
         MUL     v2,ip,v2
         ADDS    a1,a2,a1
         ADDCS   v2,v2,#0x10000
         ADDS    v1,v1,a1,LSL #16
         ADC     v2,v2,a1,LSR #16
         
         ADDS    v1,v1,a3           ;now add in c
         ADCCS   v2,v2,#0
         MOV     a1,v2              ;get quotient

         STR     v1,[a4]            ;store remainder
         LDMFD   sp!, {v1-v2}
         MOVS    pc,lr



********************************************************************


;
;  Version for Pyramid 90x and 98x computers
;  from Rod Worley, Monash University, Victoria, Australia
;
;  No muldvm() or muldvd()
;
        .text 0
        .globl _muldiv
_muldiv:
        movw    pr0,pr8             ;save a in reg 8
        movw    $0x0,pr0            ;zero reg0 so long reg 0,1 is b
        emul    pr8,pr0             ;extended multiply by a
        addw    pr2,pr1             ;add c to extended result
        addwc   $0x0,pr0         
        ediv    pr3,pr0             ;extended div by m
        movw    pr1,(pr4)           ;store remainder
        ret                         ;return qotient in pr0


************************************************************************


/* This is the transputer version, by A.H. Pepperdine   */
/* Assumes that the result will fit into a 32-bit word  */
/* The error flag will be set if                        */
/*           (a*b+c)/m >= 2**32                         */
/* ie. equivalently, if                                 */
/*           ( (a*b+c) >> 32) >= m                      */

unsigned int muldiv(unsigned int a, unsigned int b, unsigned int c,
                    unsigned int m, unsigned int * rp)
{
    unsigned int q;
    __asm
    {
      ldabc a, b, c;
      lmul  ;
      ld    m;
      ldiv  ;
      stab  q, *rp;
    }
    return q;
}

/* The base is 2**32, ie a full 32-bit unsigned integer */
/* The error flag will be set if the result will not fit*/
/* into a word, ie.                                     */
/* for muldvm that is if (a >= m)                       */
/* and for muldvd it cannot happen                      */

unsigned int muldvm(unsigned int a, unsigned int c,
                    unsigned int m, unsigned int * rp)
{
    unsigned int q;
    __asm
    {
      ldabc m, c, a;
      ldiv ;
      stab q, *rp;
    }
    return q;
}

unsigned int muldvd(unsigned int a, unsigned int b, unsigned int c,
                    unsigned int * rp)
{
    unsigned int q;
    __asm
    {
      ldabc a, b, c;
      lmul ;
      stab *rp, q;
    }
    return q;
}


*********************************************************************


/*  Now ... just to confuse you even more ....

    Blakley/Sloan 'portable' method for Modular multiplication IEEE Trans 
    Computers C-34 March 1985 pp 290-292 eliminates need for double length 
    product - but will be slow. Might suit some RISC computers with no 
    multiply/divide instructions. To speed up try completely unravelling for() 
    loops.       
    
    This method should only be used if the mr_utype data type is twice the size 
    of a "mr_hltype" data-type. This must be defined below.

    Note: DON't define MR_NOASM in mirdef.h if using this method.

    */ 

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

mr_small muldiv(a,b,c,m,rp)
mr_small a,b,c,m;
mr_small *rp;
{
    int i;
    mr_small d,q=0,r=0;
    d=m-a;
    for (i=MIRACL/4;i>0;i--)
    { /* do it bit by bit */
        r<<=1;
        if ((mr_utype)c<0) r++;
        c<<=1;
        q<<=1;
        if ((mr_utype)b<0)
        {
            if (r>=m) { r-=d; q++; }
            else        r+=a;
        }
        if (r>=m) { r-=m; q++; }
        b<<=1;
        r<<=1;
        if ((mr_utype)c<0) r++;
        c<<=1;
        q<<=1;
        if ((mr_utype)b<0)
        {
            if (r>=m) { r-=d; q++; }
            else        r+=a;
        }
        if (r>=m) { r-=m; q++; }
        b<<=1;
        r<<=1;
        if ((mr_utype)c<0) r++;
        c<<=1;
        q<<=1;
        if ((mr_utype)b<0)
        {
            if (r>=m) { r-=d; q++; }
            else        r+=a;
        }
        if (r>=m) { r-=m; q++; }
        b<<=1;
        r<<=1;
        if ((mr_utype)c<0) r++;
        c<<=1;
        q<<=1;
        if ((mr_utype)b<0)
        {
            if (r>=m) { r-=d; q++; }
            else        r+=a;
        }
        if (r>=m) { r-=m; q++; }
        b<<=1;
    }
    *rp=r;
    return q;
}

mr_small muldvm(a,c,m,rp)
mr_small a,c,m;
mr_small *rp;
{ /* modified Blakely-Sloan */
    register int i,carry;
    register mr_small q=0,r=0;
    r=a;
    for (i=MIRACL/4;i>0;i--)
    { /* do it bit by bit */
        carry=0;
        if ((mr_utype)r<0) carry=1;
        r<<=1;
        if ((mr_utype)c<0) r++;
        c<<=1;
        q<<=1;
        if (carry || r>=m) { r-=m; q++; }
        carry=0;
        if ((mr_utype)r<0) carry=1;
        r<<=1;
        if ((mr_utype)c<0) r++;
        c<<=1;
        q<<=1;
        if (carry || r>=m) { r-=m; q++; }
        carry=0;
        if ((mr_utype)r<0) carry=1;
        r<<=1;
        if ((mr_utype)c<0) r++;
        c<<=1;
        q<<=1;
        if (carry || r>=m) { r-=m; q++; }
        carry=0;
        if ((mr_utype)r<0) carry=1;
        r<<=1;
        if ((mr_utype)c<0) r++;
        c<<=1;
        q<<=1;
        if (carry || r>=m) { r-=m; q++; }
    }
    *rp=r;
    return q;
}

/* define mr_hltype as that C type that is half the size in bits of the 
   underlying type (mr_utype in mirdef.h). Perhaps short if mr_utype is long? */

#define mr_hltype short

mr_small muldvd(a,b,c,rp)
mr_small a,b,c;
mr_small *rp;
{ /* multiply by parts */
    mr_small middle,middle2;
    mr_small q,r;
    unsigned mr_hltype am,al,bm,bl;
    int hshift=(MIRACL>>1);
    am=(unsigned mr_hltype)(a>>hshift);
    al=(unsigned mr_hltype)a;
    bm=(unsigned mr_hltype)(b>>hshift);
    bl=(unsigned mr_hltype)b;
/* form partial products */
    r= (mr_small)al*bl;  
    q= (mr_small)am*bm;
    middle=(mr_small)al*bm;
    middle2=(mr_small)bl*am;
    middle+=middle2;                        /* combine them - carefully */
    if (middle<middle2) q+=((mr_small)1<<hshift);
    r+=(middle << hshift);
    if ((r>>hshift)<(unsigned mr_hltype)middle) q++;
    q+=(middle>>hshift);
    r+=c;
    if (r<c) q++;
    *rp=r;
    return q;
}   

void muldvd2(a,b,c,rp)
mr_small a,b;
mr_small *c,*rp;
{ /* multiply by parts */
    mr_small middle,middle2;
    mr_small q,r;
    unsigned mr_hltype am,al,bm,bl;
    int hshift=(MIRACL>>1);
    am=(unsigned mr_hltype)(a>>hshift);
    al=(unsigned mr_hltype)a;
    bm=(unsigned mr_hltype)(b>>hshift);
    bl=(unsigned mr_hltype)b;
/* form partial products */
    r= (mr_small)al*bl;