1 | /*
|
---|
2 | * (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
|
---|
3 | * See the copyright notice in the ACK home directory, in the file "Copyright".
|
---|
4 | *
|
---|
5 | * Author: Ceriel J.H. Jacobs
|
---|
6 | */
|
---|
7 | /* $Header: /cvsup/minix/src/lib/math/atan.c,v 1.1.1.1 2005/04/21 14:56:24 beng Exp $ */
|
---|
8 |
|
---|
9 | #include <float.h>
|
---|
10 | #include <math.h>
|
---|
11 | #include <errno.h>
|
---|
12 | #include "localmath.h"
|
---|
13 |
|
---|
14 | double
|
---|
15 | atan(double x)
|
---|
16 | {
|
---|
17 | /* Algorithm and coefficients from:
|
---|
18 | "Software manual for the elementary functions"
|
---|
19 | by W.J. Cody and W. Waite, Prentice-Hall, 1980
|
---|
20 | */
|
---|
21 |
|
---|
22 | static double p[] = {
|
---|
23 | -0.13688768894191926929e+2,
|
---|
24 | -0.20505855195861651981e+2,
|
---|
25 | -0.84946240351320683534e+1,
|
---|
26 | -0.83758299368150059274e+0
|
---|
27 | };
|
---|
28 | static double q[] = {
|
---|
29 | 0.41066306682575781263e+2,
|
---|
30 | 0.86157349597130242515e+2,
|
---|
31 | 0.59578436142597344465e+2,
|
---|
32 | 0.15024001160028576121e+2,
|
---|
33 | 1.0
|
---|
34 | };
|
---|
35 | static double a[] = {
|
---|
36 | 0.0,
|
---|
37 | 0.52359877559829887307710723554658381, /* pi/6 */
|
---|
38 | M_PI_2,
|
---|
39 | 1.04719755119659774615421446109316763 /* pi/3 */
|
---|
40 | };
|
---|
41 |
|
---|
42 | int neg = x < 0;
|
---|
43 | int n;
|
---|
44 | double g;
|
---|
45 |
|
---|
46 | if (__IsNan(x)) {
|
---|
47 | errno = EDOM;
|
---|
48 | return x;
|
---|
49 | }
|
---|
50 | if (neg) {
|
---|
51 | x = -x;
|
---|
52 | }
|
---|
53 | if (x > 1.0) {
|
---|
54 | x = 1.0/x;
|
---|
55 | n = 2;
|
---|
56 | }
|
---|
57 | else n = 0;
|
---|
58 |
|
---|
59 | if (x > 0.26794919243112270647) { /* 2-sqtr(3) */
|
---|
60 | n = n + 1;
|
---|
61 | x = (((0.73205080756887729353*x-0.5)-0.5)+x)/
|
---|
62 | (1.73205080756887729353+x);
|
---|
63 | }
|
---|
64 |
|
---|
65 | /* ??? avoid underflow ??? */
|
---|
66 |
|
---|
67 | g = x * x;
|
---|
68 | x += x * g * POLYNOM3(g, p) / POLYNOM4(g, q);
|
---|
69 | if (n > 1) x = -x;
|
---|
70 | x += a[n];
|
---|
71 | return neg ? -x : x;
|
---|
72 | }
|
---|