home *** CD-ROM | disk | FTP | other *** search
/ Atari FTP / ATARI_FTP_0693.zip / ATARI_FTP_0693 / Mint / mntlib32.zoo / _divdf3.cpp < prev    next >
C/C++ Source or Header  |  1993-06-17  |  6KB  |  261 lines
  1. |
  2. | double floating point divide routine
  3. |
  4. #ifndef __M68881__
  5.     .text
  6.     .even
  7.     .globl    __divdf3, ___divdf3
  8. #ifdef    ERROR_CHECK
  9. #include "errbase.h"
  10.     .globl    __infinitydf
  11. LC0:
  12.     .ascii "floating point division by 0\12\15\0"
  13.     .even
  14. #endif    ERROR_CHECK
  15.  
  16. __divdf3:
  17. ___divdf3:
  18. #ifdef    ERROR_CHECK
  19.     tstl    a7@(12)            | check if divisor is 0
  20.     bne    continue
  21.     tstl    a7@(16)
  22.     bne    continue
  23.  
  24.     pea    pc@(LC0)
  25.     pea    Stderr
  26.     jbsr    _fprintf
  27.     addql    #8,a7
  28.  
  29.     moveq    #Erange,d0        | set _errno to ERANGE
  30.     Emove    d0,Errno
  31.     moveml    __infinitydf,d0-d1    | return signed infinity
  32.     btst    #31,a7@(4)        | transfer sign of dividend
  33.     beq    clear            | (mjr++)
  34.     bset    #31,d0            |
  35.     rts                |
  36. clear:                    |
  37.     bclr    #31,d0            |
  38.     rts
  39. continue:
  40.  
  41. #endif    /* ERROR_CHECK */
  42. #ifndef sfp004
  43. |
  44. | written by Kai-Uwe Bloem (I5110401@dbstu1.bitnet).
  45. | Based on a 80x86 floating point packet from comp.os.minix, written by P.Housel
  46. |
  47. | Revision 1.2.4 michal 05-93 (ntomczak@vm.ucs.ualberta.ca)
  48. |   + resynchro with errno codes
  49. |   + code smoothing
  50. |
  51. | patched by Olaf Flebbe (flebbe@tat.physik.uni-tuebingen.de)
  52. |
  53. | Revision 1.2.3 olaf 4-93
  54. |   + Fixed sign for retinf, and retzero: it is in d2.w
  55. |
  56. | Revision 1.2.2 olaf 12-92
  57. |   + added support for NaN and Infinites
  58. |   + added support for -0
  59. |
  60. | Revision 1.2.1 olaf 11-92
  61. |   + prevent the tie rounding case if dividing is not exact.
  62. |      > paranoia now says: "Division appears to round correctly"
  63. |      ** requires _normdf Version 1.6.1 or later
  64. |
  65. | Revision 1.2, kub 01-90 :
  66. | added support for denormalized numbers
  67. |
  68. | Revision 1.1, kub 12-89 :
  69. | Ported over to 68k assembler
  70. |
  71. | Revision 1.0:
  72. | original 8088 code from P.S.Housel
  73.  
  74. BIAS8    =    0x3FF-1
  75.  
  76.     lea    sp@(4),a0    | pointer to parameters u and v
  77.     moveml    d2-d7,sp@-    | save registers
  78.     moveml    a0@,d4-d5/d6-d7    | d4-d5 = u, d6-d7 = v
  79.  
  80.     movel    #0x0fffff,d3
  81.     movel    d4,d0        | d0 = u.exp
  82.     andl    d3,d4        | remove exponent from u.mantissa
  83.     swap    d0
  84.     movew    d0,d2        | d2 = u.sign
  85.  
  86.     movel    d6,d1        | d1 = v.exp
  87.     andl    d3,d6        | remove exponent from v.mantissa
  88.     swap    d1
  89.     eorw    d1,d2        | d2 = u.sign ^ v.sign (in bit 15)
  90.  
  91.     moveq    #15,d3
  92.     bclr    d3,d1        | kill sign bit
  93.     bclr    d3,d0        | kill sign bit
  94.     lsrw    #4,d0
  95.     lsrw    #4,d1
  96. |
  97. |
  98. |
  99.     movew    #0x7ff,d3
  100.     cmpw    d3,d0
  101.     beq    0f        |u == NaN || u== Inf
  102.     cmpw    d3,d1
  103.     beq    1f        | v == NaN || v == Inf
  104.     tstw    d0
  105.     bne    3f        | u not zero nor denorm
  106.     movel    d5,d3
  107.     orl    d4,d3
  108.     beq    2f        | 0/ ?
  109.  
  110. 3:    tstw    d1
  111.     bne    nospec
  112.  
  113.     movel    d7,d3
  114.     orl    d6,d3
  115.     bne    nospec
  116. |    bra    retinf        | x/0 -> +/- Inf
  117. |
  118. |    Return Infinity with correct sign
  119. |
  120. retinf:    moveq    #0,d1
  121.     movel    #0xffe00000,d0
  122.     lslw    #1,d2
  123.     roxrl   #1,d0        | shift in high bit as given by d2
  124. return:    moveml    sp@+,d2-d7
  125.     rts
  126.  
  127. 0:    orl    d5,d4        | u == NaN ?
  128.     bne    retnan        | NaN/ x
  129.     cmpw    d3,d1
  130.     beq    retnan        | Inf/Inf or Inf/NaN
  131.     bra    retinf        | Inf/x | x != Inf && x != NaN
  132.  
  133. 1:    orl    d7,d6
  134.     bne    retnan        | x/NaN
  135. |    bra    retzero        | x/Inf -> +/- 0
  136. |
  137. |    Return correct signed zero
  138. |
  139. retzero:moveq    #0,d0        | zero destination
  140.     movel    d0,d1
  141.     lslw    #1,d2        | we need an extension bit
  142.     roxrl    #1,d0
  143.     bra    return
  144.  
  145. 2:    tstw    d1
  146.     bne    retzero        | 0/x ->+/- 0
  147.     orl    d5,d4
  148.     bne    retzero        | 0/x
  149. |    bra    retnan        | 0/0
  150. |
  151. |    Return NaN
  152. |
  153. retnan: moveql    #-1,d1
  154.     movel    d1,d0
  155.     lsrl    #1,d0        | 0x7fffffff -> d0
  156.     bra    return
  157. |
  158. |    End of special handling
  159. |
  160. nospec:    moveq    #20,d3
  161.     bset    d3,d4        | restore implied leading "1"
  162.     tstw    d0        | check for zero exponent - no leading "1"
  163.     bne    1f
  164.     bclr    d3,d4        | remove it
  165.     addw    #1,d0        | "normalize" exponent
  166.  
  167. 1:    bset    d3,d6        | restore implied leading "1"
  168.     tstw    d1        | check for zero exponent - no leading "1"
  169.     bne    1f
  170.     bclr    d3,d6        | remove it
  171. 0:    addw    #1,d1        | "normalize" exponent
  172.  
  173. 1:    movew    d2,a0        | save sign
  174.  
  175.     subw    d1,d0        | subtract exponents,
  176.     addw    #BIAS8-11+1,d0    |  add bias back in, account for shift
  177.     addw    #66,d0        |  add loop offset, +2 for extra rounding bits
  178.                 |   for denormalized numbers (2 implied by dbra)
  179.     movew    #24,d1        | bit number for "implied" pos (+4 for rounding)
  180.     moveq    #-1,d2        | zero the quotient
  181.     moveq    #-1,d3        |  (for speed it is a one''s complement)
  182.     subl    d7,d5        | initial subtraction,
  183.     subxl    d6,d4        | u = u - v
  184. 2:
  185.     btst    d1,d2        | divide until 1 in implied position
  186.     beq    5f
  187.  
  188.     addl    d5,d5
  189.     addxl    d4,d4
  190.     bcs    4f        | if carry is set, add, else subtract
  191.  
  192.     addxl    d3,d3        | shift quotient and set bit zero
  193.     addxl    d2,d2
  194.     subl    d7,d5        | subtract
  195.     subxl    d6,d4        | u = u - v
  196.     dbra    d0,2b        | give up if result is denormalized
  197.     bra    5f
  198. 4:
  199.     addxl    d3,d3        | shift quotient and clear bit zero
  200.     addxl    d2,d2
  201.     addl    d7,d5        | add (restore)
  202.     addxl    d6,d4        | u = u + v
  203.     dbra    d0,2b        | give up if result is denormalized
  204. 5:    subqw    #2,d0        | remove rounding offset for denormalized nums
  205.     notl    d2        | invert quotient to get it right
  206.     notl    d3
  207.  
  208.     movel   d5,d1
  209.     orl     d4,d1           | check for exact result
  210.     beq     1f
  211.     moveql  #-1,d1          | Set rounding bits for tie case
  212. 1:    movel    d2,d4        | save quotient mantissa
  213.     movel    d3,d5
  214.     movew    a0,d2        | get sign back
  215.     jmp    norm_df        | (registers on stack removed by norm_df)
  216.  
  217. #else
  218.  
  219. | double precision floating point stuff for Atari-gcc using the SFP004
  220. | developed with gas
  221. |
  222. | double precision division
  223. |
  224. | M. Ritzert (mjr at dmzrzu71)
  225. |
  226. | 4.10.1990
  227. |
  228. | no NAN checking implemented to gain compatibility with the TT-lib
  229. |
  230. | addresses of the 68881 data port. This choice is fastest when much data is
  231. | transferred between the two processors.
  232.  
  233. comm =     -6
  234. resp =    -16
  235. zahl =      0
  236.  
  237. | waiting loop ...
  238. |
  239. | wait:
  240. | ww:    cmpiw    #0x8900,a0@(resp)
  241. |     beq    ww
  242. | is coded directly by
  243. |    .long    0x0c688900, 0xfff067f8
  244.  
  245.     lea    0xfffffa50:w,a0
  246.     movew    #0x5400,a0@(comm)    | load first argument to fp0
  247.     cmpiw    #0x8900,a0@(resp)    | check
  248.     movel    a7@(4),a0@
  249.     movel    a7@(8),a0@
  250.     movew    #0x5420,a0@(comm)
  251.     .long    0x0c688900, 0xfff067f8
  252.     movel    a7@(12),a0@
  253.     movel    a7@(16),a0@
  254.     movew    #0x7400,a0@(comm)    | result to d0
  255.     .long    0x0c688900, 0xfff067f8
  256.     movel    a0@,d0
  257.     movel    a0@,d1
  258.      rts
  259. #endif    sfp004
  260. #endif /* !__M68881__ */
  261.