Return to cbrt.s CVS log

Up to [CSRG BSD Unix] / 43BSD / usr.lib / libm / VAX

Annotation of 43BSD/usr.lib/libm/VAX/cbrt.s, revision 1.1.1.1

1.1       root        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
                     77: 

This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.