Annotation of 43BSDTahoe/usr.lib/libm/tahoe/cbrt.s, revision 1.1.1.1

1.1       root        1: # Copyright (c) 1987 Regents of the University of California.
                      2: # All rights reserved.
                      3: #
                      4: # Redistribution and use in source and binary forms are permitted
                      5: # provided that the above copyright notice and this paragraph are
                      6: # duplicated in all such forms and that any documentation,
                      7: # advertising materials, and other materials related to such
                      8: # distribution and use acknowledge that the software was developed
                      9: # by the University of California, Berkeley.  The name of the
                     10: # University may not be used to endorse or promote products derived
                     11: # from this software without specific prior written permission.
                     12: # THIS SOFTWARE IS PROVIDED ``AS IS'' AND WITHOUT ANY EXPRESS OR
                     13: # IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED
                     14: # WARRANTIES OF MERCHANTIBILITY AND FITNESS FOR A PARTICULAR PURPOSE.
                     15: #
                     16: # All recipients should regard themselves as participants in an ongoing
                     17: # research project and hence should feel obligated to report their
                     18: # experiences (good or bad) with these elementary function codes, using
                     19: # the sendbug(8) program, to the authors.
                     20: #
                     21: #      @(#)cbrt.s      5.4 (Berkeley) 6/30/88
                     22: #
                     23:        .data
                     24:        .align  2
                     25: _sccsid:
                     26: .asciz "@(#)cbrt.s     5.4     (ucb.elefunt) 6/30/88"
                     27: 
                     28: # double cbrt(double arg)
                     29: # W. Kahan, 10/13/80. revised 1/13/84 for keeping sign symmetry
                     30: # Re-coded in tahoe assembly language by Z. Alex Liu (7/13/87)
                     31: # Max error less than 0.667 ULPs _if_ +,-,*,/ were all correctly rounded...
                     32:        .globl  _cbrt
                     33:        .globl  _d_cbrt
                     34:        .globl  _dcbrt_
                     35:        .text
                     36:        .align  2
                     37: _cbrt:
                     38: _d_cbrt:
                     39:        .word   0x01fc          # save r2-r8
                     40:        movl    4(fp),r0        # r0:r1 = x
                     41:        movl    8(fp),r1
                     42:        brb     1f
                     43: _dcbrt_:
                     44:        .word   0x01fc          # save r2-r8
                     45:        movl    4(fp),r8
                     46:        movl    (r8),r0
                     47:        movl    4(r8),r1        # r0:r1 = x
                     48: 
                     49: 1:     andl3   $0x7f800000,r0,r2       # biased exponent of x
                     50:        beql    return          # dcbrt(0)=0  dcbrt(res)=res. operand
                     51:        andl3   $0x80000000,r0,r8       # r8 has sign(x)
                     52:        xorl2   r8,r0           # r0 is abs(x)
                     53:        movl    r0,r2           # r2 has abs(x)
                     54:        divl2   $3,r2           # rough dcbrt with bias/3
                     55:        addl2   B,r2            # restore bias, diminish fraction
                     56:        ldf     r2              # acc = |q|=|dcbrt| to 5 bits
                     57:        mulf    r2              # acc = qq
                     58:        divf    r0              # acc = qq/|x|
                     59:        mulf    r2              # acc = qqq/|x|
                     60:        addf    C               # acc = C+qqq/|x|
                     61:        stf     r3              # r3 = s = C+qqq/|x|
                     62:        ldf     D               # acc = D
                     63:        divf    r3              # acc = D/s
                     64:        addf    E               # acc = E+D/s
                     65:        addf    r3              # acc = s+E+D/s
                     66:        stf     r3              # r3 = s+E+D/s
                     67:        ldf     F               # acc = F
                     68:        divf    r3              # acc = F/(s+E+D/s)
                     69:        addf    G               # acc = G+F/(s+E+D/s)
                     70:        mulf    r2              # acc = q*(G+F/(s+E+D/s)) = new q to 23 bits
                     71:        stf     r2              # r2 = q*(G+F/(s+E+D/s)) = new q to 23 bits
                     72:        clrl    r3              # r2:r3 = q as double float
                     73:        ldd     r2              # acc = q as double float
                     74:        muld    r2              # acc = qq exactly
                     75:        std     r4              # r4:r5 = qq exactly
                     76:        ldd     r0              # acc = |x|
                     77:        divd    r4              # acc = |x|/(q*q) rounded
                     78:        std     r0              # r0:r1 = |x|/(q*q) rounded
                     79:        subd    r2              # acc = |x|/(q*q)-q exactly
                     80:        std     r6              # r6:r7 = |x|/(q*q)-q exactly
                     81:        movl    r2,r4
                     82:        clrl    r5              # r4:r5 = q as double float
                     83:        addl2   $0x800000,r4    # r4:r5 = 2*q
                     84:        ldd     r4              # acc = 2*q
                     85:        addd    r0              # acc = 2*q+|x|/(q*q)
                     86:        std     r4              # r4:r5 = 2*q+|x|/(q*q)
                     87:        ldd     r6              # acc = |x|/(q*q)-q
                     88:        divd    r4              # acc = (|x|/(q*q)-q)/(2*q+|x|/(q*q))
                     89:        muld    r2              # acc = q*(|x|/(q*q)-q)/(2*q+|x|/(q*q))
                     90:        addd    r2              # acc = q+q*(|x|/(q*q)-q)/(2*q+|x|/(q*q))
                     91:        std     r0              # r0:r1 = |result|
                     92:        orl2    r8,r0           # restore the sign bit
                     93: return: ret                    # error less than 0.667ULPs?
                     94: 
                     95:        .data
                     96:        .align  2
                     97: B :    .long   721142941       #(86-0.03306235651)*(2^23)
                     98:        .align  2
                     99: C:     .long   0x400af8b0      #.float 0f0.5428571429  # 19/35
                    100:        .align  2
                    101: D:     .long   0xc0348ef1      #.float 0f-0.7053061224 # -864/1225
                    102:        .align  2
                    103: E:     .long   0x40b50750      #.float 0f1.414285714   # 99/70
                    104:        .align  2
                    105: F:     .long   0x40cdb6db      #.float 0f1.607142857   # 45/28
                    106:        .align  2
                    107: G:     .long   0x3fb6db6e      #.float 0f0.3571428571  # 5/14

unix.superglobalmegacorp.com

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