home *** CD-ROM | disk | FTP | other *** search
/ GEMini Atari / GEMini_Atari_CD-ROM_Walnut_Creek_December_1993.iso / files / gnu / libsrc87 / _mulsf3.cpp < prev    next >
Encoding:
Text File  |  1993-07-30  |  4.2 KB  |  203 lines
  1. | mjr: not needed on the TT
  2.  
  3. #ifndef    __M68881__
  4.  
  5.     .text
  6.     .even
  7.     .globl    __mulsf3, ___mulsf3
  8.  
  9. __mulsf3:
  10. ___mulsf3:
  11.  
  12. # ifdef    sfp004
  13.  
  14. | single precision floating point stuff for Atari-gcc using the SFP004
  15. | developed with gas
  16. |
  17. | single floating point multiplication routine
  18. |
  19. | M. Ritzert (mjr at dfg.dbp.de)
  20. |
  21. | 7. Juli 1989
  22. |
  23. | revision 1.1: adapted to the gcc lib patchlevel 58
  24. | 4.10.90
  25.  
  26. comm =     -6
  27. resp =    -16
  28. zahl =      0
  29.  
  30.     lea    0xfffffa50:w,a0
  31.     movew    #0x4400,a0@(comm)    | load first argument to fp0
  32.     cmpiw    #0x8900,a0@(resp)    | check
  33.     movel    a7@(4),a0@
  34.     movew    #0x4427,a0@(comm)
  35.     .long    0x0c688900, 0xfff067f8
  36.     movel    a7@(8),a0@
  37.     movew    #0x6400,a0@(comm)    | result to d0
  38.     .long    0x0c688900, 0xfff067f8
  39.     movel    a0@,d0
  40.     rts
  41.  
  42. # else    sfp004
  43.  
  44. | single floating point multiplication routine
  45. |
  46. | written by Kai-Uwe Bloem (I5110401@dbstu1.bitnet).
  47. | Based on a 80x86 floating point packet from comp.os.minix, written by P.Housel
  48. | patched by Olaf Flebbe (flebbe@tat.physik.uni-tuebingen.de)
  49. |
  50. | Revision 1.2.1 olaf 12-92:
  51. |   + added support for NaN and Infinites
  52. |   + added support for -0
  53. |
  54. | Revision 1.2, kub 01-90 :
  55. | added support for denormalized numbers
  56. |
  57. | Revision 1.1, kub 12-89 :
  58. | Created single float version for 68000. Code could be speed up by having
  59. | the accumulator in the 68000 register set ...
  60. |
  61. | Revision 1.0:
  62. | original 8088 code from P.S.Housel for double floats
  63.  
  64. BIAS4    =    0x7F-1
  65.  
  66.     lea    sp@(4),a0
  67.     moveml    d2-d5,sp@-
  68.     moveml    a0@,d4/d5    | d4 = v, d5 = u
  69.  
  70.  
  71.     movel    d5,d0        | d0 = u.exp
  72.     swap    d0
  73.     movew    d0,d2        | d2 = u.sign
  74.     lsrw    #7,d0
  75.     andw    #0xff,d0    | kill sign bit
  76.  
  77.     movel    d4,d1        | d1 = v.exp
  78.     swap    d1
  79.     eorw    d1,d2        | d2 = u.sign ^ v.sign (in bit 31)
  80.     lsrw    #7,d1
  81.     andw    #0xff,d1    | kill sign bit
  82.  
  83.     andl    #0x7fffff,d5    | remove exponent from u.mantissa
  84.     andl    #0x7fffff,d4    | remove exponent from v.mantissa
  85. |
  86. | Testing for NaN and Infinities
  87. |
  88.     cmpb    #0xff,d0
  89.     beq    0f
  90.     cmpb    #0xff,d1
  91.     bne    nospec
  92.     bra    1f
  93. |    first operand is special
  94. |    Nan?
  95. 0:    tstl    d5
  96.     bne    retnan
  97. |    Test v==0.
  98.     tstb    d1
  99.     bne    retinf  |  Inf  * x == Inf
  100.     tstl    d4
  101.     beq    retnan  | Inf * 0 == NaN
  102.  
  103. retinf:    tstl    d2
  104.     bpl    0f
  105.     movel    #0xff800000,d0
  106. return:    moveml    sp@+,d2-d5
  107.     rts
  108.  
  109.  
  110. 0:    movel    #0x7f800000,d0
  111.     bra    return    
  112.  
  113. retnan: movel    #0x7fffffff,d0
  114.     bra    return
  115. |
  116. | v is special
  117. |
  118. 1:    tstl    d4
  119.     bne    retnan
  120.     tstb    d0
  121.     bne    retinf
  122.     tstl    d5
  123.     beq    retnan
  124.     bra    retinf
  125. |
  126. | end of NaN and Inf.
  127. |    
  128. nospec:    subw    #8,sp        | multiplication accumulator
  129.     tstw    d0        | check for zero exponent - no leading "1"
  130.     beq    0f
  131.     orl    #0x800000,d5    | restore implied leading "1"
  132.     bra    1f
  133. 0:    addw    #1,d0        | "normalize" exponent
  134. 1:    tstl    d5
  135.     beq    retz        | multiplying zero
  136.  
  137.     tstw    d1        | check for zero exponent - no leading "1"
  138.     beq    0f
  139.     orl    #0x800000,d4    | restore implied leading "1"
  140.     bra    1f
  141. 0:    addw    #1,d1        | "normalize" exponent
  142. 1:    tstl    d4
  143.     beq    retz        | multiply by zero
  144.  
  145.     addw    d1,d0        | add exponents,
  146.     subw    #BIAS4+16-8,d0    | remove excess bias, acnt for repositioning
  147.  
  148.     clrl    sp@        | initialize 64-bit product to zero
  149.     clrl    sp@(4)
  150.  
  151. | see Knuth, Seminumerical Algorithms, section 4.3. algorithm M
  152.  
  153.     movew    d4,d3
  154.     mulu    d5,d3        | mulitply with bigit from multiplier
  155.     movel    d3,sp@(4)    | store into result
  156.  
  157.     movel    d4,d3
  158.     swap    d3
  159.     mulu    d5,d3
  160.     addl    d3,sp@(2)    | add to result
  161.  
  162.     swap    d5        | [TOP 8 BITS SHOULD BE ZERO !]
  163.  
  164.     movew    d4,d3
  165.     mulu    d5,d3        | mulitply with bigit from multiplier
  166.     addl    d3,sp@(2)    | store into result (no carry can occur here)
  167.  
  168.     movel    d4,d3
  169.     swap    d3
  170.     mulu    d5,d3
  171.     addl    d3,sp@        | add to result
  172.                 | [TOP 16 BITS SHOULD BE ZERO !]
  173.     moveml    sp@(2),d4-d5    | get the 48 valid mantissa bits
  174.     clrw    d5        | (pad to 64)
  175. 2:
  176.     cmpl    #0x0000ffff,d4    | multiply (shift) until
  177.     bhi    3f        |  1 in upper 16 result bits
  178.     cmpw    #9,d0        | give up for denormalized numbers
  179.     ble    3f
  180.     swap    d4        | (we''re getting here only when multiplying
  181.     swap    d5        |  with a denormalized number; there''s an
  182.     movew    d5,d4        |  eventual loss of 4 bits in the rounding
  183.     clrw    d5        |  byte -- what a pity 8-)
  184.     subw    #16,d0        | decrement exponent
  185.     bra    2b
  186. 3:
  187.     movel    d5,d1        | get rounding bits
  188.     roll    #8,d1
  189.     movel    d1,d3        | see if sticky bit should be set
  190.     andl    #0xffffff00,d3
  191.     beq    4f
  192.     orb    #1,d1        | set "sticky bit" if any low-order set
  193. 4:    addw    #8,sp        | remove accumulator from stack
  194.     jmp    norm_sf        | (result in d4)
  195.  
  196. retz:    clrl    d0        | save zero as result
  197.     addw    #8,sp
  198.     moveml    sp@+,d2-d5
  199.     rts            | no normalizing neccessary
  200.  
  201. # endif    sfp004
  202. #endif    __M68881__
  203.