Mercurial > hg > CbC > CbC_gcc
comparison libquadmath/math/tanq_kernel.c @ 145:1830386684a0
gcc-9.2.0
author | anatofuz |
---|---|
date | Thu, 13 Feb 2020 11:34:05 +0900 |
parents | |
children |
comparison
equal
deleted
inserted
replaced
131:84e7813d76e9 | 145:1830386684a0 |
---|---|
1 /* | |
2 * ==================================================== | |
3 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. | |
4 * | |
5 * Developed at SunPro, a Sun Microsystems, Inc. business. | |
6 * Permission to use, copy, modify, and distribute this | |
7 * software is freely granted, provided that this notice | |
8 * is preserved. | |
9 * ==================================================== | |
10 */ | |
11 | |
12 /* | |
13 Long double expansions are | |
14 Copyright (C) 2001 Stephen L. Moshier <moshier@na-net.ornl.gov> | |
15 and are incorporated herein by permission of the author. The author | |
16 reserves the right to distribute this material elsewhere under different | |
17 copying permissions. These modifications are distributed here under | |
18 the following terms: | |
19 | |
20 This library is free software; you can redistribute it and/or | |
21 modify it under the terms of the GNU Lesser General Public | |
22 License as published by the Free Software Foundation; either | |
23 version 2.1 of the License, or (at your option) any later version. | |
24 | |
25 This library is distributed in the hope that it will be useful, | |
26 but WITHOUT ANY WARRANTY; without even the implied warranty of | |
27 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU | |
28 Lesser General Public License for more details. | |
29 | |
30 You should have received a copy of the GNU Lesser General Public | |
31 License along with this library; if not, see | |
32 <http://www.gnu.org/licenses/>. */ | |
33 | |
34 /* __quadmath_kernel_tanq( x, y, k ) | |
35 * kernel tan function on [-pi/4, pi/4], pi/4 ~ 0.7854 | |
36 * Input x is assumed to be bounded by ~pi/4 in magnitude. | |
37 * Input y is the tail of x. | |
38 * Input k indicates whether tan (if k=1) or | |
39 * -1/tan (if k= -1) is returned. | |
40 * | |
41 * Algorithm | |
42 * 1. Since tan(-x) = -tan(x), we need only to consider positive x. | |
43 * 2. if x < 2^-57, return x with inexact if x!=0. | |
44 * 3. tan(x) is approximated by a rational form x + x^3 / 3 + x^5 R(x^2) | |
45 * on [0,0.67433]. | |
46 * | |
47 * Note: tan(x+y) = tan(x) + tan'(x)*y | |
48 * ~ tan(x) + (1+x*x)*y | |
49 * Therefore, for better accuracy in computing tan(x+y), let | |
50 * r = x^3 * R(x^2) | |
51 * then | |
52 * tan(x+y) = x + (x^3 / 3 + (x^2 *(r+y)+y)) | |
53 * | |
54 * 4. For x in [0.67433,pi/4], let y = pi/4 - x, then | |
55 * tan(x) = tan(pi/4-y) = (1-tan(y))/(1+tan(y)) | |
56 * = 1 - 2*(tan(y) - (tan(y)^2)/(1+tan(y))) | |
57 */ | |
58 | |
59 #include "quadmath-imp.h" | |
60 | |
61 static const __float128 | |
62 one = 1, | |
63 pio4hi = 7.8539816339744830961566084581987569936977E-1Q, | |
64 pio4lo = 2.1679525325309452561992610065108379921906E-35Q, | |
65 | |
66 /* tan x = x + x^3 / 3 + x^5 T(x^2)/U(x^2) | |
67 0 <= x <= 0.6743316650390625 | |
68 Peak relative error 8.0e-36 */ | |
69 TH = 3.333333333333333333333333333333333333333E-1Q, | |
70 T0 = -1.813014711743583437742363284336855889393E7Q, | |
71 T1 = 1.320767960008972224312740075083259247618E6Q, | |
72 T2 = -2.626775478255838182468651821863299023956E4Q, | |
73 T3 = 1.764573356488504935415411383687150199315E2Q, | |
74 T4 = -3.333267763822178690794678978979803526092E-1Q, | |
75 | |
76 U0 = -1.359761033807687578306772463253710042010E8Q, | |
77 U1 = 6.494370630656893175666729313065113194784E7Q, | |
78 U2 = -4.180787672237927475505536849168729386782E6Q, | |
79 U3 = 8.031643765106170040139966622980914621521E4Q, | |
80 U4 = -5.323131271912475695157127875560667378597E2Q; | |
81 /* 1.000000000000000000000000000000000000000E0 */ | |
82 | |
83 | |
84 __float128 | |
85 __quadmath_kernel_tanq (__float128 x, __float128 y, int iy) | |
86 { | |
87 __float128 z, r, v, w, s; | |
88 int32_t ix, sign; | |
89 ieee854_float128 u, u1; | |
90 | |
91 u.value = x; | |
92 ix = u.words32.w0 & 0x7fffffff; | |
93 if (ix < 0x3fc60000) /* x < 2**-57 */ | |
94 { | |
95 if ((int) x == 0) | |
96 { /* generate inexact */ | |
97 if ((ix | u.words32.w1 | u.words32.w2 | u.words32.w3 | |
98 | (iy + 1)) == 0) | |
99 return one / fabsq (x); | |
100 else if (iy == 1) | |
101 { | |
102 math_check_force_underflow (x); | |
103 return x; | |
104 } | |
105 else | |
106 return -one / x; | |
107 } | |
108 } | |
109 if (ix >= 0x3ffe5942) /* |x| >= 0.6743316650390625 */ | |
110 { | |
111 if ((u.words32.w0 & 0x80000000) != 0) | |
112 { | |
113 x = -x; | |
114 y = -y; | |
115 sign = -1; | |
116 } | |
117 else | |
118 sign = 1; | |
119 z = pio4hi - x; | |
120 w = pio4lo - y; | |
121 x = z + w; | |
122 y = 0.0; | |
123 } | |
124 z = x * x; | |
125 r = T0 + z * (T1 + z * (T2 + z * (T3 + z * T4))); | |
126 v = U0 + z * (U1 + z * (U2 + z * (U3 + z * (U4 + z)))); | |
127 r = r / v; | |
128 | |
129 s = z * x; | |
130 r = y + z * (s * r + y); | |
131 r += TH * s; | |
132 w = x + r; | |
133 if (ix >= 0x3ffe5942) | |
134 { | |
135 v = (__float128) iy; | |
136 w = (v - 2.0 * (x - (w * w / (w + v) - r))); | |
137 /* SIGN is set for arguments that reach this code, but not | |
138 otherwise, resulting in warnings that it may be used | |
139 uninitialized although in the cases where it is used it has | |
140 always been set. */ | |
141 | |
142 | |
143 if (sign < 0) | |
144 w = -w; | |
145 | |
146 return w; | |
147 } | |
148 if (iy == 1) | |
149 return w; | |
150 else | |
151 { /* if allow error up to 2 ulp, | |
152 simply return -1.0/(x+r) here */ | |
153 /* compute -1.0/(x+r) accurately */ | |
154 u1.value = w; | |
155 u1.words32.w2 = 0; | |
156 u1.words32.w3 = 0; | |
157 v = r - (u1.value - x); /* u1+v = r+x */ | |
158 z = -1.0 / w; | |
159 u.value = z; | |
160 u.words32.w2 = 0; | |
161 u.words32.w3 = 0; | |
162 s = 1.0 + u.value * u1.value; | |
163 return u.value + z * (s + u.value * v); | |
164 } | |
165 } |