1 |
|
|
/* $OpenBSD: e_atan2l.c,v 1.2 2016/09/12 19:47:02 guenther 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 */ |
4 |
|
|
/* |
5 |
|
|
* ==================================================== |
6 |
|
|
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. |
7 |
|
|
* |
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 |
11 |
|
|
* is preserved. |
12 |
|
|
* ==================================================== |
13 |
|
|
* |
14 |
|
|
*/ |
15 |
|
|
|
16 |
|
|
/* |
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>. |
20 |
|
|
*/ |
21 |
|
|
|
22 |
|
|
#include <float.h> |
23 |
|
|
#include <math.h> |
24 |
|
|
|
25 |
|
|
#include "invtrig.h" |
26 |
|
|
#include "math.h" |
27 |
|
|
#include "math_private.h" |
28 |
|
|
|
29 |
|
|
#ifdef EXT_IMPLICIT_NBIT |
30 |
|
|
#define LDBL_NBIT 0 |
31 |
|
|
#else /* EXT_IMPLICIT_NBIT */ |
32 |
|
|
#define LDBL_NBIT 0x80000000 |
33 |
|
|
#endif /* EXT_IMPLICIT_NBIT */ |
34 |
|
|
|
35 |
|
|
static volatile long double |
36 |
|
|
tiny = 1.0e-300; |
37 |
|
|
static const long double |
38 |
|
|
zero = 0.0; |
39 |
|
|
|
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) |
46 |
|
|
#else |
47 |
|
|
static const long double |
48 |
|
|
pi = 3.14159265358979323846264338327950280e+00L; |
49 |
|
|
#endif |
50 |
|
|
|
51 |
|
|
long double |
52 |
|
|
atan2l(long double y, long double x) |
53 |
|
|
{ |
54 |
|
|
union { |
55 |
|
|
long double e; |
56 |
|
|
struct ieee_ext bits; |
57 |
|
|
} ux, uy; |
58 |
|
|
long double z; |
59 |
|
|
int32_t k,m; |
60 |
|
|
int16_t exptx, expsignx, expty, expsigny; |
61 |
|
|
|
62 |
|
|
uy.e = y; |
63 |
|
|
expsigny = (uy.bits.ext_sign << 15) | uy.bits.ext_exp; |
64 |
|
|
expty = expsigny & 0x7fff; |
65 |
|
|
ux.e = x; |
66 |
|
|
expsignx = (ux.bits.ext_sign << 15) | ux.bits.ext_exp; |
67 |
|
|
exptx = expsignx & 0x7fff; |
68 |
|
|
|
69 |
|
|
if ((exptx==BIAS+LDBL_MAX_EXP && |
70 |
|
|
((ux.bits.ext_frach&~LDBL_NBIT) |
71 |
|
|
#ifdef EXT_FRACHMBITS |
72 |
|
|
| ux.bits.ext_frachm |
73 |
|
|
#endif /* EXT_FRACHMBITS */ |
74 |
|
|
#ifdef EXT_FRACLMBITS |
75 |
|
|
| ux.bits.ext_fraclm |
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) |
80 |
|
|
#ifdef EXT_FRACHMBITS |
81 |
|
|
| uy.bits.ext_frachm |
82 |
|
|
#endif /* EXT_FRACHMBITS */ |
83 |
|
|
#ifdef EXT_FRACLMBITS |
84 |
|
|
| uy.bits.ext_fraclm |
85 |
|
|
#endif /* EXT_FRACLMBITS */ |
86 |
|
|
| uy.bits.ext_fracl)!=0)) /* y is NaN */ |
87 |
|
|
return x+y; |
88 |
|
|
if (expsignx==BIAS && ((ux.bits.ext_frach&~LDBL_NBIT) |
89 |
|
|
#ifdef EXT_FRACHMBITS |
90 |
|
|
| ux.bits.ext_frachm |
91 |
|
|
#endif /* EXT_FRACHMBITS */ |
92 |
|
|
#ifdef EXT_FRACLMBITS |
93 |
|
|
| ux.bits.ext_fraclm |
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) */ |
98 |
|
|
|
99 |
|
|
/* when y = 0 */ |
100 |
|
|
if(expty==0 && ((uy.bits.ext_frach&~LDBL_NBIT) |
101 |
|
|
#ifdef EXT_FRACHMBITS |
102 |
|
|
| uy.bits.ext_frachm |
103 |
|
|
#endif /* EXT_FRACHMBITS */ |
104 |
|
|
#ifdef EXT_FRACLMBITS |
105 |
|
|
| uy.bits.ext_fraclm |
106 |
|
|
#endif /* EXT_FRACLMBITS */ |
107 |
|
|
| uy.bits.ext_fracl)==0) { |
108 |
|
|
switch(m) { |
109 |
|
|
case 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 */ |
113 |
|
|
} |
114 |
|
|
} |
115 |
|
|
/* when x = 0 */ |
116 |
|
|
if(exptx==0 && ((ux.bits.ext_frach&~LDBL_NBIT) |
117 |
|
|
#ifdef EXT_FRACHMBITS |
118 |
|
|
| ux.bits.ext_frachm |
119 |
|
|
#endif /* EXT_FRACHMBITS */ |
120 |
|
|
#ifdef EXT_FRACLMBITS |
121 |
|
|
| ux.bits.ext_fraclm |
122 |
|
|
#endif /* EXT_FRACLMBITS */ |
123 |
|
|
| ux.bits.ext_fracl)==0) |
124 |
|
|
return (expsigny<0)? -pio2_hi-tiny: pio2_hi+tiny; |
125 |
|
|
|
126 |
|
|
/* when x is INF */ |
127 |
|
|
if(exptx==BIAS+LDBL_MAX_EXP) { |
128 |
|
|
if(expty==BIAS+LDBL_MAX_EXP) { |
129 |
|
|
switch(m) { |
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)*/ |
134 |
|
|
} |
135 |
|
|
} else { |
136 |
|
|
switch(m) { |
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) */ |
141 |
|
|
} |
142 |
|
|
} |
143 |
|
|
} |
144 |
|
|
/* when y is INF */ |
145 |
|
|
if(expty==BIAS+LDBL_MAX_EXP) |
146 |
|
|
return (expsigny<0)? -pio2_hi-tiny: pio2_hi+tiny; |
147 |
|
|
|
148 |
|
|
/* compute y/x */ |
149 |
|
|
k = expty-exptx; |
150 |
|
|
if(k > LDBL_MANT_DIG+2) { /* |y/x| huge */ |
151 |
|
|
z=pio2_hi+pio2_lo; |
152 |
|
|
m&=1; |
153 |
|
|
} |
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 */ |
156 |
|
|
switch (m) { |
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(-,-) */ |
162 |
|
|
} |
163 |
|
|
} |
164 |
|
|
DEF_STD(atan2l); |