|
|
1.1 root 1: /* @(#)atan.c 4.1/4.2 10/31/84 CCI-CPG */
2:
3:
4: /*
5: * Double-precision arctangent function recoded with
6: * new algorithm polynomial at CCI-CPG.
7: *
8: * atan returns the value of the arctangent of its
9: * argument in the range [-pi/2,pi/2].
10: *
11: * atan2 returns the arctangent of arg1/arg2
12: * in the range [-pi,pi].
13: *
14: * There are no error returns.
15: *
16: * New version by Les Powers (3/31/85)
17: */
18:
19: #include "atand.h"
20:
21:
22: static double pio2 = 1.570796326794896619231e0;
23: static double pio4 = 0.785398163397448309615e0;
24: static double ap0 = 1.0;
25: static double ap1 = -1.0/3.0;
26: static double ap2 = 1.0/5.0;
27: static double ap3 = -1.0/7.0;
28:
29:
30: /*
31: * atan reduces its argument
32: * and returns the arctan.
33: */
34: double
35: atan(x)
36: double x;
37: {
38: double t,s;
39: unsigned ta;
40: union {
41: double d;
42: struct { unsigned i0,i1; } i;
43: struct { unsigned char b0,b1,b2,b3; } b;
44: } u;
45: if ( x < 0 )
46: return( -atan(-x) );
47: u.d = x;
48: if (u.i.i0 >= 0x41000000 ) { /* if (x >= 2.0) */
49: return( pio2 - atan( 1.0/x ) );
50: }
51: if (u.i.i0 >= 0x40000000 ) { /* if (x >= 0.5) */
52: u.i.i0 = (u.i.i0 & 0x40ff8000) | 0x00008000;
53: u.i.i1 = 0;
54: t = (x-u.d)/(x*u.d+1.0);
55: s = t*t;
56: return( t*(ap0+s*(ap1+s*ap2)) + at1[u.b.b1] );
57: }
58: if (u.i.i0 >= 0x3d000000) { /* if (x >= 1/128) */
59: u.i.i0 += 0x03800000;
60: ta = u.d;
61: u.d = ta;
62: u.i.i0 -= 0x03800000;
63: t = (x-u.d)/(x*u.d+1.0);
64: }
65: else {
66: ta = 0;
67: t = x;
68: }
69: s = t*t;
70: return( t*(ap0+s*(ap1+s*(ap2+s*ap3))) + at0[ta] );
71: }
72:
73:
74: /*
75: * atan2 discovers what quadrant the angle
76: * is in and calls atan.
77: */
78: double
79: atan2(y,x)
80: double y,x;
81: {
82: double atan();
83:
84: if (y > 0)
85: if ( x >= y )
86: return( atan(y/x));
87: else if ( x >= 0 )
88: return( pio2 - atan(x/y));
89: else if (-x <= y )
90: return( pio2 + atan(-x/y));
91: else
92: return( pio2 + pio2 - atan(-y/x));
93: else if (y < 0)
94: if ( x >= -y )
95: return( -atan(-y/x));
96: else if ( x >= 0 )
97: return(-pio2 + atan(x/(-y)));
98: else if (-x <= -y )
99: return(-pio2 - atan(x/y));
100: else
101: return(-pio2 - pio2 + atan(y/x));
102: else
103: if ( x > 0 )
104: return( 0 );
105: else if ( x < 0 )
106: return( pio2 + pio2 );
107: else
108: return( pio2 );
109: }
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.