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     ! 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.