|
|
1.1 root 1: # double cbrt(arg)
2: # double arg
3: # no error exits
4: #method: range reduction to [1/8,1], poly appox, newtons method
5: # J F Jarvis, August 10,1978
6: .globl _cbrt
7: .text
8: .align 1
9: _cbrt:
10: .word 0x00c0
11: bispsw $0xe0
12: clrl r3
13: movd 4(ap),r4
14: jgtr range
15: jeql retz
16: mnegd r4,r4 # |arg| in r0,r1
17: movl $0x100,r3 # sign bit of result
18: range:
19: extzv $7,$8,r4,r6
20: insv $128,$7,$8,r4 # 0.5<= frac: r4,r5 <1.0
21: clrl r7
22: ediv $3,r6,r6,r7 # r6= expnt/3; r7= expnt%3
23: addb2 $86,r6
24: bisl2 r3,r6 # sign,exponent of result
25: polyf r4,$3,pcoef # initial estimate is Hart&Cheney CBRT 0642
26: # D=4.1
27: muld3 r0,r0,r2 # Newtons method, iteration 1, H&C 6.1.10
28: divd3 r2,r4,r2
29: subd3 r2,r0,r2
30: muld2 third,r2
31: subd2 r2,r0 # D=8.2
32: muld3 r0,r0,r2 # iteration 2
33: divd3 r2,r4,r2
34: subd3 r2,r0,r2
35: muld2 third,r2
36: subd2 r2,r0
37: muld2 hc[r7],r0 # set range
38: insv r6,$7,$9,r0 # set sign,exponent
39: ret
40: retz:
41: clrd r0
42: ret
43: .data
44: .align 2
45: third: .double 0d0.33333333333333333333e+0
46: hc:
47: .double 0d1.25992104989487316476e+0
48: .double 0d1.58740105196819947475e+0
49: .double 0d1.0e+0
50: pcoef:
51: .float 0f0.1467073818e+0
52: .float 0f-0.5173964673e+0
53: .float 0f0.9319858515e+0
54: .float 0f0.4387762363e+0
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.