home *** CD-ROM | disk | FTP | other *** search
/ Geek Gadgets 1 / ADE-1.bin / ade-dist / libm-5.4-src.tgz / tar.out / bsd / libm / tahoe / cbrt.s < prev    next >
Text File  |  1996-09-28  |  4KB  |  119 lines
  1. # Copyright (c) 1987 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. #    @(#)cbrt.s    5.5 (Berkeley) 10/9/90
  33. #
  34.     .data
  35.     .align    2
  36. _sccsid:
  37. .asciz    "@(#)cbrt.s    5.5    (ucb.elefunt) 10/9/90"
  38.  
  39. # double cbrt(double arg)
  40. # W. Kahan, 10/13/80. revised 1/13/84 for keeping sign symmetry
  41. # Re-coded in tahoe assembly language by Z. Alex Liu (7/13/87)
  42. # Max error less than 0.667 ULPs _if_ +,-,*,/ were all correctly rounded...
  43.     .globl    _cbrt
  44.     .globl    _d_cbrt
  45.     .globl    _dcbrt_
  46.     .text
  47.     .align    2
  48. _cbrt:
  49. _d_cbrt:
  50.     .word    0x01fc        # save r2-r8
  51.     movl    4(fp),r0    # r0:r1 = x
  52.     movl    8(fp),r1
  53.     brb    1f
  54. _dcbrt_:
  55.     .word    0x01fc        # save r2-r8
  56.     movl    4(fp),r8
  57.     movl    (r8),r0
  58.     movl    4(r8),r1    # r0:r1 = x
  59.  
  60. 1:    andl3    $0x7f800000,r0,r2    # biased exponent of x
  61.     beql    return        # dcbrt(0)=0  dcbrt(res)=res. operand
  62.     andl3    $0x80000000,r0,r8    # r8 has sign(x)
  63.     xorl2    r8,r0        # r0 is abs(x)
  64.     movl    r0,r2        # r2 has abs(x)
  65.     divl2    $3,r2        # rough dcbrt with bias/3
  66.     addl2    B,r2        # restore bias, diminish fraction
  67.     ldf    r2        # acc = |q|=|dcbrt| to 5 bits
  68.     mulf    r2        # acc = qq
  69.     divf    r0        # acc = qq/|x|
  70.     mulf    r2        # acc = qqq/|x|
  71.     addf    C        # acc = C+qqq/|x|
  72.     stf    r3        # r3 = s = C+qqq/|x|
  73.     ldf    D        # acc = D
  74.     divf    r3        # acc = D/s
  75.     addf    E        # acc = E+D/s
  76.     addf    r3        # acc = s+E+D/s
  77.     stf    r3        # r3 = s+E+D/s
  78.     ldf    F        # acc = F
  79.     divf    r3        # acc = F/(s+E+D/s)
  80.     addf    G        # acc = G+F/(s+E+D/s)
  81.     mulf    r2        # acc = q*(G+F/(s+E+D/s)) = new q to 23 bits
  82.     stf    r2        # r2 = q*(G+F/(s+E+D/s)) = new q to 23 bits
  83.     clrl    r3        # r2:r3 = q as double float
  84.     ldd    r2        # acc = q as double float
  85.     muld    r2        # acc = qq exactly
  86.     std    r4        # r4:r5 = qq exactly
  87.     ldd    r0        # acc = |x|
  88.     divd    r4        # acc = |x|/(q*q) rounded
  89.     std    r0        # r0:r1 = |x|/(q*q) rounded
  90.     subd    r2        # acc = |x|/(q*q)-q exactly
  91.     std    r6        # r6:r7 = |x|/(q*q)-q exactly
  92.     movl    r2,r4
  93.     clrl    r5        # r4:r5 = q as double float
  94.     addl2    $0x800000,r4    # r4:r5 = 2*q
  95.     ldd    r4        # acc = 2*q
  96.     addd    r0        # acc = 2*q+|x|/(q*q)
  97.     std    r4        # r4:r5 = 2*q+|x|/(q*q)
  98.     ldd    r6        # acc = |x|/(q*q)-q
  99.     divd    r4        # acc = (|x|/(q*q)-q)/(2*q+|x|/(q*q))
  100.     muld    r2        # acc = q*(|x|/(q*q)-q)/(2*q+|x|/(q*q))
  101.     addd    r2        # acc = q+q*(|x|/(q*q)-q)/(2*q+|x|/(q*q))
  102.     std    r0        # r0:r1 = |result|
  103.     orl2    r8,r0        # restore the sign bit
  104. return: ret            # error less than 0.667ULPs?
  105.  
  106.     .data
  107.     .align    2
  108. B :    .long    721142941    #(86-0.03306235651)*(2^23)
  109.     .align    2
  110. C:    .long    0x400af8b0    #.float    0f0.5428571429    # 19/35
  111.     .align    2
  112. D:    .long    0xc0348ef1    #.float    0f-0.7053061224    # -864/1225
  113.     .align    2
  114. E:    .long    0x40b50750    #.float    0f1.414285714    # 99/70
  115.     .align    2
  116. F:    .long    0x40cdb6db    #.float    0f1.607142857    # 45/28
  117.     .align    2
  118. G:    .long    0x3fb6db6e    #.float    0f0.3571428571    # 5/14
  119.