home *** CD-ROM | disk | FTP | other *** search
/ GEMini Atari / GEMini_Atari_CD-ROM_Walnut_Creek_December_1993.iso / zip / gnu / mntlib16.lzh / MNTLIB16 / _MULDF3.S < prev    next >
Text File  |  1993-07-29  |  3KB  |  136 lines
  1. | double floating point multiplication routine
  2. |
  3. | written by Kai-Uwe Bloem (I5110401@dbstu1.bitnet).
  4. | Based on a 80x86 floating point packet from comp.os.minix, written by P.Housel
  5. |
  6. |
  7. | Revision 1.2, kub 01-90 :
  8. | added support for denormalized numbers
  9. |
  10. | Revision 1.1, kub 12-89 :
  11. | Ported over to 68k assembler
  12. |
  13. | Revision 1.0:
  14. | original 8088 code from P.S.Housel
  15.  
  16. BIAS8    =    0x3FF-1
  17.  
  18.     .text
  19.     .even
  20.     .globl    __muldf3, ___muldf3
  21.  
  22. __muldf3:
  23. ___muldf3:
  24.     lea    sp@(4),a0
  25.     moveml    d2-d7,sp@-
  26.     moveml    a0@,d4-d5/d6-d7 | d4-d5 = v, d6-d7 = u
  27.     subw    #16,sp        | multiplication accumulator
  28.  
  29.     movel    d6,d0        | d0 = u.exp
  30.     swap    d0
  31.     movew    d0,d2        | d2 = u.sign
  32.     lsrw    #4,d0
  33.     andw    #0x07ff,d0    | kill sign bit
  34.  
  35.     movel    d4,d1        | d1 = v.exp
  36.     swap    d1
  37.     eorw    d1,d2        | d2 = u.sign ^ v.sign (in bit 31)
  38.     lsrw    #4,d1
  39.     andw    #0x07ff,d1    | kill sign bit
  40.  
  41.     andl    #0x0fffff,d6    | remove exponent from u.mantissa
  42.     tstw    d0        | check for zero exponent - no leading "1"
  43.     beq    0f
  44.     orl    #0x100000,d6    | restore implied leading "1"
  45.     bra    1f
  46. 0:    addw    #1,d0        | "normalize" exponent
  47. 1:    movel    d6,d3
  48.     orl    d7,d3
  49.     beq    retz        | multiplying by zero
  50.  
  51.     andl    #0x0fffff,d4    | remove exponent from v.mantissa
  52.     tstw    d1        | check for zero exponent - no leading "1"
  53.     beq    0f
  54.     orl    #0x100000,d4    | restore implied leading "1"
  55.     bra    1f
  56. 0:    addw    #1,d1        | "normalize" exponent
  57. 1:    movel    d4,d3
  58.     orl    d5,d3
  59.     beq    retz        | multiplying by zero
  60.  
  61.     addw    d1,d0        | add exponents,
  62.     subw    #BIAS8+16-11,d0    | remove excess bias, acnt for repositioning
  63.  
  64.     clrl    sp@        | initialize 128-bit product to zero
  65.     clrl    sp@(4)
  66.     clrl    sp@(8)
  67.     clrl    sp@(12)
  68.     lea    sp@(8),a1    | accumulator pointer
  69.  
  70. | see Knuth, Seminumerical Algorithms, section 4.3. algorithm M
  71.  
  72.     swap    d2
  73.     movew    #4-1,d2
  74. 1:
  75.     movew    d5,d3
  76.     mulu    d7,d3        | mulitply with bigit from multiplier
  77.     addl    d3,a1@(4)    | store into result
  78.     movew    d4,d3
  79.     mulu    d7,d3
  80.     movel    a1@,d1        | add to result
  81.     addxl    d3,d1
  82.     movel    d1,a1@
  83.     roxlw    a1@-        | rotate carry in
  84.  
  85.     movel    d5,d3
  86.     swap    d3
  87.     mulu    d7,d3
  88.     addl    d3,a1@(4)    | add to result
  89.     movel    d4,d3
  90.     swap    d3
  91.     mulu    d7,d3
  92.     movel    a1@,d1        | add to result
  93.     addxl    d3,d1
  94.     movel    d1,a1@
  95.  
  96.     movew    d6,d7
  97.     swap    d6
  98.     swap    d7
  99.     dbra    d2,1b
  100.  
  101.     swap    d2        | [TOP 16 BITS SHOULD BE ZERO !]
  102.  
  103.     moveml    sp@(2),d4-d7    | get the 112 valid bits
  104.     clrw    d7        | (pad to 128)
  105. 2:
  106.     cmpl    #0x0000ffff,d4    | multiply (shift) until
  107.     bhi    3f        |  1 in upper 16 result bits
  108.     cmpw    #9,d0        | give up for denormalized numbers
  109.     ble    3f
  110.     swap    d4        | (we're getting here only when multiplying
  111.     swap    d5        |  with a denormalized number; there's an
  112.     swap    d6        |  eventual loss of 4 bits in the rounding
  113.     swap    d7        |  byte -- what a pity 8-)
  114.     movew    d5,d4
  115.     movew    d6,d5
  116.     movew    d7,d6
  117.     clrw    d7
  118.     subw    #16,d0        | decrement exponent
  119.     bra    2b
  120. 3:
  121.     movel    d6,d1        | get rounding bits
  122.     roll    #8,d1
  123.     movel    d1,d3        | see if sticky bit should be set
  124.     orl    d7,d3        | (lower 16 bits of d7 are guaranteed to be 0)
  125.     andl    #0xffffff00,d3
  126.     beq    4f
  127.     orb    #1,d1        | set "sticky bit" if any low-order set
  128. 4:    addw    #16,sp        | remove accumulator from stack
  129.     jmp    norm_df        | (result in d4/d5)
  130.  
  131. retz:    clrl    d0        | save zero as result
  132.     clrl    d1
  133.     addw    #16,sp
  134.     moveml    sp@+,d2-d7
  135.     rts            | no normalizing neccessary
  136.