1 /* $OpenBSD: e_atan2l.c,v 1.1 2008/12/09 20:00:35 martynas Exp $ */
2 /* @(#)e_atan2.c 1.3 95/01/18 */
3 /* FreeBSD: head/lib/msun/src/e_atan2.c 176451 2008-02-22 02:30:36Z das */
5 * ====================================================
6 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
8 * Developed at SunSoft, a Sun Microsystems, Inc. business.
9 * Permission to use, copy, modify, and distribute this
10 * software is freely granted, provided that this notice
12 * ====================================================
17 * See comments in e_atan2.c.
18 * Converted to long double by David Schultz <das@FreeBSD.ORG>.
19 * Adapted for OpenBSD by Martynas Venckus <martynas@openbsd.org>.
27 #include "math_private.h"
29 #ifdef EXT_IMPLICIT_NBIT
31 #else /* EXT_IMPLICIT_NBIT */
32 #define LDBL_NBIT 0x80000000
33 #endif /* EXT_IMPLICIT_NBIT */
35 static volatile long double
37 static const long double
40 #if defined(__amd64__) || defined(__i386__)
41 /* XXX Work around the fact that gcc truncates long double constants on i386 */
42 static volatile double
43 pi1
= 3.14159265358979311600e+00, /* 0x1.921fb54442d18p+1 */
44 pi2
= 1.22514845490862001043e-16; /* 0x1.1a80000000000p-53 */
45 #define pi ((long double)pi1 + pi2)
47 static const long double
48 pi
= 3.14159265358979323846264338327950280e+00L;
52 atan2l(long double y
, long double x
)
60 int16_t exptx
, expsignx
, expty
, expsigny
;
63 expsigny
= (uy
.bits
.ext_sign
<< 15) | uy
.bits
.ext_exp
;
64 expty
= expsigny
& 0x7fff;
66 expsignx
= (ux
.bits
.ext_sign
<< 15) | ux
.bits
.ext_exp
;
67 exptx
= expsignx
& 0x7fff;
69 if ((exptx
==BIAS
+LDBL_MAX_EXP
&&
70 ((ux
.bits
.ext_frach
&~LDBL_NBIT
)
73 #endif /* EXT_FRACHMBITS */
76 #endif /* EXT_FRACLMBITS */
77 | ux
.bits
.ext_fracl
)!=0) || /* x is NaN */
78 (expty
==BIAS
+LDBL_MAX_EXP
&&
79 ((uy
.bits
.ext_frach
&~LDBL_NBIT
)
82 #endif /* EXT_FRACHMBITS */
85 #endif /* EXT_FRACLMBITS */
86 | uy
.bits
.ext_fracl
)!=0)) /* y is NaN */
88 if (expsignx
==BIAS
&& ((ux
.bits
.ext_frach
&~LDBL_NBIT
)
91 #endif /* EXT_FRACHMBITS */
94 #endif /* EXT_FRACLMBITS */
95 | ux
.bits
.ext_fracl
)==0)
96 return atanl(y
); /* x=1.0 */
97 m
= ((expsigny
>>15)&1)|((expsignx
>>14)&2); /* 2*sign(x)+sign(y) */
100 if(expty
==0 && ((uy
.bits
.ext_frach
&~LDBL_NBIT
)
101 #ifdef EXT_FRACHMBITS
103 #endif /* EXT_FRACHMBITS */
104 #ifdef EXT_FRACLMBITS
106 #endif /* EXT_FRACLMBITS */
107 | uy
.bits
.ext_fracl
)==0) {
110 case 1: return y
; /* atan(+-0,+anything)=+-0 */
111 case 2: return pi
+tiny
;/* atan(+0,-anything) = pi */
112 case 3: return -pi
-tiny
;/* atan(-0,-anything) =-pi */
116 if(exptx
==0 && ((ux
.bits
.ext_frach
&~LDBL_NBIT
)
117 #ifdef EXT_FRACHMBITS
119 #endif /* EXT_FRACHMBITS */
120 #ifdef EXT_FRACLMBITS
122 #endif /* EXT_FRACLMBITS */
123 | ux
.bits
.ext_fracl
)==0)
124 return (expsigny
<0)? -pio2_hi
-tiny
: pio2_hi
+tiny
;
127 if(exptx
==BIAS
+LDBL_MAX_EXP
) {
128 if(expty
==BIAS
+LDBL_MAX_EXP
) {
130 case 0: return pio2_hi
*0.5+tiny
;/* atan(+INF,+INF) */
131 case 1: return -pio2_hi
*0.5-tiny
;/* atan(-INF,+INF) */
132 case 2: return 1.5*pio2_hi
+tiny
;/*atan(+INF,-INF)*/
133 case 3: return -1.5*pio2_hi
-tiny
;/*atan(-INF,-INF)*/
137 case 0: return zero
; /* atan(+...,+INF) */
138 case 1: return -zero
; /* atan(-...,+INF) */
139 case 2: return pi
+tiny
; /* atan(+...,-INF) */
140 case 3: return -pi
-tiny
; /* atan(-...,-INF) */
145 if(expty
==BIAS
+LDBL_MAX_EXP
)
146 return (expsigny
<0)? -pio2_hi
-tiny
: pio2_hi
+tiny
;
150 if(k
> LDBL_MANT_DIG
+2) { /* |y/x| huge */
154 else if(expsignx
<0&&k
<-LDBL_MANT_DIG
-2) z
=0.0; /* |y/x| tiny, x<0 */
155 else z
=atanl(fabsl(y
/x
)); /* safe to do y/x */
157 case 0: return z
; /* atan(+,+) */
158 case 1: return -z
; /* atan(-,+) */
159 case 2: return pi
-(z
-pi_lo
);/* atan(+,-) */
160 default: /* case 3 */
161 return (z
-pi_lo
)-pi
;/* atan(-,-) */