[9] | 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 | }
|
---|