home *** CD-ROM | disk | FTP | other *** search
/ The C Users' Group Library 1994 August / wc-cdrom-cusersgrouplibrary-1994-08.iso / vol_200 / 247_01 / bnmuldv.any < prev    next >
Text File  |  1989-04-19  |  12KB  |  378 lines
  1. /*
  2.  *  MIRACL routine muldiv
  3.  *  bnmuldv.c
  4.  * 
  5.  *  calculates (a*b+c)/m and (a*b+c)%m as quickly as
  6.  *  possible. Should ideally be written in assembly
  7.  *  language of target machine for speed and/or to allow
  8.  *  largest possible MAXBASE. See mirdef.h. The problem is
  9.  *  to avoid overflow in the calculation of the intermediate
  10.  *  product a*b+c. Note that this routine needs only work
  11.  *  for UNSIGNED parameters 
  12.  *
  13.  *  First are shown two straight-forward portable approaches, 
  14.  *  which needs a long data type capable of holding
  15.  *  the product of two int types (i.e. twice as long).
  16.  *  The second method uses the 'ldiv' function which
  17.  *  is defined in the new ANSI standard for C. If you are
  18.  *  using a C compiler which complies with this standard, this
  19.  *  should be faster, as it calculates quotient and remainder
  20.  *  simultaneously (somebody up there heard me!).
  21.  *
  22.  *  This is followed by various other assembly language
  23.  *  implementations for popular processors, computers and 
  24.  *  compilers.
  25.  *
  26.  *  Since parameter passing and returning is time-consuming,
  27.  *  this routine should be generated 'inline', if compiler
  28.  *  allows it.
  29.  *
  30.  
  31. int muldiv(a,b,c,m,rp)
  32. int a,b,c,m,*rp;
  33. {
  34.     long p;
  35.     int q;
  36.     p=(long)a*b+c;
  37.     q=p/m;
  38.     *rp=p-(long)q*m;
  39.     return q;
  40. }
  41.  
  42. #include <stdlib.h>
  43.  
  44. int muldiv(a,b,c,m,rp)
  45. int a,b,c,m,*rp;
  46. {
  47.     ldiv_t r;
  48.     r=ldiv((long)a*b+c,(long)m);
  49.     *rp=(int)r.rem;
  50.     return (int)r.quot;
  51. }
  52.  
  53. ;
  54. ;   VAX11 version for Dec C compiler
  55. ;   with 32 bit int using 64-bit quadword
  56. ;   for the intermediate product. 
  57. ;
  58.     .entry muldiv,0
  59.     subl   #4,sp
  60.     emul   4(ap),8(ap),12(ap),r0   ;a*b+c
  61.     ediv   16(ap),r0,r0,@20(ap)    ;quo. in r0, rem. in *rp
  62.     ret
  63.     .end
  64.  
  65. ;
  66. ;  IBM-PC versions - small memory model only
  67. ;  Easily modified for other memory models
  68. ;
  69. ;  Microsoft C compiler V4.0+
  70. ;  Written for MS macro-assembler
  71. ;
  72.         ASSUME CS:_TEXT
  73. _TEXT   SEGMENT BYTE PUBLIC 'CODE'
  74.  
  75.         PUBLIC _muldiv
  76. _muldiv PROC NEAR
  77.         push bp                 ;standard C linkage
  78.         mov  bp,sp
  79.  
  80.         mov  ax,[bp+4]          ;get a
  81.         mul  WORD PTR [bp+6]    ;multiply by b
  82.         add  ax,[bp+8]          ;add c to low word
  83.         adc  dx,0h              ;add carry to high word
  84.         div  WORD PTR [bp+10]   ;divide by m
  85.         mov  bx,[bp+12]         ;get address for remainder
  86.         mov  [bx],dx            ;store remainder
  87.  
  88.         pop  bp                 ;standard C return
  89.         ret                     ;quotient in ax
  90.  
  91. _muldiv endP
  92.  
  93. _TEXT   ENDS
  94. END
  95.  
  96. /*
  97.  *  Turbo C compiler V1.0. Uses inline assembly feature
  98.  *  Generates code identical to above Microsoft version
  99.  */
  100. int muldiv(a,b,c,m,rp)
  101. int a,b,c,m,*rp;
  102. {
  103.     asm mov  ax,a             ;/* get a                     */
  104.     asm mul  WORD PTR b       ;/* multiply by b             */
  105.     asm add  ax,c             ;/* add c to low word         */
  106.     asm adc  dx,0h            ;/* add carry to high word    */
  107.     asm div  WORD PTR m       ;/* divide by m               */
  108.     asm mov  bx,rp            ;/* get address for remainder */
  109.     asm mov  [bx],dx          ;/* store remainder           */
  110.     return (_AX);
  111. }
  112. /*    Replace last two asm lines when using large data memory models */
  113. /*    asm les  bx, DWORD PTR rp          ; get address for remainder */
  114. /*    asm mov  WORD PTR es:[bx],dx       ; store remainder           */
  115.  
  116. ;
  117. ;  Zorland C compiler Version 1.1 (Lattice C compatible)
  118. ;  Note that later versions of Zortech C (and C++) use the MS version
  119. ;  Written for MS macro-assembler
  120. ; 
  121. PGROUP  GROUP PROG
  122. PROG    SEGMENT BYTE PUBLIC 'PROG'
  123.         ASSUME CS:PGROUP
  124.         PUBLIC muldiv
  125. muldiv  PROC NEAR
  126.         push bp                 ;standard C linkage
  127.         mov  bp,sp
  128.  
  129.         mov  ax,[bp+4]          ;get a
  130.         mul  WORD PTR [bp+6]    ;multiply by b
  131.         add  ax,[bp+8]          ;add c to low word
  132.         adc  dx,0h              ;add carry to high word
  133.         div  WORD PTR [bp+10]   ;divide by m
  134.         mov  bx,[bp+12]         ;get address for remainder
  135.         mov  [bx],dx            ;store remainder
  136.  
  137.         pop  bp                 ;standard C return
  138.         ret                     ;quotient in ax
  139.  
  140. muldiv  endP
  141.  
  142. PROG    ENDS
  143. END
  144.  
  145. ;
  146. ;   Aztec C compiler V3.4
  147. ;
  148. CODESEG SEGMENT BYTE PUBLIC 'CODESEG'
  149.         ASSUME CS:CODESEG
  150.         PUBLIC muldiv_
  151. muldiv_ PROC NEAR
  152.         push bp                 ;standard C linkage
  153.         mov  bp,sp
  154.  
  155.         mov  ax,[bp+4]          ;get a
  156.         mul  WORD PTR [bp+6]    ;multiply by b
  157.         add  ax,[bp+8]          ;add c to low word
  158.         adc  dx,0h              ;add carry to high word
  159.         div  WORD PTR [bp+10]   ;divide by m
  160.         mov  bx,[bp+12]         ;get address for remainder
  161.         mov  [bx],dx            ;store remainder
  162.  
  163.         pop  bp                 ;standard C return
  164.         ret                     ;quotient in ax
  165.  
  166. muldiv_ endP
  167.  
  168. CODESEG ENDS
  169. END
  170.  
  171. /
  172. /  Mark Williams C compiler V2.0
  173. /
  174. .globl muldiv_
  175.  
  176. muldiv_:
  177.         push bp                 /standard C linkage
  178.         mov  bp,sp
  179.  
  180.         mov  ax,4(bp)           /get a
  181.         mul  6(bp)              /multiply by b
  182.         add  ax,8(bp)           /add c to low word
  183.         adc  dx,$0x00           /add carry to high word
  184.         div  10(bp)             /divide by m
  185.         mov  bx,12(bp)          /get address for remainder
  186.         mov  (bx),dx            /store remainder
  187.  
  188.         pop  bp                 /standard C return
  189.         ret                     /quotient in ax
  190.  
  191. ;
  192. ;  IBM-PC-8087 for Microsoft C compiler V4.0+
  193. ;  with 'small' defined as 'long' (32 bit). See mirdef.h
  194. ;  This allows IBM-PC to look like 32-bit computer
  195. ;  (which it isn't). To make use of this option:
  196. ;
  197. ;  (1) Must have 8087 Maths Co-processor (for speed and to hold 64-bit
  198. ;      intermediate product).
  199. ;
  200. ;  (2) Must use 'ANSI' enhanced type C compiler, e.g. Microsoft V3.0+
  201. ;      and must use header 'miracl.h' which declares function
  202. ;      parameter types. Changes will be necessary in 'mirdef.h'
  203. ;
  204. ;      Note: many compilation warnings may be generated - ignore them.  
  205. ;
  206. ;  Note: This is NOT, in most cases, faster, but it does allow
  207. ;        very high precision calculations, e.g. 1000!
  208. ;
  209.         ASSUME CS:_TEXT
  210. _TEXT   SEGMENT BYTE PUBLIC 'CODE'
  211.  
  212.         PUBLIC _MULDIV
  213. _muldiv PROC NEAR
  214.         push si                 ;standard C linkage
  215.         push bp          
  216.         mov  bp,sp
  217.  
  218.         finit                   ;initialise 8087
  219.         fild  DWORD PTR [bp+6]  ;get a
  220.         fimul DWORD PTR [bp+0ah];multiply by b
  221.         fiadd DWORD PTR [bp+0eh];add c
  222.         fild  DWORD PTR [bp+12h];get m
  223.         fld st(1)               ;duplicate a*b+c on stack
  224.         fprem                   ;get remainder
  225.         fist  DWORD PTR [bp+0ah];store remainder in b
  226.         fsubr st,st(2)          ;subtract rem from total
  227.         fdiv st,st(1)           ;divide by m
  228.         fist  DWORD PTR [bp+6]  ;store quotient in a
  229.         wait
  230.  
  231.         mov  si,[bp+22]         ;get address for remainder
  232.         mov  ax,[bp+10]
  233.         mov  dx,[bp+12]         ;get remainder
  234.         mov  [si],ax
  235.         mov  [si+2],dx          ;store remainder
  236.         mov  ax,[bp+6]
  237.         mov  dx,[bp+8]          ;get quotient in dx:ax
  238.  
  239.         pop  bp                 ;standard C return
  240.         pop  si
  241.         ret
  242.  
  243. _muldiv endP
  244.  
  245. _TEXT   ENDS
  246. END
  247.  
  248. ;
  249. ;  IBM-80386 version - for Microsoft C V5.0+
  250. ;  Written for MS macro-assembler V5.0+ by Andrej Sauer 
  251. ;  Same comments apply as above (except for 8087 requirement)   
  252. ;
  253.         .386
  254.         ASSUME CS:_TEXT
  255. _TEXT   SEGMENT USE16 PUBLIC 'CODE'
  256.  
  257.         PUBLIC _MULDIV
  258. _muldiv PROC  NEAR
  259.         push  bp                     ;standard C linkage
  260.         mov   ebp,esp
  261.  
  262.         mov   eax,[ebp+4]            ;get a
  263.         mul   DWORD PTR [ebp+8]      ;multiply by b
  264.         add   eax,DWORD PTR [ebp+12] ;add c to low word
  265.         adc   edx,0h                 ;add carry to high word
  266.         div   DWORD PTR [ebp+16]     ;divide by m
  267.         movzx ebx,WORD PTR [ebp+20]  ;get address for remainder, zero high
  268.         mov   [ebx],edx              ;store remainder
  269.         shld  edx,eax,16             ;shift higher half of quotient
  270.                                      ;into lower half of edx
  271.  
  272.         pop   bp                     ;standard C return
  273.         ret                          ;quotient: high bits in dx, lows in ax
  274.  
  275. _muldiv endP
  276.  
  277. _TEXT   ENDS
  278. END
  279.  
  280. int muldiv(a,b,c,m,rp)
  281. int a,b,c,m,*rp;
  282. {
  283.     asm
  284.     {
  285.     ;
  286.     ;  MACintosh version for Megamax C compiler
  287.     ;  with 16-bit int, 68000 processor
  288.     ;
  289.         move   a(A6),D1      ;get a
  290.         mulu   b(A6),D1      ;multiply by b
  291.         move   c(A6),D0      ;get c
  292.         ext.l  D0            ;extend to longword
  293.         add.l  D0,D1         ;D1 contains a*b+c
  294.         divu   m(A6),D1      ;divide by m
  295.         move   D1,D0         ;return with quotient in D0
  296.         swap   D1            ;get remainder
  297.         move.l rp(A6),A0     ;get address for remainder
  298.         move   D1,(A0)       ;store remainder
  299.     }
  300. }
  301.  
  302. int muldiv(a,b,c,m,rp)
  303. int a,b,c,m,*rp;
  304. {
  305.     asm
  306.     {
  307.     ;
  308.     ; 32016 processor version for BBC Master Scientific
  309.     ; with 32-bit int, by Dudley Long, Rutherford-Appleton Labs.
  310.     ;
  311.         movd   a,0           ;move a to R0
  312.         meid   b,0           ;multiply  by b, result extended
  313.         addd   c,0           ;add c to extended number in R0 & R1
  314.         addcd  #0,1
  315.         deid   m,0           ;divide by m
  316.         movd   0,0(rp)       ;remainder to *rp
  317.         movd   1,0           ;quotient returned in R0
  318.     }
  319. }
  320.  
  321. ;
  322. ; Acorn ARM Risc version (32-bit) for Archimedes micro
  323. ; unsigned only
  324. ; Wingpass Macro Assembler
  325. ;
  326. .INCLUDE "A.REGNAMES"
  327.  
  328. .AREA C$$code, .CODE, .READONLY
  329.  
  330. muldiv::
  331.          MOV     ip, sp             ;standard linkage
  332.          STMFD   sp!, {v1-v4}
  333.  
  334.          CMPS    a2,#0x40000000     ;check for b=MAXBASE
  335.          MOVEQ   v3,a1,LSL #30      ;this idea is quicker because
  336.          MOVEQ   v4,a1,LSR #2       ;of ARM barrel shifting capability
  337.          BEQ     addin
  338.          MOV     v1,a1,LSR #16      ;do it the hard way
  339.          MOV     v2,a2,LSR #16
  340.          BIC     a1,a1,v1,LSL #16
  341.          BIC     a2,a2,v2,LSL #16
  342.          MUL     v3,a1,a2           ;form partial products of a*b
  343.          MUL     v4,v1,v2
  344.          SUB     v1,v1,a1
  345.          SUB     v2,a2,v2
  346.          MLA     v1,v2,v1,v3        ;look - only 3 MULs!
  347.          ADD     v1,v1,v4
  348.          ADDS    v3,v3,v1,LSL #16
  349.          ADC     v4,v4,v1,LSR #16
  350. addin:
  351.          ADDS    v3,v3,a3           ;now add in c
  352.          ADCCS   v4,v4,#0
  353.  
  354.          CMPS    a4,#0x40000000     ;check for m=MAXBASE
  355.          MOVEQ   a1,v3,LSR #30
  356.          ADDEQ   a1,a1,v4,LSL #2
  357.          BICEQ   v4,v3,#0xC0000000
  358.          BEQ     leave
  359.          MOV     a1,#0              ;do long division by m
  360.  
  361. divlp:
  362.  
  363. .REPEAT 32                          ;2xfaster than a loop!
  364.          MOVS    v3,v3,ASL #1       ;get next bit into carry
  365.          ADC     v4,v4,v4           ;accumulate remainder
  366.          CMPS    v4,a4
  367.          SUBCS   v4,v4,a4
  368.          ADC     a1,a1,a1           ;accumulate quotient
  369. .ENDREPEAT
  370.  
  371. leave:
  372.          LDR     v3,[ip]
  373.          STR     v4,[v3]       ;store remainder
  374.          LDMFD   sp!, {v1-v4}
  375.          MOVS pc,lr
  376. */
  377.  
  378.