home *** CD-ROM | disk | FTP | other *** search
/ Fresh Fish 4 / FreshFish_May-June1994.bin / bbs / gnu / libm-src.lha / src / amiga / libm / national / support.s < prev   
Encoding:
Text File  |  1993-09-23  |  10.6 KB  |  450 lines
  1. ; Copyright (c) 1985 Regents of the University of California.
  2. ; All rights reserved.
  3. ;
  4. ; Redistribution and use in source and binary forms, with or without
  5. ; modification, are permitted provided that the following conditions
  6. ; are met:
  7. ; 1. Redistributions of source code must retain the above copyright
  8. ;    notice, this list of conditions and the following disclaimer.
  9. ; 2. Redistributions in binary form must reproduce the above copyright
  10. ;    notice, this list of conditions and the following disclaimer in the
  11. ;    documentation and/or other materials provided with the distribution.
  12. ; 3. All advertising materials mentioning features or use of this software
  13. ;    must display the following acknowledgement:
  14. ;    This product includes software developed by the University of
  15. ;    California, Berkeley and its contributors.
  16. ; 4. Neither the name of the University nor the names of its contributors
  17. ;    may be used to endorse or promote products derived from this software
  18. ;    without specific prior written permission.
  19. ;
  20. ; THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
  21. ; ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  22. ; IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  23. ; ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
  24. ; FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
  25. ; DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
  26. ; OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
  27. ; HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
  28. ; LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
  29. ; OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
  30. ; SUCH DAMAGE.
  31. ;
  32. ;    @(#)support.s    5.4 (Berkeley) 10/9/90
  33. ;
  34.  
  35. ; IEEE recommended functions
  36. ; 
  37. ; double copysign(x,y)
  38. ; double x,y;
  39. ; IEEE 754 recommended function, return x*sign(y)
  40. ; Coded by K.C. Ng in National 32k assembler, 11/9/85.
  41. ;
  42.     .vers    2
  43.     .text
  44.     .align    2
  45.     .globl    _copysign
  46. _copysign:
  47.     movl    4(sp),f0
  48.     movd    8(sp),r0
  49.     movd    16(sp),r1
  50.     xord    r0,r1
  51.     andd    0x80000000,r1
  52.     cmpqd    0,r1
  53.     beq    end
  54.     negl    f0,f0
  55. end:    ret    0
  56.  
  57. ; 
  58. ; double logb(x)
  59. ; double x;
  60. ; IEEE p854 recommended function, return the exponent of x (return float(N) 
  61. ; such that 1 <= x*2**-N < 2, even for subnormal number.
  62. ; Coded by K.C. Ng in National 32k assembler, 11/9/85.
  63. ; Note: subnormal number (if implemented) will be taken care of. 
  64. ;
  65.     .vers    2
  66.     .text
  67.     .align    2
  68.     .globl    _logb
  69. _logb:
  70. ;
  71. ; extract the exponent of x
  72. ; glossaries:    r0 = high part of x
  73. ;        r1 = unbias exponent of x
  74. ;        r2 = 20 (first exponent bit position)
  75. ;
  76.     movd    8(sp),r0
  77.     movd    20,r2
  78.     extd    r2,r0,r1,11    ; extract the exponent of x
  79.     cmpqd    0,r1        ; if exponent bits = 0, goto L3
  80.     beq    L3
  81.     cmpd    0x7ff,r1
  82.     beq    L2        ; if exponent bits = 0x7ff, goto L2
  83. L1:    subd    1023,r1        ; unbias the exponent
  84.     movdl    r1,f0        ; convert the exponent to floating value
  85.     ret    0
  86. ;
  87. ; x is INF or NaN, simply return x
  88. ;
  89. L2:
  90.     movl    4(sp),f0    ; logb(+inf)=+inf, logb(NaN)=NaN
  91.     ret    0
  92. ;
  93. ; x is 0 or subnormal
  94. ;
  95. L3:
  96.     movl    4(sp),f0
  97.     cmpl    0f0,f0
  98.     beq    L5        ; x is 0 , goto L5 (return -inf)
  99. ;
  100. ; Now x is subnormal
  101. ;
  102.     mull    L64,f0        ; scale up f0 with 2**64
  103.     movl    f0,tos
  104.     movd    tos,r0
  105.     movd    tos,r0        ; now r0 = new high part of x
  106.     extd    r2,r0,r1,11    ; extract the exponent of x to r1
  107.     subd    1087,r1        ; unbias the exponent with correction 
  108.     movdl    r1,f0        ; convert the exponent to floating value
  109.     ret    0
  110. ;
  111. ; x is 0, return logb(0)= -INF
  112. ;
  113. L5:
  114.     movl    0f1.0e300,f0
  115.     mull    0f-1.0e300,f0    ; multiply two big numbers to get -INF
  116.     ret    0
  117. ; 
  118. ; double rint(x)
  119. ; double x;
  120. ; ... delivers integer nearest x in direction of prevailing rounding
  121. ; ... mode
  122. ; Coded by K.C. Ng in National 32k assembler, 11/9/85.
  123. ; Note: subnormal number (if implemented) will be taken care of. 
  124. ;
  125.     .vers    2
  126.     .text
  127.     .align    2
  128.     .globl    _rint
  129. _rint:
  130. ;
  131.     movd    8(sp),r0
  132.     movd    20,r2
  133.     extd    r2,r0,r1,11    ; extract the exponent of x
  134.     cmpd    0x433,r1
  135.     ble    itself
  136.     movl    L52,f2        ; f2 = L = 2**52
  137.     cmpqd    0,r0
  138.     ble    L1
  139.     negl    f2,f2        ; f2 = s = copysign(L,x)
  140. L1:    addl    f2,f0        ; f0 = x + s
  141.     subl    f2,f0        ; f0 = f0 - s
  142.     ret    0
  143. itself:    movl    4(sp),f0
  144.     ret    0
  145. L52:    .double    0x0,0x43300000    ; L52=2**52
  146. ; 
  147. ; int finite(x)
  148. ; double x;
  149. ; IEEE 754 recommended function, return 0 if x is NaN or INF, else 0
  150. ; Coded by K.C. Ng in National 32k assembler, 11/9/85.
  151. ;
  152.     .vers    2
  153.     .text
  154.     .align    2
  155.     .globl    _finite
  156. _finite:
  157.     movd    4(sp),r1
  158.     andd    0x800fffff,r1
  159.     cmpd    0x7ff00000,r1
  160.     sned    r0        ; r0=0 if exponent(x) = 0x7ff
  161.     ret    0
  162. ; 
  163. ; double scalb(x,N)
  164. ; double x; int N;
  165. ; IEEE 754 recommended function, return x*2**N by adjusting 
  166. ; exponent of x.
  167. ; Coded by K.C. Ng in National 32k assembler, 11/9/85. 
  168. ; Note: subnormal number (if implemented) will be taken care of 
  169. ;
  170.     .vers    2
  171.     .text
  172.     .align    2
  173.     .globl    _scalb
  174. _scalb:
  175. ;
  176. ; if x=0 return 0
  177. ;
  178.     movl    4(sp),f0
  179.     cmpl    0f0,f0
  180.     beq    end        ; scalb(0,N) is x itself
  181. ;
  182. ; extract the exponent of x
  183. ; glossaries:    r0 = high part of x, 
  184. ;        r1 = unbias exponent of x,
  185. ;        r2 = 20 (first exponent bit position).
  186. ;
  187.     movd    8(sp),r0    ; r0 = high part of x
  188.     movd    20,r2        ; r2 = 20
  189.     extd    r2,r0,r1,11    ; extract the exponent of x in r1
  190.     cmpd    0x7ff,r1    
  191. ;
  192. ; if exponent of x is 0x7ff, then x is NaN or INF; simply return x  
  193. ;
  194.     beq    end        
  195.     cmpqd    0,r1
  196. ;
  197. ; if exponent of x is zero, then x is subnormal; goto L19
  198. ;
  199.     beq    L19        
  200.     addd    12(sp),r1    ; r1 = (exponent of x) + N
  201.     bfs    inof        ; if integer overflows, goto inof
  202.     cmpqd    0,r1        ; if new exponent <= 0, goto underflow
  203.     bge    underflow
  204.     cmpd    2047,r1        ; if new exponent >= 2047 goto overflow
  205.     ble    overflow
  206.     insd    r2,r1,r0,11    ; insert the new exponent 
  207.     movd    r0,tos
  208.     movd    8(sp),tos
  209.     movl    tos,f0        ; return x*2**N
  210. end:    ret    0
  211. inof:    bcs    underflow    ; negative int overflow if Carry bit is set
  212. overflow:
  213.     andd    0x80000000,r0    ; keep the sign of x
  214.     ord    0x7fe00000,r0    ; set x to a huge number
  215.     movd    r0,tos
  216.     movqd    0,tos
  217.     movl    tos,f0
  218.     mull    0f1.0e300,f0    ; multiply two huge number to get overflow
  219.     ret    0
  220. underflow:
  221.     addd    64,r1        ; add 64 to exonent to see if it is subnormal
  222.     cmpqd    0,r1
  223.     bge    zero        ; underflow to zero
  224.     insd    r2,r1,r0,11    ; insert the new exponent 
  225.     movd    r0,tos
  226.     movd    8(sp),tos
  227.     movl    tos,f0
  228.     mull    L30,f0        ; readjust x by multiply it with 2**-64
  229.     ret    0
  230. zero:    andd    0x80000000,r0    ; keep the sign of x
  231.     ord    0x00100000,r0    ; set x to a tiny number
  232.     movd    r0,tos
  233.     movqd    0,tos
  234.     movl    tos,f0
  235.     mull    0f1.0e-300,f0    ; underflow to 0  by multipling two tiny nos.
  236.     ret    0
  237. L19:        ; subnormal number
  238.     mull    L32,f0        ; scale up x by 2**64
  239.     movl    f0,tos
  240.     movd    tos,r0
  241.     movd    tos,r0        ; get the high part of new x
  242.     extd    r2,r0,r1,11    ; extract the exponent of x in r1
  243.     addd    12(sp),r1    ; exponent of x + N
  244.     subd    64,r1        ; adjust it by subtracting 64
  245.     cmpqd    0,r1
  246.     bge    underflow
  247.     cmpd    2047,r1
  248.     ble    overflow
  249.     insd    r2,r1,r0,11    ; insert back the incremented exponent 
  250.     movd    r0,tos
  251.     movd    8(sp),tos
  252.     movl    tos,f0
  253. end:    ret    0
  254. L30:    .double    0x0,0x3bf00000    ; floating point 2**-64
  255. L32:    .double    0x0,0x43f00000    ; floating point 2**64
  256. ; 
  257. ; double drem(x,y)
  258. ; double x,y;
  259. ; IEEE double remainder function, return x-n*y, where n=x/y rounded to 
  260. ; nearest integer (half way case goes to even). Result exact.
  261. ; Coded by K.C. Ng in National 32k assembly, 11/19/85.
  262. ;
  263.     .vers    2
  264.     .text
  265.     .align    2
  266.     .globl    _drem
  267. _drem:
  268. ;
  269. ; glossaries:    
  270. ;        r2 = high part of x
  271. ;        r3 = exponent of x
  272. ;        r4 = high part of y
  273. ;        r5 = exponent of y
  274. ;        r6 = sign of x
  275. ;        r7 = constant 0x7ff00000
  276. ;
  277. ;  16(fp) : y
  278. ;   8(fp) : x
  279. ; -12(fp) : adjustment on y when y is subnormal
  280. ; -16(fp) : fsr
  281. ; -20(fp) : nx
  282. ; -28(fp) : t
  283. ; -36(fp) : t1
  284. ; -40(fp) : nf
  285. ; 
  286. ;
  287.     enter    [r3,r4,r5,r6,r7],40
  288.     movl    f6,tos
  289.     movl    f4,tos
  290.     movl    0f0,-12(fp)
  291.     movd    0,-20(fp)
  292.     movd    0,-40(fp)
  293.     movd    0x7ff00000,r7    ; initialize r7=0x7ff00000
  294.     movd    12(fp),r2    ; r2 = high(x)
  295.     movd    r2,r3
  296.     andd    r7,r3        ; r3 = xexp
  297.     cmpd    r7,r3
  298. ; if x is NaN or INF goto L1
  299.     beq    L1        
  300.     movd    20(fp),r4
  301.     bicd    [31],r4        ; r4 = high part of |y|
  302.     movd    r4,20(fp)    ; y = |y|
  303.     movd    r4,r5
  304.     andd    r7,r5        ; r5 = yexp
  305.     cmpd    r7,r5
  306.     beq    L2        ; if y is NaN or INF goto L2
  307.     cmpd    0x04000000,r5    ; 
  308.     bgt    L3        ; if y is tiny goto L3
  309. ;
  310. ; now y != 0 , x is finite
  311. ;
  312. L10:
  313.     movd    r2,r6
  314.     andd    0x80000000,r6    ; r6 = sign(x)
  315.     bicd    [31],r2        ; x <- |x|
  316.     sfsr    r1
  317.     movd    r1,-16(fp)    ; save fsr in -16(fp)
  318.     bicd    [5],r1
  319.     lfsr    r1        ; disable inexact interupt
  320.     movd    16(fp),r0    ; r0 = low part of y
  321.     movd    r0,r1        ; r1 = r0 = low part of y
  322.     andd    0xf8000000,r1    ; mask off the lsb 27 bits of y
  323.  
  324.     movd    r2,12(fp)    ; update x to |x|
  325.     movd    r0,-28(fp)    ; 
  326.     movd    r4,-24(fp)    ; t  = y
  327.     movd    r4,-32(fp)    ; 
  328.     movd    r1,-36(fp)    ; t1 = y with trialing 27 zeros
  329.     movd    0x01900000,r1    ; r1 = 25 in exponent field
  330. LOOP:
  331.     movl    8(fp),f0    ; f0 = x
  332.     movl    16(fp),f2    ; f2 = y
  333.     cmpl    f0,f2
  334.     ble    fnad        ; goto fnad (final adjustment) if x <= y
  335.     movd    r4,-32(fp)
  336.     movd    r3,r0
  337.     subd    r5,r0        ; xexp - yexp
  338.     subd    r1,r0        ; r0 = xexp - yexp - m25
  339.     cmpqd    0,r0        ; r0 > 0 ?
  340.     bge    1f
  341.     addd    r4,r0        ; scale up (high) y
  342.     movd    r0,-24(fp)    ; scale up t
  343.     movl    -28(fp),f2    ; t
  344.     movd    r0,-32(fp)    ; scale up t1
  345. 1:
  346.     movl    -36(fp),f4    ; t1
  347.     movl    f0,f6
  348.     divl    f2,f6        ; f6 = x/t
  349.     floorld    f6,r0        ; r0 = [x/t]
  350.     movdl    r0,f6        ; f6 = n
  351.     subl    f4,f2        ; t = t - t1 (tail of t1)
  352.     mull    f6,f4        ; f4 = n*t1    ...exact
  353.     subl    f4,f0        ; x = x - n*t1
  354.     mull    f6,f2        ; n*(t-t1)    ...exact
  355.     subl    f2,f0        ; x = x - n*(t-t1)
  356. ; update xexp
  357.     movl    f0,8(fp)
  358.     movd    12(fp),r3    
  359.     andd    r7,r3
  360.     jump    LOOP
  361. fnad:
  362.     cmpqd    0,-20(fp)    ; 0 = nx?
  363.     beq    final
  364.     mull    -12(fp),8(fp)    ; scale up x the same amount as y
  365.     movd    0,-20(fp)
  366.     movd    12(fp),r2
  367.     movd    r2,r3
  368.     andd    r7,r3        ; update exponent of x
  369.     jump    LOOP
  370.  
  371. final:
  372.     movl    16(fp),f2    ; f2 = y (f0=x, r0=n)
  373.     subd    0x100000,r4    ; high y /2
  374.     movd    r4,-24(fp)
  375.     movl    -28(fp),f4    ; f4 = y/2
  376.     cmpl    f0,f4        ; x > y/2 ?
  377.     bgt    1f
  378.     bne    2f
  379.     andd    1,r0        ; n is odd or even
  380.     cmpqd    0,r0
  381.     beq    2f
  382. 1:
  383.     subl    f2,f0        ; x = x - y
  384. 2:
  385.     cmpqd    0,-40(fp)
  386.     beq    3f
  387.     divl    -12(fp),f0    ; scale down the answer
  388. 3:
  389.     movl    f0,tos
  390.     xord    r6,tos
  391.     movl    tos,f0
  392.     movd    -16(fp),r0
  393.     lfsr    r0        ; restore the fsr
  394.  
  395. end:    movl    tos,f4
  396.     movl    tos,f6
  397.     exit    [r3,r4,r5,r6,r7]
  398.     ret    0
  399. ;
  400. ; y is NaN or INF
  401. ;
  402. L2:    
  403.     movd    16(fp),r0    ; r0 = low part of y
  404.     andd    0xfffff,r4    ; r4 = high part of y & 0x000fffff
  405.     ord    r4,r0
  406.     cmpqd    0,r0
  407.     beq    L4
  408.     movl    16(fp),f0    ; y is NaN, return y
  409.     jump    end
  410. L4:    movl    8(fp),f0    ; y is inf, return x
  411.     jump    end
  412. ;
  413. ; exponent of y is less than 64, y may be zero or subnormal
  414. ;
  415. L3:
  416.     movl    16(fp),f0
  417.     cmpl    0f0,f0
  418.     bne    L5
  419.     divl    f0,f0        ; y is 0, return NaN by doing 0/0
  420.     jump    end
  421. ;
  422. ; subnormal y or tiny y
  423. ;
  424. L5:    
  425.     movd    0x04000000,-20(fp)    ; nx = 64 in exponent field
  426.     movl    L64,f2
  427.     movl    f2,-12(fp)
  428.     mull    f2,f0
  429.     cmpl    f0,LTINY
  430.     bgt    L6
  431.     mull    f2,f0
  432.     addd    0x04000000,-20(fp)    ; nx = nx + 64 in exponent field
  433.     mull    f2,-12(fp)
  434. L6:
  435.     movd    -20(fp),-40(fp)
  436.     movl    f0,16(fp)
  437.     movd    20(fp),r4
  438.     movd    r4,r5
  439.     andd    r7,r5        ; exponent of new y
  440.     jump    L10
  441. ;
  442. ; x is NaN or INF, return x-x
  443. ;
  444. L1:
  445.     movl    8(fp),f0
  446.     subl    f0,f0        ; if x is INF, then INF-INF is NaN
  447.     ret    0
  448. L64:    .double 0x0,0x43f00000    ; L64 = 2**64
  449. LTINY:    .double 0x0,0x04000000    ; LTINY = 2**-959
  450.