overby.me
this repo has no description
1/* This is an implementation of atan. It is written in standard C except:
2
3 double is expected be an IEEE 754 double-precision implementation.
4
5 "volatile" is used to attempt to force certain floating-point
6 operations to occur at run time (to generate exceptions that might not
7 be generated if the operations are performed at compile time).
8*/
9
10
11// Include math.h to ensure we match the declarations.
12#include <math.h>
13
14
15/* Declare certain constants volatile to force the compiler to access them
16 when we reference them. This in turn forces arithmetic operations on them
17 to be performed at run time (or as if at run time). We use such operations
18 to generate exceptions such as underflow or inexact.
19*/
20static volatile const double Tiny = 0x1p-1022;
21
22
23#if defined __STDC__ && 199901L <= __STDC_VERSION__ && !defined __GNUC__
24 // GCC does not currently support FENV_ACCESS. Maybe someday.
25 #pragma STDC FENV_ACCESS ON
26#endif
27
28
29/* double atan(double x).
30
31 (This routine appears below, following subroutines.)
32
33 Notes:
34
35 Citations in parentheses below indicate the source of a requirement.
36
37 "C" stands for ISO/IEC 9899:TC2.
38
39 The Open Group specification (IEEE Std 1003.1, 2004 edition) adds no
40 requirements since it defers to C and requires errno behavior only if
41 we choose to support it by arranging for "math_errhandling &
42 MATH_ERRNO" to be non-zero, which we do not.
43
44 Return value:
45
46 For arctangent of +/- zero, return zero with same sign (C F.9 12 and
47 F.9.1.3).
48
49 For arctangent of +/- infinity, return +/- pi/2 (C F.9.1.3).
50
51 For a NaN, return the same NaN (C F.9 11 and 13). (If the NaN is a
52 signalling NaN, we return the "same" NaN quieted.)
53
54 Otherwise:
55
56 If the rounding mode is round-to-nearest, return arctangent(x)
57 faithfully rounded. This is not proven but seems likely.
58 Generally, the largest source of errors is the evaluation of the
59 polynomial using double precision. Some analysis might bound this
60 and prove faithful rounding. The largest observed error is .814
61 ULP.
62
63 Return a value in [-pi/2, +pi/2] (C 7.12.4.3 3).
64
65 Not implemented: In other rounding modes, return arctangent(x)
66 possibly with slightly worse error, not necessarily honoring the
67 rounding mode (Ali Sazegari narrowing C F.9 10).
68
69 Exceptions:
70
71 Raise underflow for a denormal result (C F.9 7 and Draft Standard for
72 Floating-Point Arithmetic P754 Draft 1.2.5 9.5). If the input is the
73 smallest normal, underflow may or may not be raised. This is stricter
74 than the older 754 standard.
75
76 May or may not raise inexact, even if the result is exact (C F.9 8).
77
78 Raise invalid if the input is a signalling NaN (C 5.2.4.2.2 3, in spite
79 of C 4.2.1), but not if the input is a quiet NaN (C F.9 11).
80
81 May not raise exceptions otherwise (C F.9 9).
82
83 Properties:
84
85 Not proven: Monotonic.
86*/
87
88
89// Return arctangent(x) given that 2 < x, with the same properties as atan.
90static double Tail(double x)
91{
92 {
93 static const double HalfPi = 0x3.243f6a8885a308d313198a2e037ap-1;
94
95 // For large x, generate inexact and return pi/2.
96 if (0x1p53 <= x)
97 return HalfPi + Tiny;
98 if (isnan(x))
99 return x - x;
100 }
101
102 static const double p03 = -0x1.5555555554A51p-2;
103 static const double p05 = +0x1.999999989EBCAp-3;
104 static const double p07 = -0x1.249248E1422E3p-3;
105 static const double p09 = +0x1.C71C5EDFED480p-4;
106 static const double p11 = -0x1.745B7F2D72663p-4;
107 static const double p13 = +0x1.3AFD7A0E6EB75p-4;
108 static const double p15 = -0x1.104146B1A1AE8p-4;
109 static const double p17 = +0x1.D78252FA69C1Cp-5;
110 static const double p19 = -0x1.81D33E401836Dp-5;
111 static const double p21 = +0x1.007733E06CEB3p-5;
112 static const double p23 = -0x1.83DAFDA7BD3FDp-7;
113
114 static const double p000 = +0x1.921FB54442D18p0;
115 static const double p001 = +0x1.1A62633145C07p-54;
116
117 double y = 1/x;
118
119 // Square y.
120 double y2 = y * y;
121
122 return p001 - ((((((((((((
123 + p23) * y2
124 + p21) * y2
125 + p19) * y2
126 + p17) * y2
127 + p15) * y2
128 + p13) * y2
129 + p11) * y2
130 + p09) * y2
131 + p07) * y2
132 + p05) * y2
133 + p03) * y2 * y + y) + p000;
134}
135
136
137/* Return arctangent(x) given that 0x1p-27 < |x| <= 1/2, with the same
138 properties as atan.
139*/
140static double atani0(double x)
141{
142 static const double p03 = -0x1.555555555551Bp-2;
143 static const double p05 = +0x1.99999999918D8p-3;
144 static const double p07 = -0x1.2492492179CA3p-3;
145 static const double p09 = +0x1.C71C7096C2725p-4;
146 static const double p11 = -0x1.745CF51795B21p-4;
147 static const double p13 = +0x1.3B113F18AC049p-4;
148 static const double p15 = -0x1.10F31279EC05Dp-4;
149 static const double p17 = +0x1.DFE7B9674AE37p-5;
150 static const double p19 = -0x1.A38CF590469ECp-5;
151 static const double p21 = +0x1.56CDB5D887934p-5;
152 static const double p23 = -0x1.C0EB85F543412p-6;
153 static const double p25 = +0x1.4A9F5C4724056p-7;
154
155 // Square x.
156 double x2 = x * x;
157
158 return ((((((((((((
159 + p25) * x2
160 + p23) * x2
161 + p21) * x2
162 + p19) * x2
163 + p17) * x2
164 + p15) * x2
165 + p13) * x2
166 + p11) * x2
167 + p09) * x2
168 + p07) * x2
169 + p05) * x2
170 + p03) * x2 * x + x;
171}
172
173
174/* Return arctangent(x) given that 1/2 < x <= 3/4, with the same properties as
175 atan.
176*/
177static double atani1(double x)
178{
179 static const double p00 = +0x1.1E00BABDEFED0p-1;
180 static const double p01 = +0x1.702E05C0B8155p-1;
181 static const double p02 = -0x1.4AF2B78215A1Bp-2;
182 static const double p03 = +0x1.5D0B7E9E69054p-6;
183 static const double p04 = +0x1.A1247CA5D9475p-4;
184 static const double p05 = -0x1.519E110F61B54p-4;
185 static const double p06 = +0x1.A759263F377F2p-7;
186 static const double p07 = +0x1.094966BE2B531p-5;
187 static const double p08 = -0x1.09BC0AB7F914Cp-5;
188 static const double p09 = +0x1.FF3B7C531AA4Ap-8;
189 static const double p10 = +0x1.950E69DCDD967p-7;
190 static const double p11 = -0x1.D88D31ABC3AE5p-7;
191 static const double p12 = +0x1.10F3E20F6A2E2p-8;
192
193 double y = x - 0x1.4000000000027p-1;
194
195 return ((((((((((((
196 + p12) * y
197 + p11) * y
198 + p10) * y
199 + p09) * y
200 + p08) * y
201 + p07) * y
202 + p06) * y
203 + p05) * y
204 + p04) * y
205 + p03) * y
206 + p02) * y
207 + p01) * y
208 + p00;
209}
210
211
212/* Return arctangent(x) given that 3/4 < x <= 1, with the same properties as
213 atan.
214*/
215static double atani2(double x)
216{
217 static const double p00 = +0x1.700A7C580EA7Ep-01;
218 static const double p01 = +0x1.21FB781196AC3p-01;
219 static const double p02 = -0x1.1F6A8499714A2p-02;
220 static const double p03 = +0x1.41B15E5E8DCD0p-04;
221 static const double p04 = +0x1.59BC93F81895Ap-06;
222 static const double p05 = -0x1.63B543EFFA4EFp-05;
223 static const double p06 = +0x1.C90E92AC8D86Cp-06;
224 static const double p07 = -0x1.91F7E2A7A338Fp-08;
225 static const double p08 = -0x1.AC1645739E676p-08;
226 static const double p09 = +0x1.152311B180E6Cp-07;
227 static const double p10 = -0x1.265EF51B17DB7p-08;
228 static const double p11 = +0x1.CA7CDE5DE9BD7p-14;
229
230 double y = x - 0x1.c0000000f4213p-1;
231
232 return (((((((((((
233 + p11) * y
234 + p10) * y
235 + p09) * y
236 + p08) * y
237 + p07) * y
238 + p06) * y
239 + p05) * y
240 + p04) * y
241 + p03) * y
242 + p02) * y
243 + p01) * y
244 + p00;
245}
246
247
248/* Return arctangent(x) given that 1 < x <= 4/3, with the same properties as
249 atan.
250*/
251static double atani3(double x)
252{
253 static const double p00 = +0x1.B96E5A78C5C40p-01;
254 static const double p01 = +0x1.B1B1B1B1B1B3Dp-02;
255 static const double p02 = -0x1.AC97826D58470p-03;
256 static const double p03 = +0x1.3FD2B9F586A67p-04;
257 static const double p04 = -0x1.BC317394714B7p-07;
258 static const double p05 = -0x1.2B01FC60CC37Ap-07;
259 static const double p06 = +0x1.73A9328786665p-07;
260 static const double p07 = -0x1.C0B993A09CE31p-08;
261 static const double p08 = +0x1.2FCDACDD6E5B5p-09;
262 static const double p09 = +0x1.CBD49DA316282p-13;
263 static const double p10 = -0x1.0120E602F6336p-10;
264 static const double p11 = +0x1.A89224FF69018p-11;
265 static const double p12 = -0x1.883D8959134B3p-12;
266
267 double y = x - 0x1.2aaaaaaaaaa96p0;
268
269 return ((((((((((((
270 + p12) * y
271 + p11) * y
272 + p10) * y
273 + p09) * y
274 + p08) * y
275 + p07) * y
276 + p06) * y
277 + p05) * y
278 + p04) * y
279 + p03) * y
280 + p02) * y
281 + p01) * y
282 + p00;
283}
284
285
286/* Return arctangent(x) given that 4/3 < x <= 5/3, with the same properties as
287 atan.
288*/
289static double atani4(double x)
290{
291 static const double p00 = +0x1.F730BD281F69Dp-01;
292 static const double p01 = +0x1.3B13B13B13B0Cp-02;
293 static const double p02 = -0x1.22D719C06115Ep-03;
294 static const double p03 = +0x1.C963C83985742p-05;
295 static const double p04 = -0x1.135A0938EC462p-06;
296 static const double p05 = +0x1.13A254D6E5B7Cp-09;
297 static const double p06 = +0x1.DFAA5E77B7375p-10;
298 static const double p07 = -0x1.F4AC1342182D2p-10;
299 static const double p08 = +0x1.25BAD4D85CBE1p-10;
300 static const double p09 = -0x1.E4EEF429EB680p-12;
301 static const double p10 = +0x1.B4E30D1BA3819p-14;
302 static const double p11 = +0x1.0280537F097F3p-15;
303
304 double y = x - 0x1.8000000000003p0;
305
306 return (((((((((((
307 + p11) * y
308 + p10) * y
309 + p09) * y
310 + p08) * y
311 + p07) * y
312 + p06) * y
313 + p05) * y
314 + p04) * y
315 + p03) * y
316 + p02) * y
317 + p01) * y
318 + p00;
319}
320
321
322/* Return arctangent(x) given that 5/3 < x <= 2, with the same properties as
323 atan.
324*/
325static double atani5(double x)
326{
327 static const double p00 = +0x1.124A85750FB5Cp+00;
328 static const double p01 = +0x1.D59AE78C11C49p-03;
329 static const double p02 = -0x1.8AD3C44F10DC3p-04;
330 static const double p03 = +0x1.2B090AAD5F9DCp-05;
331 static const double p04 = -0x1.881EC3D15241Fp-07;
332 static const double p05 = +0x1.8CB82A74E0699p-09;
333 static const double p06 = -0x1.3182219E21362p-12;
334 static const double p07 = -0x1.2B9AD13DB35A8p-12;
335 static const double p08 = +0x1.10F884EAC0E0Ap-12;
336 static const double p09 = -0x1.3045B70E93129p-13;
337 static const double p10 = +0x1.00B6A460AC05Dp-14;
338
339 double y = x - 0x1.d555555461337p0;
340
341 return ((((((((((
342 + p10) * y
343 + p09) * y
344 + p08) * y
345 + p07) * y
346 + p06) * y
347 + p05) * y
348 + p04) * y
349 + p03) * y
350 + p02) * y
351 + p01) * y
352 + p00;
353}
354
355
356// See documentation above.
357double atan(double x)
358{
359 if (x < 0)
360 if (x < -1)
361 if (x < -5/3.)
362 if (x < -2)
363 return -Tail(-x);
364 else
365 return -atani5(-x);
366 else
367 if (x < -4/3.)
368 return -atani4(-x);
369 else
370 return -atani3(-x);
371 else
372 if (x < -.5)
373 if (x < -.75)
374 return -atani2(-x);
375 else
376 return -atani1(-x);
377 else
378 if (x < -0x1.d12ed0af1a27fp-27)
379 return atani0(x);
380 else
381 if (x <= -0x1p-1022)
382 // Generate inexact and return x.
383 return (Tiny + 1) * x;
384 else
385 if (x == 0)
386 return x;
387 else
388 // Generate underflow and return x.
389 return x*Tiny + x;
390 else
391 if (x <= +1)
392 if (x <= +.5)
393 if (x <= +0x1.d12ed0af1a27fp-27)
394 if (x < +0x1p-1022)
395 if (x == 0)
396 return x;
397 else
398 // Generate underflow and return x.
399 return x*Tiny + x;
400 else
401 // Generate inexact and return x.
402 return (Tiny + 1) * x;
403 else
404 return atani0(x);
405 else
406 if (x <= +.75)
407 return +atani1(+x);
408 else
409 return +atani2(+x);
410 else
411 if (x <= +5/3.)
412 if (x <= +4/3.)
413 return +atani3(+x);
414 else
415 return +atani4(+x);
416 else
417 if (x <= +2)
418 return +atani5(+x);
419 else
420 return +Tail(+x);
421}