home *** CD-ROM | disk | FTP | other *** search
/ GEMini Atari / GEMini_Atari_CD-ROM_Walnut_Creek_December_1993.iso / files / gnu / libsrc87 / _divsf3.cpp < prev    next >
Encoding:
C/C++ Source or Header  |  1993-07-30  |  5.6 KB  |  258 lines
  1. |
  2. | single floating point divide routine
  3. |
  4. #ifndef __M68881__
  5.     .text
  6.     .even
  7.     .globl    __divsf3, ___divsf3
  8. #ifdef    ERROR_CHECK
  9. #include "errbase.h"
  10.     .globl    __infinitysf
  11.  
  12. LC0:
  13.     .ascii "floating point division by 0\12\15\0"
  14.     .even
  15. #endif    ERROR_CHECK
  16.  
  17. __divsf3:
  18. ___divsf3:
  19.  
  20. #ifdef    ERROR_CHECK
  21.     tstl    a7@(8)            | check if divisor is 0
  22.     bne    no_exception
  23.  
  24.     pea    pc@(LC0)
  25.     pea    Stderr
  26.     jbsr    _fprintf    |
  27.     addql    #8,a7        |
  28.                     | set _errno to ERANGE
  29.     moveq    #ERANGE,d0
  30.     Emove    d0,Errno
  31.     movel    __infinitysf,d0        | 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.     rts
  38.  
  39. no_exception:
  40. #endif    ERROR_CHECK
  41.  
  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. | patched by Olaf Flebbe (flebbe@tat.physik.uni-tuebingen.de)
  47. |
  48. | Revision 1.2.2 olaf 12-92
  49. |   + added support for NaN and Infinites
  50. |   + added support for -0
  51. |
  52. | Revision 1.2.1 olaf 11-92
  53. |   + prevent the tie rounding case if dividing is not exact.
  54. |      > paranoia now says: "Division appears to round correctly"
  55. |      ** requires _normsf Version 1.4.2 or later
  56. |
  57. | Revision 1.2, kub 01-90 :
  58. | added support for denormalized numbers
  59. |
  60. | Revision 1.1, kub 12-89 :
  61. | Created single float version for 68000
  62. |
  63. | Revision 1.0:
  64. | original 8088 code from P.S.Housel for double floats
  65.  
  66. BIAS4    =    0x7F-1
  67.  
  68.     lea    sp@(4),a0    | pointer to parameters u and v
  69.     moveml    d2-d5,sp@-    | save registers
  70.     moveml    a0@,d4/d5    | d4 = u, d5 = v
  71.  
  72.     movel    d4,d0        | d0 = u.exp
  73.     swap    d0
  74.     movew    d0,d2        | d2 = u.sign
  75.     lsrw    #7,d0
  76.     andw    #0xff,d0    | kill sign bit
  77.  
  78.     movel    d5,d1        | d1 = v.exp
  79.     swap    d1
  80.     eorw    d1,d2        | d2 = u.sign ^ v.sign (in bit 31)
  81.     lsrw    #7,d1
  82.     andw    #0xff,d1    | kill sign bit
  83.  
  84.     andl    #0x7fffff,d4    | remove exponent from u.mantissa
  85.     andl    #0x7fffff,d5    | remove exponent from v.mantissa
  86. |
  87. |
  88. |
  89.     cmpb    #0xff,d0    
  90.     beq    0f        |u == NaN || u== Inf
  91.     cmpb    #0xff,d1
  92.     beq    1f        | v == NaN || v == Inf
  93.     tstb    d0
  94.     bne    3f        | u not zero nor denorm
  95.     tstl    d4
  96.     beq    2f        | 0/ ?
  97.     
  98. 3:    tstw    d1
  99.     bne    nospec
  100.  
  101.     tstl    d5
  102.     bne    nospec
  103.     bra    retinf        | x/0 -> +/- Inf
  104.  
  105. 0:    tstl    d4        | u == NaN ?
  106.     bne    retnan        | NaN/ x
  107.     cmpb    #0xff,d1    
  108.     beq    retnan        | Inf/Inf or Inf/NaN 
  109.     bra    retinf        | Inf/x | x != Inf && x != NaN
  110.  
  111. 1:    tstl    d5
  112.     bne    retnan        | x/NaN
  113.     bra    retzero        | x/Inf -> +/- 0
  114.  
  115. 2:    tstw    d1
  116.     bne    retzero        | 0/x ->+/- 0
  117.     tstl    d4
  118.     bne    retzero        | 0/x 
  119.     bra    retnan        | 0/0
  120. |
  121. |    Return Infinity with correct sign
  122. |    
  123. retinf:    tstl    d2
  124.     bpl    0f
  125.     movel    #0xff800000,d0
  126. return:    moveml    sp@+,d2-d5
  127.     rts
  128.  
  129. 0:    movel    #0x7f800000,d0
  130.     bra    return    
  131. |
  132. |    Return NaN
  133. |
  134. retnan: movel    #0x7fffffff,d0
  135.     bra    return
  136. |
  137. |    Return correct signed zero
  138. |
  139. retzero:clrl    d0        | zero destination
  140.     tstl    d2
  141.     bmi    return
  142.     bset    #31,d0
  143.     bra    return
  144. |
  145. |    End of special handling
  146. |    
  147. nospec:    tstw    d0        | check for zero exponent - no leading "1"
  148.     beq    0f
  149.     orl    #0x800000,d4    | restore implied leading "1"
  150.     bra    1f
  151. 0:    addw    #1,d0        | "normalize" exponent
  152. 1:
  153. # ifndef ERROR_CHECK
  154.     tstl    d4
  155.     beq    retz        | dividing zero
  156. # endif    ERROR_CHECK
  157.  
  158.     tstw    d1        | check for zero exponent - no leading "1"
  159.     beq    0f
  160.     orl    #0x800000,d5    | restore implied leading "1"
  161.     bra    1f
  162. 0:    addw    #1,d1        | "normalize" exponent
  163. 1:    tstl    d5
  164.     beq    divz        | divide by zero
  165.  
  166.     subw    d1,d0        | subtract exponents,
  167.     addw    #BIAS4-8+1,d0    |  add bias back in, account for shift
  168.     addw    #34,d0        |  add loop offset, +2 for extra rounding bits
  169.                 |   for denormalized numbers (2 implied by dbra)
  170.     movew    #27,d1        | bit number for "implied" pos (+4 for rounding)
  171.     movel    #-1,d3        |  zero quotient (for speed a one''s complement)
  172.     subl    d5,d4        | initial subtraction, u = u - v
  173. 2:
  174.     btst    d1,d3        | divide until 1 in implied position
  175.     beq    5f
  176.  
  177.     addl    d4,d4
  178.     bcs    4f        | if carry is set, add, else subtract
  179.  
  180.     addxl    d3,d3        | shift quotient and set bit zero
  181.     subl    d5,d4        | subtract, u = u - v
  182.     dbra    d0,2b        | give up if result is denormalized
  183.     bra    5f
  184. 4:
  185.     addxl    d3,d3        | shift quotient and clear bit zero
  186.     addl    d5,d4        | add (restore), u = u + v
  187.     dbra    d0,2b        | give up if result is denormalized
  188. 5:    subw    #2,d0        | remove rounding offset for denormalized nums
  189.     notl    d3        | invert quotient to get it right
  190.  
  191.     clrl    d1        | zero rounding bits
  192.     tstl     d4        | check for exact result
  193.     beq    1f
  194.     moveql    #-1,d1        | prevent tie case
  195. 1:    movel    d3,d4        | save quotient mantissa
  196.     jmp    norm_sf        | (registers on stack removed by norm_sf)
  197.  
  198. # ifndef ERROR_CHECK
  199. retz:    clrl    d0        | zero destination
  200.     moveml    sp@+,d2-d5
  201.     rts            | no normalization needed
  202.  
  203. divz:    movel    __infinitysf,d0    | return infinity value
  204.     moveml    sp@+,d2-d5    | should really cause trap ?!?
  205.     btst    #31,sp@(4)    | transfer sign of dividend
  206.     beq    clear        | (mjr++)
  207.     bset    #31,d0        |
  208.     rts            |
  209. clear:                |
  210.     bclr    #31,d0        |
  211.     rts
  212. # endif    ERROR_CHECK
  213. #else
  214.  
  215. | single precision floating point stuff for Atari-gcc using the SFP004
  216. | or compatible boards with a memory mapped 68881 
  217. | developed with gas
  218. |
  219. |  single floating point divide routine
  220. |
  221. | M. Ritzert (mjr at dmzrzu71)
  222. |            (ritzert@dfg.dbp.de)
  223. | 4.10.1990
  224. |
  225. | +_infinitysf returned instead of a NAN
  226. | the DOMAIN exception is not supported yet. In case of an exception
  227. | _errno is always set to ERANGE
  228.  
  229. | addresses of the 68881 data port. This choice is fastest when much data is
  230. | transferred between the two processors.
  231.  
  232. comm =     -6
  233. resp =    -16
  234. zahl =      0
  235.  
  236. | waiting loop ...
  237. |
  238. | wait:
  239. | ww:    cmpiw    #0x8900,a0@(resp)
  240. |     beq    ww
  241. | is coded directly by
  242. |    .long    0x0c688900, 0xfff067f8
  243.  
  244.     lea    0xfffffa50:w,a0
  245.     movew    #0x4400,a0@(comm)    | load first argument to fp0
  246.     cmpiw    #0x8900,a0@(resp)    | check
  247.     movel    a7@(4),a0@
  248.     movew    #0x4424,a0@(comm)
  249.     .long    0x0c688900, 0xfff067f8
  250.     movel    a7@(8),a0@
  251.     movew    #0x6400,a0@(comm)    | result to d0
  252.     .long    0x0c688900, 0xfff067f8
  253.     movel    a0@,d0            | REMARK: 0/0 returns a NAN
  254.     rts                | if ERROR_CHECK is disabled
  255.  
  256. #endif    sfp004
  257. #endif /* !__M68881__ */
  258.