|
|
1.1 root 1: #
2: # Copyright (c) 1985 Regents of the University of California.
3: #
4: # Use and reproduction of this software are granted in accordance with
5: # the terms and conditions specified in the Berkeley Software License
6: # Agreement (in particular, this entails acknowledgement of the programs'
7: # source, and inclusion of this notice) with the additional understanding
8: # that all recipients should regard themselves as participants in an
9: # ongoing research project and hence should feel obligated to report
10: # their experiences (good or bad) with these elementary function codes,
11: # using "sendbug 4bsd-bugs@BERKELEY", to the authors.
12: #
13:
14: # @(#)cabs.s 1.2 (Berkeley) 8/21/85
15:
16: # double precision complex absolute value
17: # CABS by W. Kahan, 9/7/80.
18: # Revised for reserved operands by E. LeBlanc, 8/18/82
19: # argument for complex absolute value by reference, *4(ap)
20: # argument for cabs and hypot (C fcns) by value, 4(ap)
21: # output is in r0:r1 (error less than 0.86 ulps)
22:
23: .text
24: .align 1
25: .globl _cabs
26: .globl _hypot
27: .globl _z_abs
28: .globl libm$cdabs_r6
29: .globl libm$dsqrt_r5
30:
31: # entry for c functions cabs and hypot
32: _cabs:
33: _hypot:
34: .word 0x807c # save r2-r6, enable floating overflow
35: movq 4(ap),r0 # r0:1 = x
36: movq 12(ap),r2 # r2:3 = y
37: jmp cabs2
38: # entry for Fortran use, call by: d = abs(z)
39: _z_abs:
40: .word 0x807c # save r2-r6, enable floating overflow
41: movl 4(ap),r2 # indirect addressing is necessary here
42: movq (r2)+,r0 # r0:1 = x
43: movq (r2),r2 # r2:3 = y
44:
45: cabs2:
46: bicw3 $0x7f,r0,r4 # r4 has signed biased exp of x
47: cmpw $0x8000,r4
48: jeql return # x is a reserved operand, so return it
49: bicw3 $0x7f,r2,r5 # r5 has signed biased exp of y
50: cmpw $0x8000,r5
51: jneq cont # y isn't a reserved operand
52: movq r2,r0 # return y if it's reserved
53: ret
54:
55: cont:
56: bsbb regs_set # r0:1 = dsqrt(x^2+y^2)/2^r6
57: addw2 r6,r0 # unscaled cdabs in r0:1
58: jvc return # unless it overflows
59: subw2 $0x80,r0 # halve r0 to get meaningful overflow
60: addd2 r0,r0 # overflow; r0 is half of true abs value
61: return:
62: ret
63:
64: libm$cdabs_r6: # ENTRY POINT for cdsqrt
65: # calculates a scaled (factor in r6)
66: # complex absolute value
67:
68: movq (r4)+,r0 # r0:r1 = x via indirect addressing
69: movq (r4),r2 # r2:r3 = y via indirect addressing
70:
71: bicw3 $0x7f,r0,r5 # r5 has signed biased exp of x
72: cmpw $0x8000,r5
73: jeql cdreserved # x is a reserved operand
74: bicw3 $0x7f,r2,r5 # r5 has signed biased exp of y
75: cmpw $0x8000,r5
76: jneq regs_set # y isn't a reserved operand either?
77:
78: cdreserved:
79: movl *4(ap),r4 # r4 -> (u,v), if x or y is reserved
80: movq r0,(r4)+ # copy u and v as is and return
81: movq r2,(r4) # (again addressing is indirect)
82: ret
83:
84: regs_set:
85: bicw2 $0x8000,r0 # r0:r1 = dabs(x)
86: bicw2 $0x8000,r2 # r2:r3 = dabs(y)
87: cmpw r0,r2
88: jgeq ordered
89: movq r0,r4
90: movq r2,r0
91: movq r4,r2 # force y's exp <= x's exp
92: ordered:
93: bicw3 $0x7f,r0,r6 # r6 = exponent(x) + bias(129)
94: jeql retsb # if x = y = 0 then cdabs(x,y) = 0
95: subw2 $0x4780,r6 # r6 = exponent(x) - 14
96: subw2 r6,r0 # 2^14 <= scaled x < 2^15
97: bitw $0xff80,r2
98: jeql retsb # if y = 0 return dabs(x)
99: subw2 r6,r2
100: cmpw $0x3780,r2 # if scaled y < 2^-18
101: jgtr retsb # return dabs(x)
102: emodd r0,$0,r0,r4,r0 # r4 + r0:1 = scaled x^2
103: emodd r2,$0,r2,r5,r2 # r5 + r2:3 = scaled y^2
104: addd2 r2,r0
105: addl2 r5,r4
106: cvtld r4,r2
107: addd2 r2,r0 # r0:1 = scaled x^2 + y^2
108: jmp libm$dsqrt_r5 # r0:1 = dsqrt(x^2+y^2)/2^r6
109: retsb:
110: rsb # error < 0.86 ulp
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.