Mercurial > hg > CbC > CbC_gcc

comparison libquadmath/math/expq.c @ 111:04ced10e8804

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
gcc 7
author kono
date Fri, 27 Oct 2017 22:46:09 +0900
parents 561a7518be6b
children 1830386684a0
comparison
equal deleted inserted replaced
68:561a7518be6b 111:04ced10e8804
1 /* Quad-precision floating point e^x. 1 /* Quad-precision floating point e^x.
2 Copyright (C) 1999 Free Software Foundation, Inc. 2 Copyright (C) 1999-2017 Free Software Foundation, Inc.
3 This file is part of the GNU C Library. 3 This file is part of the GNU C Library.
4 Contributed by Jakub Jelinek <jj@ultra.linux.cz> 4 Contributed by Jakub Jelinek <jj@ultra.linux.cz>
5 Partly based on double-precision code 5 Partly based on double-precision code
6 by Geoffrey Keating <geoffk@ozemail.com.au> 6 by Geoffrey Keating <geoffk@ozemail.com.au>
7 7
28 # define USE_FENV_H 28 # define USE_FENV_H
29 # endif 29 # endif
30 #endif 30 #endif
31 31
32 32
33 /* __expl_table basically consists of four tables, T_EXPL_ARG{1,2} and 33 /* __expq_table basically consists of four tables, T_EXPL_ARG{1,2} and
34 T_EXPL_RES{1,2}. All tables use positive and negative indexes, the 0 points 34 T_EXPL_RES{1,2}. All tables use positive and negative indexes, the 0 points
35 are marked by T_EXPL_* defines. 35 are marked by T_EXPL_* defines.
36 For ARG1 and RES1 tables lets B be 89 and S 256.0, for ARG2 and RES2 B is 65 36 For ARG1 and RES1 tables lets B be 89 and S 256.0, for ARG2 and RES2 B is 65
37 and S 32768.0. 37 and S 32768.0.
38 These table have the property that, for all integers -B <= i <= B 38 These table have the property that, for all integers -B <= i <= B
39 expl(__expl_table[T_EXPL_ARGN+2*i]+__expl_table[T_EXPL_ARGN+2*i+1]+r) == 39 expl(__expq_table[T_EXPL_ARGN+2*i]+__expq_table[T_EXPL_ARGN+2*i+1]+r) ==
40 __expl_table[T_EXPL_RESN+i], __expl_table[T_EXPL_RESN+i] is some exact number 40 __expq_table[T_EXPL_RESN+i], __expq_table[T_EXPL_RESN+i] is some exact number
41 with the low 58 bits of the mantissa 0, 41 with the low 58 bits of the mantissa 0,
42 __expl_table[T_EXPL_ARGN+2*i] == i/S+s 42 __expq_table[T_EXPL_ARGN+2*i] == i/S+s
43 where absl(s) <= 2^-54 and absl(r) <= 2^-212. */ 43 where absl(s) <= 2^-54 and absl(r) <= 2^-212. */
44 44
45 static const __float128 __expl_table [] = { 45 static const __float128 __expq_table [] = {
46 -3.47656250000000000584188889839535373E-01Q, /* bffd640000000000002b1b04213cf000 */ 46 -3.47656250000000000584188889839535373E-01Q, /* bffd640000000000002b1b04213cf000 */
47 6.90417668990715641167244540876988960E-32Q, /* 3f97667c3fdb588a6ae1af8748357a17 */ 47 6.90417668990715641167244540876988960E-32Q, /* 3f97667c3fdb588a6ae1af8748357a17 */
48 -3.43749999999999981853132895957607418E-01Q, /* bffd5ffffffffffffac4ff5f4050b000 */ 48 -3.43749999999999981853132895957607418E-01Q, /* bffd5ffffffffffffac4ff5f4050b000 */
49 -7.16021898043268093462818380603370350E-33Q, /* bf94296c8219427edc1431ac2498583e */ 49 -7.16021898043268093462818380603370350E-33Q, /* bf94296c8219427edc1431ac2498583e */
50 -3.39843750000000013418643523138766329E-01Q, /* bffd5c000000000003de1f027a30e000 */ 50 -3.39843750000000013418643523138766329E-01Q, /* bffd5c000000000003de1f027a30e000 */
1073 1073
1074 /* 32768 */ 1074 /* 32768 */
1075 #define TWO15 C[11] 1075 #define TWO15 C[11]
1076 32768.0Q, 1076 32768.0Q,
1077 1077
1078 /* Chebyshev polynom coeficients for (exp(x)-1)/x */ 1078 /* Chebyshev polynom coefficients for (exp(x)-1)/x */
1079 #define P1 C[12] 1079 #define P1 C[12]
1080 #define P2 C[13] 1080 #define P2 C[13]
1081 #define P3 C[14] 1081 #define P3 C[14]
1082 #define P4 C[15] 1082 #define P4 C[15]
1083 #define P5 C[16] 1083 #define P5 C[16]
1120 t -= THREEp103; 1120 t -= THREEp103;
1121 1121
1122 /* Compute tval1 = t. */ 1122 /* Compute tval1 = t. */
1123 tval1 = (int) (t * TWO8); 1123 tval1 = (int) (t * TWO8);
1124 1124
1125 x -= __expl_table[T_EXPL_ARG1+2*tval1]; 1125 x -= __expq_table[T_EXPL_ARG1+2*tval1];
1126 xl -= __expl_table[T_EXPL_ARG1+2*tval1+1]; 1126 xl -= __expq_table[T_EXPL_ARG1+2*tval1+1];
1127 1127
1128 /* Calculate t/32768. */ 1128 /* Calculate t/32768. */
1129 t = x + THREEp96; 1129 t = x + THREEp96;
1130 t -= THREEp96; 1130 t -= THREEp96;
1131 1131
1132 /* Compute tval2 = t. */ 1132 /* Compute tval2 = t. */
1133 tval2 = (int) (t * TWO15); 1133 tval2 = (int) (t * TWO15);
1134 1134
1135 x -= __expl_table[T_EXPL_ARG2+2*tval2]; 1135 x -= __expq_table[T_EXPL_ARG2+2*tval2];
1136 xl -= __expl_table[T_EXPL_ARG2+2*tval2+1]; 1136 xl -= __expq_table[T_EXPL_ARG2+2*tval2+1];
1137 1137
1138 x = x + xl; 1138 x = x + xl;
1139 1139
1140 /* Compute ex2 = 2^n_0 e^(argtable[tval1]) e^(argtable[tval2]). */ 1140 /* Compute ex2 = 2^n_0 e^(argtable[tval1]) e^(argtable[tval2]). */
1141 ex2_u.value = __expl_table[T_EXPL_RES1 + tval1] 1141 ex2_u.value = __expq_table[T_EXPL_RES1 + tval1]
1142 * __expl_table[T_EXPL_RES2 + tval2]; 1142 * __expq_table[T_EXPL_RES2 + tval2];
1143 n_i = (int)n; 1143 n_i = (int)n;
1144 /* 'unsafe' is 1 iff n_1 != 0. */ 1144 /* 'unsafe' is 1 iff n_1 != 0. */
1145 unsafe = abs(n_i) >= -FLT128_MIN_EXP - 1; 1145 unsafe = abs(n_i) >= 15000;
1146 ex2_u.ieee.exponent += n_i >> unsafe; 1146 ex2_u.ieee.exponent += n_i >> unsafe;
1147 1147
1148 /* Compute scale = 2^n_1. */ 1148 /* Compute scale = 2^n_1. */
1149 scale_u.value = 1.0Q; 1149 scale_u.value = 1.0Q;
1150 scale_u.ieee.exponent += n_i - (n_i >> unsafe); 1150 scale_u.ieee.exponent += n_i - (n_i >> unsafe);
1177 #endif 1177 #endif
1178 #endif 1178 #endif
1179 ex3_u.d = (result - ex2_u.d) - x22 * ex2_u.d; 1179 ex3_u.d = (result - ex2_u.d) - x22 * ex2_u.d;
1180 ex2_u.d = result; 1180 ex2_u.d = result;
1181 ex3_u.ieee.exponent += LDBL_MANT_DIG + 15 + IEEE854_LONG_DOUBLE_BIAS 1181 ex3_u.ieee.exponent += LDBL_MANT_DIG + 15 + IEEE854_LONG_DOUBLE_BIAS
1182 - ex2_u.ieee.exponent; 1182 - ex2_u.ieee.exponent;
1183 n_i = abs (ex3_u.d); 1183 n_i = abs (ex3_u.d);
1184 n_i = (n_i + 1) / 2; 1184 n_i = (n_i + 1) / 2;
1185 #ifdef USE_FENV_H 1185 #ifdef USE_FENV_H
1186 fesetenv (&oldenv); 1186 fesetenv (&oldenv);
1187 #ifdef FE_TONEAREST 1187 #ifdef FE_TONEAREST
1194 } 1194 }
1195 */ 1195 */
1196 if (!unsafe) 1196 if (!unsafe)
1197 return result; 1197 return result;
1198 else 1198 else
1199 return result * scale_u.value; 1199 {
1200 result *= scale_u.value;
1201 math_check_force_underflow_nonneg (result);
1202 return result;
1203 }
1200 } 1204 }
1201 /* Exceptional cases: */ 1205 /* Exceptional cases: */
1202 else if (__builtin_isless (x, himark)) 1206 else if (__builtin_isless (x, himark))
1203 { 1207 {
1204 if (isinfq (x)) 1208 if (isinfq (x))