home *** CD-ROM | disk | FTP | other *** search
/ InfoMagic Source Code 1993 July / THE_SOURCE_CODE_CD_ROM.iso / bsd_srcs / lib / libm / vax / cbrt.s < prev    next >
Encoding:
Text File  |  1991-04-12  |  3.7 KB  |  101 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. #    @(#)cbrt.s    5.4 (Berkeley) 10/9/90
  33. #
  34.     .data
  35.     .align    2
  36. _sccsid:
  37. .asciz    "@(#)cbrt.s    1.1 (Berkeley) 5/23/85; 5.4 (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. # error check by E LeBlanc, 8/18/82
  42. # Revised and tested by K.C. Ng, 5/2/85  
  43. # Max error less than 0.667 ulps (unit in the last places)
  44.     .globl    _cbrt
  45.     .globl    _d_cbrt
  46.     .globl  _dcbrt_
  47.     .text
  48.     .align    1
  49.  
  50. _cbrt:
  51. _d_cbrt:
  52.     .word    0x00fc        # save r2 to r7
  53.     movq    4(ap),r0    # r0 = argument x
  54.     jmp     dcbrt2
  55. _dcbrt_:
  56.     .word    0x00fc        # save r2 to r7
  57.     movq    *4(ap),r0    # r0 = argument x
  58.  
  59. dcbrt2:    bicw3    $0x807f,r0,r2    # biased exponent of x
  60.     jeql    return        # dcbrt(0)=0  dcbrt(res)=res. operand
  61.     bicw3    $0x7fff,r0,ap    # ap has sign(x)
  62.     xorw2    ap,r0        # r0 is abs(x)
  63.     movl    r0,r2        # r2 has abs(x)
  64.     rotl    $16,r2,r2    # r2 = |x| with bits unscrambled
  65.     divl2    $3,r2        # rough dcbrt with bias/3
  66.     addl2    B,r2        # restore bias, diminish fraction
  67.     rotl    $16,r2,r2    # r2=|q|=|dcbrt| to 5 bits
  68.     mulf3    r2,r2,r3    # r3 =qq
  69.     divf2    r0,r3        # r3 = qq/x
  70.     mulf2    r2,r3
  71.     addf2    C,r3        # r3 = s = C + qqq/x
  72.     divf3    r3,D,r4        # r4 = D/s
  73.     addf2    E,r4
  74.     addf2    r4,r3        # r3 = s + E + D/s
  75.     divf3    r3,F,r3        # r3 = F / (s + E + D/s)
  76.     addf2    G,r3        # r3 = G + F / (s + E + D/s)
  77.     mulf2    r3,r2        # r2 = qr3 = new q to 23 bits
  78.     clrl    r3        # r2:r3 = q as double float
  79.     muld3    r2,r2,r4    # r4:r5 = qq exactly
  80.     divd2    r4,r0        # r0:r1 = x/(q*q) rounded
  81.     subd3    r2,r0,r6    # r6:r7 = x/(q*q) - q exactly
  82.     movq    r2,r4        # r4:r5 = q
  83.     addw2    $0x80,r4    # r4:r5 = 2 * q
  84.     addd2    r0,r4        # r4:r5 = 2*q + x/(q*q)
  85.     divd2    r4,r6        # r6:r7 = (x/(q*q)-q)/(2*q+x/(q*q))
  86.     muld2    r2,r6        # r6:r7 = q*(x/(q*q)-q)/(2*q+x/(q*q))
  87.     addd3    r6,r2,r0    # r0:r1 = q + r6:r7
  88.     bisw2    ap,r0        # restore the sign bit
  89. return:
  90.     ret            # error less than 0.667 ulps
  91.  
  92. .data
  93. .align    2
  94. B :    .long         721142941        # (86-0.03306235651)*(2^23)
  95. C :    .float        0f0.5428571429        # 19/35
  96. D :    .float        0f-0.7053061224        # -864/1225
  97. E :    .float        0f1.414285714        # 99/70
  98. F :    .float        0f1.607142857        # 45/28
  99. G :    .float        0f0.3571428571        # 5/14
  100.  
  101.