home *** CD-ROM | disk | FTP | other *** search
/ GEMini Atari / GEMini_Atari_CD-ROM_Walnut_Creek_December_1993.iso / files / gnu / libsrc87 / _muldf3.cpp < prev    next >
Encoding:
Text File  |  1993-07-30  |  4.8 KB  |  244 lines
  1. | mjr: not needed on the TT
  2.  
  3. #ifndef    __M68881__
  4. # ifdef    sfp004
  5.  
  6. | double precision floating point stuff for Atari-gcc using the SFP004
  7. | developed with gas
  8. |
  9. | double precision multiplication
  10. |
  11. | M. Ritzert (mjr at dfg.dbp.de)
  12. |
  13. | 4.10.1990
  14. |
  15. | no NAN checking implemented since the 68881 treats this situation "correctly",
  16. | i.e. according to IEEE
  17.  
  18. | addresses of the 68881 data port. This choice is fastest when much data is
  19. | transferred between the two processors.
  20.  
  21. comm =     -6
  22. resp =    -16
  23. zahl =      0
  24.  
  25. | waiting loop ...
  26. |
  27. | wait:
  28. | ww:    cmpiw    #0x8900,a0@(resp)
  29. |     beq    ww
  30. | is coded directly by
  31. |    .long    0x0c688900, 0xfff067f8
  32.  
  33.     .text
  34.     .even
  35.     .globl    __muldf3, ___muldf3
  36.  
  37. __muldf3:
  38. ___muldf3:
  39.     lea    0xfffa50,a0
  40.     movew    #0x5400,a0@(comm)    | load first argument to fp0
  41.     cmpiw    #0x8900,a0@(resp)    | check
  42.     movel    a7@(4),a0@
  43.     movel    a7@(8),a0@
  44.     movew    #0x5423,a0@(comm)
  45.     .long    0x0c688900, 0xfff067f8
  46.     movel    a7@(12),a0@
  47.     movel    a7@(16),a0@
  48.     movew    #0x7400,a0@(comm)    | result to d0/d1
  49.     .long    0x0c688900, 0xfff067f8
  50.     movel    a0@,d0
  51.     movel    a0@,d1
  52.     rts
  53.  
  54. # else    sfp004
  55.  
  56. | double floating point multiplication routine
  57. |
  58. | written by Kai-Uwe Bloem (I5110401@dbstu1.bitnet).
  59. | Based on a 80x86 floating point packet from comp.os.minix, written by P.Housel
  60. | patched by Olaf Flebbe (flebbe@tat.physik.uni-tuebingen.de)
  61. |
  62. | Revision 1.2.1 olaf 12-92:
  63. |   + added support for NaN and Infinites
  64. |   + added support for -0
  65. |
  66. | Revision 1.2, kub 01-90 :
  67. | added support for denormalized numbers
  68. |
  69. | Revision 1.1, kub 12-89 :
  70. | Ported over to 68k assembler
  71. |
  72. | Revision 1.0:
  73. | original 8088 code from P.S.Housel
  74.  
  75. BIAS8    =    0x3FF-1
  76.  
  77.     .text
  78.     .even
  79.     .globl    __muldf3, ___muldf3
  80.  
  81. __muldf3:
  82. ___muldf3:
  83.     lea    sp@(4),a0
  84.     moveml    d2-d7,sp@-
  85.     moveml    a0@,d4-d5/d6-d7 | d4-d5 = v, d6-d7 = u
  86.  
  87.     movel    d6,d0        | d0 = u.exp
  88.     swap    d0
  89.     movew    d0,d2        | d2 = u.sign
  90.     lsrw    #4,d0
  91.     andw    #0x07ff,d0    | kill sign bit
  92.  
  93.     movel    d4,d1        | d1 = v.exp
  94.     swap    d1
  95.     eorw    d1,d2        | d2 = u.sign ^ v.sign (in bit 31)
  96.     lsrw    #4,d1
  97.     andw    #0x07ff,d1    | kill sign bit
  98.  
  99.     andl    #0x0fffff,d6    | remove exponent from u.mantissa
  100.     andl    #0x0fffff,d4    | remove exponent from v.mantissa
  101. |
  102. | Testing for NaN and Infinities
  103. |
  104.     cmpw    #0x7ff,d0
  105.     beq    0f
  106.     cmpw    #0x7ff,d1
  107.     bne    nospec
  108.     bra    1f
  109. |    first operand is special
  110. |    Nan?
  111. 0:    orl    d7,d6
  112.     bne    retnan
  113. |    Test v==0.
  114.     tstw    d1
  115.     bne    retinf  |  Inf  * x == Inf
  116.     orl    d5,d4
  117.     beq    retnan  | Inf * 0 == NaN
  118. retinf:    
  119.     clrl    d1
  120.     tstl    d2
  121.     bpl    0f
  122.     movel    #0xfff00000,d0
  123. return:    moveml    sp@+,d2-d7
  124.     rts
  125.  
  126.  
  127. 0:    movel    #0x7ff00000,d0
  128.     bra    return    
  129.  
  130. retnan: movel    #0x7fffffff,d0
  131.     moveql    #-1,d1
  132.     bra    return
  133. |
  134. | v is special
  135. |
  136. 1:    orl    d5,d4
  137.     bne    retnan
  138.     tstw    d0
  139.     bne    retinf
  140.     orl    d7,d6
  141.     beq    retnan
  142.     bra    retinf
  143. |
  144. | end of NaN and Inf.
  145. |    
  146. nospec:    subw    #16,sp        | multiplication accumulator
  147.  
  148.     tstw    d0        | check for zero exponent - no leading "1"
  149.     beq    0f
  150.     orl    #0x100000,d6    | restore implied leading "1"
  151.     bra    1f
  152. 0:    addw    #1,d0        | "normalize" exponent
  153. 1:    movel    d6,d3
  154.     orl    d7,d3
  155.     beq    retz        | multiplying by zero
  156.  
  157.     tstw    d1        | check for zero exponent - no leading "1"
  158.     beq    0f
  159.     orl    #0x100000,d4    | restore implied leading "1"
  160.     bra    1f
  161. 0:    addw    #1,d1        | "normalize" exponent
  162. 1:    movel    d4,d3
  163.     orl    d5,d3
  164.     beq    retz        | multiplying by zero
  165.  
  166.     addw    d1,d0        | add exponents,
  167.     subw    #BIAS8+16-11,d0    | remove excess bias, acnt for repositioning
  168.  
  169.     clrl    sp@        | initialize 128-bit product to zero
  170.     clrl    sp@(4)
  171.     clrl    sp@(8)
  172.     clrl    sp@(12)
  173.     lea    sp@(8),a1    | accumulator pointer
  174.  
  175. | see Knuth, Seminumerical Algorithms, section 4.3. algorithm M
  176.  
  177.     swap    d2
  178.     movew    #4-1,d2
  179. 1:
  180.     movew    d5,d3
  181.     mulu    d7,d3        | mulitply with bigit from multiplier
  182.     addl    d3,a1@(4)    | store into result
  183.     movew    d4,d3
  184.     mulu    d7,d3
  185.     movel    a1@,d1        | add to result
  186.     addxl    d3,d1
  187.     movel    d1,a1@
  188.     roxlw    a1@-        | rotate carry in
  189.  
  190.     movel    d5,d3
  191.     swap    d3
  192.     mulu    d7,d3
  193.     addl    d3,a1@(4)    | add to result
  194.     movel    d4,d3
  195.     swap    d3
  196.     mulu    d7,d3
  197.     movel    a1@,d1        | add to result
  198.     addxl    d3,d1
  199.     movel    d1,a1@
  200.  
  201.     movew    d6,d7
  202.     swap    d6
  203.     swap    d7
  204.     dbra    d2,1b
  205.  
  206.     swap    d2        | [TOP 16 BITS SHOULD BE ZERO !]
  207.  
  208.     moveml    sp@(2),d4-d7    | get the 112 valid bits
  209.     clrw    d7        | (pad to 128)
  210. 2:
  211.     cmpl    #0x0000ffff,d4    | multiply (shift) until
  212.     bhi    3f        |  1 in upper 16 result bits
  213.     cmpw    #9,d0        | give up for denormalized numbers
  214.     ble    3f
  215.     swap    d4        | (we''re getting here only when multiplying
  216.     swap    d5        |  with a denormalized number; there''s an
  217.     swap    d6        |  eventual loss of 4 bits in the rounding
  218.     swap    d7        |  byte -- what a pity 8-)
  219.     movew    d5,d4
  220.     movew    d6,d5
  221.     movew    d7,d6
  222.     clrw    d7
  223.     subw    #16,d0        | decrement exponent
  224.     bra    2b
  225. 3:
  226.     movel    d6,d1        | get rounding bits
  227.     roll    #8,d1
  228.     movel    d1,d3        | see if sticky bit should be set
  229.     orl    d7,d3        | (lower 16 bits of d7 are guaranteed to be 0)
  230.     andl    #0xffffff00,d3
  231.     beq    4f
  232.     orb    #1,d1        | set "sticky bit" if any low-order set
  233. 4:    addw    #16,sp        | remove accumulator from stack
  234.     jmp    norm_df        | (result in d4/d5)
  235.  
  236. retz:    clrl    d0        | save zero as result
  237.     clrl    d1
  238.     addw    #16,sp
  239.     moveml    sp@+,d2-d7
  240.     rts            | no normalizing neccessary
  241.  
  242. # endif    sfp004
  243. #endif    __M68881__
  244.