|
|
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: bispsw $0xe0
16: movd 4(ap),r0
17: bsbb dsqrt_r5
18: ret
19: dsqrt_r5:
20: movd r0,r4
21: jleq nonpos #argument is not positive
22: movzwl r4,r2
23: ashl $-1,r2,r0
24: addw2 $0x203c,r0 #r0 has magic initial appx
25:
26: # Do two steps of Heron's rule
27:
28: divf3 r0,r4,r2
29: addf2 r2,r0
30: subw2 $0x80,r0
31:
32: divf3 r0,r4,r2
33: addf2 r2,r0
34: subw2 $0x80,r0
35:
36:
37: # Scale argument and approximation to prevent over/underflow
38: # NOTE: The following four steps would not be necessary if underflow
39: # were gentle.
40:
41: bicw3 $0xffff807f,r4,r1
42: subw2 $0x4080,r1 # r1 contains scaling factor
43: subw2 r1,r4
44: movl r0,r2
45: subw2 r1,r2
46:
47: # Cubic step
48:
49: clrl r1
50: clrl r3
51: muld2 r0,r2
52: subd2 r2,r4
53: addw2 $0x100,r2
54: addd2 r4,r2
55: muld2 r0,r4
56: divd2 r2,r4
57: addw2 $0x80,r4
58: addd2 r4,r0
59: rsb
60: nonpos:
61: jneq negarg
62: clrd r0 #argument is zero
63: rsb
64: negarg:
65: movl $EDOM,_errno
66: mnegd r4,-(sp)
67: calls $2,_sqrt
68: mnegd r0,r0 # returns -sqrt(-arg)
69: ret
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.