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)
34 extern 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
43 static const double TBL[] = {
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
239 static 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
247 static const float two24 = 16777216.0;
248
249 float
250 atan2f(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 }
|
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 /*
23 * Copyright 2011 Nexenta Systems, Inc. All rights reserved.
24 */
25
26 /*
27 * Copyright 2005 Sun Microsystems, Inc. All rights reserved.
28 * Use is subject to license terms.
29 */
30
31 #pragma weak __atan2f = atan2f
32
33 #include "libm.h"
34
35 #if defined(__i386) && !defined(__amd64)
36 extern int __swapRP(int);
37 #endif
38
39 /*
40 * For i = 0, ..., 192, let x[i] be the double precision number whose
41 * high order 32 bits are 0x3f900000 + (i << 16) and whose low order
42 * 32 bits are zero. Then TBL[i] := atan(x[i]) to double precision.
43 */
44
45 static const double TBL[] = {
228 1.54807296595325550e+00,
229 1.54906061995310385e+00,
230 1.54996600675867957e+00,
231 1.55079899282174605e+00,
232 1.55156792769518947e+00,
233 1.55227992472688747e+00,
234 1.55294108165534417e+00,
235 1.55355665560036682e+00,
236 1.55413120308095598e+00,
237 1.55466869295126031e+00,
238 1.55517259817441977e+00,
239 };
240
241 static const double
242 pio4 = 7.8539816339744827900e-01,
243 pio2 = 1.5707963267948965580e+00,
244 negpi = -3.1415926535897931160e+00,
245 q1 = -3.3333333333296428046e-01,
246 q2 = 1.9999999186853752618e-01,
247 zero = 0.0;
248 static const float two24 = 16777216.0;
249
250 float
251 atan2f(float fy, float fx)
252 {
253 double a, t, s, dbase;
254 float x, y, base;
255 int i, k, hx, hy, ix, iy, sign;
256
257 #if defined(__i386) && !defined(__amd64)
258 int rp;
259 #endif
260
261 iy = *(int *)&fy;
262 ix = *(int *)&fx;
263 hy = iy & ~0x80000000;
264 hx = ix & ~0x80000000;
265
266 sign = 0;
267
268 if (hy > hx) {
269 x = fy;
270 y = fx;
271 i = hx;
272 hx = hy;
273 hy = i;
274
275 if (iy < 0) {
276 x = -x;
277 sign = 1;
278 }
279
280 if (ix < 0) {
281 y = -y;
282 a = pio2;
283 } else {
284 a = -pio2;
285 sign = 1 - sign;
286 }
287 } else {
288 y = fy;
289 x = fx;
290
291 if (iy < 0) {
292 y = -y;
293 sign = 1;
294 }
295
296 if (ix < 0) {
297 x = -x;
298 a = negpi;
299 sign = 1 - sign;
300 } else {
301 a = zero;
302 }
303 }
304
305 if (hx >= 0x7f800000 || hx - hy >= 0x0c800000) {
306 if (hx >= 0x7f800000) {
307 if (hx > 0x7f800000) /* nan */
308 return (x * y);
309 else if (hy >= 0x7f800000)
310 a += pio4;
311 } else if ((int)a == 0) {
312 a = (double)y / x;
313 }
314
315 return ((float)((sign) ? -a : a));
316 }
317
318 if (hy < 0x00800000) {
319 if (hy == 0)
320 return ((float)((sign) ? -a : a));
321
322 /* scale subnormal y */
323 y *= two24;
324 x *= two24;
325 hy = *(int *)&y;
326 hx = *(int *)&x;
327 }
328
329 #if defined(__i386) && !defined(__amd64)
330 rp = __swapRP(fp_extended);
331 #endif
332 k = (hy - hx + 0x3f800000) & 0xfff80000;
333
334 if (k >= 0x3c800000) { /* |y/x| >= 1/64 */
335 *(int *)&base = k;
336 k = (k - 0x3c800000) >> 19;
337 a += TBL[k];
338 } else {
339 /*
340 * For some reason this is faster on USIII than just
341 * doing t = y/x in this case.
342 */
343 *(int *)&base = 0;
344 }
345
346 dbase = (double)base;
347 t = (y - x * dbase) / (x + y * dbase);
348 s = t * t;
349 a = (a + t) + t * s * (q1 + s * q2);
350 #if defined(__i386) && !defined(__amd64)
351 if (rp != fp_extended)
352 (void) __swapRP(rp);
353 #endif
354 return ((float)((sign) ? -a : a));
355 }
|