1/*
2 * CDDL HEADER START
3 *
4 * The contents of this file are subject to the terms of the
5 * Common Development and Distribution License (the "License").
6 * You may not use this file except in compliance with the License.
7 *
8 * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
9 * or http://www.opensolaris.org/os/licensing.
10 * See the License for the specific language governing permissions
11 * and limitations under the License.
12 *
13 * When distributing Covered Code, include this CDDL HEADER in each
14 * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
15 * If applicable, add the following below this CDDL HEADER, with the
16 * fields enclosed by brackets "[]" replaced with your own identifying
17 * information: Portions Copyright [yyyy] [name of copyright owner]
18 *
19 * CDDL HEADER END
20 */
21/*
22 * Copyright 2011 Nexenta Systems, Inc.  All rights reserved.
23 */
24/*
25 * Copyright 2005 Sun Microsystems, Inc.  All rights reserved.
26 * Use is subject to license terms.
27 */
28
29#pragma weak __atan2f = atan2f
30
31#include "libm.h"
32
33#if defined(__i386) && !defined(__amd64)
34extern int __swapRP(int);
35#endif
36
37/*
38 * For i = 0, ..., 192, let x[i] be the double precision number whose
39 * high order 32 bits are 0x3f900000 + (i << 16) and whose low order
40 * 32 bits are zero.  Then TBL[i] := atan(x[i]) to double precision.
41 */
42
43static const double TBL[] = {
44	1.56237286204768313e-02,
45	1.66000375562312640e-02,
46	1.75763148444955872e-02,
47	1.85525586258889763e-02,
48	1.95287670414137082e-02,
49	2.05049382324763683e-02,
50	2.14810703409090559e-02,
51	2.24571615089905717e-02,
52	2.34332098794675855e-02,
53	2.44092135955758099e-02,
54	2.53851708010611396e-02,
55	2.63610796402007873e-02,
56	2.73369382578244127e-02,
57	2.83127447993351995e-02,
58	2.92884974107309737e-02,
59	3.02641942386252458e-02,
60	3.12398334302682774e-02,
61	3.31909314971115949e-02,
62	3.51417768027967800e-02,
63	3.70923545503918164e-02,
64	3.90426499551669928e-02,
65	4.09926482452637811e-02,
66	4.29423346623621707e-02,
67	4.48916944623464972e-02,
68	4.68407129159696539e-02,
69	4.87893753095156174e-02,
70	5.07376669454602178e-02,
71	5.26855731431300420e-02,
72	5.46330792393594777e-02,
73	5.65801705891457105e-02,
74	5.85268325663017702e-02,
75	6.04730505641073168e-02,
76	6.24188099959573500e-02,
77	6.63088949198234884e-02,
78	7.01969710718705203e-02,
79	7.40829225490337306e-02,
80	7.79666338315423008e-02,
81	8.18479898030765457e-02,
82	8.57268757707448092e-02,
83	8.96031774848717461e-02,
84	9.34767811585894698e-02,
85	9.73475734872236709e-02,
86	1.01215441667466668e-01,
87	1.05080273416329528e-01,
88	1.08941956989865793e-01,
89	1.12800381201659389e-01,
90	1.16655435441069349e-01,
91	1.20507009691224562e-01,
92	1.24354994546761438e-01,
93	1.32039761614638762e-01,
94	1.39708874289163648e-01,
95	1.47361481088651630e-01,
96	1.54996741923940973e-01,
97	1.62613828597948568e-01,
98	1.70211925285474408e-01,
99	1.77790228992676075e-01,
100	1.85347949995694761e-01,
101	1.92884312257974672e-01,
102	2.00398553825878512e-01,
103	2.07889927202262986e-01,
104	2.15357699697738048e-01,
105	2.22801153759394521e-01,
106	2.30219587276843718e-01,
107	2.37612313865471242e-01,
108	2.44978663126864143e-01,
109	2.59629629408257512e-01,
110	2.74167451119658789e-01,
111	2.88587361894077410e-01,
112	3.02884868374971417e-01,
113	3.17055753209147029e-01,
114	3.31096076704132103e-01,
115	3.45002177207105132e-01,
116	3.58770670270572245e-01,
117	3.72398446676754202e-01,
118	3.85882669398073752e-01,
119	3.99220769575252543e-01,
120	4.12410441597387323e-01,
121	4.25449637370042266e-01,
122	4.38336559857957830e-01,
123	4.51069655988523499e-01,
124	4.63647609000806094e-01,
125	4.88333951056405535e-01,
126	5.12389460310737732e-01,
127	5.35811237960463704e-01,
128	5.58599315343562441e-01,
129	5.80756353567670414e-01,
130	6.02287346134964152e-01,
131	6.23199329934065904e-01,
132	6.43501108793284371e-01,
133	6.63202992706093286e-01,
134	6.82316554874748071e-01,
135	7.00854407884450192e-01,
136	7.18829999621624527e-01,
137	7.36257428981428097e-01,
138	7.53151280962194414e-01,
139	7.69526480405658297e-01,
140	7.85398163397448279e-01,
141	8.15691923316223422e-01,
142	8.44153986113171051e-01,
143	8.70903457075652976e-01,
144	8.96055384571343927e-01,
145	9.19719605350416858e-01,
146	9.42000040379463610e-01,
147	9.62994330680936206e-01,
148	9.82793723247329054e-01,
149	1.00148313569423464e+00,
150	1.01914134426634972e+00,
151	1.03584125300880014e+00,
152	1.05165021254837376e+00,
153	1.06663036531574362e+00,
154	1.08083900054116833e+00,
155	1.09432890732118993e+00,
156	1.10714871779409041e+00,
157	1.13095374397916038e+00,
158	1.15257199721566761e+00,
159	1.17227388112847630e+00,
160	1.19028994968253166e+00,
161	1.20681737028525249e+00,
162	1.22202532321098967e+00,
163	1.23605948947808186e+00,
164	1.24904577239825443e+00,
165	1.26109338225244039e+00,
166	1.27229739520871732e+00,
167	1.28274087974427076e+00,
168	1.29249666778978534e+00,
169	1.30162883400919616e+00,
170	1.31019393504755555e+00,
171	1.31824205101683711e+00,
172	1.32581766366803255e+00,
173	1.33970565959899957e+00,
174	1.35212738092095464e+00,
175	1.36330010035969384e+00,
176	1.37340076694501589e+00,
177	1.38257482149012589e+00,
178	1.39094282700241845e+00,
179	1.39860551227195762e+00,
180	1.40564764938026987e+00,
181	1.41214106460849531e+00,
182	1.41814699839963154e+00,
183	1.42371797140649403e+00,
184	1.42889927219073276e+00,
185	1.43373015248470903e+00,
186	1.43824479449822262e+00,
187	1.44247309910910193e+00,
188	1.44644133224813509e+00,
189	1.45368758222803240e+00,
190	1.46013910562100091e+00,
191	1.46591938806466282e+00,
192	1.47112767430373470e+00,
193	1.47584462045214027e+00,
194	1.48013643959415142e+00,
195	1.48405798811891154e+00,
196	1.48765509490645531e+00,
197	1.49096634108265924e+00,
198	1.49402443552511865e+00,
199	1.49685728913695626e+00,
200	1.49948886200960629e+00,
201	1.50193983749385196e+00,
202	1.50422816301907281e+00,
203	1.50636948736934317e+00,
204	1.50837751679893928e+00,
205	1.51204050407917401e+00,
206	1.51529782154917969e+00,
207	1.51821326518395483e+00,
208	1.52083793107295384e+00,
209	1.52321322351791322e+00,
210	1.52537304737331958e+00,
211	1.52734543140336587e+00,
212	1.52915374769630819e+00,
213	1.53081763967160667e+00,
214	1.53235373677370856e+00,
215	1.53377621092096650e+00,
216	1.53509721411557254e+00,
217	1.53632722579538861e+00,
218	1.53747533091664934e+00,
219	1.53854944435964280e+00,
220	1.53955649336462841e+00,
221	1.54139303859089161e+00,
222	1.54302569020147562e+00,
223	1.54448660954197448e+00,
224	1.54580153317597646e+00,
225	1.54699130060982659e+00,
226	1.54807296595325550e+00,
227	1.54906061995310385e+00,
228	1.54996600675867957e+00,
229	1.55079899282174605e+00,
230	1.55156792769518947e+00,
231	1.55227992472688747e+00,
232	1.55294108165534417e+00,
233	1.55355665560036682e+00,
234	1.55413120308095598e+00,
235	1.55466869295126031e+00,
236	1.55517259817441977e+00,
237};
238
239static const double
240	pio4	=  7.8539816339744827900e-01,
241	pio2	=  1.5707963267948965580e+00,
242	negpi	= -3.1415926535897931160e+00,
243	q1	= -3.3333333333296428046e-01,
244	q2	=  1.9999999186853752618e-01,
245	zero	=  0.0;
246
247static const float two24 = 16777216.0;
248
249float
250atan2f(float fy, float fx)
251{
252	double	a, t, s, dbase;
253	float	x, y, base;
254	int	i, k, hx, hy, ix, iy, sign;
255#if defined(__i386) && !defined(__amd64)
256	int	rp;
257#endif
258
259	iy = *(int *)&fy;
260	ix = *(int *)&fx;
261	hy = iy & ~0x80000000;
262	hx = ix & ~0x80000000;
263
264	sign = 0;
265	if (hy > hx) {
266		x = fy;
267		y = fx;
268		i = hx;
269		hx = hy;
270		hy = i;
271		if (iy < 0) {
272			x = -x;
273			sign = 1;
274		}
275		if (ix < 0) {
276			y = -y;
277			a = pio2;
278		} else {
279			a = -pio2;
280			sign = 1 - sign;
281		}
282	} else {
283		y = fy;
284		x = fx;
285		if (iy < 0) {
286			y = -y;
287			sign = 1;
288		}
289		if (ix < 0) {
290			x = -x;
291			a = negpi;
292			sign = 1 - sign;
293		} else {
294			a = zero;
295		}
296	}
297
298	if (hx >= 0x7f800000 || hx - hy >= 0x0c800000) {
299		if (hx >= 0x7f800000) {
300			if (hx > 0x7f800000) /* nan */
301				return (x * y);
302			else if (hy >= 0x7f800000)
303				a += pio4;
304		} else if ((int)a == 0) {
305			a = (double)y / x;
306		}
307		return ((float)((sign)? -a : a));
308	}
309
310	if (hy < 0x00800000) {
311		if (hy == 0)
312			return ((float)((sign)? -a : a));
313		/* scale subnormal y */
314		y *= two24;
315		x *= two24;
316		hy = *(int *)&y;
317		hx = *(int *)&x;
318	}
319
320#if defined(__i386) && !defined(__amd64)
321	rp = __swapRP(fp_extended);
322#endif
323	k = (hy - hx + 0x3f800000) & 0xfff80000;
324	if (k >= 0x3c800000) {	/* |y/x| >= 1/64 */
325		*(int *)&base = k;
326		k = (k - 0x3c800000) >> 19;
327		a += TBL[k];
328	} else {
329		/*
330		 * For some reason this is faster on USIII than just
331		 * doing t = y/x in this case.
332		 */
333		*(int *)&base = 0;
334	}
335	dbase = (double)base;
336	t = (y - x * dbase) / (x + y * dbase);
337	s = t * t;
338	a = (a + t) + t * s * (q1 + s * q2);
339#if defined(__i386) && !defined(__amd64)
340	if (rp != fp_extended)
341		(void) __swapRP(rp);
342#endif
343	return ((float)((sign)? -a : a));
344}
345