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 2006 Sun Microsystems, Inc. All rights reserved.
26 * Use is subject to license terms.
27 */
28
29 #pragma weak __atan2 = atan2
30
31 #include "libm.h"
32
33 /*
34 * Let t(0) = 1 and for i = 1, ..., 160, let t(i) be the slope of
35 * the line bisecting the conical hull of the set of points (x,y)
36 * where x and y are positive normal floating point numbers and
37 * the high order words hx and hy of their binary representations
38 * satisfy |hx - hy - i * 0x8000| <= 0x4000. Then:
39 *
40 * TBL[4*i+2] is t(i) rounded to 21 significant bits (i.e., the
41 * low order word is zero), and
42 *
43 * TBL[4*i] + TBL[4*i+1] form a doubled-double approximation to
389 -3.33333333333327571893331786354179101074860633009e-0001,
390 +1.99999999942671624230086497610394721817438631379e-0001,
391 -1.42856965565428636896183013324727205980484158356e-0001,
392 +1.10894981496317081405107718475040168084164825641e-0001,
393 };
394
395 #define zero C[0]
396 #define twom3 C[1]
397 #define two110 C[2]
398 #define pio4 C[3]
399 #define pio2 C[4]
400 #define pio2_lo C[5]
401 #define mpi C[6]
402 #define mpi_lo C[7]
403 #define p1 C[8]
404 #define p2 C[9]
405 #define p3 C[10]
406 #define p4 C[11]
407
408 double
409 atan2(double oy, double ox) {
410 double ah, al, t, xh, x, y, z;
411 int i, k, hx, hy, sx, sy;
412 #ifndef lint
413 volatile int inexact __unused;
414 #endif
415
416 hy = ((int *)&oy)[HIWORD];
417 sy = hy & 0x80000000;
418 hy &= ~0x80000000;
419
420 hx = ((int *)&ox)[HIWORD];
421 sx = hx & 0x80000000;
422 hx &= ~0x80000000;
423
424 if (hy > hx || (hy == hx && ((unsigned *)&oy)[LOWORD] >
425 ((unsigned *)&ox)[LOWORD])) {
426 i = hx;
427 hx = hy;
428 hy = i;
429 x = fabs(oy);
430 y = fabs(ox);
431 if (sx) {
432 ah = pio2;
433 al = pio2_lo;
434 } else {
435 ah = -pio2;
436 al = -pio2_lo;
437 sy ^= 0x80000000;
438 }
439 } else {
440 x = fabs(ox);
441 y = fabs(oy);
442 if (sx) {
443 ah = mpi;
444 al = mpi_lo;
445 sy ^= 0x80000000;
446 } else {
447 ah = al = zero;
448 }
449 }
450
451 if (hx >= 0x7fe00000 || hx - hy >= 0x03600000) {
452 if (hx >= 0x7ff00000) {
453 if (((hx ^ 0x7ff00000) | ((int *)&x)[LOWORD]) != 0)
454 return (ox * oy);
455 if (hy >= 0x7ff00000)
456 ah += pio4;
457 #ifndef lint
458 inexact = (int)ah; /* inexact if ah != 0 */
459 #endif
460 return ((sy)? -ah : ah);
461 }
462 if (hx - hy >= 0x03600000) {
463 if ((int)ah == 0)
464 ah = y / x;
465 return ((sy)? -ah : ah);
466 }
467 y *= twom3;
468 x *= twom3;
469 hy -= 0x00300000;
470 hx -= 0x00300000;
471 } else if (hy < 0x00100000) {
472 if ((hy | ((int *)&y)[LOWORD]) == 0) {
473 if ((hx | ((int *)&x)[LOWORD]) == 0)
474 return (_SVID_libm_err(ox, oy, 3));
475 #ifndef lint
476 inexact = (int)ah; /* inexact if ah != 0 */
477 #endif
478 return ((sy)? -ah : ah);
479 }
480 y *= two110;
481 x *= two110;
482 hy = ((int *)&y)[HIWORD];
483 hx = ((int *)&x)[HIWORD];
484 }
485
486 k = (((hx - hy) + 0x00004000) >> 13) & ~0x3;
487 if (k > 644)
488 k = 644;
489 ah += TBL[k];
490 al += TBL[k+1];
491 t = TBL[k+2];
492
493 xh = x;
494 ((int *)&xh)[LOWORD] = 0;
495 z = ((y - t * xh) - t * (x - xh)) / (x + y * t);
496 x = z * z;
497 t = ah + (z + (al + (z * x) * (p1 + x * (p2 + x * (p3 + x * p4)))));
498 return ((sy)? -t : t);
499 }
|
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 2006 Sun Microsystems, Inc. All rights reserved.
28 * Use is subject to license terms.
29 */
30
31 #pragma weak __atan2 = atan2
32
33 #include "libm.h"
34
35 /*
36 * Let t(0) = 1 and for i = 1, ..., 160, let t(i) be the slope of
37 * the line bisecting the conical hull of the set of points (x,y)
38 * where x and y are positive normal floating point numbers and
39 * the high order words hx and hy of their binary representations
40 * satisfy |hx - hy - i * 0x8000| <= 0x4000. Then:
41 *
42 * TBL[4*i+2] is t(i) rounded to 21 significant bits (i.e., the
43 * low order word is zero), and
44 *
45 * TBL[4*i] + TBL[4*i+1] form a doubled-double approximation to
391 -3.33333333333327571893331786354179101074860633009e-0001,
392 +1.99999999942671624230086497610394721817438631379e-0001,
393 -1.42856965565428636896183013324727205980484158356e-0001,
394 +1.10894981496317081405107718475040168084164825641e-0001,
395 };
396
397 #define zero C[0]
398 #define twom3 C[1]
399 #define two110 C[2]
400 #define pio4 C[3]
401 #define pio2 C[4]
402 #define pio2_lo C[5]
403 #define mpi C[6]
404 #define mpi_lo C[7]
405 #define p1 C[8]
406 #define p2 C[9]
407 #define p3 C[10]
408 #define p4 C[11]
409
410 double
411 atan2(double oy, double ox)
412 {
413 double ah, al, t, xh, x, y, z;
414 int i, k, hx, hy, sx, sy;
415
416 #ifndef lint
417 volatile int inexact __unused;
418 #endif
419
420 hy = ((int *)&oy)[HIWORD];
421 sy = hy & 0x80000000;
422 hy &= ~0x80000000;
423
424 hx = ((int *)&ox)[HIWORD];
425 sx = hx & 0x80000000;
426 hx &= ~0x80000000;
427
428 if (hy > hx || (hy == hx && ((unsigned *)&oy)[LOWORD] >
429 ((unsigned *)&ox)[LOWORD])) {
430 i = hx;
431 hx = hy;
432 hy = i;
433 x = fabs(oy);
434 y = fabs(ox);
435
436 if (sx) {
437 ah = pio2;
438 al = pio2_lo;
439 } else {
440 ah = -pio2;
441 al = -pio2_lo;
442 sy ^= 0x80000000;
443 }
444 } else {
445 x = fabs(ox);
446 y = fabs(oy);
447
448 if (sx) {
449 ah = mpi;
450 al = mpi_lo;
451 sy ^= 0x80000000;
452 } else {
453 ah = al = zero;
454 }
455 }
456
457 if (hx >= 0x7fe00000 || hx - hy >= 0x03600000) {
458 if (hx >= 0x7ff00000) {
459 if (((hx ^ 0x7ff00000) | ((int *)&x)[LOWORD]) != 0)
460 return (ox * oy);
461
462 if (hy >= 0x7ff00000)
463 ah += pio4;
464
465 #ifndef lint
466 inexact = (int)ah; /* inexact if ah != 0 */
467 #endif
468 return ((sy) ? -ah : ah);
469 }
470
471 if (hx - hy >= 0x03600000) {
472 if ((int)ah == 0)
473 ah = y / x;
474
475 return ((sy) ? -ah : ah);
476 }
477
478 y *= twom3;
479 x *= twom3;
480 hy -= 0x00300000;
481 hx -= 0x00300000;
482 } else if (hy < 0x00100000) {
483 if ((hy | ((int *)&y)[LOWORD]) == 0) {
484 if ((hx | ((int *)&x)[LOWORD]) == 0)
485 return (_SVID_libm_err(ox, oy, 3));
486
487 #ifndef lint
488 inexact = (int)ah; /* inexact if ah != 0 */
489 #endif
490 return ((sy) ? -ah : ah);
491 }
492
493 y *= two110;
494 x *= two110;
495 hy = ((int *)&y)[HIWORD];
496 hx = ((int *)&x)[HIWORD];
497 }
498
499 k = (((hx - hy) + 0x00004000) >> 13) & ~0x3;
500
501 if (k > 644)
502 k = 644;
503
504 ah += TBL[k];
505 al += TBL[k + 1];
506 t = TBL[k + 2];
507
508 xh = x;
509 ((int *)&xh)[LOWORD] = 0;
510 z = ((y - t * xh) - t * (x - xh)) / (x + y * t);
511 x = z * z;
512 t = ah + (z + (al + (z * x) * (p1 + x * (p2 + x * (p3 + x * p4)))));
513 return ((sy) ? -t : t);
514 }
|