amath  1.8.5
Simple command line calculator
atan.c
Go to the documentation of this file.
1 /*-
2  * Copyright (c) 2014-2018 Carsten Sonne Larsen <cs@innolan.net>
3  * All rights reserved.
4  *
5  * Redistribution and use in source and binary forms, with or without
6  * modification, are permitted provided that the following conditions
7  * are met:
8  * 1. Redistributions of source code must retain the above copyright
9  * notice, this list of conditions and the following disclaimer.
10  * 2. Redistributions in binary form must reproduce the above copyright
11  * notice, this list of conditions and the following disclaimer in the
12  * documentation and/or other materials provided with the distribution.
13  *
14  * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
15  * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
16  * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
17  * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
18  * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
19  * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
20  * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
21  * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
22  * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
23  * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
24  *
25  * Project homepage:
26  * https://amath.innolan.net
27  *
28  * The original source code can be obtained from:
29  * http://www.netlib.org/fdlibm/s_atan.c
30  *
31  * =================================================================
32  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
33  *
34  * Developed at SunSoft, a Sun Microsystems, Inc. business.
35  * Permission to use, copy, modify, and distribute this
36  * software is freely granted, provided that this notice
37  * is preserved.
38  * =================================================================
39  */
40 
41 /**
42  * @file atan.c
43  * @brief Inverse tangent function
44  */
45 
46 #include "prim.h"
47 
48 static const double atanhi[] = {
49  4.63647609000806093515e-01, /* atan(0.5)hi 0x3FDDAC67, 0x0561BB4F */
50  7.85398163397448278999e-01, /* atan(1.0)hi 0x3FE921FB, 0x54442D18 */
51  9.82793723247329054082e-01, /* atan(1.5)hi 0x3FEF730B, 0xD281F69B */
52  1.57079632679489655800e+00, /* atan(inf)hi 0x3FF921FB, 0x54442D18 */
53 };
54 
55 static const double atanlo[] = {
56  2.26987774529616870924e-17, /* atan(0.5)lo 0x3C7A2B7F, 0x222F65E2 */
57  3.06161699786838301793e-17, /* atan(1.0)lo 0x3C81A626, 0x33145C07 */
58  1.39033110312309984516e-17, /* atan(1.5)lo 0x3C700788, 0x7AF0CBBD */
59  6.12323399573676603587e-17, /* atan(inf)lo 0x3C91A626, 0x33145C07 */
60 };
61 
62 static const double aT[] = {
63  3.33333333333329318027e-01, /* 0x3FD55555, 0x5555550D */
64  -1.99999999998764832476e-01, /* 0xBFC99999, 0x9998EBC4 */
65  1.42857142725034663711e-01, /* 0x3FC24924, 0x920083FF */
66  -1.11111104054623557880e-01, /* 0xBFBC71C6, 0xFE231671 */
67  9.09088713343650656196e-02, /* 0x3FB745CD, 0xC54C206E */
68  -7.69187620504482999495e-02, /* 0xBFB3B0F2, 0xAF749A6D */
69  6.66107313738753120669e-02, /* 0x3FB10D66, 0xA0D03D51 */
70  -5.83357013379057348645e-02, /* 0xBFADDE2D, 0x52DEFD9A */
71  4.97687799461593236017e-02, /* 0x3FA97B4B, 0x24760DEB */
72  -3.65315727442169155270e-02, /* 0xBFA2B444, 0x2C6A6C2F */
73  1.62858201153657823623e-02, /* 0x3F90AD3A, 0xE322DA11 */
74 };
75 
76 static const double
77  one = 1.0,
78  huge = 1.0e300;
79 
80 /**
81  * @brief Inverse tangent function
82  * @details
83  * <pre>
84  * Method
85  * 1. Reduce x to positive by atan(x) = -atan(-x).
86  * 2. According to the integer k=4t+0.25 chopped, t=x, the argument
87  * is further reduced to one of the following intervals and the
88  * arctangent of t is evaluated by the corresponding formula:
89  *
90  * [0,7/16] atan(x) = t-t^3*(a1+t^2*(a2+...(a10+t^2*a11)...)
91  * [7/16,11/16] atan(x) = atan(1/2) + atan( (t-0.5)/(1+t/2) )
92  * [11/16.19/16] atan(x) = atan( 1 ) + atan( (t-1)/(1+t) )
93  * [19/16,39/16] atan(x) = atan(3/2) + atan( (t-1.5)/(1+1.5t) )
94  * [39/16,INF] atan(x) = atan(INF) + atan( -1/t )
95  *
96  * Constants
97  * The hexadecimal values are the intended ones for the following
98  * constants. The decimal values may be used, provided that the
99  * compiler will convert from decimal to binary accurately enough
100  * to produce the hexadecimal values shown.
101  * </pre>
102  */
103 double atan(double x)
104 {
105  double w, s1, s2, z;
106  int32_t ix, hx, id;
107 
108  GET_HIGH_WORD(hx, x);
109  ix = hx & 0x7fffffff;
110  if (ix >= 0x44100000)
111  { /* if |x| >= 2^66 */
112  uint32_t low;
113 
114  GET_LOW_WORD(low, x);
115  if (ix > 0x7ff00000 || (ix == 0x7ff00000 && (low != 0)))
116  return NAN;
117  if (hx > 0)
118  return atanhi[3] + atanlo[3];
119  else
120  return -atanhi[3] - atanlo[3];
121  }
122  if (ix < 0x3fdc0000)
123  { /* |x| < 0.4375 */
124  if (ix < 0x3e200000)
125  { /* |x| < 2^-29 */
126  if (huge + x > one)
127  return x; /* raise inexact */
128  }
129  id = -1;
130  }
131  else
132  {
133  x = fabs(x);
134  if (ix < 0x3ff30000)
135  { /* |x| < 1.1875 */
136  if (ix < 0x3fe60000)
137  { /* 7/16 <=|x|<11/16 */
138  id = 0;
139  x = (2.0 * x - one) / (2.0 + x);
140  }
141  else
142  { /* 11/16<=|x|< 19/16 */
143  id = 1;
144  x = (x - one) / (x + one);
145  }
146  }
147  else
148  {
149  if (ix < 0x40038000)
150  { /* |x| < 2.4375 */
151  id = 2;
152  x = (x - 1.5) / (one + 1.5 * x);
153  }
154  else
155  { /* 2.4375 <= |x| < 2^66 */
156  id = 3;
157  x = -1.0 / x;
158  }
159  }
160  }
161  /* end of argument reduction */
162  z = x * x;
163  w = z * z;
164  /* break sum from i=0 to 10 aT[i]z**(i+1) into odd and even poly */
165  s1 = z * (aT[0] + w * (aT[2] + w * (aT[4] + w * (aT[6] + w * (aT[8] + w * aT[10])))));
166  s2 = w * (aT[1] + w * (aT[3] + w * (aT[5] + w * (aT[7] + w * aT[9]))));
167  if (id < 0)
168  {
169  return x - x * (s1 + s2);
170  }
171 
172  z = atanhi[id] - ((x * (s1 + s2) - atanlo[id]) - x);
173  return (hx < 0) ? -z : z;
174 }
#define GET_HIGH_WORD(i, d)
Get the more significant 32 bit int from a double.
Definition: prim.h:167
static const double huge
Definition: atan.c:78
#define GET_LOW_WORD(i, d)
Get the less significant 32 bit int from a double.
Definition: prim.h:177
static const double aT[]
Definition: atan.c:62
static const double atanhi[]
Definition: atan.c:48
#define NAN
Definition: mathr.h:53
double fabs(double x)
Returns the absolute value of x.
Definition: fabs.c:51
static const double one
Definition: atan.c:77
double atan(double x)
Inverse tangent function.
Definition: atan.c:103
static const double atanlo[]
Definition: atan.c:55