![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to cabs.s CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [CSRG BSD Unix] / 43BSDReno / lib / libm / vax|
Return to cabs.s CVS log | Up to [CSRG BSD Unix] / 43BSDReno / lib / libm / vax |
1.1 root 1: # Copyright (c) 1985 Regents of the University of California.
2: # All rights reserved.
3: #
4: # Redistribution and use in source and binary forms are permitted
5: # provided that the above copyright notice and this paragraph are
6: # duplicated in all such forms and that any documentation,
7: # advertising materials, and other materials related to such
8: # distribution and use acknowledge that the software was developed
9: # by the University of California, Berkeley. The name of the
10: # University may not be used to endorse or promote products derived
11: # from this software without specific prior written permission.
12: # THIS SOFTWARE IS PROVIDED ``AS IS'' AND WITHOUT ANY EXPRESS OR
13: # IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED
14: # WARRANTIES OF MERCHANTIBILITY AND FITNESS FOR A PARTICULAR PURPOSE.
15: #
16: # All recipients should regard themselves as participants in an ongoing
17: # research project and hence should feel obligated to report their
18: # experiences (good or bad) with these elementary function codes, using
19: # the sendbug(8) program, to the authors.
20: #
21: # @(#)cabs.s 5.3 (Berkeley) 6/30/88
22: #
23: .data
24: .align 2
25: _sccsid:
26: .asciz "@(#)cabs.s 1.2 (Berkeley) 8/21/85; 5.3 (ucb.elefunt) 6/30/88"
27:
28: # double precision complex absolute value
29: # CABS by W. Kahan, 9/7/80.
30: # Revised for reserved operands by E. LeBlanc, 8/18/82
31: # argument for complex absolute value by reference, *4(ap)
32: # argument for cabs and hypot (C fcns) by value, 4(ap)
33: # output is in r0:r1 (error less than 0.86 ulps)
34:
35: .text
36: .align 1
37: .globl _cabs
38: .globl _hypot
39: .globl _z_abs
40: .globl libm$cdabs_r6
41: .globl libm$dsqrt_r5
42:
43: # entry for c functions cabs and hypot
44: _cabs:
45: _hypot:
46: .word 0x807c # save r2-r6, enable floating overflow
47: movq 4(ap),r0 # r0:1 = x
48: movq 12(ap),r2 # r2:3 = y
49: jmp cabs2
50: # entry for Fortran use, call by: d = abs(z)
51: _z_abs:
52: .word 0x807c # save r2-r6, enable floating overflow
53: movl 4(ap),r2 # indirect addressing is necessary here
54: movq (r2)+,r0 # r0:1 = x
55: movq (r2),r2 # r2:3 = y
56:
57: cabs2:
58: bicw3 $0x7f,r0,r4 # r4 has signed biased exp of x
59: cmpw $0x8000,r4
60: jeql return # x is a reserved operand, so return it
61: bicw3 $0x7f,r2,r5 # r5 has signed biased exp of y
62: cmpw $0x8000,r5
63: jneq cont # y isn't a reserved operand
64: movq r2,r0 # return y if it's reserved
65: ret
66:
67: cont:
68: bsbb regs_set # r0:1 = dsqrt(x^2+y^2)/2^r6
69: addw2 r6,r0 # unscaled cdabs in r0:1
70: jvc return # unless it overflows
71: subw2 $0x80,r0 # halve r0 to get meaningful overflow
72: addd2 r0,r0 # overflow; r0 is half of true abs value
73: return:
74: ret
75:
76: libm$cdabs_r6: # ENTRY POINT for cdsqrt
77: # calculates a scaled (factor in r6)
78: # complex absolute value
79:
80: movq (r4)+,r0 # r0:r1 = x via indirect addressing
81: movq (r4),r2 # r2:r3 = y via indirect addressing
82:
83: bicw3 $0x7f,r0,r5 # r5 has signed biased exp of x
84: cmpw $0x8000,r5
85: jeql cdreserved # x is a reserved operand
86: bicw3 $0x7f,r2,r5 # r5 has signed biased exp of y
87: cmpw $0x8000,r5
88: jneq regs_set # y isn't a reserved operand either?
89:
90: cdreserved:
91: movl *4(ap),r4 # r4 -> (u,v), if x or y is reserved
92: movq r0,(r4)+ # copy u and v as is and return
93: movq r2,(r4) # (again addressing is indirect)
94: ret
95:
96: regs_set:
97: bicw2 $0x8000,r0 # r0:r1 = dabs(x)
98: bicw2 $0x8000,r2 # r2:r3 = dabs(y)
99: cmpw r0,r2
100: jgeq ordered
101: movq r0,r4
102: movq r2,r0
103: movq r4,r2 # force y's exp <= x's exp
104: ordered:
105: bicw3 $0x7f,r0,r6 # r6 = exponent(x) + bias(129)
106: jeql retsb # if x = y = 0 then cdabs(x,y) = 0
107: subw2 $0x4780,r6 # r6 = exponent(x) - 14
108: subw2 r6,r0 # 2^14 <= scaled x < 2^15
109: bitw $0xff80,r2
110: jeql retsb # if y = 0 return dabs(x)
111: subw2 r6,r2
112: cmpw $0x3780,r2 # if scaled y < 2^-18
113: jgtr retsb # return dabs(x)
114: emodd r0,$0,r0,r4,r0 # r4 + r0:1 = scaled x^2
115: emodd r2,$0,r2,r5,r2 # r5 + r2:3 = scaled y^2
116: addd2 r2,r0
117: addl2 r5,r4
118: cvtld r4,r2
119: addd2 r2,r0 # r0:1 = scaled x^2 + y^2
120: jmp libm$dsqrt_r5 # r0:1 = dsqrt(x^2+y^2)/2^r6
121: retsb:
122: rsb # error < 0.86 ulp
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.