home *** CD-ROM | disk | FTP | other *** search
/ Geek Gadgets 1 / ADE-1.bin / ade-dist / libm-5.4-src.tgz / tar.out / bsd / libm / national / support.s < prev   
Text File  |  1996-09-28  |  11KB  |  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. ; double copysign(x,y)
  37. ; double x,y;
  38. ; IEEE 754 recommended function, return x*sign(y)
  39. ; Coded by K.C. Ng in National 32k assembler, 11/9/85.
  40. ;
  41.     .vers    2
  42.     .text
  43.     .align    2
  44.     .globl    _copysign
  45. _copysign:
  46.     movl    4(sp),f0
  47.     movd    8(sp),r0
  48.     movd    16(sp),r1
  49.     xord    r0,r1
  50.     andd    0x80000000,r1
  51.     cmpqd    0,r1
  52.     beq    end
  53.     negl    f0,f0
  54. end:    ret    0
  55.  
  56. ; double logb(x)
  57. ; double x;
  58. ; IEEE p854 recommended function, return the exponent of x (return float(N) 
  59. ; such that 1 <= x*2**-N < 2, even for subnormal number.
  60. ; Coded by K.C. Ng in National 32k assembler, 11/9/85.
  61. ; Note: subnormal number (if implemented) will be taken care of. 
  62. ;
  63.     .vers    2
  64.     .text
  65.     .align    2
  66.     .globl    _logb
  67. _logb:
  68. ;
  69. ; extract the exponent of x
  70. ; glossaries:    r0 = high part of x
  71. ;        r1 = unbias exponent of x
  72. ;        r2 = 20 (first exponent bit position)
  73. ;
  74.     movd    8(sp),r0
  75.     movd    20,r2
  76.     extd    r2,r0,r1,11    ; extract the exponent of x
  77.     cmpqd    0,r1        ; if exponent bits = 0, goto L3
  78.     beq    L3
  79.     cmpd    0x7ff,r1
  80.     beq    L2        ; if exponent bits = 0x7ff, goto L2
  81. L1:    subd    1023,r1        ; unbias the exponent
  82.     movdl    r1,f0        ; convert the exponent to floating value
  83.     ret    0
  84. ;
  85. ; x is INF or NaN, simply return x
  86. ;
  87. L2:
  88.     movl    4(sp),f0    ; logb(+inf)=+inf, logb(NaN)=NaN
  89.     ret    0
  90. ;
  91. ; x is 0 or subnormal
  92. ;
  93. L3:
  94.     movl    4(sp),f0
  95.     cmpl    0f0,f0
  96.     beq    L5        ; x is 0 , goto L5 (return -inf)
  97. ;
  98. ; Now x is subnormal
  99. ;
  100.     mull    L64,f0        ; scale up f0 with 2**64
  101.     movl    f0,tos
  102.     movd    tos,r0
  103.     movd    tos,r0        ; now r0 = new high part of x
  104.     extd    r2,r0,r1,11    ; extract the exponent of x to r1
  105.     subd    1087,r1        ; unbias the exponent with correction 
  106.     movdl    r1,f0        ; convert the exponent to floating value
  107.     ret    0
  108. ;
  109. ; x is 0, return logb(0)= -INF
  110. ;
  111. L5:
  112.     movl    0f1.0e300,f0
  113.     mull    0f-1.0e300,f0    ; multiply two big numbers to get -INF
  114.     ret    0
  115. ; double rint(x)
  116. ; double x;
  117. ; ... delivers integer nearest x in direction of prevailing rounding
  118. ; ... mode
  119. ; Coded by K.C. Ng in National 32k assembler, 11/9/85.
  120. ; Note: subnormal number (if implemented) will be taken care of. 
  121. ;
  122.     .vers    2
  123.     .text
  124.     .align    2
  125.     .globl    _rint
  126. _rint:
  127. ;
  128.     movd    8(sp),r0
  129.     movd    20,r2
  130.     extd    r2,r0,r1,11    ; extract the exponent of x
  131.     cmpd    0x433,r1
  132.     ble    itself
  133.     movl    L52,f2        ; f2 = L = 2**52
  134.     cmpqd    0,r0
  135.     ble    L1
  136.     negl    f2,f2        ; f2 = s = copysign(L,x)
  137. L1:    addl    f2,f0        ; f0 = x + s
  138.     subl    f2,f0        ; f0 = f0 - s
  139.     ret    0
  140. itself:    movl    4(sp),f0
  141.     ret    0
  142. L52:    .double    0x0,0x43300000    ; L52=2**52
  143. ; int finite(x)
  144. ; double x;
  145. ; IEEE 754 recommended function, return 0 if x is NaN or INF, else 0
  146. ; Coded by K.C. Ng in National 32k assembler, 11/9/85.
  147. ;
  148.     .vers    2
  149.     .text
  150.     .align    2
  151.     .globl    _finite
  152. _finite:
  153.     movd    4(sp),r1
  154.     andd    0x800fffff,r1
  155.     cmpd    0x7ff00000,r1
  156.     sned    r0        ; r0=0 if exponent(x) = 0x7ff
  157.     ret    0
  158. ; double scalb(x,N)
  159. ; double x; int N;
  160. ; IEEE 754 recommended function, return x*2**N by adjusting 
  161. ; exponent of x.
  162. ; Coded by K.C. Ng in National 32k assembler, 11/9/85. 
  163. ; Note: subnormal number (if implemented) will be taken care of 
  164. ;
  165.     .vers    2
  166.     .text
  167.     .align    2
  168.     .globl    _scalb
  169. _scalb:
  170. ;
  171. ; if x=0 return 0
  172. ;
  173.     movl    4(sp),f0
  174.     cmpl    0f0,f0
  175.     beq    end        ; scalb(0,N) is x itself
  176. ;
  177. ; extract the exponent of x
  178. ; glossaries:    r0 = high part of x, 
  179. ;        r1 = unbias exponent of x,
  180. ;        r2 = 20 (first exponent bit position).
  181. ;
  182.     movd    8(sp),r0    ; r0 = high part of x
  183.     movd    20,r2        ; r2 = 20
  184.     extd    r2,r0,r1,11    ; extract the exponent of x in r1
  185.     cmpd    0x7ff,r1    
  186. ;
  187. ; if exponent of x is 0x7ff, then x is NaN or INF; simply return x  
  188. ;
  189.     beq    end        
  190.     cmpqd    0,r1
  191. ;
  192. ; if exponent of x is zero, then x is subnormal; goto L19
  193. ;
  194.     beq    L19        
  195.     addd    12(sp),r1    ; r1 = (exponent of x) + N
  196.     bfs    inof        ; if integer overflows, goto inof
  197.     cmpqd    0,r1        ; if new exponent <= 0, goto underflow
  198.     bge    underflow
  199.     cmpd    2047,r1        ; if new exponent >= 2047 goto overflow
  200.     ble    overflow
  201.     insd    r2,r1,r0,11    ; insert the new exponent 
  202.     movd    r0,tos
  203.     movd    8(sp),tos
  204.     movl    tos,f0        ; return x*2**N
  205. end:    ret    0
  206. inof:    bcs    underflow    ; negative int overflow if Carry bit is set
  207. overflow:
  208.     andd    0x80000000,r0    ; keep the sign of x
  209.     ord    0x7fe00000,r0    ; set x to a huge number
  210.     movd    r0,tos
  211.     movqd    0,tos
  212.     movl    tos,f0
  213.     mull    0f1.0e300,f0    ; multiply two huge number to get overflow
  214.     ret    0
  215. underflow:
  216.     addd    64,r1        ; add 64 to exonent to see if it is subnormal
  217.     cmpqd    0,r1
  218.     bge    zero        ; underflow to zero
  219.     insd    r2,r1,r0,11    ; insert the new exponent 
  220.     movd    r0,tos
  221.     movd    8(sp),tos
  222.     movl    tos,f0
  223.     mull    L30,f0        ; readjust x by multiply it with 2**-64
  224.     ret    0
  225. zero:    andd    0x80000000,r0    ; keep the sign of x
  226.     ord    0x00100000,r0    ; set x to a tiny number
  227.     movd    r0,tos
  228.     movqd    0,tos
  229.     movl    tos,f0
  230.     mull    0f1.0e-300,f0    ; underflow to 0  by multipling two tiny nos.
  231.     ret    0
  232. L19:        ; subnormal number
  233.     mull    L32,f0        ; scale up x by 2**64
  234.     movl    f0,tos
  235.     movd    tos,r0
  236.     movd    tos,r0        ; get the high part of new x
  237.     extd    r2,r0,r1,11    ; extract the exponent of x in r1
  238.     addd    12(sp),r1    ; exponent of x + N
  239.     subd    64,r1        ; adjust it by subtracting 64
  240.     cmpqd    0,r1
  241.     bge    underflow
  242.     cmpd    2047,r1
  243.     ble    overflow
  244.     insd    r2,r1,r0,11    ; insert back the incremented exponent 
  245.     movd    r0,tos
  246.     movd    8(sp),tos
  247.     movl    tos,f0
  248. end:    ret    0
  249. L30:    .double    0x0,0x3bf00000    ; floating point 2**-64
  250. L32:    .double    0x0,0x43f00000    ; floating point 2**64
  251. ; double drem(x,y)
  252. ; double x,y;
  253. ; IEEE double remainder function, return x-n*y, where n=x/y rounded to 
  254. ; nearest integer (half way case goes to even). Result exact.
  255. ; Coded by K.C. Ng in National 32k assembly, 11/19/85.
  256. ;
  257.     .vers    2
  258.     .text
  259.     .align    2
  260.     .globl    _drem
  261. _drem:
  262. ;
  263. ; glossaries:    
  264. ;        r2 = high part of x
  265. ;        r3 = exponent of x
  266. ;        r4 = high part of y
  267. ;        r5 = exponent of y
  268. ;        r6 = sign of x
  269. ;        r7 = constant 0x7ff00000
  270. ;
  271. ;  16(fp) : y
  272. ;   8(fp) : x
  273. ; -12(fp) : adjustment on y when y is subnormal
  274. ; -16(fp) : fsr
  275. ; -20(fp) : nx
  276. ; -28(fp) : t
  277. ; -36(fp) : t1
  278. ; -40(fp) : nf
  279. ;
  280.     enter    [r3,r4,r5,r6,r7],40
  281.     movl    f6,tos
  282.     movl    f4,tos
  283.     movl    0f0,-12(fp)
  284.     movd    0,-20(fp)
  285.     movd    0,-40(fp)
  286.     movd    0x7ff00000,r7    ; initialize r7=0x7ff00000
  287.     movd    12(fp),r2    ; r2 = high(x)
  288.     movd    r2,r3
  289.     andd    r7,r3        ; r3 = xexp
  290.     cmpd    r7,r3
  291. ; if x is NaN or INF goto L1
  292.     beq    L1        
  293.     movd    20(fp),r4
  294.     bicd    [31],r4        ; r4 = high part of |y|
  295.     movd    r4,20(fp)    ; y = |y|
  296.     movd    r4,r5
  297.     andd    r7,r5        ; r5 = yexp
  298.     cmpd    r7,r5
  299.     beq    L2        ; if y is NaN or INF goto L2
  300.     cmpd    0x04000000,r5    ; 
  301.     bgt    L3        ; if y is tiny goto L3
  302. ;
  303. ; now y != 0 , x is finite
  304. ;
  305. L10:
  306.     movd    r2,r6
  307.     andd    0x80000000,r6    ; r6 = sign(x)
  308.     bicd    [31],r2        ; x <- |x|
  309.     sfsr    r1
  310.     movd    r1,-16(fp)    ; save fsr in -16(fp)
  311.     bicd    [5],r1
  312.     lfsr    r1        ; disable inexact interupt
  313.     movd    16(fp),r0    ; r0 = low part of y
  314.     movd    r0,r1        ; r1 = r0 = low part of y
  315.     andd    0xf8000000,r1    ; mask off the lsb 27 bits of y
  316.  
  317.     movd    r2,12(fp)    ; update x to |x|
  318.     movd    r0,-28(fp)    ; 
  319.     movd    r4,-24(fp)    ; t  = y
  320.     movd    r4,-32(fp)    ; 
  321.     movd    r1,-36(fp)    ; t1 = y with trialing 27 zeros
  322.     movd    0x01900000,r1    ; r1 = 25 in exponent field
  323. LOOP:
  324.     movl    8(fp),f0    ; f0 = x
  325.     movl    16(fp),f2    ; f2 = y
  326.     cmpl    f0,f2
  327.     ble    fnad        ; goto fnad (final adjustment) if x <= y
  328.     movd    r4,-32(fp)
  329.     movd    r3,r0
  330.     subd    r5,r0        ; xexp - yexp
  331.     subd    r1,r0        ; r0 = xexp - yexp - m25
  332.     cmpqd    0,r0        ; r0 > 0 ?
  333.     bge    1f
  334.     addd    r4,r0        ; scale up (high) y
  335.     movd    r0,-24(fp)    ; scale up t
  336.     movl    -28(fp),f2    ; t
  337.     movd    r0,-32(fp)    ; scale up t1
  338. 1:
  339.     movl    -36(fp),f4    ; t1
  340.     movl    f0,f6
  341.     divl    f2,f6        ; f6 = x/t
  342.     floorld    f6,r0        ; r0 = [x/t]
  343.     movdl    r0,f6        ; f6 = n
  344.     subl    f4,f2        ; t = t - t1 (tail of t1)
  345.     mull    f6,f4        ; f4 = n*t1    ...exact
  346.     subl    f4,f0        ; x = x - n*t1
  347.     mull    f6,f2        ; n*(t-t1)    ...exact
  348.     subl    f2,f0        ; x = x - n*(t-t1)
  349. ; update xexp
  350.     movl    f0,8(fp)
  351.     movd    12(fp),r3    
  352.     andd    r7,r3
  353.     jump    LOOP
  354. fnad:
  355.     cmpqd    0,-20(fp)    ; 0 = nx?
  356.     beq    final
  357.     mull    -12(fp),8(fp)    ; scale up x the same amount as y
  358.     movd    0,-20(fp)
  359.     movd    12(fp),r2
  360.     movd    r2,r3
  361.     andd    r7,r3        ; update exponent of x
  362.     jump    LOOP
  363.  
  364. final:
  365.     movl    16(fp),f2    ; f2 = y (f0=x, r0=n)
  366.     subd    0x100000,r4    ; high y /2
  367.     movd    r4,-24(fp)
  368.     movl    -28(fp),f4    ; f4 = y/2
  369.     cmpl    f0,f4        ; x > y/2 ?
  370.     bgt    1f
  371.     bne    2f
  372.     andd    1,r0        ; n is odd or even
  373.     cmpqd    0,r0
  374.     beq    2f
  375. 1:
  376.     subl    f2,f0        ; x = x - y
  377. 2:
  378.     cmpqd    0,-40(fp)
  379.     beq    3f
  380.     divl    -12(fp),f0    ; scale down the answer
  381. 3:
  382.     movl    f0,tos
  383.     xord    r6,tos
  384.     movl    tos,f0
  385.     movd    -16(fp),r0
  386.     lfsr    r0        ; restore the fsr
  387.  
  388. end:    movl    tos,f4
  389.     movl    tos,f6
  390.     exit    [r3,r4,r5,r6,r7]
  391.     ret    0
  392. ;
  393. ; y is NaN or INF
  394. ;
  395. L2:    
  396.     movd    16(fp),r0    ; r0 = low part of y
  397.     andd    0xfffff,r4    ; r4 = high part of y & 0x000fffff
  398.     ord    r4,r0
  399.     cmpqd    0,r0
  400.     beq    L4
  401.     movl    16(fp),f0    ; y is NaN, return y
  402.     jump    end
  403. L4:    movl    8(fp),f0    ; y is inf, return x
  404.     jump    end
  405. ;
  406. ; exponent of y is less than 64, y may be zero or subnormal
  407. ;
  408. L3:
  409.     movl    16(fp),f0
  410.     cmpl    0f0,f0
  411.     bne    L5
  412.     divl    f0,f0        ; y is 0, return NaN by doing 0/0
  413.     jump    end
  414. ;
  415. ; subnormal y or tiny y
  416. ;
  417. L5:    
  418.     movd    0x04000000,-20(fp)    ; nx = 64 in exponent field
  419.     movl    L64,f2
  420.     movl    f2,-12(fp)
  421.     mull    f2,f0
  422.     cmpl    f0,LTINY
  423.     bgt    L6
  424.     mull    f2,f0
  425.     addd    0x04000000,-20(fp)    ; nx = nx + 64 in exponent field
  426.     mull    f2,-12(fp)
  427. L6:
  428.     movd    -20(fp),-40(fp)
  429.     movl    f0,16(fp)
  430.     movd    20(fp),r4
  431.     movd    r4,r5
  432.     andd    r7,r5        ; exponent of new y
  433.     jump    L10
  434. ;
  435. ; x is NaN or INF, return x-x
  436. ;
  437. L1:
  438.     movl    8(fp),f0
  439.     subl    f0,f0        ; if x is INF, then INF-INF is NaN
  440.     ret    0
  441. L64:    .double 0x0,0x43f00000    ; L64 = 2**64
  442. LTINY:    .double 0x0,0x04000000    ; LTINY = 2**-959
  443.