Commit | Line | Data |
---|---|---|
24a603be | 1 | # Copyright (c) 1985 Regents of the University of California. |
fe5be67b KB |
2 | # All rights reserved. |
3 | # | |
4 | # Redistribution and use in source and binary forms are permitted | |
5 | # provided that this notice is preserved and that due credit is given | |
6 | # to the University of California at Berkeley. The name of the University | |
7 | # may not be used to endorse or promote products derived from this | |
8 | # software without specific prior written permission. This software | |
9 | # is provided ``as is'' without express or implied warranty. | |
10 | # | |
11 | # All recipients should regard themselves as participants in an ongoing | |
12 | # research project and hence should feel obligated to report their | |
13 | # experiences (good or bad) with these elementary function codes, using | |
14 | # the sendbug(8) program, to the authors. | |
15 | # | |
16 | # @(#)cabs.s 5.2 (Berkeley) %G% | |
24a603be | 17 | # |
7c0a3811 GK |
18 | .data |
19 | .align 2 | |
20 | _sccsid: | |
fe5be67b | 21 | .asciz "@(#)cabs.s 1.2 (Berkeley) 8/21/85; 5.2 (ucb.elefunt) %G%" |
24a603be ZAL |
22 | |
23 | # double precision complex absolute value | |
24 | # CABS by W. Kahan, 9/7/80. | |
25 | # Revised for reserved operands by E. LeBlanc, 8/18/82 | |
26 | # argument for complex absolute value by reference, *4(ap) | |
27 | # argument for cabs and hypot (C fcns) by value, 4(ap) | |
28 | # output is in r0:r1 (error less than 0.86 ulps) | |
29 | ||
30 | .text | |
31 | .align 1 | |
32 | .globl _cabs | |
33 | .globl _hypot | |
34 | .globl _z_abs | |
35 | .globl libm$cdabs_r6 | |
36 | .globl libm$dsqrt_r5 | |
37 | ||
38 | # entry for c functions cabs and hypot | |
39 | _cabs: | |
40 | _hypot: | |
41 | .word 0x807c # save r2-r6, enable floating overflow | |
42 | movq 4(ap),r0 # r0:1 = x | |
43 | movq 12(ap),r2 # r2:3 = y | |
44 | jmp cabs2 | |
45 | # entry for Fortran use, call by: d = abs(z) | |
46 | _z_abs: | |
47 | .word 0x807c # save r2-r6, enable floating overflow | |
48 | movl 4(ap),r2 # indirect addressing is necessary here | |
49 | movq (r2)+,r0 # r0:1 = x | |
50 | movq (r2),r2 # r2:3 = y | |
51 | ||
52 | cabs2: | |
53 | bicw3 $0x7f,r0,r4 # r4 has signed biased exp of x | |
54 | cmpw $0x8000,r4 | |
55 | jeql return # x is a reserved operand, so return it | |
56 | bicw3 $0x7f,r2,r5 # r5 has signed biased exp of y | |
57 | cmpw $0x8000,r5 | |
58 | jneq cont # y isn't a reserved operand | |
59 | movq r2,r0 # return y if it's reserved | |
60 | ret | |
61 | ||
62 | cont: | |
63 | bsbb regs_set # r0:1 = dsqrt(x^2+y^2)/2^r6 | |
64 | addw2 r6,r0 # unscaled cdabs in r0:1 | |
65 | jvc return # unless it overflows | |
66 | subw2 $0x80,r0 # halve r0 to get meaningful overflow | |
67 | addd2 r0,r0 # overflow; r0 is half of true abs value | |
68 | return: | |
69 | ret | |
70 | ||
71 | libm$cdabs_r6: # ENTRY POINT for cdsqrt | |
72 | # calculates a scaled (factor in r6) | |
73 | # complex absolute value | |
74 | ||
75 | movq (r4)+,r0 # r0:r1 = x via indirect addressing | |
76 | movq (r4),r2 # r2:r3 = y via indirect addressing | |
77 | ||
78 | bicw3 $0x7f,r0,r5 # r5 has signed biased exp of x | |
79 | cmpw $0x8000,r5 | |
80 | jeql cdreserved # x is a reserved operand | |
81 | bicw3 $0x7f,r2,r5 # r5 has signed biased exp of y | |
82 | cmpw $0x8000,r5 | |
83 | jneq regs_set # y isn't a reserved operand either? | |
84 | ||
85 | cdreserved: | |
86 | movl *4(ap),r4 # r4 -> (u,v), if x or y is reserved | |
87 | movq r0,(r4)+ # copy u and v as is and return | |
88 | movq r2,(r4) # (again addressing is indirect) | |
89 | ret | |
90 | ||
91 | regs_set: | |
92 | bicw2 $0x8000,r0 # r0:r1 = dabs(x) | |
93 | bicw2 $0x8000,r2 # r2:r3 = dabs(y) | |
94 | cmpw r0,r2 | |
95 | jgeq ordered | |
96 | movq r0,r4 | |
97 | movq r2,r0 | |
98 | movq r4,r2 # force y's exp <= x's exp | |
99 | ordered: | |
100 | bicw3 $0x7f,r0,r6 # r6 = exponent(x) + bias(129) | |
101 | jeql retsb # if x = y = 0 then cdabs(x,y) = 0 | |
102 | subw2 $0x4780,r6 # r6 = exponent(x) - 14 | |
103 | subw2 r6,r0 # 2^14 <= scaled x < 2^15 | |
104 | bitw $0xff80,r2 | |
105 | jeql retsb # if y = 0 return dabs(x) | |
106 | subw2 r6,r2 | |
107 | cmpw $0x3780,r2 # if scaled y < 2^-18 | |
108 | jgtr retsb # return dabs(x) | |
109 | emodd r0,$0,r0,r4,r0 # r4 + r0:1 = scaled x^2 | |
110 | emodd r2,$0,r2,r5,r2 # r5 + r2:3 = scaled y^2 | |
111 | addd2 r2,r0 | |
112 | addl2 r5,r4 | |
113 | cvtld r4,r2 | |
114 | addd2 r2,r0 # r0:1 = scaled x^2 + y^2 | |
115 | jmp libm$dsqrt_r5 # r0:1 = dsqrt(x^2+y^2)/2^r6 | |
116 | retsb: | |
117 | rsb # error < 0.86 ulp |