amath  1.8.5
Simple command line calculator
ktan.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/k_tan.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 ktan.c
43  * @brief Kernel tan function
44  */
45 
46 #include "prim.h"
47 
48 static const double xxx[] = {
49  3.33333333333334091986e-01, /* 3FD55555, 55555563 */
50  1.33333333333201242699e-01, /* 3FC11111, 1110FE7A */
51  5.39682539762260521377e-02, /* 3FABA1BA, 1BB341FE */
52  2.18694882948595424599e-02, /* 3F9664F4, 8406D637 */
53  8.86323982359930005737e-03, /* 3F8226E3, E96E8493 */
54  3.59207910759131235356e-03, /* 3F6D6D22, C9560328 */
55  1.45620945432529025516e-03, /* 3F57DBC8, FEE08315 */
56  5.88041240820264096874e-04, /* 3F4344D8, F2F26501 */
57  2.46463134818469906812e-04, /* 3F3026F7, 1A8D1068 */
58  7.81794442939557092300e-05, /* 3F147E88, A03792A6 */
59  7.14072491382608190305e-05, /* 3F12B80F, 32F0A7E9 */
60  -1.85586374855275456654e-05, /* BEF375CB, DB605373 */
61  2.59073051863633712884e-05, /* 3EFB2A70, 74BF7AD4 */
62  /* one */ 1.00000000000000000000e+00, /* 3FF00000, 00000000 */
63  /* pio4 */ 7.85398163397448278999e-01, /* 3FE921FB, 54442D18 */
64  /* pio4lo */ 3.06161699786838301793e-17 /* 3C81A626, 33145C07 */
65 };
66 
67 #define one xxx[13]
68 #define pio4 xxx[14]
69 #define pio4lo xxx[15]
70 #define T xxx
71 
72 /**
73  * @brief Kernel tan function
74  * @details
75  * <pre>
76  * Kernel tan function on [-pi/4, pi/4], pi/4 ~ 0.7854
77  * Input x is assumed to be bounded by ~pi/4 in magnitude.
78  * Input y is the tail of x.
79  * Input k indicates whether tan (if k = 1) or -1/tan (if k = -1) is returned.
80  *
81  * Algorithm
82  * 1. Since tan(-x) = -tan(x), we need only to consider positive x.
83  * 2. if x < 2^-28 (hx<0x3e300000 0), return x with inexact if x!=0.
84  * 3. tan(x) is approximated by a odd polynomial of degree 27 on
85  * [0,0.67434]
86  * 3 27
87  * tan(x) ~ x + T1*x + ... + T13*x
88  * where
89  *
90  * |tan(x) 2 4 26 | -59.2
91  * |----- - (1+T1*x +T2*x +.... +T13*x )| <= 2
92  * | x |
93  *
94  * Note: tan(x+y) = tan(x) + tan'(x)*y
95  * ~ tan(x) + (1+x*x)*y
96  * Therefore, for better accuracy in computing tan(x+y), let
97  * 3 2 2 2 2
98  * r = x *(T2+x *(T3+x *(...+x *(T12+x *T13))))
99  * then
100  * 3 2
101  * tan(x+y) = x + (T1*x + (x *(r+y)+y))
102  *
103  * 4. For x in [0.67434,pi/4], let y = pi/4 - x, then
104  * tan(x) = tan(pi/4-y) = (1-tan(y))/(1+tan(y))
105  * = 1 - 2*(tan(y) - (tan(y)^2)/(1+tan(y)))
106  * </pre>
107  */
108 double __kernel_tan(double x, double y, int iy)
109 {
110  double z, r, v, w, s;
111  double a, t;
112  int32_t ix, hx;
113  uint32_t low;
114 
115  // high word of x
116  GET_HIGH_WORD(hx, x);
117 
118  // high word of |x|
119  ix = hx & 0x7FFFFFFF;
120 
121  // x < 2**-28
122  if (ix < 0x3E300000)
123  {
124  // generate inexact
125  if ((int32_t)x == 0)
126  {
127  GET_LOW_WORD(low, x);
128 
129  if (((ix | low) | (iy + 1)) == 0)
130  {
131  return one / fabs(x);
132  }
133 
134  if (iy == 1)
135  {
136  return x;
137  }
138 
139  // compute -1 / (x+y) carefully
140  z = w = x + y;
141  SET_LOW_WORD(z, 0);
142  v = y - (z - x);
143  t = a = -one / w;
144  SET_LOW_WORD(t, 0);
145  s = one + t * z;
146  return t + a * (s + t * v);
147  }
148  }
149 
150  // |x| >= 0.6744
151  if (ix >= 0x3FE59428)
152  {
153  if (hx < 0)
154  {
155  x = -x;
156  y = -y;
157  }
158  z = pio4 - x;
159  w = pio4lo - y;
160  x = z + w;
161  y = 0.0;
162  }
163 
164  z = x * x;
165  w = z * z;
166 
167  /*
168  * Break x^5*(T[1]+x^2*T[2]+...) into
169  * x^5(T[1]+x^4*T[3]+...+x^20*T[11]) +
170  * x^5(x^2*(T[2]+x^4*T[4]+...+x^22*[T12]))
171  */
172  r = T[1] + w * (T[3] + w * (T[5] + w * (T[7] + w * (T[9] + w * T[11]))));
173  v = z * (T[2] + w * (T[4] + w * (T[6] + w * (T[8] + w * (T[10] + w * T[12])))));
174  s = z * x;
175  r = y + z * (s * (r + v) + y);
176  r += T[0] * s;
177  w = x + r;
178 
179  if (ix >= 0x3FE59428)
180  {
181  v = (double)iy;
182  return (double)(1 - ((hx >> 30) & 2)) *
183  (v - 2.0 * (x - (w * w / (w + v) - r)));
184  }
185 
186  if (iy == 1)
187  {
188  return w;
189  }
190 
191  // compute -1.0 / (x+r) accurately
192  z = w;
193  SET_LOW_WORD(z, 0);
194  v = r - (z - x); // z+v = r+x
195  t = a = -1.0 / w; // a = -1.0/w
196  SET_LOW_WORD(t, 0);
197  s = 1.0 + t * z;
198  return t + a * (s + t * v);
199 }
#define pio4
Definition: ktan.c:68
#define GET_HIGH_WORD(i, d)
Get the more significant 32 bit int from a double.
Definition: prim.h:167
#define pio4lo
Definition: ktan.c:69
#define GET_LOW_WORD(i, d)
Get the less significant 32 bit int from a double.
Definition: prim.h:177
static const double xxx[]
Definition: ktan.c:48
double fabs(double x)
Returns the absolute value of x.
Definition: fabs.c:51
double __kernel_tan(double x, double y, int iy)
Kernel tan function.
Definition: ktan.c:108
#define T
Definition: ktan.c:70
#define SET_LOW_WORD(d, v)
Set the less significant 32 bits of a double from an int.
Definition: prim.h:209
#define one
Definition: ktan.c:67