|
|
1.1 root 1: # Copyright (c) 1987 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.4 (Berkeley) 6/30/88
22: #
23: .data
24: .align 2
25: _sccsid:
26: .asciz "@(#)cabs.s 5.4 5.4 (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(fp)
32: # argument for cabs and hypot (C fcns) by value, 4(fp)
33: # output is in r0:r1
34:
35: .text
36: .align 2
37: .globl _cabs
38: .globl _hypot
39: .globl _z_abs
40:
41: # entry for c functions cabs and hypot
42: _cabs:
43: _hypot:
44: .word 0x807c # save r2-r6, enable floating overflow
45: movl 16(fp),r3
46: movl 12(fp),r2 # r2:3 = y
47: movl 8(fp),r1
48: movl 4(fp),r0 # r0:1 = x
49: brb 1f
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(fp),r4 # indirect addressing is necessary here
54: movl 12(r4),r3 #
55: movl 8(r4),r2 # r2:3 = y
56: movl 4(r4),r1 #
57: movl (r4),r0 # r0:1 = x
58: 1: andl3 $0xff800000,r0,r4 # r4 has signed biased exp of x
59: cmpl $0x80000000,r4
60: beql 2f # x is a reserved operand, so return it
61: andl3 $0xff800000,r2,r5 # r5 has signed biased exp of y
62: cmpl $0x80000000,r5
63: bneq 3f # y isn't a reserved operand
64: movl r3,r1
65: movl r2,r0 # return y if it's reserved
66: 2: ret
67:
68: 3: callf $4,regs_set # r0:1 = dsqrt(x^2+y^2)/2^r6
69: addl2 r6,r0 # unscaled cdabs in r0:1
70: jvc 2b # unless it overflows
71: subl2 $0x800000,r0 # halve r0 to get meaningful overflow
72: ldd r0
73: addd r0 # overflow; r0 is half of true abs value
74: ret
75:
76: regs_set:
77: .word 0x0000
78: andl2 $0x7fffffff,r0 # r0:r1 = dabs(x)
79: andl2 $0x7fffffff,r2 # r2:r3 = dabs(y)
80: cmpl r0,r2
81: bgeq 4f
82: movl r1,r5
83: movl r0,r4
84: movl r3,r1
85: movl r2,r0
86: movl r5,r3
87: movl r4,r2 # force y's exp <= x's exp
88: 4: andl3 $0xff800000,r0,r6 # r6 = exponent(x) + bias(129)
89: beql 5f # if x = y = 0 then cdabs(x,y) = 0
90: subl2 $0x47800000,r6 # r6 = exponent(x) - 14
91: subl2 r6,r0 # 2^14 <= scaled x < 2^15
92: bitl $0xff800000,r2
93: beql 5f # if y = 0 return dabs(x)
94: subl2 r6,r2
95: cmpl $0x37800000,r2 # if scaled y < 2^-18
96: bgtr 5f # return dabs(x)
97: ldd r0
98: muld r0
99: std r0 # r0:1 = scaled x^2
100: ldd r2
101: muld r2 # acc = scaled y^2
102: addd r0
103: std r0
104: pushl r1
105: pushl r0
106: callf $12,_sqrt # r0:1 = dsqrt(x^2+y^2)/2^r6
107: 5: ret
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.