![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to cbrt.s CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [CSRG BSD Unix] / 43BSDReno / lib / libm / vax|
Return to cbrt.s CVS log | Up to [CSRG BSD Unix] / 43BSDReno / lib / libm / vax |
1.1 root 1: # Copyright (c) 1985 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.3 (Berkeley) 6/30/88
22: #
23: .data
24: .align 2
25: _sccsid:
26: .asciz "@(#)cbrt.s 1.1 (Berkeley) 5/23/85; 5.3 (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: # error check by E LeBlanc, 8/18/82
31: # Revised and tested by K.C. Ng, 5/2/85
32: # Max error less than 0.667 ulps (unit in the last places)
33: .globl _cbrt
34: .globl _d_cbrt
35: .globl _dcbrt_
36: .text
37: .align 1
38:
39: _cbrt:
40: _d_cbrt:
41: .word 0x00fc # save r2 to r7
42: movq 4(ap),r0 # r0 = argument x
43: jmp dcbrt2
44: _dcbrt_:
45: .word 0x00fc # save r2 to r7
46: movq *4(ap),r0 # r0 = argument x
47:
48: dcbrt2: bicw3 $0x807f,r0,r2 # biased exponent of x
49: jeql return # dcbrt(0)=0 dcbrt(res)=res. operand
50: bicw3 $0x7fff,r0,ap # ap has sign(x)
51: xorw2 ap,r0 # r0 is abs(x)
52: movl r0,r2 # r2 has abs(x)
53: rotl $16,r2,r2 # r2 = |x| with bits unscrambled
54: divl2 $3,r2 # rough dcbrt with bias/3
55: addl2 B,r2 # restore bias, diminish fraction
56: rotl $16,r2,r2 # r2=|q|=|dcbrt| to 5 bits
57: mulf3 r2,r2,r3 # r3 =qq
58: divf2 r0,r3 # r3 = qq/x
59: mulf2 r2,r3
60: addf2 C,r3 # r3 = s = C + qqq/x
61: divf3 r3,D,r4 # r4 = D/s
62: addf2 E,r4
63: addf2 r4,r3 # r3 = s + E + D/s
64: divf3 r3,F,r3 # r3 = F / (s + E + D/s)
65: addf2 G,r3 # r3 = G + F / (s + E + D/s)
66: mulf2 r3,r2 # r2 = qr3 = new q to 23 bits
67: clrl r3 # r2:r3 = q as double float
68: muld3 r2,r2,r4 # r4:r5 = qq exactly
69: divd2 r4,r0 # r0:r1 = x/(q*q) rounded
70: subd3 r2,r0,r6 # r6:r7 = x/(q*q) - q exactly
71: movq r2,r4 # r4:r5 = q
72: addw2 $0x80,r4 # r4:r5 = 2 * q
73: addd2 r0,r4 # r4:r5 = 2*q + x/(q*q)
74: divd2 r4,r6 # r6:r7 = (x/(q*q)-q)/(2*q+x/(q*q))
75: muld2 r2,r6 # r6:r7 = q*(x/(q*q)-q)/(2*q+x/(q*q))
76: addd3 r6,r2,r0 # r0:r1 = q + r6:r7
77: bisw2 ap,r0 # restore the sign bit
78: return:
79: ret # error less than 0.667 ulps
80:
81: .data
82: .align 2
83: B : .long 721142941 # (86-0.03306235651)*(2^23)
84: C : .float 0f0.5428571429 # 19/35
85: D : .float 0f-0.7053061224 # -864/1225
86: E : .float 0f1.414285714 # 99/70
87: F : .float 0f1.607142857 # 45/28
88: G : .float 0f0.3571428571 # 5/14
89:
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.