# Python code coverage for Modules/_math.c

# | count | content |
---|---|---|

1 | n/a | /* Definitions of some C99 math library functions, for those platforms |

2 | n/a | that don't implement these functions already. */ |

3 | n/a | |

4 | n/a | #include "Python.h" |

5 | n/a | #include <float.h> |

6 | n/a | #include "_math.h" |

7 | n/a | |

8 | n/a | /* The following copyright notice applies to the original |

9 | n/a | implementations of acosh, asinh and atanh. */ |

10 | n/a | |

11 | n/a | /* |

12 | n/a | * ==================================================== |

13 | n/a | * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. |

14 | n/a | * |

15 | n/a | * Developed at SunPro, a Sun Microsystems, Inc. business. |

16 | n/a | * Permission to use, copy, modify, and distribute this |

17 | n/a | * software is freely granted, provided that this notice |

18 | n/a | * is preserved. |

19 | n/a | * ==================================================== |

20 | n/a | */ |

21 | n/a | |

22 | n/a | #if !defined(HAVE_ACOSH) || !defined(HAVE_ASINH) |

23 | n/a | static const double ln2 = 6.93147180559945286227E-01; |

24 | n/a | static const double two_pow_p28 = 268435456.0; /* 2**28 */ |

25 | n/a | #endif |

26 | n/a | #if !defined(HAVE_ASINH) || !defined(HAVE_ATANH) |

27 | n/a | static const double two_pow_m28 = 3.7252902984619141E-09; /* 2**-28 */ |

28 | n/a | #endif |

29 | n/a | #if !defined(HAVE_ATANH) && !defined(Py_NAN) |

30 | n/a | static const double zero = 0.0; |

31 | n/a | #endif |

32 | n/a | |

33 | n/a | |

34 | n/a | #ifndef HAVE_ACOSH |

35 | n/a | /* acosh(x) |

36 | n/a | * Method : |

37 | n/a | * Based on |

38 | n/a | * acosh(x) = log [ x + sqrt(x*x-1) ] |

39 | n/a | * we have |

40 | n/a | * acosh(x) := log(x)+ln2, if x is large; else |

41 | n/a | * acosh(x) := log(2x-1/(sqrt(x*x-1)+x)) if x>2; else |

42 | n/a | * acosh(x) := log1p(t+sqrt(2.0*t+t*t)); where t=x-1. |

43 | n/a | * |

44 | n/a | * Special cases: |

45 | n/a | * acosh(x) is NaN with signal if x<1. |

46 | n/a | * acosh(NaN) is NaN without signal. |

47 | n/a | */ |

48 | n/a | |

49 | n/a | double |

50 | n/a | _Py_acosh(double x) |

51 | n/a | { |

52 | n/a | if (Py_IS_NAN(x)) { |

53 | n/a | return x+x; |

54 | n/a | } |

55 | n/a | if (x < 1.) { /* x < 1; return a signaling NaN */ |

56 | n/a | errno = EDOM; |

57 | n/a | #ifdef Py_NAN |

58 | n/a | return Py_NAN; |

59 | n/a | #else |

60 | n/a | return (x-x)/(x-x); |

61 | n/a | #endif |

62 | n/a | } |

63 | n/a | else if (x >= two_pow_p28) { /* x > 2**28 */ |

64 | n/a | if (Py_IS_INFINITY(x)) { |

65 | n/a | return x+x; |

66 | n/a | } |

67 | n/a | else { |

68 | n/a | return log(x) + ln2; /* acosh(huge)=log(2x) */ |

69 | n/a | } |

70 | n/a | } |

71 | n/a | else if (x == 1.) { |

72 | n/a | return 0.0; /* acosh(1) = 0 */ |

73 | n/a | } |

74 | n/a | else if (x > 2.) { /* 2 < x < 2**28 */ |

75 | n/a | double t = x * x; |

76 | n/a | return log(2.0 * x - 1.0 / (x + sqrt(t - 1.0))); |

77 | n/a | } |

78 | n/a | else { /* 1 < x <= 2 */ |

79 | n/a | double t = x - 1.0; |

80 | n/a | return m_log1p(t + sqrt(2.0 * t + t * t)); |

81 | n/a | } |

82 | n/a | } |

83 | n/a | #endif /* HAVE_ACOSH */ |

84 | n/a | |

85 | n/a | |

86 | n/a | #ifndef HAVE_ASINH |

