1/* expm1q.c 2 * 3 * Exponential function, minus 1 4 * 128-bit long double precision 5 * 6 * 7 * 8 * SYNOPSIS: 9 * 10 * long double x, y, expm1q(); 11 * 12 * y = expm1q( x ); 13 * 14 * 15 * 16 * DESCRIPTION: 17 * 18 * Returns e (2.71828...) raised to the x power, minus one. 19 * 20 * Range reduction is accomplished by separating the argument 21 * into an integer k and fraction f such that 22 * 23 * x k f 24 * e = 2 e. 25 * 26 * An expansion x + .5 x^2 + x^3 R(x) approximates exp(f) - 1 27 * in the basic range [-0.5 ln 2, 0.5 ln 2]. 28 * 29 * 30 * ACCURACY: 31 * 32 * Relative error: 33 * arithmetic domain # trials peak rms 34 * IEEE -79,+MAXLOG 100,000 1.7e-34 4.5e-35 35 * 36 */ 37 38/* Copyright 2001 by Stephen L. Moshier 39 40 This library is free software; you can redistribute it and/or 41 modify it under the terms of the GNU Lesser General Public 42 License as published by the Free Software Foundation; either 43 version 2.1 of the License, or (at your option) any later version. 44 45 This library is distributed in the hope that it will be useful, 46 but WITHOUT ANY WARRANTY; without even the implied warranty of 47 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 48 Lesser General Public License for more details. 49 50 You should have received a copy of the GNU Lesser General Public 51 License along with this library; if not, see 52 <http://www.gnu.org/licenses/>. */ 53 54#include "quadmath-imp.h" 55 56/* exp(x) - 1 = x + 0.5 x^2 + x^3 P(x)/Q(x) 57 -.5 ln 2 < x < .5 ln 2 58 Theoretical peak relative error = 8.1e-36 */ 59 60static const __float128 61 P0 = 2.943520915569954073888921213330863757240E8Q, 62 P1 = -5.722847283900608941516165725053359168840E7Q, 63 P2 = 8.944630806357575461578107295909719817253E6Q, 64 P3 = -7.212432713558031519943281748462837065308E5Q, 65 P4 = 4.578962475841642634225390068461943438441E4Q, 66 P5 = -1.716772506388927649032068540558788106762E3Q, 67 P6 = 4.401308817383362136048032038528753151144E1Q, 68 P7 = -4.888737542888633647784737721812546636240E-1Q, 69 Q0 = 1.766112549341972444333352727998584753865E9Q, 70 Q1 = -7.848989743695296475743081255027098295771E8Q, 71 Q2 = 1.615869009634292424463780387327037251069E8Q, 72 Q3 = -2.019684072836541751428967854947019415698E7Q, 73 Q4 = 1.682912729190313538934190635536631941751E6Q, 74 Q5 = -9.615511549171441430850103489315371768998E4Q, 75 Q6 = 3.697714952261803935521187272204485251835E3Q, 76 Q7 = -8.802340681794263968892934703309274564037E1Q, 77 /* Q8 = 1.000000000000000000000000000000000000000E0 */ 78/* C1 + C2 = ln 2 */ 79 80 C1 = 6.93145751953125E-1Q, 81 C2 = 1.428606820309417232121458176568075500134E-6Q, 82/* ln 2^-114 */ 83 minarg = -7.9018778583833765273564461846232128760607E1Q, big = 1e4932Q; 84 85 86__float128 87expm1q (__float128 x) 88{ 89 __float128 px, qx, xx; 90 int32_t ix, sign; 91 ieee854_float128 u; 92 int k; 93 94 /* Detect infinity and NaN. */ 95 u.value = x; 96 ix = u.words32.w0; 97 sign = ix & 0x80000000; 98 ix &= 0x7fffffff; 99 if (!sign && ix >= 0x40060000) 100 { 101 /* If num is positive and exp >= 6 use plain exp. */ 102 return expq (x); 103 } 104 if (ix >= 0x7fff0000) 105 { 106 /* Infinity (which must be negative infinity). */ 107 if (((ix & 0xffff) | u.words32.w1 | u.words32.w2 | u.words32.w3) == 0) 108 return -1; 109 /* NaN. Invalid exception if signaling. */ 110 return x + x; 111 } 112 113 /* expm1(+- 0) = +- 0. */ 114 if ((ix == 0) && (u.words32.w1 | u.words32.w2 | u.words32.w3) == 0) 115 return x; 116 117 /* Minimum value. */ 118 if (x < minarg) 119 return (4.0/big - 1); 120 121 /* Avoid internal underflow when result does not underflow, while 122 ensuring underflow (without returning a zero of the wrong sign) 123 when the result does underflow. */ 124 if (fabsq (x) < 0x1p-113Q) 125 { 126 math_check_force_underflow (x); 127 return x; 128 } 129 130 /* Express x = ln 2 (k + remainder), remainder not exceeding 1/2. */ 131 xx = C1 + C2; /* ln 2. */ 132 px = floorq (0.5 + x / xx); 133 k = px; 134 /* remainder times ln 2 */ 135 x -= px * C1; 136 x -= px * C2; 137 138 /* Approximate exp(remainder ln 2). */ 139 px = (((((((P7 * x 140 + P6) * x 141 + P5) * x + P4) * x + P3) * x + P2) * x + P1) * x + P0) * x; 142 143 qx = (((((((x 144 + Q7) * x 145 + Q6) * x + Q5) * x + Q4) * x + Q3) * x + Q2) * x + Q1) * x + Q0; 146 147 xx = x * x; 148 qx = x + (0.5 * xx + xx * px / qx); 149 150 /* exp(x) = exp(k ln 2) exp(remainder ln 2) = 2^k exp(remainder ln 2). 151 152 We have qx = exp(remainder ln 2) - 1, so 153 exp(x) - 1 = 2^k (qx + 1) - 1 154 = 2^k qx + 2^k - 1. */ 155 156 px = ldexpq (1, k); 157 x = px * qx + (px - 1.0); 158 return x; 159} 160