1 /*
2 * IBM Accurate Mathematical Library
3 * written by International Business Machines Corp.
4 * Copyright (C) 2001-2023 Free Software Foundation, Inc.
5 *
6 * This program is free software; you can redistribute it and/or modify
7 * it under the terms of the GNU Lesser General Public License as published by
8 * the Free Software Foundation; either version 2.1 of the License, or
9 * (at your option) any later version.
10 *
11 * This program is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU Lesser General Public License for more details.
15 *
16 * You should have received a copy of the GNU Lesser General Public License
17 * along with this program; if not, see <https://www.gnu.org/licenses/>.
18 */
19 /************************************************************************/
20 /* MODULE_NAME: atnat.c */
21 /* */
22 /* FUNCTIONS: uatan */
23 /* signArctan */
24 /* */
25 /* FILES NEEDED: dla.h endian.h mydefs.h atnat.h */
26 /* uatan.tbl */
27 /* */
28 /************************************************************************/
29
30 #include <dla.h>
31 #include "mydefs.h"
32 #include "uatan.tbl"
33 #include "atnat.h"
34 #include <fenv.h>
35 #include <float.h>
36 #include <libm-alias-double.h>
37 #include <math.h>
38 #include <fenv_private.h>
39 #include <math-underflow.h>
40
41 #define TWO52 0x1.0p52
42
43 /* Fix the sign of y and return */
44 static double
45 __signArctan (double x, double y)
46 {
47 return copysign (y, x);
48 }
49
50 /* atan with max ULP of ~0.523 based on random sampling. */
51 double
52 __atan (double x)
53 {
54 double cor, t1, t2, t3, u,
55 v, w, ww, y, yy, z;
56 int i, ux, dx;
57 mynumber num;
58
59 num.d = x;
60 ux = num.i[HIGH_HALF];
61 dx = num.i[LOW_HALF];
62
63 /* x=NaN */
64 if (((ux & 0x7ff00000) == 0x7ff00000)
65 && (((ux & 0x000fffff) | dx) != 0x00000000))
66 return x + x;
67
68 /* Regular values of x, including denormals +-0 and +-INF */
69 SET_RESTORE_ROUND (FE_TONEAREST);
70 u = (x < 0) ? -x : x;
71 if (u < C)
72 {
73 if (u < B)
74 {
75 if (u < A)
76 {
77 math_check_force_underflow_nonneg (u);
78 return x;
79 }
80 else
81 { /* A <= u < B */
82 v = x * x;
83 yy = d11.d + v * d13.d;
84 yy = d9.d + v * yy;
85 yy = d7.d + v * yy;
86 yy = d5.d + v * yy;
87 yy = d3.d + v * yy;
88 yy *= x * v;
89
90 y = x + yy;
91 /* Max ULP is 0.511. */
92 return y;
93 }
94 }
95 else
96 { /* B <= u < C */
97 i = (TWO52 + 256 * u) - TWO52;
98 i -= 16;
99 z = u - cij[i][0].d;
100 yy = cij[i][5].d + z * cij[i][6].d;
101 yy = cij[i][4].d + z * yy;
102 yy = cij[i][3].d + z * yy;
103 yy = cij[i][2].d + z * yy;
104 yy *= z;
105
106 t1 = cij[i][1].d;
107 y = t1 + yy;
108 /* Max ULP is 0.56. */
109 return __signArctan (x, y);
110 }
111 }
112 else
113 {
114 if (u < D)
115 { /* C <= u < D */
116 w = 1 / u;
117 EMULV (w, u, t1, t2);
118 ww = w * ((1 - t1) - t2);
119 i = (TWO52 + 256 * w) - TWO52;
120 i -= 16;
121 z = (w - cij[i][0].d) + ww;
122
123 yy = cij[i][5].d + z * cij[i][6].d;
124 yy = cij[i][4].d + z * yy;
125 yy = cij[i][3].d + z * yy;
126 yy = cij[i][2].d + z * yy;
127 yy = HPI1 - z * yy;
128
129 t1 = HPI - cij[i][1].d;
130 y = t1 + yy;
131 /* Max ULP is 0.503. */
132 return __signArctan (x, y);
133 }
134 else
135 {
136 if (u < E)
137 { /* D <= u < E */
138 w = 1 / u;
139 v = w * w;
140 EMULV (w, u, t1, t2);
141
142 yy = d11.d + v * d13.d;
143 yy = d9.d + v * yy;
144 yy = d7.d + v * yy;
145 yy = d5.d + v * yy;
146 yy = d3.d + v * yy;
147 yy *= w * v;
148
149 ww = w * ((1 - t1) - t2);
150 ESUB (HPI, w, t3, cor);
151 yy = ((HPI1 + cor) - ww) - yy;
152 y = t3 + yy;
153 /* Max ULP is 0.5003. */
154 return __signArctan (x, y);
155 }
156 else
157 {
158 /* u >= E */
159 if (x > 0)
160 return HPI;
161 else
162 return MHPI;
163 }
164 }
165 }
166 }
167
168 #ifndef __atan
169 libm_alias_double (__atan, atan)
170 #endif