87 | n/a | /* asinh(x) |

88 | n/a | * Method : |

89 | n/a | * Based on |

90 | n/a | * asinh(x) = sign(x) * log [ |x| + sqrt(x*x+1) ] |

91 | n/a | * we have |

92 | n/a | * asinh(x) := x if 1+x*x=1, |

93 | n/a | * := sign(x)*(log(x)+ln2)) for large |x|, else |

94 | n/a | * := sign(x)*log(2|x|+1/(|x|+sqrt(x*x+1))) if|x|>2, else |

95 | n/a | * := sign(x)*log1p(|x| + x^2/(1 + sqrt(1+x^2))) |

96 | n/a | */ |

97 | n/a | |

98 | n/a | double |

99 | n/a | _Py_asinh(double x) |

100 | n/a | { |

101 | n/a | double w; |

102 | n/a | double absx = fabs(x); |

103 | n/a | |

104 | n/a | if (Py_IS_NAN(x) || Py_IS_INFINITY(x)) { |

105 | n/a | return x+x; |

106 | n/a | } |

107 | n/a | if (absx < two_pow_m28) { /* |x| < 2**-28 */ |

108 | n/a | return x; /* return x inexact except 0 */ |

109 | n/a | } |

110 | n/a | if (absx > two_pow_p28) { /* |x| > 2**28 */ |

111 | n/a | w = log(absx) + ln2; |

112 | n/a | } |

113 | n/a | else if (absx > 2.0) { /* 2 < |x| < 2**28 */ |

114 | n/a | w = log(2.0 * absx + 1.0 / (sqrt(x * x + 1.0) + absx)); |

115 | n/a | } |

116 | n/a | else { /* 2**-28 <= |x| < 2= */ |

117 | n/a | double t = x*x; |

118 | n/a | w = m_log1p(absx + t / (1.0 + sqrt(1.0 + t))); |

119 | n/a | } |

120 | n/a | return copysign(w, x); |

121 | n/a | |

122 | n/a | } |

123 | n/a | #endif /* HAVE_ASINH */ |

124 | n/a | |

125 | n/a | |

126 | n/a | #ifndef HAVE_ATANH |

127 | n/a | /* atanh(x) |

128 | n/a | * Method : |

129 | n/a | * 1.Reduced x to positive by atanh(-x) = -atanh(x) |

130 | n/a | * 2.For x>=0.5 |

131 | n/a | * 1 2x x |

132 | n/a | * atanh(x) = --- * log(1 + -------) = 0.5 * log1p(2 * -------) |

133 | n/a | * 2 1 - x 1 - x |

134 | n/a | * |

135 | n/a | * For x<0.5 |

136 | n/a | * atanh(x) = 0.5*log1p(2x+2x*x/(1-x)) |

137 | n/a | * |

138 | n/a | * Special cases: |

139 | n/a | * atanh(x) is NaN if |x| >= 1 with signal; |

140 | n/a | * atanh(NaN) is that NaN with no signal; |

141 | n/a | * |

142 | n/a | */ |

143 | n/a | |

144 | n/a | double |

145 | n/a | _Py_atanh(double x) |

146 | n/a | { |

147 | n/a | double absx; |

148 | n/a | double t; |

149 | n/a | |

150 | n/a | if (Py_IS_NAN(x)) { |

151 | n/a | return x+x; |

152 | n/a | } |

153 | n/a | absx = fabs(x); |

154 | n/a | if (absx >= 1.) { /* |x| >= 1 */ |

155 | n/a | errno = EDOM; |

156 | n/a | #ifdef Py_NAN |

157 | n/a | return Py_NAN; |

158 | n/a | #else |

159 | n/a | return x / zero; |

160 | n/a | #endif |

161 | n/a | } |

162 | n/a | if (absx < two_pow_m28) { /* |x| < 2**-28 */ |

163 | n/a | return x; |

164 | n/a | } |

165 | n/a | if (absx < 0.5) { /* |x| < 0.5 */ |

166 | n/a | t = absx+absx; |

167 | n/a | t = 0.5 * m_log1p(t + t*absx / (1.0 - absx)); |

168 | n/a | } |

169 | n/a | else { /* 0.5 <= |x| <= 1.0 */ |

170 | n/a | t = 0.5 * m_log1p((absx + absx) / (1.0 - absx)); |

171 | n/a | } |

172 | n/a | return copysign(t, x); |

173 | n/a | } |

174 | n/a | #endif /* HAVE_ATANH */ |

175 | n/a | |

176 | n/a | |

177 | n/a | #ifndef HAVE_EXPM1 |

