![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to atan2.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 atan2.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: # @(#)atan2.s 5.3 (Berkeley) 6/30/88
22: #
23: .data
24: .align 2
25: _sccsid:
26: .asciz "@(#)atan2.s 1.2 (Berkeley) 8/21/85; 5.3 (ucb.elefunt) 6/30/88"
27:
28: # ATAN2(Y,X)
29: # RETURN ARG (X+iY)
30: # VAX D FORMAT (56 BITS PRECISION)
31: # CODED IN VAX ASSEMBLY LANGUAGE BY K.C. NG, 4/16/85;
32: #
33: #
34: # Method :
35: # 1. Reduce y to positive by atan2(y,x)=-atan2(-y,x).
36: # 2. Reduce x to positive by (if x and y are unexceptional):
37: # ARG (x+iy) = arctan(y/x) ... if x > 0,
38: # ARG (x+iy) = pi - arctan[y/(-x)] ... if x < 0,
39: # 3. According to the integer k=4t+0.25 truncated , t=y/x, the argument
40: # is further reduced to one of the following intervals and the
41: # arctangent of y/x is evaluated by the corresponding formula:
42: #
43: # [0,7/16] atan(y/x) = t - t^3*(a1+t^2*(a2+...(a10+t^2*a11)...)
44: # [7/16,11/16] atan(y/x) = atan(1/2) + atan( (y-x/2)/(x+y/2) )
45: # [11/16.19/16] atan(y/x) = atan( 1 ) + atan( (y-x)/(x+y) )
46: # [19/16,39/16] atan(y/x) = atan(3/2) + atan( (y-1.5x)/(x+1.5y) )
47: # [39/16,INF] atan(y/x) = atan(INF) + atan( -x/y )
48: #
49: # Special cases:
50: # Notations: atan2(y,x) == ARG (x+iy) == ARG(x,y).
51: #
52: # ARG( NAN , (anything) ) is NaN;
53: # ARG( (anything), NaN ) is NaN;
54: # ARG(+(anything but NaN), +-0) is +-0 ;
55: # ARG(-(anything but NaN), +-0) is +-PI ;
56: # ARG( 0, +-(anything but 0 and NaN) ) is +-PI/2;
57: # ARG( +INF,+-(anything but INF and NaN) ) is +-0 ;
58: # ARG( -INF,+-(anything but INF and NaN) ) is +-PI;
59: # ARG( +INF,+-INF ) is +-PI/4 ;
60: # ARG( -INF,+-INF ) is +-3PI/4;
61: # ARG( (anything but,0,NaN, and INF),+-INF ) is +-PI/2;
62: #
63: # Accuracy:
64: # atan2(y,x) returns the exact ARG(x+iy) nearly rounded.
65: #
66: .text
67: .align 1
68: .globl _atan2
69: _atan2 :
70: .word 0x0ff4
71: movq 4(ap),r2 # r2 = y
72: movq 12(ap),r4 # r4 = x
73: bicw3 $0x7f,r2,r0
74: bicw3 $0x7f,r4,r1
75: cmpw r0,$0x8000 # y is the reserved operand
76: jeql resop
77: cmpw r1,$0x8000 # x is the reserved operand
78: jeql resop
79: subl2 $8,sp
80: bicw3 $0x7fff,r2,-4(fp) # copy y sign bit to -4(fp)
81: bicw3 $0x7fff,r4,-8(fp) # copy x sign bit to -8(fp)
82: cmpd r4,$0x4080 # x = 1.0 ?
83: bneq xnot1
84: movq r2,r0
85: bicw2 $0x8000,r0 # t = |y|
86: movq r0,r2 # y = |y|
87: brb begin
88: xnot1:
89: bicw3 $0x807f,r2,r11 # yexp
90: jeql yeq0 # if y=0 goto yeq0
91: bicw3 $0x807f,r4,r10 # xexp
92: jeql pio2 # if x=0 goto pio2
93: subw2 r10,r11 # k = yexp - xexp
94: cmpw r11,$0x2000 # k >= 64 (exp) ?
95: jgeq pio2 # atan2 = +-pi/2
96: divd3 r4,r2,r0 # t = y/x never overflow
97: bicw2 $0x8000,r0 # t > 0
98: bicw2 $0xff80,r2 # clear the exponent of y
99: bicw2 $0xff80,r4 # clear the exponent of x
100: bisw2 $0x4080,r2 # normalize y to [1,2)
101: bisw2 $0x4080,r4 # normalize x to [1,2)
102: subw2 r11,r4 # scale x so that yexp-xexp=k
103: begin:
104: cmpw r0,$0x411c # t : 39/16
105: jgeq L50
106: addl3 $0x180,r0,r10 # 8*t
107: cvtrfl r10,r10 # [8*t] rounded to int
108: ashl $-1,r10,r10 # [8*t]/2
109: casel r10,$0,$4
110: L1:
111: .word L20-L1
112: .word L20-L1
113: .word L30-L1
114: .word L40-L1
115: .word L40-L1
116: L10:
117: movq $0xb4d9940f985e407b,r6 # Hi=.98279372324732906796d0
118: movq $0x21b1879a3bc2a2fc,r8 # Lo=-.17092002525602665777d-17
119: subd3 r4,r2,r0 # y-x
120: addw2 $0x80,r0 # 2(y-x)
121: subd2 r4,r0 # 2(y-x)-x
122: addw2 $0x80,r4 # 2x
123: movq r2,r10
124: addw2 $0x80,r10 # 2y
125: addd2 r10,r2 # 3y
126: addd2 r4,r2 # 3y+2x
127: divd2 r2,r0 # (2y-3x)/(2x+3y)
128: brw L60
129: L20:
130: cmpw r0,$0x3280 # t : 2**(-28)
131: jlss L80
132: clrq r6 # Hi=r6=0, Lo=r8=0
133: clrq r8
134: brw L60
135: L30:
136: movq $0xda7b2b0d63383fed,r6 # Hi=.46364760900080611433d0
137: movq $0xf0ea17b2bf912295,r8 # Lo=.10147340032515978826d-17
138: movq r2,r0
139: addw2 $0x80,r0 # 2y
140: subd2 r4,r0 # 2y-x
141: addw2 $0x80,r4 # 2x
142: addd2 r2,r4 # 2x+y
143: divd2 r4,r0 # (2y-x)/(2x+y)
144: brb L60
145: L50:
146: movq $0x68c2a2210fda40c9,r6 # Hi=1.5707963267948966135d1
147: movq $0x06e0145c26332326,r8 # Lo=.22517417741562176079d-17
148: cmpw r0,$0x5100 # y : 2**57
149: bgeq L90
150: divd3 r2,r4,r0
151: bisw2 $0x8000,r0 # -x/y
152: brb L60
153: L40:
154: movq $0x68c2a2210fda4049,r6 # Hi=.78539816339744830676d0
155: movq $0x06e0145c263322a6,r8 # Lo=.11258708870781088040d-17
156: subd3 r4,r2,r0 # y-x
157: addd2 r4,r2 # y+x
158: divd2 r2,r0 # (y-x)/(y+x)
159: L60:
160: movq r0,r10
161: muld2 r0,r0
162: polyd r0,$12,ptable
163: muld2 r10,r0
164: subd2 r0,r8
165: addd3 r8,r10,r0
166: addd2 r6,r0
167: L80:
168: movw -8(fp),r2
169: bneq pim
170: bisw2 -4(fp),r0 # return sign(y)*r0
171: ret
172: L90: # x >= 2**25
173: movq r6,r0
174: brb L80
175: pim:
176: subd3 r0,$0x68c2a2210fda4149,r0 # pi-t
177: bisw2 -4(fp),r0
178: ret
179: yeq0:
180: movw -8(fp),r2
181: beql zero # if sign(x)=1 return pi
182: movq $0x68c2a2210fda4149,r0 # pi=3.1415926535897932270d1
183: ret
184: zero:
185: clrq r0 # return 0
186: ret
187: pio2:
188: movq $0x68c2a2210fda40c9,r0 # pi/2=1.5707963267948966135d1
189: bisw2 -4(fp),r0 # return sign(y)*pi/2
190: ret
191: resop:
192: movq $0x8000,r0 # propagate the reserved operand
193: ret
194: .align 2
195: ptable:
196: .quad 0xb50f5ce96e7abd60
197: .quad 0x51e44a42c1073e02
198: .quad 0x3487e3289643be35
199: .quad 0xdb62066dffba3e54
200: .quad 0xcf8e2d5199abbe70
201: .quad 0x26f39cb884883e88
202: .quad 0x135117d18998be9d
203: .quad 0x602ce9742e883eba
204: .quad 0xa35ad0be8e38bee3
205: .quad 0xffac922249243f12
206: .quad 0x7f14ccccccccbf4c
207: .quad 0xaa8faaaaaaaa3faa
208: .quad 0x0000000000000000
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.