| blob | a6d1fe8fd83b541d2dedb681454955a1a32918c3 |
1 /* $NetBSD: n_atan2.S,v 1.9 2014/10/10 20:58:09 martin 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 WEAK_ALIAS(_atan2l, _atan2)
88 #endif
90 ENTRY(_atan2, 0x0fc0)
91 movq 4(%ap),%r2 # %r2 = y
92 movq 12(%ap),%r4 # %r4 = x
93 bicw3 $0x7f,%r2,%r0
94 bicw3 $0x7f,%r4,%r1
95 cmpw %r0,$0x8000 # y is the reserved operand
96 jeql resop
97 cmpw %r1,$0x8000 # x is the reserved operand
98 jeql resop
99 subl2 $8,%sp
100 bicw3 $0x7fff,%r2,-4(%fp) # copy y sign bit to -4(%fp)
101 bicw3 $0x7fff,%r4,-8(%fp) # copy x sign bit to -8(%fp)
102 cmpd %r4,$0x4080 # x = 1.0 ?
103 bneq xnot1
104 movq %r2,%r0
105 bicw2 $0x8000,%r0 # t = |y|
106 movq %r0,%r2 # y = |y|
107 jbr begin
108 xnot1:
109 bicw3 $0x807f,%r2,%r11 # yexp
110 jeql yeq0 # if y=0 goto yeq0
111 bicw3 $0x807f,%r4,%r10 # xexp
112 jeql pio2 # if x=0 goto pio2
113 subw2 %r10,%r11 # k = yexp - xexp
114 cmpw %r11,$0x2000 # k >= 64 (exp) ?
115 jgeq pio2 # atan2 = +-pi/2
116 divd3 %r4,%r2,%r0 # t = y/x never overflow
117 bicw2 $0x8000,%r0 # t > 0
118 bicw2 $0xff80,%r2 # clear the exponent of y
119 bicw2 $0xff80,%r4 # clear the exponent of x
120 bisw2 $0x4080,%r2 # normalize y to [1,2)
121 bisw2 $0x4080,%r4 # normalize x to [1,2)
122 subw2 %r11,%r4 # scale x so that yexp-xexp=k
123 begin:
124 cmpw %r0,$0x411c # t : 39/16
125 jgeq L50
126 addl3 $0x180,%r0,%r10 # 8*t
127 cvtrfl %r10,%r10 # [8*t] rounded to int
128 ashl $-1,%r10,%r10 # [8*t]/2
129 casel %r10,$0,$4
130 L1:
131 .word L20-L1
132 .word L20-L1
133 .word L30-L1
134 .word L40-L1
135 .word L40-L1
136 L10:
137 movq $0xb4d9940f985e407b,%r6 # Hi=.98279372324732906796d0
138 movq $0x21b1879a3bc2a2fc,%r8 # Lo=-.17092002525602665777d-17
139 subd3 %r4,%r2,%r0 # y-x
140 addw2 $0x80,%r0 # 2(y-x)
141 subd2 %r4,%r0 # 2(y-x)-x
142 addw2 $0x80,%r4 # 2x
143 movq %r2,%r10
144 addw2 $0x80,%r10 # 2y
145 addd2 %r10,%r2 # 3y
146 addd2 %r4,%r2 # 3y+2x
147 divd2 %r2,%r0 # (2y-3x)/(2x+3y)
148 jbr L60
149 L20:
150 cmpw %r0,$0x3280 # t : 2**(-28)
151 jlss L80
152 clrq %r6 # Hi=%r6=0, Lo=%r8=0
153 clrq %r8
154 jbr L60
155 L30:
156 movq $0xda7b2b0d63383fed,%r6 # Hi=.46364760900080611433d0
157 movq $0xf0ea17b2bf912295,%r8 # Lo=.10147340032515978826d-17
158 movq %r2,%r0
159 addw2 $0x80,%r0 # 2y
160 subd2 %r4,%r0 # 2y-x
161 addw2 $0x80,%r4 # 2x
162 addd2 %r2,%r4 # 2x+y
163 divd2 %r4,%r0 # (2y-x)/(2x+y)
164 jbr L60
165 L50:
166 movq $0x68c2a2210fda40c9,%r6 # Hi=1.5707963267948966135d1
167 movq $0x06e0145c26332326,%r8 # Lo=.22517417741562176079d-17
168 cmpw %r0,$0x5100 # y : 2**57
169 bgeq L90
170 divd3 %r2,%r4,%r0
171 bisw2 $0x8000,%r0 # -x/y
172 jbr L60
173 L40:
174 movq $0x68c2a2210fda4049,%r6 # Hi=.78539816339744830676d0
175 movq $0x06e0145c263322a6,%r8 # Lo=.11258708870781088040d-17
176 subd3 %r4,%r2,%r0 # y-x
177 addd2 %r4,%r2 # y+x
178 divd2 %r2,%r0 # (y-x)/(y+x)
179 L60:
180 movq %r0,%r10
181 muld2 %r0,%r0
182 polyd %r0,$12,ptable
183 muld2 %r10,%r0
184 subd2 %r0,%r8
185 addd3 %r8,%r10,%r0
186 addd2 %r6,%r0
187 L80:
188 movw -8(%fp),%r2
189 bneq pim
190 bisw2 -4(%fp),%r0 # return sign(y)*%r0
191 ret
192 L90: # x >= 2**25
193 movq %r6,%r0
194 jbr L80
195 pim:
196 subd3 %r0,$0x68c2a2210fda4149,%r0 # pi-t
197 bisw2 -4(%fp),%r0
198 ret
199 yeq0:
200 movw -8(%fp),%r2
201 beql zero # if sign(x)=1 return pi
202 movq $0x68c2a2210fda4149,%r0 # pi=3.1415926535897932270d1
203 ret
204 zero:
205 clrq %r0 # return 0
206 ret
207 pio2:
208 movq $0x68c2a2210fda40c9,%r0 # pi/2=1.5707963267948966135d1
209 bisw2 -4(%fp),%r0 # return sign(y)*pi/2
210 ret
211 resop:
212 movq $0x8000,%r0 # propagate the reserved operand
213 ret
215 _ALIGN_TEXT
216 ptable:
217 .quad 0xb50f5ce96e7abd60
218 .quad 0x51e44a42c1073e02
219 .quad 0x3487e3289643be35
220 .quad 0xdb62066dffba3e54
221 .quad 0xcf8e2d5199abbe70
222 .quad 0x26f39cb884883e88
223 .quad 0x135117d18998be9d
224 .quad 0x602ce9742e883eba
225 .quad 0xa35ad0be8e38bee3
226 .quad 0xffac922249243f12
227 .quad 0x7f14ccccccccbf4c
228 .quad 0xaa8faaaaaaaa3faa
229 .quad 0x0000000000000000