1 | /* expm1l.c |
---|---|

2 | * |

3 | * Exponential function, minus 1 |

4 | * 128-bit __float128 precision |

5 | * |

6 | * |

7 | * |

8 | * SYNOPSIS: |

9 | * |

10 | * __float128 x, y, expm1l(); |

11 | * |

12 | * y = expm1l( 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, write to the Free Software |

52 | Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */ |

53 | |

54 | |

55 | |

56 | #include <errno.h> |

57 | #include "quadmath-imp.h" |

58 | |

59 | /* exp(x) - 1 = x + 0.5 x^2 + x^3 P(x)/Q(x) |

60 | -.5 ln 2 < x < .5 ln 2 |

61 | Theoretical peak relative error = 8.1e-36 */ |

62 | |

63 | static const __float128 |

64 | P0 = 2.943520915569954073888921213330863757240E8Q, |

65 | P1 = -5.722847283900608941516165725053359168840E7Q, |

66 | P2 = 8.944630806357575461578107295909719817253E6Q, |

67 | P3 = -7.212432713558031519943281748462837065308E5Q, |

68 | P4 = 4.578962475841642634225390068461943438441E4Q, |

69 | P5 = -1.716772506388927649032068540558788106762E3Q, |

70 | P6 = 4.401308817383362136048032038528753151144E1Q, |

71 | P7 = -4.888737542888633647784737721812546636240E-1Q, |

72 | Q0 = 1.766112549341972444333352727998584753865E9Q, |

73 | Q1 = -7.848989743695296475743081255027098295771E8Q, |

74 | Q2 = 1.615869009634292424463780387327037251069E8Q, |

75 | Q3 = -2.019684072836541751428967854947019415698E7Q, |

76 | Q4 = 1.682912729190313538934190635536631941751E6Q, |

77 | Q5 = -9.615511549171441430850103489315371768998E4Q, |

78 | Q6 = 3.697714952261803935521187272204485251835E3Q, |

79 | Q7 = -8.802340681794263968892934703309274564037E1Q, |

80 | /* Q8 = 1.000000000000000000000000000000000000000E0 */ |

81 | /* C1 + C2 = ln 2 */ |

82 | |

83 | C1 = 6.93145751953125E-1Q, |

84 | C2 = 1.428606820309417232121458176568075500134E-6Q, |

85 | /* ln 2^-114 */ |

86 | minarg = -7.9018778583833765273564461846232128760607E1Q; |

87 | |

88 | |

89 | __float128 |

90 | expm1q (__float128 x) |

91 | { |

92 | __float128 px, qx, xx; |

93 | int32_t ix, sign; |

94 | ieee854_float128 u; |

95 | int k; |

96 | |

97 | /* Detect infinity and NaN. */ |

98 | u.value = x; |

99 | ix = u.words32.w0; |

100 | sign = ix & 0x80000000; |

101 | ix &= 0x7fffffff; |

102 | if (!sign && ix >= 0x40060000) |

103 | { |

104 | /* If num is positive and exp >= 6 use plain exp. */ |

105 | return expq (x); |

106 | } |

107 | if (ix >= 0x7fff0000) |

108 | { |

109 | /* Infinity (which must be negative infinity). */ |

110 | if (((ix & 0xffff) | u.words32.w1 | u.words32.w2 | u.words32.w3) == 0) |

111 | return -1.0Q; |

112 | /* NaN. Invalid exception if signaling. */ |

113 | return x + x; |

114 | } |

115 | |

116 | /* expm1(+- 0) = +- 0. */ |

117 | if ((ix == 0) && (u.words32.w1 | u.words32.w2 | u.words32.w3) == 0) |

118 | return x; |

119 | |

120 | /* Minimum value. */ |

121 | if (x < minarg) |

122 | return (4.0/HUGE_VALQ - 1.0Q); |

123 | |

124 | /* Avoid internal underflow when result does not underflow, while |

125 | ensuring underflow (without returning a zero of the wrong sign) |

126 | when the result does underflow. */ |

127 | if (fabsq (x) < 0x1p-113Q) |

128 | { |

129 | math_check_force_underflow (x); |

130 | return x; |

131 | } |

132 | |

133 | /* Express x = ln 2 (k + remainder), remainder not exceeding 1/2. */ |

134 | xx = C1 + C2; /* ln 2. */ |

135 | px = floorq (0.5 + x / xx); |

136 | k = px; |

137 | /* remainder times ln 2 */ |

138 | x -= px * C1; |

139 | x -= px * C2; |

140 | |

141 | /* Approximate exp(remainder ln 2). */ |

142 | px = (((((((P7 * x |

143 | + P6) * x |

144 | + P5) * x + P4) * x + P3) * x + P2) * x + P1) * x + P0) * x; |

145 | |

146 | qx = (((((((x |

147 | + Q7) * x |

148 | + Q6) * x + Q5) * x + Q4) * x + Q3) * x + Q2) * x + Q1) * x + Q0; |

149 | |

150 | xx = x * x; |

151 | qx = x + (0.5 * xx + xx * px / qx); |

152 | |

153 | /* exp(x) = exp(k ln 2) exp(remainder ln 2) = 2^k exp(remainder ln 2). |

154 | |

155 | We have qx = exp(remainder ln 2) - 1, so |

156 | exp(x) - 1 = 2^k (qx + 1) - 1 |

157 | = 2^k qx + 2^k - 1. */ |

158 | |

159 | px = ldexpq (1.0Q, k); |

160 | x = px * qx + (px - 1.0); |

161 | return x; |

162 | } |

163 |