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