Mercurial > hg > CbC > CbC_gcc

comparison libquadmath/math/log1pq.c @ 145:1830386684a0

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
gcc-9.2.0
author anatofuz
date Thu, 13 Feb 2020 11:34:05 +0900
parents 04ced10e8804
children
comparison
equal deleted inserted replaced
131:84e7813d76e9 145:1830386684a0
1 /* log1pl.c 1 /* log1pq.c
2 * 2 *
3 * Relative error logarithm 3 * Relative error logarithm
4 * Natural logarithm of 1+x for __float128 precision 4 * Natural logarithm of 1+x, 128-bit long double precision
5 * 5 *
6 * 6 *
7 * 7 *
8 * SYNOPSIS: 8 * SYNOPSIS:
9 * 9 *
10 * __float128 x, y, log1pl(); 10 * long double x, y, log1pq();
11 * 11 *
12 * y = log1pq( x ); 12 * y = log1pq( x );
13 * 13 *
14 * 14 *
15 * 15 *
47 but WITHOUT ANY WARRANTY; without even the implied warranty of 47 but WITHOUT ANY WARRANTY; without even the implied warranty of
48 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 48 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
49 Lesser General Public License for more details. 49 Lesser General Public License for more details.
50 50
51 You should have received a copy of the GNU Lesser General Public 51 You should have received a copy of the GNU Lesser General Public
52 License along with this library; if not, write to the Free Software 52 License along with this library; if not, see
53 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */ 53 <http://www.gnu.org/licenses/>. */
54
55 54
56 #include "quadmath-imp.h" 55 #include "quadmath-imp.h"
57 56
58 /* Coefficients for log(1+x) = x - x^2 / 2 + x^3 P(x)/Q(x) 57 /* Coefficients for log(1+x) = x - x^2 / 2 + x^3 P(x)/Q(x)
59 * 1/sqrt(2) <= 1+x < sqrt(2) 58 * 1/sqrt(2) <= 1+x < sqrt(2)
72 P4 = 2.854829159639697837788887080758954924001E5Q, 71 P4 = 2.854829159639697837788887080758954924001E5Q,
73 P3 = 3.007007295140399532324943111654767187848E5Q, 72 P3 = 3.007007295140399532324943111654767187848E5Q,
74 P2 = 2.014652742082537582487669938141683759923E5Q, 73 P2 = 2.014652742082537582487669938141683759923E5Q,
75 P1 = 7.771154681358524243729929227226708890930E4Q, 74 P1 = 7.771154681358524243729929227226708890930E4Q,
76 P0 = 1.313572404063446165910279910527789794488E4Q, 75 P0 = 1.313572404063446165910279910527789794488E4Q,
77 /* Q12 = 1.000000000000000000000000000000000000000E0Q, */ 76 /* Q12 = 1.000000000000000000000000000000000000000E0L, */
78 Q11 = 4.839208193348159620282142911143429644326E1Q, 77 Q11 = 4.839208193348159620282142911143429644326E1Q,
79 Q10 = 9.104928120962988414618126155557301584078E2Q, 78 Q10 = 9.104928120962988414618126155557301584078E2Q,
80 Q9 = 9.147150349299596453976674231612674085381E3Q, 79 Q9 = 9.147150349299596453976674231612674085381E3Q,
81 Q8 = 5.605842085972455027590989944010492125825E4Q, 80 Q8 = 5.605842085972455027590989944010492125825E4Q,
82 Q7 = 2.248234257620569139969141618556349415120E5Q, 81 Q7 = 2.248234257620569139969141618556349415120E5Q,
99 R4 = 8.057002716646055371965756206836056074715E1Q, 98 R4 = 8.057002716646055371965756206836056074715E1Q,
100 R3 = -2.024301798136027039250415126250455056397E3Q, 99 R3 = -2.024301798136027039250415126250455056397E3Q,
101 R2 = 2.048819892795278657810231591630928516206E4Q, 100 R2 = 2.048819892795278657810231591630928516206E4Q,
102 R1 = -8.977257995689735303686582344659576526998E4Q, 101 R1 = -8.977257995689735303686582344659576526998E4Q,
103 R0 = 1.418134209872192732479751274970992665513E5Q, 102 R0 = 1.418134209872192732479751274970992665513E5Q,
104 /* S6 = 1.000000000000000000000000000000000000000E0Q, */ 103 /* S6 = 1.000000000000000000000000000000000000000E0L, */
105 S5 = -1.186359407982897997337150403816839480438E2Q, 104 S5 = -1.186359407982897997337150403816839480438E2Q,
106 S4 = 3.998526750980007367835804959888064681098E3Q, 105 S4 = 3.998526750980007367835804959888064681098E3Q,
107 S3 = -5.748542087379434595104154610899551484314E4Q, 106 S3 = -5.748542087379434595104154610899551484314E4Q,
108 S2 = 4.001557694070773974936904547424676279307E5Q, 107 S2 = 4.001557694070773974936904547424676279307E5Q,
109 S1 = -1.332535117259762928288745111081235577029E6Q, 108 S1 = -1.332535117259762928288745111081235577029E6Q,
112 /* C1 + C2 = ln 2 */ 111 /* C1 + C2 = ln 2 */
113 static const __float128 C1 = 6.93145751953125E-1Q; 112 static const __float128 C1 = 6.93145751953125E-1Q;
114 static const __float128 C2 = 1.428606820309417232121458176568075500134E-6Q; 113 static const __float128 C2 = 1.428606820309417232121458176568075500134E-6Q;
115 114
116 static const __float128 sqrth = 0.7071067811865475244008443621048490392848Q; 115 static const __float128 sqrth = 0.7071067811865475244008443621048490392848Q;
117 static const __float128 zero = 0.0Q; 116 /* ln (2^16384 * (1 - 2^-113)) */
118 117 static const __float128 zero = 0;
119 118
120 __float128 119 __float128
121 log1pq (__float128 xm1) 120 log1pq (__float128 xm1)
122 { 121 {
123 __float128 x, y, z, r, s; 122 __float128 x, y, z, r, s;
138 137
139 if ((hx & 0x7fffffff) < 0x3f8e0000) 138 if ((hx & 0x7fffffff) < 0x3f8e0000)
140 { 139 {
141 math_check_force_underflow (xm1); 140 math_check_force_underflow (xm1);
142 if ((int) xm1 == 0) 141 if ((int) xm1 == 0)
143 return xm1; 142 return xm1;
144 } 143 }
145 144
146 if (xm1 >= 0x1p113Q) 145 if (xm1 >= 0x1p113Q)
147 x = xm1; 146 x = xm1;
148 else 147 else
149 x = xm1 + 1.0Q; 148 x = xm1 + 1;
150 149
151 /* log1p(-1) = -inf */ 150 /* log1p(-1) = -inf */
152 if (x <= 0.0Q) 151 if (x <= 0)
153 { 152 {
154 if (x == 0.0Q) 153 if (x == 0)
155 return (-1.0Q / zero); /* log1p(-1) = -inf */ 154 return (-1 / zero); /* log1p(-1) = -inf */
156 else 155 else
157 return (zero / (x - x)); 156 return (zero / (x - x));
158 } 157 }
159 158
160 /* Separate mantissa from exponent. */ 159 /* Separate mantissa from exponent. */
205 204
206 if (x < sqrth) 205 if (x < sqrth)
207 { 206 {
208 e -= 1; 207 e -= 1;
209 if (e != 0) 208 if (e != 0)
210 x = 2.0Q * x - 1.0Q; /* 2x - 1 */ 209 x = 2 * x - 1; /* 2x - 1 */
211 else 210 else
212 x = xm1; 211 x = xm1;
213 } 212 }
214 else 213 else
215 { 214 {
216 if (e != 0) 215 if (e != 0)
217 x = x - 1.0Q; 216 x = x - 1;
218 else 217 else
219 x = xm1; 218 x = xm1;
220 } 219 }
221 z = x * x; 220 z = x * x;
222 r = (((((((((((P12 * x 221 r = (((((((((((P12 * x