4.3BSD - /usr/src/usr.lib/libm/VAX/cbrt.s

   1: #
   2: # Copyright (c) 1985 Regents of the University of California.
   3: #
   4: # Use and reproduction of this software are granted  in  accordance  with
   5: # the terms and conditions specified in  the  Berkeley  Software  License
   6: # Agreement (in particular, this entails acknowledgement of the programs'
   7: # source, and inclusion of this notice) with the additional understanding
   8: # that  all  recipients  should regard themselves as participants  in  an
   9: # ongoing  research  project and hence should  feel  obligated  to report
  10: # their  experiences (good or bad) with these elementary function  codes,
  11: # using "sendbug 4bsd-bugs@BERKELEY", to the authors.
  12: #
  13: 
  14: # @(#)cbrt.s	1.1 (Berkeley) 5/23/85
  15: 
  16: # double cbrt(double arg)
  17: # W. Kahan, 10/13/80. revised 1/13/84 for keeping sign symmetry
  18: # error check by E LeBlanc, 8/18/82
  19: # Revised and tested by K.C. Ng, 5/2/85
  20: # Max error less than 0.667 ulps (unit in the last places)
  21:         .globl  _cbrt
  22:         .globl  _d_cbrt
  23:         .globl  _dcbrt_
  24:         .text
  25:         .align  1
  26: 
  27: _cbrt:
  28: _d_cbrt:
  29:         .word   0x00fc          # save r2 to r7
  30:         movq    4(ap),r0        # r0 = argument x
  31:         jmp     dcbrt2
  32: _dcbrt_:
  33:         .word   0x00fc          # save r2 to r7
  34:         movq    *4(ap),r0       # r0 = argument x
  35: 
  36: dcbrt2: bicw3   $0x807f,r0,r2   # biased exponent of x
  37:         jeql    return          # dcbrt(0)=0  dcbrt(res)=res. operand
  38:         bicw3   $0x7fff,r0,ap   # ap has sign(x)
  39:         xorw2   ap,r0           # r0 is abs(x)
  40:         movl    r0,r2           # r2 has abs(x)
  41:         rotl    $16,r2,r2       # r2 = |x| with bits unscrambled
  42:         divl2   $3,r2           # rough dcbrt with bias/3
  43:         addl2   B,r2            # restore bias, diminish fraction
  44:         rotl    $16,r2,r2       # r2=|q|=|dcbrt| to 5 bits
  45:         mulf3   r2,r2,r3        # r3 =qq
  46:         divf2   r0,r3           # r3 = qq/x
  47:         mulf2   r2,r3
  48:         addf2   C,r3            # r3 = s = C + qqq/x
  49:         divf3   r3,D,r4         # r4 = D/s
  50:         addf2   E,r4
  51:         addf2   r4,r3           # r3 = s + E + D/s
  52:         divf3   r3,F,r3         # r3 = F / (s + E + D/s)
  53:         addf2   G,r3            # r3 = G + F / (s + E + D/s)
  54:         mulf2   r3,r2           # r2 = qr3 = new q to 23 bits
  55:         clrl    r3              # r2:r3 = q as double float
  56:         muld3   r2,r2,r4        # r4:r5 = qq exactly
  57:         divd2   r4,r0           # r0:r1 = x/(q*q) rounded
  58:         subd3   r2,r0,r6        # r6:r7 = x/(q*q) - q exactly
  59:         movq    r2,r4           # r4:r5 = q
  60:         addw2   $0x80,r4        # r4:r5 = 2 * q
  61:         addd2   r0,r4           # r4:r5 = 2*q + x/(q*q)
  62:         divd2   r4,r6           # r6:r7 = (x/(q*q)-q)/(2*q+x/(q*q))
  63:         muld2   r2,r6           # r6:r7 = q*(x/(q*q)-q)/(2*q+x/(q*q))
  64:         addd3   r6,r2,r0        # r0:r1 = q + r6:r7
  65:         bisw2   ap,r0           # restore the sign bit
  66: return:
  67:         ret                     # error less than 0.667 ulps
  68: 
  69: .data
  70: .align  2
  71: B :     .long            721142941              # (86-0.03306235651)*(2^23)
  72: C :     .float          0f0.5428571429          # 19/35
  73: D :     .float          0f-0.7053061224         # -864/1225
  74: E :     .float          0f1.414285714           # 99/70
  75: F :     .float          0f1.607142857           # 45/28
  76: G :     .float          0f0.3571428571          # 5/14

Defined functions

Defined variables