home *** CD-ROM | disk | FTP | other *** search
/ InfoMagic Source Code 1993 July / THE_SOURCE_CODE_CD_ROM.iso / bsd_srcs / lib / libm / tahoe / support.s < prev    next >
Encoding:
Text File  |  1991-04-12  |  7.2 KB  |  240 lines
  1. /*
  2.  * Copyright (c) 1987 Regents of the University of California.
  3.  * All rights reserved.
  4.  *
  5.  * Redistribution and use in source and binary forms, with or without
  6.  * modification, are permitted provided that the following conditions
  7.  * are met:
  8.  * 1. Redistributions of source code must retain the above copyright
  9.  *    notice, this list of conditions and the following disclaimer.
  10.  * 2. Redistributions in binary form must reproduce the above copyright
  11.  *    notice, this list of conditions and the following disclaimer in the
  12.  *    documentation and/or other materials provided with the distribution.
  13.  * 3. All advertising materials mentioning features or use of this software
  14.  *    must display the following acknowledgement:
  15.  *    This product includes software developed by the University of
  16.  *    California, Berkeley and its contributors.
  17.  * 4. Neither the name of the University nor the names of its contributors
  18.  *    may be used to endorse or promote products derived from this software
  19.  *    without specific prior written permission.
  20.  *
  21.  * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
  22.  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  23.  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  24.  * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
  25.  * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
  26.  * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
  27.  * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
  28.  * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
  29.  * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
  30.  * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
  31.  * SUCH DAMAGE.
  32.  */
  33.     .data
  34.     .align    2
  35. _sccsid:
  36.     .asciz    "@(#)support.s    5.6    (ucb.elefunt)    10/9/90"
  37. /*
  38.  * copysign(x,y),
  39.  * logb(x),
  40.  * scalb(x,N),
  41.  * finite(x),
  42.  * drem(x,y),
  43.  * Coded in vax assembly language by K. C. Ng 4/9/85.
  44.  * Re-coded in tahoe assembly language by Z. Alex Liu 7/13/87.
  45.  */
  46. /*
  47.  * double copysign(x,y)
  48.  * double x,y;
  49.  */
  50.     .globl    _copysign
  51.     .text
  52.     .align    2
  53. _copysign:
  54.     .word    0x0004            # save r2
  55.     movl    8(fp),r1
  56.     movl    4(fp),r0        # r0:r1 = x
  57.     andl3    $0x7f800000,r0,r2    # r2 = biased exponent of x
  58.     beql    1f            # if 0 or reserved op then return x
  59.     andl3    $0x80000000,12(fp),r2    # r2 = sign bit of y at bit-31
  60.     andl2    $0x7fffffff,r0        # replace x by |x|
  61.     orl2    r2,r0            # copy the sign bit of y to x
  62. 1:    ret
  63. /*
  64.  * double logb(x)
  65.  * double x;
  66.  */
  67.     .globl    _logb
  68.     .text
  69.     .align    2
  70. _logb:
  71.     .word    0x0000            # save nothing
  72.     andl3    $0x7f800000,4(fp),r0    # r0[b23:b30] = biased exponent of x
  73.     beql    1f
  74.     shrl    $23,r0,r0        # r0[b0:b7] = biased exponent of x
  75.     subl2    $129,r0            # r0 = unbiased exponent of x
  76.     cvld    r0            # acc = unbiased exponent of x (double)
  77.     std    r0            # r0 =  unbiased exponent of x (double)
  78.     ret
  79. 1:    movl    8(fp),r1        # 8(fp) must be moved first
  80.     movl    4(fp),r0        # r0:r1 = x (zero or reserved op)
  81.     blss    2f            # simply return if reserved op
  82.     movl    $0xfe000000,r1
  83.     movl    $0xcfffffff,r0        # -2147483647.0
  84. 2:    ret
  85. /*
  86.  * long finite(x)
  87.  * double x;
  88.  */
  89.     .globl    _finite
  90.     .text
  91.     .align    2
  92. _finite:
  93.     .word    0x0000            # save nothing
  94.     andl3    $0xff800000,4(fp),r0    # r0 = sign of x & its biased exponent
  95.     cmpl    r0,$0x80000000        # is x a reserved op?
  96.     beql    1f            # if so, return FALSE (0)
  97.     movl    $1,r0            # else return TRUE (1)
  98.     ret
  99. 1:    clrl    r0
  100.     ret
  101. /*
  102.  * double scalb(x,N)
  103.  * double x; int N;
  104.  */
  105.     .globl    _scalb
  106.     .set    ERANGE,34
  107.     .text
  108.     .align    2
  109. _scalb:
  110.     .word    0x000c            # save r2-r3
  111.     movl    8(fp),r1
  112.     movl    4(fp),r0        # r0:r1 = x (-128 <= Ex <= 126)
  113.     andl3    $0x7f800000,r0,r3    # r3[b23:b30] = biased exponent of x
  114.     beql    1f            # is x a 0 or a reserved operand?
  115.     movl    12(fp),r2        # r2 = N
  116.     cmpl    r2,$0xff        # if N >= 255
  117.     bgeq    2f            # then the result must overflow
  118.     cmpl    r2,$-0xff        # if N <= -255
  119.     bleq    3f            # then the result must underflow
  120.     shrl    $23,r3,r3        # r3[b0:b7] = biased exponent of x
  121.     addl2    r2,r3            # r3 = biased exponent of the result
  122.     bleq    3f            # if <= 0 then the result underflows
  123.     cmpl    r3,$0x100        # if >= 256 then the result overflows
  124.     bgeq    2f
  125.     shll    $23,r3,r3        # r3[b23:b30] = biased exponent of res.
  126.     andl2    $0x807fffff,r0
  127.     orl2    r3,r0            # r0:r1 = x*2^N
  128. 1:    ret
  129. 2:    pushl    $ERANGE            # if the result would overflow
  130.     callf    $8,_infnan        # and _infnan returns
  131.     andl3    $0x80000000,4(fp),r2    # get the sign of input arg
  132.     orl2    r2,r0            # re-attach the sign to r0:r1
  133.     ret
  134. 3:    clrl    r1            # if the result would underflow
  135.     clrl    r0            # then return 0
  136.     ret
  137. /*
  138.  * double drem(x,y)
  139.  * double x,y;
  140.  * Returns x-n*y where n=[x/y] rounded (to even in the half way case).
  141.  */
  142.     .globl    _drem
  143.     .set    EDOM,33
  144.     .text
  145.     .align    2
  146. _drem:
  147.     .word    0x1ffc            # save r2-r12
  148.     movl    16(fp),r3
  149.     movl    12(fp),r2        # r2:r3 = y
  150.     movl    8(fp),r1
  151.     movl    4(fp),r0        # r0:r1 = x
  152.     andl3    $0xff800000,r0,r4
  153.     cmpl    r4,$0x80000000        # is x a reserved operand?
  154.     beql    1f            # if yes then propagate x and return
  155.     andl3    $0xff800000,r2,r4
  156.     cmpl    r4,$0x80000000        # is y a reserved operand?
  157.     bneq    2f
  158.     movl    r3,r1
  159.     movl    r2,r0            # if yes then propagate y and return
  160. 1:    ret
  161.  
  162. 2:    tstl    r4            # is y a 0?
  163.     bneq    3f
  164.     pushl    $EDOM            # if so then generate reserved op fault
  165.     callf    $8,_infnan
  166.     ret
  167.  
  168. 3:    andl2    $0x7fffffff,r2        # r2:r3 = y <- |y|
  169.     clrl    r12            # r12 = nx := 0
  170.     cmpl    r2,$0x1c800000        # Ey ? 57 
  171.     bgtr    4f            # if Ey > 57 goto 4
  172.     addl2    $0x1c800000,r2        # scale up y by 2**57
  173.     movl    $0x1c800000,r12        # r12[b23:b30] = nx = 57
  174. 4:    pushl    r12            # pushed onto stack: nf := nx
  175.     andl3    $0x80000000,r0,-(sp)    # pushed onto stack: sign of x
  176.     andl2    $0x7fffffff,r0        # r0:r1 = x <- |x|
  177.     movl    r3,r11            # r10:r11 = y1 = y w/ last 27 bits 0
  178.     andl3    $0xf8000000,r10,r11    # clear last 27 bits of y1
  179.  
  180. Loop:    cmpd2    r0,r2            # x ? y
  181.     bleq    6f            # if x <= y goto 6
  182.  /*                     # begin argument reduction */
  183.     movl    r3,r5
  184.     movl    r2,r4            # r4:r5 = t = y
  185.     movl    r11,r7
  186.     movl    r10,r6            # r6:r7 = t1 = y1
  187.     andl3    $0x7f800000,r0,r8    # r8[b23:b30] = Ex:biased exponent of x
  188.     andl3    $0x7f800000,r2,r9    # r9[b23:b30] = Ey:biased exponent of y
  189.     subl2    r9,r8            # r8[b23:b30] = Ex-Ey
  190.     subl2    $0x0c800000,r8        # r8[b23:b30] = k = Ex-Ey-25
  191.     blss    5f            # if k < 0 goto 5
  192.     addl2    r8,r4            # t += k    
  193.     addl2    r8,r6            # t1 += k, scale up t and t1
  194. 5:    ldd    r0            # acc = x
  195.     divd    r4            # acc = x/t
  196.     cvdl    r8            # r8 = n = [x/t] truncated
  197.     cvld    r8            # acc = dble(n)
  198.     std    r8            # r8:r9 = dble(n)
  199.     ldd    r4            # acc = t
  200.     subd    r6            # acc = t-t1
  201.     muld    r8            # acc = n*(t-t1)
  202.     std    r4            # r4:r5 = n*(t-t1)
  203.     ldd    r6            # acc = t1
  204.     muld    r8            # acc = n*t1
  205.     subd    r0            # acc = n*t1-x
  206.     negd                # acc = x-n*t1
  207.     subd    r4            # acc = (x-n*t1)-n*(t-t1)
  208.     std    r0            # r0:r1 = (x-n*t1)-n*(t-t1)
  209.     brb    Loop
  210.  
  211. 6:    movl    r12,r6            # r6 = nx
  212.     beql    7f            # if nx == 0 goto 7
  213.     addl2    r6,r0            # x <- x*2**57:scale x up by nx
  214.     clrl    r12            # clear nx
  215.     brb    Loop
  216.  
  217. 7:    movl    r3,r5
  218.     movl    r2,r4            # r4:r5 = y
  219.     subl2    $0x800000,r4        # r4:r5 = y/2
  220.     cmpd2    r0,r4            # x ? y/2
  221.     blss    9f            # if x < y/2 goto 9
  222.     bgtr    8f            # if x > y/2 goto 8
  223.     ldd    r8            # acc = dble(n)
  224.     cvdl    r8            # r8 = ifix(dble(n))
  225.     bbc    $0,r8,9f        # if the last bit is zero, goto 9
  226. 8:    ldd    r0            # acc = x
  227.     subd    r2            # acc = x-y
  228.     std    r0            # r0:r1 = x-y
  229. 9:    xorl2    (sp)+,r0        # x^sign (exclusive or)
  230.     movl    (sp)+,r6        # r6 = nf
  231.     andl3    $0x7f800000,r0,r8    # r8 = biased exponent of x
  232.     andl2    $0x807fffff,r0        # r0 = x w/ exponent zapped
  233.     subl2    r6,r8            # r8 = Ex-nf
  234.     bgtr    0f            # if Ex-nf > 0 goto 0
  235.     clrl    r8            # clear r8
  236.     clrl    r0
  237.     clrl    r1            # x underflows to zero
  238. 0:    orl2    r8,r0            # put r8 into x's exponent field
  239.     ret
  240.