| blob | accb40df2d600cae23fb83a33d1a181f7a34d0ce |
1 /* $NetBSD: n_atan2.S,v 1.7 2008/03/20 16:41:26 mhitch Exp $ */
2 /*
3 * Copyright (c) 1985, 1993
4 * The Regents of the University of California. All rights reserved.
5 *
6 * Redistribution and use in source and binary forms, with or without
7 * modification, are permitted provided that the following conditions
8 * are met:
9 * 1. Redistributions of source code must retain the above copyright
10 * notice, this list of conditions and the following disclaimer.
11 * 2. Redistributions in binary form must reproduce the above copyright
12 * notice, this list of conditions and the following disclaimer in the
13 * documentation and/or other materials provided with the distribution.
14 * 3. Neither the name of the University nor the names of its contributors
15 * may be used to endorse or promote products derived from this software
16 * without specific prior written permission.
17 *
18 * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
19 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
20 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
21 * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
22 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
23 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
24 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
25 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
26 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
27 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
28 * SUCH DAMAGE.
29 *
30 * @(#)atan2.s 8.1 (Berkeley) 6/4/93
31 */
33 #include <machine/asm.h>
35 /*
36 * ATAN2(Y,X)
37 * RETURN ARG (X+iY)
38 * VAX D FORMAT (56 BITS PRECISION)
39 * CODED IN VAX ASSEMBLY LANGUAGE BY K.C. NG, 4/16/85;
40 *
41 *
42 * Method :
43 * 1. Reduce y to positive by atan2(y,x)=-atan2(-y,x).
44 * 2. Reduce x to positive by (if x and y are unexceptional):
45 * ARG (x+iy) = arctan(y/x) ... if x > 0,
46 * ARG (x+iy) = pi - arctan[y/(-x)] ... if x < 0,
47 * 3. According to the integer k=4t+0.25 truncated , t=y/x, the argument
48 * is further reduced to one of the following intervals and the
49 * arctangent of y/x is evaluated by the corresponding formula:
50 *
51 * [0,7/16] atan(y/x) = t - t^3*(a1+t^2*(a2+...(a10+t^2*a11)...)
52 * [7/16,11/16] atan(y/x) = atan(1/2) + atan( (y-x/2)/(x+y/2) )
53 * [11/16.19/16] atan(y/x) = atan( 1 ) + atan( (y-x)/(x+y) )
54 * [19/16,39/16] atan(y/x) = atan(3/2) + atan( (y-1.5x)/(x+1.5y) )
55 * [39/16,INF] atan(y/x) = atan(INF) + atan( -x/y )
56 *
57 * Special cases:
58 * Notations: atan2(y,x) == ARG (x+iy) == ARG(x,y).
59 *
60 * ARG( NAN , (anything) ) is NaN;
61 * ARG( (anything), NaN ) is NaN;
62 * ARG(+(anything but NaN), +-0) is +-0 ;
63 * ARG(-(anything but NaN), +-0) is +-PI ;
64 * ARG( 0, +-(anything but 0 and NaN) ) is +-PI/2;
65 * ARG( +INF,+-(anything but INF and NaN) ) is +-0 ;
66 * ARG( -INF,+-(anything but INF and NaN) ) is +-PI;
67 * ARG( +INF,+-INF ) is +-PI/4 ;
68 * ARG( -INF,+-INF ) is +-3PI/4;
69 * ARG( (anything but,0,NaN, and INF),+-INF ) is +-PI/2;
70 *
71 * Accuracy:
72 * atan2(y,x) returns the exact ARG(x+iy) nearly rounded.
73 */
75 #ifdef WEAK_ALIAS
76 WEAK_ALIAS(atan2f, _atan2f)
77 #endif
79 ENTRY(_atan2f, 0)
80 cvtfd 4(%ap),-(%sp)
81 calls $2,_C_LABEL(_atan2)
82 cvtdf %r0,%r0
83 ret
85 #ifdef WEAK_ALIAS
86 WEAK_ALIAS(atan2, _atan2)
87 #endif
89 ENTRY(_atan2, 0x0fc0)
90 movq 4(%ap),%r2 # %r2 = y
91 movq 12(%ap),%r4 # %r4 = x
92 bicw3 $0x7f,%r2,%r0
93 bicw3 $0x7f,%r4,%r1
94 cmpw %r0,$0x8000 # y is the reserved operand
95 jeql resop
96 cmpw %r1,$0x8000 # x is the reserved operand
97 jeql resop
98 subl2 $8,%sp
99 bicw3 $0x7fff,%r2,-4(%fp) # copy y sign bit to -4(%fp)
100 bicw3 $0x7fff,%r4,-8(%fp) # copy x sign bit to -8(%fp)
101 cmpd %r4,$0x4080 # x = 1.0 ?
102 bneq xnot1
103 movq %r2,%r0
104 bicw2 $0x8000,%r0 # t = |y|
105 movq %r0,%r2 # y = |y|
106 jbr begin
107 xnot1:
108 bicw3 $0x807f,%r2,%r11 # yexp
109 jeql yeq0 # if y=0 goto yeq0
110 bicw3 $0x807f,%r4,%r10 # xexp
111 jeql pio2 # if x=0 goto pio2
112 subw2 %r10,%r11 # k = yexp - xexp
113 cmpw %r11,$0x2000 # k >= 64 (exp) ?
114 jgeq pio2 # atan2 = +-pi/2
115 divd3 %r4,%r2,%r0 # t = y/x never overflow
116 bicw2 $0x8000,%r0 # t > 0
117 bicw2 $0xff80,%r2 # clear the exponent of y
118 bicw2 $0xff80,%r4 # clear the exponent of x
119 bisw2 $0x4080,%r2 # normalize y to [1,2)
120 bisw2 $0x4080,%r4 # normalize x to [1,2)
121 subw2 %r11,%r4 # scale x so that yexp-xexp=k
122 begin:
123 cmpw %r0,$0x411c # t : 39/16
124 jgeq L50
125 addl3 $0x180,%r0,%r10 # 8*t
126 cvtrfl %r10,%r10 # [8*t] rounded to int
127 ashl $-1,%r10,%r10 # [8*t]/2
128 casel %r10,$0,$4
129 L1:
130 .word L20-L1
131 .word L20-L1
132 .word L30-L1
133 .word L40-L1
134 .word L40-L1
135 L10:
136 movq $0xb4d9940f985e407b,%r6 # Hi=.98279372324732906796d0
137 movq $0x21b1879a3bc2a2fc,%r8 # Lo=-.17092002525602665777d-17
138 subd3 %r4,%r2,%r0 # y-x
139 addw2 $0x80,%r0 # 2(y-x)
140 subd2 %r4,%r0 # 2(y-x)-x
141 addw2 $0x80,%r4 # 2x
142 movq %r2,%r10
143 addw2 $0x80,%r10 # 2y
144 addd2 %r10,%r2 # 3y
145 addd2 %r4,%r2 # 3y+2x
146 divd2 %r2,%r0 # (2y-3x)/(2x+3y)
147 jbr L60
148 L20:
149 cmpw %r0,$0x3280 # t : 2**(-28)
150 jlss L80
151 clrq %r6 # Hi=%r6=0, Lo=%r8=0
152 clrq %r8
153 jbr L60
154 L30:
155 movq $0xda7b2b0d63383fed,%r6 # Hi=.46364760900080611433d0
156 movq $0xf0ea17b2bf912295,%r8 # Lo=.10147340032515978826d-17
157 movq %r2,%r0
158 addw2 $0x80,%r0 # 2y
159 subd2 %r4,%r0 # 2y-x
160 addw2 $0x80,%r4 # 2x
161 addd2 %r2,%r4 # 2x+y
162 divd2 %r4,%r0 # (2y-x)/(2x+y)
163 jbr L60
164 L50:
165 movq $0x68c2a2210fda40c9,%r6 # Hi=1.5707963267948966135d1
166 movq $0x06e0145c26332326,%r8 # Lo=.22517417741562176079d-17
167 cmpw %r0,$0x5100 # y : 2**57
168 bgeq L90
169 divd3 %r2,%r4,%r0
170 bisw2 $0x8000,%r0 # -x/y
171 jbr L60
172 L40:
173 movq $0x68c2a2210fda4049,%r6 # Hi=.78539816339744830676d0
174 movq $0x06e0145c263322a6,%r8 # Lo=.11258708870781088040d-17
175 subd3 %r4,%r2,%r0 # y-x
176 addd2 %r4,%r2 # y+x
177 divd2 %r2,%r0 # (y-x)/(y+x)
178 L60:
179 movq %r0,%r10
180 muld2 %r0,%r0
181 polyd %r0,$12,ptable
182 muld2 %r10,%r0
183 subd2 %r0,%r8
184 addd3 %r8,%r10,%r0
185 addd2 %r6,%r0
186 L80:
187 movw -8(%fp),%r2
188 bneq pim
189 bisw2 -4(%fp),%r0 # return sign(y)*%r0
190 ret
191 L90: # x >= 2**25
192 movq %r6,%r0
193 jbr L80
194 pim:
195 subd3 %r0,$0x68c2a2210fda4149,%r0 # pi-t
196 bisw2 -4(%fp),%r0
197 ret
198 yeq0:
199 movw -8(%fp),%r2
200 beql zero # if sign(x)=1 return pi
201 movq $0x68c2a2210fda4149,%r0 # pi=3.1415926535897932270d1
202 ret
203 zero:
204 clrq %r0 # return 0
205 ret
206 pio2:
207 movq $0x68c2a2210fda40c9,%r0 # pi/2=1.5707963267948966135d1
208 bisw2 -4(%fp),%r0 # return sign(y)*pi/2
209 ret
210 resop:
211 movq $0x8000,%r0 # propagate the reserved operand
212 ret
214 _ALIGN_TEXT
215 ptable:
216 .quad 0xb50f5ce96e7abd60
217 .quad 0x51e44a42c1073e02
218 .quad 0x3487e3289643be35
219 .quad 0xdb62066dffba3e54
220 .quad 0xcf8e2d5199abbe70
221 .quad 0x26f39cb884883e88
222 .quad 0x135117d18998be9d
223 .quad 0x602ce9742e883eba
224 .quad 0xa35ad0be8e38bee3
225 .quad 0xffac922249243f12
226 .quad 0x7f14ccccccccbf4c
227 .quad 0xaa8faaaaaaaa3faa
228 .quad 0x0000000000000000