178 | n/a | /* Mathematically, expm1(x) = exp(x) - 1. The expm1 function is designed |

179 | n/a | to avoid the significant loss of precision that arises from direct |

180 | n/a | evaluation of the expression exp(x) - 1, for x near 0. */ |

181 | n/a | |

182 | n/a | double |

183 | n/a | _Py_expm1(double x) |

184 | n/a | { |

185 | n/a | /* For abs(x) >= log(2), it's safe to evaluate exp(x) - 1 directly; this |

186 | n/a | also works fine for infinities and nans. |

187 | n/a | |

188 | n/a | For smaller x, we can use a method due to Kahan that achieves close to |

189 | n/a | full accuracy. |

190 | n/a | */ |

191 | n/a | |

192 | n/a | if (fabs(x) < 0.7) { |

193 | n/a | double u; |

194 | n/a | u = exp(x); |

195 | n/a | if (u == 1.0) |

196 | n/a | return x; |

197 | n/a | else |

198 | n/a | return (u - 1.0) * x / log(u); |

199 | n/a | } |

200 | n/a | else |

201 | n/a | return exp(x) - 1.0; |

202 | n/a | } |

203 | n/a | #endif /* HAVE_EXPM1 */ |

204 | n/a | |

205 | n/a | |

206 | n/a | /* log1p(x) = log(1+x). The log1p function is designed to avoid the |

207 | n/a | significant loss of precision that arises from direct evaluation when x is |

208 | n/a | small. */ |

209 | n/a | |

210 | n/a | double |

211 | n/a | _Py_log1p(double x) |

212 | n/a | { |

213 | n/a | #ifdef HAVE_LOG1P |

214 | n/a | /* Some platforms supply a log1p function but don't respect the sign of |

215 | n/a | zero: log1p(-0.0) gives 0.0 instead of the correct result of -0.0. |

216 | n/a | |

217 | n/a | To save fiddling with configure tests and platform checks, we handle the |

218 | n/a | special case of zero input directly on all platforms. |

219 | n/a | */ |

220 | n/a | if (x == 0.0) { |

221 | n/a | return x; |

222 | n/a | } |

223 | n/a | else { |

224 | n/a | return log1p(x); |

225 | n/a | } |

226 | n/a | #else |

227 | n/a | /* For x small, we use the following approach. Let y be the nearest float |

228 | n/a | to 1+x, then |

229 | n/a | |

230 | n/a | 1+x = y * (1 - (y-1-x)/y) |

231 | n/a | |

232 | n/a | so log(1+x) = log(y) + log(1-(y-1-x)/y). Since (y-1-x)/y is tiny, the |

233 | n/a | second term is well approximated by (y-1-x)/y. If abs(x) >= |

234 | n/a | DBL_EPSILON/2 or the rounding-mode is some form of round-to-nearest |

235 | n/a | then y-1-x will be exactly representable, and is computed exactly by |

236 | n/a | (y-1)-x. |

237 | n/a | |

238 | n/a | If abs(x) < DBL_EPSILON/2 and the rounding mode is not known to be |

239 | n/a | round-to-nearest then this method is slightly dangerous: 1+x could be |

240 | n/a | rounded up to 1+DBL_EPSILON instead of down to 1, and in that case |

241 | n/a | y-1-x will not be exactly representable any more and the result can be |

242 | n/a | off by many ulps. But this is easily fixed: for a floating-point |

243 | n/a | number |x| < DBL_EPSILON/2., the closest floating-point number to |

244 | n/a | log(1+x) is exactly x. |

245 | n/a | */ |

246 | n/a | |

247 | n/a | double y; |

248 | n/a | if (fabs(x) < DBL_EPSILON / 2.) { |

249 | n/a | return x; |

250 | n/a | } |

251 | n/a | else if (-0.5 <= x && x <= 1.) { |

252 | n/a | /* WARNING: it's possible that an overeager compiler |

253 | n/a | will incorrectly optimize the following two lines |

254 | n/a | to the equivalent of "return log(1.+x)". If this |

255 | n/a | happens, then results from log1p will be inaccurate |

256 | n/a | for small x. */ |

257 | n/a | y = 1.+x; |

258 | n/a | return log(y) - ((y - 1.) - x) / y; |

259 | n/a | } |

260 | n/a | else { |

261 | n/a | /* NaNs and infinities should end up here */ |

262 | n/a | return log(1.+x); |

263 | n/a | } |

264 | n/a | #endif /* ifdef HAVE_LOG1P */ |

265 | n/a | } |

266 | n/a |