|
|
1.1 root 1: # double sqrt(arg):revised July 18,1980
2: # double arg
3: # if(arg<0.0) { _errno=EDOM; return(-sqrt(-arg)) }
4: # W. Kahan's magic sqrt
5: # coded by Heidi Stettner
6:
7: .set EDOM,98
8: .text
9: .align 1
10: .globl _sqrt
11: .globl dsqrt_r5
12: .globl _errno
13: _sqrt:
14: .word 0x003c # save r5,r4,r3,r2
15: movd 4(ap),r0
16: bsbb dsqrt_r5
17: ret
18: dsqrt_r5:
19: movd r0,r4
20: jleq nonpos #argument is not positive
21: movzwl r4,r2
22: ashl $-1,r2,r0
23: addw2 $0x203c,r0 #r0 has magic initial appx
24:
25: # Do two steps of Heron's rule
26:
27: divf3 r0,r4,r2
28: addf2 r2,r0
29: subw2 $0x80,r0
30:
31: divf3 r0,r4,r2
32: addf2 r2,r0
33: subw2 $0x80,r0
34:
35:
36: # Scale argument and approximation to prevent over/underflow
37: # NOTE: The following four steps would not be necessary if underflow
38: # were gentle.
39:
40: bicw3 $0xffff807f,r4,r1
41: subw2 $0x4080,r1 # r1 contains scaling factor
42: subw2 r1,r4
43: movl r0,r2
44: subw2 r1,r2
45:
46: # Cubic step
47:
48: clrl r1
49: clrl r3
50: muld2 r0,r2
51: subd2 r2,r4
52: addw2 $0x100,r2
53: addd2 r4,r2
54: muld2 r0,r4
55: divd2 r2,r4
56: addw2 $0x80,r4
57: addd2 r4,r0
58: rsb
59: nonpos:
60: jneq negarg
61: clrd r0 #argument is zero
62: rsb
63: negarg:
64: movl $EDOM,_errno
65: mnegd r4,-(sp)
66: calls $2,_sqrt
67: mnegd r0,r0 # returns -sqrt(-arg)
68: ret
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.