# Python code coverage for Modules/mathmodule.c

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

1 | n/a | /* Math module -- standard C math library functions, pi and e */ |

2 | n/a | |

3 | n/a | /* Here are some comments from Tim Peters, extracted from the |

4 | n/a | discussion attached to http://bugs.python.org/issue1640. They |

5 | n/a | describe the general aims of the math module with respect to |

6 | n/a | special values, IEEE-754 floating-point exceptions, and Python |

7 | n/a | exceptions. |

8 | n/a | |

9 | n/a | These are the "spirit of 754" rules: |

10 | n/a | |

11 | n/a | 1. If the mathematical result is a real number, but of magnitude too |

12 | n/a | large to approximate by a machine float, overflow is signaled and the |

13 | n/a | result is an infinity (with the appropriate sign). |

14 | n/a | |

15 | n/a | 2. If the mathematical result is a real number, but of magnitude too |

16 | n/a | small to approximate by a machine float, underflow is signaled and the |

17 | n/a | result is a zero (with the appropriate sign). |

18 | n/a | |

19 | n/a | 3. At a singularity (a value x such that the limit of f(y) as y |

20 | n/a | approaches x exists and is an infinity), "divide by zero" is signaled |

21 | n/a | and the result is an infinity (with the appropriate sign). This is |

22 | n/a | complicated a little by that the left-side and right-side limits may |

23 | n/a | not be the same; e.g., 1/x approaches +inf or -inf as x approaches 0 |

24 | n/a | from the positive or negative directions. In that specific case, the |

25 | n/a | sign of the zero determines the result of 1/0. |

26 | n/a | |

27 | n/a | 4. At a point where a function has no defined result in the extended |

28 | n/a | reals (i.e., the reals plus an infinity or two), invalid operation is |

29 | n/a | signaled and a NaN is returned. |

30 | n/a | |

31 | n/a | And these are what Python has historically /tried/ to do (but not |

32 | n/a | always successfully, as platform libm behavior varies a lot): |

33 | n/a | |

34 | n/a | For #1, raise OverflowError. |

35 | n/a | |

36 | n/a | For #2, return a zero (with the appropriate sign if that happens by |

37 | n/a | accident ;-)). |

38 | n/a | |

39 | n/a | For #3 and #4, raise ValueError. It may have made sense to raise |

40 | n/a | Python's ZeroDivisionError in #3, but historically that's only been |

41 | n/a | raised for division by zero and mod by zero. |

42 | n/a | |

43 | n/a | */ |

44 | n/a | |

45 | n/a | /* |

46 | n/a | In general, on an IEEE-754 platform the aim is to follow the C99 |

47 | n/a | standard, including Annex 'F', whenever possible. Where the |

48 | n/a | standard recommends raising the 'divide-by-zero' or 'invalid' |

49 | n/a | floating-point exceptions, Python should raise a ValueError. Where |

50 | n/a | the standard recommends raising 'overflow', Python should raise an |

51 | n/a | OverflowError. In all other circumstances a value should be |

52 | n/a | returned. |

53 | n/a | */ |

54 | n/a | |

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

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

57 | n/a | |

58 | n/a | #include "clinic/mathmodule.c.h" |

59 | n/a | |

60 | n/a | /*[clinic input] |

61 | n/a | module math |

62 | n/a | [clinic start generated code]*/ |

63 | n/a | /*[clinic end generated code: output=da39a3ee5e6b4b0d input=76bc7002685dd942]*/ |

64 | n/a | |

65 | n/a | |

66 | n/a | /* |

67 | n/a | sin(pi*x), giving accurate results for all finite x (especially x |

68 | n/a | integral or close to an integer). This is here for use in the |

69 | n/a | reflection formula for the gamma function. It conforms to IEEE |

70 | n/a | 754-2008 for finite arguments, but not for infinities or nans. |

71 | n/a | */ |

72 | n/a | |

73 | n/a | static const double pi = 3.141592653589793238462643383279502884197; |

74 | n/a | static const double sqrtpi = 1.772453850905516027298167483341145182798; |

75 | n/a | static const double logpi = 1.144729885849400174143427351353058711647; |

76 | n/a | |

77 | n/a | static double |

78 | n/a | sinpi(double x) |

79 | n/a | { |

80 | n/a | double y, r; |

81 | n/a | int n; |

82 | n/a | /* this function should only ever be called for finite arguments */ |

83 | n/a | assert(Py_IS_FINITE(x)); |

84 | n/a | y = fmod(fabs(x), 2.0); |

85 | n/a | n = (int)round(2.0*y); |

86 | n/a | assert(0 <= n && n <= 4); |

87 | n/a | switch (n) { |

88 | n/a | case 0: |

89 | n/a | r = sin(pi*y); |

90 | n/a | break; |

91 | n/a | case 1: |

92 | n/a | r = cos(pi*(y-0.5)); |

93 | n/a | break; |

94 | n/a | case 2: |

95 | n/a | /* N.B. -sin(pi*(y-1.0)) is *not* equivalent: it would give |

96 | n/a | -0.0 instead of 0.0 when y == 1.0. */ |

97 | n/a | r = sin(pi*(1.0-y)); |

98 | n/a | break; |

99 | n/a | case 3: |

100 | n/a | r = -cos(pi*(y-1.5)); |

101 | n/a | break; |

102 | n/a | case 4: |

103 | n/a | r = sin(pi*(y-2.0)); |

104 | n/a | break; |

105 | n/a | default: |

106 | n/a | assert(0); /* should never get here */ |

107 | n/a | r = -1.23e200; /* silence gcc warning */ |

108 | n/a | } |

109 | n/a | return copysign(1.0, x)*r; |

110 | n/a | } |

111 | n/a | |

112 | n/a | /* Implementation of the real gamma function. In extensive but non-exhaustive |

113 | n/a | random tests, this function proved accurate to within <= 10 ulps across the |

114 | n/a | entire float domain. Note that accuracy may depend on the quality of the |

115 | n/a | system math functions, the pow function in particular. Special cases |

116 | n/a | follow C99 annex F. The parameters and method are tailored to platforms |

117 | n/a | whose double format is the IEEE 754 binary64 format. |

118 | n/a | |

119 | n/a | Method: for x > 0.0 we use the Lanczos approximation with parameters N=13 |

120 | n/a | and g=6.024680040776729583740234375; these parameters are amongst those |

121 | n/a | used by the Boost library. Following Boost (again), we re-express the |

122 | n/a | Lanczos sum as a rational function, and compute it that way. The |

123 | n/a | coefficients below were computed independently using MPFR, and have been |

124 | n/a | double-checked against the coefficients in the Boost source code. |

125 | n/a | |

126 | n/a | For x < 0.0 we use the reflection formula. |

127 | n/a | |

128 | n/a | There's one minor tweak that deserves explanation: Lanczos' formula for |

129 | n/a | Gamma(x) involves computing pow(x+g-0.5, x-0.5) / exp(x+g-0.5). For many x |

130 | n/a | values, x+g-0.5 can be represented exactly. However, in cases where it |

131 | n/a | can't be represented exactly the small error in x+g-0.5 can be magnified |

132 | n/a | significantly by the pow and exp calls, especially for large x. A cheap |

133 | n/a | correction is to multiply by (1 + e*g/(x+g-0.5)), where e is the error |

134 | n/a | involved in the computation of x+g-0.5 (that is, e = computed value of |

135 | n/a | x+g-0.5 - exact value of x+g-0.5). Here's the proof: |

136 | n/a | |

137 | n/a | Correction factor |

138 | n/a | ----------------- |

139 | n/a | Write x+g-0.5 = y-e, where y is exactly representable as an IEEE 754 |

140 | n/a | double, and e is tiny. Then: |

141 | n/a | |

142 | n/a | pow(x+g-0.5,x-0.5)/exp(x+g-0.5) = pow(y-e, x-0.5)/exp(y-e) |

143 | n/a | = pow(y, x-0.5)/exp(y) * C, |

144 | n/a | |

145 | n/a | where the correction_factor C is given by |

146 | n/a | |

147 | n/a | C = pow(1-e/y, x-0.5) * exp(e) |

148 | n/a | |

149 | n/a | Since e is tiny, pow(1-e/y, x-0.5) ~ 1-(x-0.5)*e/y, and exp(x) ~ 1+e, so: |

150 | n/a | |

151 | n/a | C ~ (1-(x-0.5)*e/y) * (1+e) ~ 1 + e*(y-(x-0.5))/y |

152 | n/a | |

153 | n/a | But y-(x-0.5) = g+e, and g+e ~ g. So we get C ~ 1 + e*g/y, and |

154 | n/a | |

155 | n/a | pow(x+g-0.5,x-0.5)/exp(x+g-0.5) ~ pow(y, x-0.5)/exp(y) * (1 + e*g/y), |

156 | n/a | |

157 | n/a | Note that for accuracy, when computing r*C it's better to do |

158 | n/a | |

159 | n/a | r + e*g/y*r; |

160 | n/a | |

161 | n/a | than |

162 | n/a | |

163 | n/a | r * (1 + e*g/y); |

164 | n/a | |

165 | n/a | since the addition in the latter throws away most of the bits of |

166 | n/a | information in e*g/y. |

167 | n/a | */ |

168 | n/a | |

169 | n/a | #define LANCZOS_N 13 |

170 | n/a | static const double lanczos_g = 6.024680040776729583740234375; |

171 | n/a | static const double lanczos_g_minus_half = 5.524680040776729583740234375; |

172 | n/a | static const double lanczos_num_coeffs[LANCZOS_N] = { |

173 | n/a | 23531376880.410759688572007674451636754734846804940, |

174 | n/a | 42919803642.649098768957899047001988850926355848959, |

175 | n/a | 35711959237.355668049440185451547166705960488635843, |

176 | n/a | 17921034426.037209699919755754458931112671403265390, |

177 | n/a | 6039542586.3520280050642916443072979210699388420708, |

178 | n/a | 1439720407.3117216736632230727949123939715485786772, |

179 | n/a | 248874557.86205415651146038641322942321632125127801, |

180 | n/a | 31426415.585400194380614231628318205362874684987640, |

181 | n/a | 2876370.6289353724412254090516208496135991145378768, |

182 | n/a | 186056.26539522349504029498971604569928220784236328, |

183 | n/a | 8071.6720023658162106380029022722506138218516325024, |

184 | n/a | 210.82427775157934587250973392071336271166969580291, |

185 | n/a | 2.5066282746310002701649081771338373386264310793408 |

186 | n/a | }; |

187 | n/a | |

188 | n/a | /* denominator is x*(x+1)*...*(x+LANCZOS_N-2) */ |

189 | n/a | static const double lanczos_den_coeffs[LANCZOS_N] = { |

190 | n/a | 0.0, 39916800.0, 120543840.0, 150917976.0, 105258076.0, 45995730.0, |

191 | n/a | 13339535.0, 2637558.0, 357423.0, 32670.0, 1925.0, 66.0, 1.0}; |

192 | n/a | |

193 | n/a | /* gamma values for small positive integers, 1 though NGAMMA_INTEGRAL */ |

194 | n/a | #define NGAMMA_INTEGRAL 23 |

195 | n/a | static const double gamma_integral[NGAMMA_INTEGRAL] = { |

196 | n/a | 1.0, 1.0, 2.0, 6.0, 24.0, 120.0, 720.0, 5040.0, 40320.0, 362880.0, |

197 | n/a | 3628800.0, 39916800.0, 479001600.0, 6227020800.0, 87178291200.0, |

198 | n/a | 1307674368000.0, 20922789888000.0, 355687428096000.0, |

199 | n/a | 6402373705728000.0, 121645100408832000.0, 2432902008176640000.0, |

200 | n/a | 51090942171709440000.0, 1124000727777607680000.0, |

201 | n/a | }; |

202 | n/a | |

203 | n/a | /* Lanczos' sum L_g(x), for positive x */ |

204 | n/a | |

205 | n/a | static double |

206 | n/a | lanczos_sum(double x) |

207 | n/a | { |

208 | n/a | double num = 0.0, den = 0.0; |

209 | n/a | int i; |

210 | n/a | assert(x > 0.0); |

211 | n/a | /* evaluate the rational function lanczos_sum(x). For large |

212 | n/a | x, the obvious algorithm risks overflow, so we instead |

213 | n/a | rescale the denominator and numerator of the rational |

214 | n/a | function by x**(1-LANCZOS_N) and treat this as a |

215 | n/a | rational function in 1/x. This also reduces the error for |

216 | n/a | larger x values. The choice of cutoff point (5.0 below) is |

217 | n/a | somewhat arbitrary; in tests, smaller cutoff values than |

218 | n/a | this resulted in lower accuracy. */ |

219 | n/a | if (x < 5.0) { |

220 | n/a | for (i = LANCZOS_N; --i >= 0; ) { |

221 | n/a | num = num * x + lanczos_num_coeffs[i]; |

222 | n/a | den = den * x + lanczos_den_coeffs[i]; |

223 | n/a | } |

224 | n/a | } |

225 | n/a | else { |

226 | n/a | for (i = 0; i < LANCZOS_N; i++) { |

227 | n/a | num = num / x + lanczos_num_coeffs[i]; |

228 | n/a | den = den / x + lanczos_den_coeffs[i]; |

229 | n/a | } |

230 | n/a | } |

231 | n/a | return num/den; |

232 | n/a | } |

233 | n/a | |

234 | n/a | /* Constant for +infinity, generated in the same way as float('inf'). */ |

235 | n/a | |

236 | n/a | static double |

237 | n/a | m_inf(void) |

238 | n/a | { |

239 | n/a | #ifndef PY_NO_SHORT_FLOAT_REPR |

240 | n/a | return _Py_dg_infinity(0); |

241 | n/a | #else |

242 | n/a | return Py_HUGE_VAL; |

243 | n/a | #endif |

244 | n/a | } |

245 | n/a | |

246 | n/a | /* Constant nan value, generated in the same way as float('nan'). */ |

247 | n/a | /* We don't currently assume that Py_NAN is defined everywhere. */ |

248 | n/a | |

249 | n/a | #if !defined(PY_NO_SHORT_FLOAT_REPR) || defined(Py_NAN) |

250 | n/a | |

251 | n/a | static double |

252 | n/a | m_nan(void) |

253 | n/a | { |

254 | n/a | #ifndef PY_NO_SHORT_FLOAT_REPR |

255 | n/a | return _Py_dg_stdnan(0); |

256 | n/a | #else |

257 | n/a | return Py_NAN; |

258 | n/a | #endif |

259 | n/a | } |

260 | n/a | |

261 | n/a | #endif |

262 | n/a | |

263 | n/a | static double |

264 | n/a | m_tgamma(double x) |

265 | n/a | { |

266 | n/a | double absx, r, y, z, sqrtpow; |

267 | n/a | |

268 | n/a | /* special cases */ |

269 | n/a | if (!Py_IS_FINITE(x)) { |

270 | n/a | if (Py_IS_NAN(x) || x > 0.0) |

271 | n/a | return x; /* tgamma(nan) = nan, tgamma(inf) = inf */ |

272 | n/a | else { |

273 | n/a | errno = EDOM; |

274 | n/a | return Py_NAN; /* tgamma(-inf) = nan, invalid */ |

275 | n/a | } |

276 | n/a | } |

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

278 | n/a | errno = EDOM; |

279 | n/a | /* tgamma(+-0.0) = +-inf, divide-by-zero */ |

280 | n/a | return copysign(Py_HUGE_VAL, x); |

281 | n/a | } |

282 | n/a | |

283 | n/a | /* integer arguments */ |

284 | n/a | if (x == floor(x)) { |

285 | n/a | if (x < 0.0) { |

286 | n/a | errno = EDOM; /* tgamma(n) = nan, invalid for */ |

287 | n/a | return Py_NAN; /* negative integers n */ |

288 | n/a | } |

289 | n/a | if (x <= NGAMMA_INTEGRAL) |

290 | n/a | return gamma_integral[(int)x - 1]; |

291 | n/a | } |

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

293 | n/a | |

294 | n/a | /* tiny arguments: tgamma(x) ~ 1/x for x near 0 */ |

295 | n/a | if (absx < 1e-20) { |

296 | n/a | r = 1.0/x; |

297 | n/a | if (Py_IS_INFINITY(r)) |

298 | n/a | errno = ERANGE; |

299 | n/a | return r; |

300 | n/a | } |

301 | n/a | |

302 | n/a | /* large arguments: assuming IEEE 754 doubles, tgamma(x) overflows for |

303 | n/a | x > 200, and underflows to +-0.0 for x < -200, not a negative |

304 | n/a | integer. */ |

305 | n/a | if (absx > 200.0) { |

306 | n/a | if (x < 0.0) { |

307 | n/a | return 0.0/sinpi(x); |

308 | n/a | } |

309 | n/a | else { |

310 | n/a | errno = ERANGE; |

311 | n/a | return Py_HUGE_VAL; |

312 | n/a | } |

313 | n/a | } |

314 | n/a | |

315 | n/a | y = absx + lanczos_g_minus_half; |

316 | n/a | /* compute error in sum */ |

317 | n/a | if (absx > lanczos_g_minus_half) { |

318 | n/a | /* note: the correction can be foiled by an optimizing |

319 | n/a | compiler that (incorrectly) thinks that an expression like |

320 | n/a | a + b - a - b can be optimized to 0.0. This shouldn't |

321 | n/a | happen in a standards-conforming compiler. */ |

322 | n/a | double q = y - absx; |

323 | n/a | z = q - lanczos_g_minus_half; |

324 | n/a | } |

325 | n/a | else { |

326 | n/a | double q = y - lanczos_g_minus_half; |

327 | n/a | z = q - absx; |

328 | n/a | } |

329 | n/a | z = z * lanczos_g / y; |

330 | n/a | if (x < 0.0) { |

331 | n/a | r = -pi / sinpi(absx) / absx * exp(y) / lanczos_sum(absx); |

332 | n/a | r -= z * r; |

333 | n/a | if (absx < 140.0) { |

334 | n/a | r /= pow(y, absx - 0.5); |

335 | n/a | } |

336 | n/a | else { |

337 | n/a | sqrtpow = pow(y, absx / 2.0 - 0.25); |

338 | n/a | r /= sqrtpow; |

339 | n/a | r /= sqrtpow; |

340 | n/a | } |

341 | n/a | } |

342 | n/a | else { |

343 | n/a | r = lanczos_sum(absx) / exp(y); |

344 | n/a | r += z * r; |

345 | n/a | if (absx < 140.0) { |

346 | n/a | r *= pow(y, absx - 0.5); |

347 | n/a | } |

348 | n/a | else { |

349 | n/a | sqrtpow = pow(y, absx / 2.0 - 0.25); |

350 | n/a | r *= sqrtpow; |

351 | n/a | r *= sqrtpow; |

352 | n/a | } |

353 | n/a | } |

354 | n/a | if (Py_IS_INFINITY(r)) |

355 | n/a | errno = ERANGE; |

356 | n/a | return r; |

357 | n/a | } |

358 | n/a | |

359 | n/a | /* |

360 | n/a | lgamma: natural log of the absolute value of the Gamma function. |

361 | n/a | For large arguments, Lanczos' formula works extremely well here. |

362 | n/a | */ |

363 | n/a | |

364 | n/a | static double |

365 | n/a | m_lgamma(double x) |

366 | n/a | { |

367 | n/a | double r, absx; |

368 | n/a | |

369 | n/a | /* special cases */ |

370 | n/a | if (!Py_IS_FINITE(x)) { |

371 | n/a | if (Py_IS_NAN(x)) |

372 | n/a | return x; /* lgamma(nan) = nan */ |

373 | n/a | else |

374 | n/a | return Py_HUGE_VAL; /* lgamma(+-inf) = +inf */ |

375 | n/a | } |

376 | n/a | |

377 | n/a | /* integer arguments */ |

378 | n/a | if (x == floor(x) && x <= 2.0) { |

379 | n/a | if (x <= 0.0) { |

380 | n/a | errno = EDOM; /* lgamma(n) = inf, divide-by-zero for */ |

381 | n/a | return Py_HUGE_VAL; /* integers n <= 0 */ |

382 | n/a | } |

383 | n/a | else { |

384 | n/a | return 0.0; /* lgamma(1) = lgamma(2) = 0.0 */ |

385 | n/a | } |

386 | n/a | } |

387 | n/a | |

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

389 | n/a | /* tiny arguments: lgamma(x) ~ -log(fabs(x)) for small x */ |

390 | n/a | if (absx < 1e-20) |

391 | n/a | return -log(absx); |

392 | n/a | |

393 | n/a | /* Lanczos' formula. We could save a fraction of a ulp in accuracy by |

394 | n/a | having a second set of numerator coefficients for lanczos_sum that |

395 | n/a | absorbed the exp(-lanczos_g) term, and throwing out the lanczos_g |

396 | n/a | subtraction below; it's probably not worth it. */ |

397 | n/a | r = log(lanczos_sum(absx)) - lanczos_g; |

398 | n/a | r += (absx - 0.5) * (log(absx + lanczos_g - 0.5) - 1); |

399 | n/a | if (x < 0.0) |

400 | n/a | /* Use reflection formula to get value for negative x. */ |

401 | n/a | r = logpi - log(fabs(sinpi(absx))) - log(absx) - r; |

402 | n/a | if (Py_IS_INFINITY(r)) |

403 | n/a | errno = ERANGE; |

404 | n/a | return r; |

405 | n/a | } |

406 | n/a | |

407 | n/a | /* |

408 | n/a | Implementations of the error function erf(x) and the complementary error |

409 | n/a | function erfc(x). |

410 | n/a | |

411 | n/a | Method: we use a series approximation for erf for small x, and a continued |

412 | n/a | fraction approximation for erfc(x) for larger x; |

413 | n/a | combined with the relations erf(-x) = -erf(x) and erfc(x) = 1.0 - erf(x), |

414 | n/a | this gives us erf(x) and erfc(x) for all x. |

415 | n/a | |

416 | n/a | The series expansion used is: |

417 | n/a | |

418 | n/a | erf(x) = x*exp(-x*x)/sqrt(pi) * [ |

419 | n/a | 2/1 + 4/3 x**2 + 8/15 x**4 + 16/105 x**6 + ...] |

420 | n/a | |

421 | n/a | The coefficient of x**(2k-2) here is 4**k*factorial(k)/factorial(2*k). |

422 | n/a | This series converges well for smallish x, but slowly for larger x. |

423 | n/a | |

424 | n/a | The continued fraction expansion used is: |

425 | n/a | |

426 | n/a | erfc(x) = x*exp(-x*x)/sqrt(pi) * [1/(0.5 + x**2 -) 0.5/(2.5 + x**2 - ) |

427 | n/a | 3.0/(4.5 + x**2 - ) 7.5/(6.5 + x**2 - ) ...] |

428 | n/a | |

429 | n/a | after the first term, the general term has the form: |

430 | n/a | |

431 | n/a | k*(k-0.5)/(2*k+0.5 + x**2 - ...). |

432 | n/a | |

433 | n/a | This expansion converges fast for larger x, but convergence becomes |

434 | n/a | infinitely slow as x approaches 0.0. The (somewhat naive) continued |

435 | n/a | fraction evaluation algorithm used below also risks overflow for large x; |

436 | n/a | but for large x, erfc(x) == 0.0 to within machine precision. (For |

437 | n/a | example, erfc(30.0) is approximately 2.56e-393). |

438 | n/a | |

439 | n/a | Parameters: use series expansion for abs(x) < ERF_SERIES_CUTOFF and |

440 | n/a | continued fraction expansion for ERF_SERIES_CUTOFF <= abs(x) < |

441 | n/a | ERFC_CONTFRAC_CUTOFF. ERFC_SERIES_TERMS and ERFC_CONTFRAC_TERMS are the |

442 | n/a | numbers of terms to use for the relevant expansions. */ |

443 | n/a | |

444 | n/a | #define ERF_SERIES_CUTOFF 1.5 |

445 | n/a | #define ERF_SERIES_TERMS 25 |

446 | n/a | #define ERFC_CONTFRAC_CUTOFF 30.0 |

447 | n/a | #define ERFC_CONTFRAC_TERMS 50 |

448 | n/a | |

449 | n/a | /* |

450 | n/a | Error function, via power series. |

451 | n/a | |

452 | n/a | Given a finite float x, return an approximation to erf(x). |

453 | n/a | Converges reasonably fast for small x. |

454 | n/a | */ |

455 | n/a | |

456 | n/a | static double |

457 | n/a | m_erf_series(double x) |

458 | n/a | { |

459 | n/a | double x2, acc, fk, result; |

460 | n/a | int i, saved_errno; |

461 | n/a | |

462 | n/a | x2 = x * x; |

463 | n/a | acc = 0.0; |

464 | n/a | fk = (double)ERF_SERIES_TERMS + 0.5; |

465 | n/a | for (i = 0; i < ERF_SERIES_TERMS; i++) { |

466 | n/a | acc = 2.0 + x2 * acc / fk; |

467 | n/a | fk -= 1.0; |

468 | n/a | } |

469 | n/a | /* Make sure the exp call doesn't affect errno; |

470 | n/a | see m_erfc_contfrac for more. */ |

471 | n/a | saved_errno = errno; |

472 | n/a | result = acc * x * exp(-x2) / sqrtpi; |

473 | n/a | errno = saved_errno; |

474 | n/a | return result; |

475 | n/a | } |

476 | n/a | |

477 | n/a | /* |

478 | n/a | Complementary error function, via continued fraction expansion. |

479 | n/a | |

480 | n/a | Given a positive float x, return an approximation to erfc(x). Converges |

481 | n/a | reasonably fast for x large (say, x > 2.0), and should be safe from |

482 | n/a | overflow if x and nterms are not too large. On an IEEE 754 machine, with x |

483 | n/a | <= 30.0, we're safe up to nterms = 100. For x >= 30.0, erfc(x) is smaller |

484 | n/a | than the smallest representable nonzero float. */ |

485 | n/a | |

486 | n/a | static double |

487 | n/a | m_erfc_contfrac(double x) |

488 | n/a | { |

489 | n/a | double x2, a, da, p, p_last, q, q_last, b, result; |

490 | n/a | int i, saved_errno; |

491 | n/a | |

492 | n/a | if (x >= ERFC_CONTFRAC_CUTOFF) |

493 | n/a | return 0.0; |

494 | n/a | |

495 | n/a | x2 = x*x; |

496 | n/a | a = 0.0; |

497 | n/a | da = 0.5; |

498 | n/a | p = 1.0; p_last = 0.0; |

499 | n/a | q = da + x2; q_last = 1.0; |

500 | n/a | for (i = 0; i < ERFC_CONTFRAC_TERMS; i++) { |

501 | n/a | double temp; |

502 | n/a | a += da; |

503 | n/a | da += 2.0; |

504 | n/a | b = da + x2; |

505 | n/a | temp = p; p = b*p - a*p_last; p_last = temp; |

506 | n/a | temp = q; q = b*q - a*q_last; q_last = temp; |

507 | n/a | } |

508 | n/a | /* Issue #8986: On some platforms, exp sets errno on underflow to zero; |

509 | n/a | save the current errno value so that we can restore it later. */ |

510 | n/a | saved_errno = errno; |

511 | n/a | result = p / q * x * exp(-x2) / sqrtpi; |

512 | n/a | errno = saved_errno; |

513 | n/a | return result; |

514 | n/a | } |

515 | n/a | |

516 | n/a | /* Error function erf(x), for general x */ |

517 | n/a | |

518 | n/a | static double |

519 | n/a | m_erf(double x) |

520 | n/a | { |

521 | n/a | double absx, cf; |

522 | n/a | |

523 | n/a | if (Py_IS_NAN(x)) |

524 | n/a | return x; |

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

526 | n/a | if (absx < ERF_SERIES_CUTOFF) |

527 | n/a | return m_erf_series(x); |

528 | n/a | else { |

529 | n/a | cf = m_erfc_contfrac(absx); |

530 | n/a | return x > 0.0 ? 1.0 - cf : cf - 1.0; |

531 | n/a | } |

532 | n/a | } |

533 | n/a | |

534 | n/a | /* Complementary error function erfc(x), for general x. */ |

535 | n/a | |

536 | n/a | static double |

537 | n/a | m_erfc(double x) |

538 | n/a | { |

539 | n/a | double absx, cf; |

540 | n/a | |

541 | n/a | if (Py_IS_NAN(x)) |

542 | n/a | return x; |

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

544 | n/a | if (absx < ERF_SERIES_CUTOFF) |

545 | n/a | return 1.0 - m_erf_series(x); |

546 | n/a | else { |

547 | n/a | cf = m_erfc_contfrac(absx); |

548 | n/a | return x > 0.0 ? cf : 2.0 - cf; |

549 | n/a | } |

550 | n/a | } |

551 | n/a | |

552 | n/a | /* |

553 | n/a | wrapper for atan2 that deals directly with special cases before |

554 | n/a | delegating to the platform libm for the remaining cases. This |

555 | n/a | is necessary to get consistent behaviour across platforms. |

556 | n/a | Windows, FreeBSD and alpha Tru64 are amongst platforms that don't |

557 | n/a | always follow C99. |

558 | n/a | */ |

559 | n/a | |

560 | n/a | static double |

561 | n/a | m_atan2(double y, double x) |

562 | n/a | { |

563 | n/a | if (Py_IS_NAN(x) || Py_IS_NAN(y)) |

564 | n/a | return Py_NAN; |

565 | n/a | if (Py_IS_INFINITY(y)) { |

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

567 | n/a | if (copysign(1., x) == 1.) |

568 | n/a | /* atan2(+-inf, +inf) == +-pi/4 */ |

569 | n/a | return copysign(0.25*Py_MATH_PI, y); |

570 | n/a | else |

571 | n/a | /* atan2(+-inf, -inf) == +-pi*3/4 */ |

572 | n/a | return copysign(0.75*Py_MATH_PI, y); |

573 | n/a | } |

574 | n/a | /* atan2(+-inf, x) == +-pi/2 for finite x */ |

575 | n/a | return copysign(0.5*Py_MATH_PI, y); |

576 | n/a | } |

577 | n/a | if (Py_IS_INFINITY(x) || y == 0.) { |

578 | n/a | if (copysign(1., x) == 1.) |

579 | n/a | /* atan2(+-y, +inf) = atan2(+-0, +x) = +-0. */ |

580 | n/a | return copysign(0., y); |

581 | n/a | else |

582 | n/a | /* atan2(+-y, -inf) = atan2(+-0., -x) = +-pi. */ |

583 | n/a | return copysign(Py_MATH_PI, y); |

584 | n/a | } |

585 | n/a | return atan2(y, x); |

586 | n/a | } |

587 | n/a | |

588 | n/a | /* |

589 | n/a | Various platforms (Solaris, OpenBSD) do nonstandard things for log(0), |

590 | n/a | log(-ve), log(NaN). Here are wrappers for log and log10 that deal with |

591 | n/a | special values directly, passing positive non-special values through to |

592 | n/a | the system log/log10. |

593 | n/a | */ |

594 | n/a | |

595 | n/a | static double |

596 | n/a | m_log(double x) |

597 | n/a | { |

598 | n/a | if (Py_IS_FINITE(x)) { |

599 | n/a | if (x > 0.0) |

600 | n/a | return log(x); |

601 | n/a | errno = EDOM; |

602 | n/a | if (x == 0.0) |

603 | n/a | return -Py_HUGE_VAL; /* log(0) = -inf */ |

604 | n/a | else |

605 | n/a | return Py_NAN; /* log(-ve) = nan */ |

606 | n/a | } |

607 | n/a | else if (Py_IS_NAN(x)) |

608 | n/a | return x; /* log(nan) = nan */ |

609 | n/a | else if (x > 0.0) |

610 | n/a | return x; /* log(inf) = inf */ |

611 | n/a | else { |

612 | n/a | errno = EDOM; |

613 | n/a | return Py_NAN; /* log(-inf) = nan */ |

614 | n/a | } |

615 | n/a | } |

616 | n/a | |

617 | n/a | /* |

618 | n/a | log2: log to base 2. |

619 | n/a | |

620 | n/a | Uses an algorithm that should: |

621 | n/a | |

622 | n/a | (a) produce exact results for powers of 2, and |

623 | n/a | (b) give a monotonic log2 (for positive finite floats), |

624 | n/a | assuming that the system log is monotonic. |

625 | n/a | */ |

626 | n/a | |

627 | n/a | static double |

628 | n/a | m_log2(double x) |

629 | n/a | { |

630 | n/a | if (!Py_IS_FINITE(x)) { |

631 | n/a | if (Py_IS_NAN(x)) |

632 | n/a | return x; /* log2(nan) = nan */ |

633 | n/a | else if (x > 0.0) |

634 | n/a | return x; /* log2(+inf) = +inf */ |

635 | n/a | else { |

636 | n/a | errno = EDOM; |

637 | n/a | return Py_NAN; /* log2(-inf) = nan, invalid-operation */ |

638 | n/a | } |

639 | n/a | } |

640 | n/a | |

641 | n/a | if (x > 0.0) { |

642 | n/a | #ifdef HAVE_LOG2 |

643 | n/a | return log2(x); |

644 | n/a | #else |

645 | n/a | double m; |

646 | n/a | int e; |

647 | n/a | m = frexp(x, &e); |

648 | n/a | /* We want log2(m * 2**e) == log(m) / log(2) + e. Care is needed when |

649 | n/a | * x is just greater than 1.0: in that case e is 1, log(m) is negative, |

650 | n/a | * and we get significant cancellation error from the addition of |

651 | n/a | * log(m) / log(2) to e. The slight rewrite of the expression below |

652 | n/a | * avoids this problem. |

653 | n/a | */ |

654 | n/a | if (x >= 1.0) { |

655 | n/a | return log(2.0 * m) / log(2.0) + (e - 1); |

656 | n/a | } |

657 | n/a | else { |

658 | n/a | return log(m) / log(2.0) + e; |

659 | n/a | } |

660 | n/a | #endif |

661 | n/a | } |

662 | n/a | else if (x == 0.0) { |

663 | n/a | errno = EDOM; |

664 | n/a | return -Py_HUGE_VAL; /* log2(0) = -inf, divide-by-zero */ |

665 | n/a | } |

666 | n/a | else { |

667 | n/a | errno = EDOM; |

668 | n/a | return Py_NAN; /* log2(-inf) = nan, invalid-operation */ |

669 | n/a | } |

670 | n/a | } |

671 | n/a | |

672 | n/a | static double |

673 | n/a | m_log10(double x) |

674 | n/a | { |

675 | n/a | if (Py_IS_FINITE(x)) { |

676 | n/a | if (x > 0.0) |

677 | n/a | return log10(x); |

678 | n/a | errno = EDOM; |

679 | n/a | if (x == 0.0) |

680 | n/a | return -Py_HUGE_VAL; /* log10(0) = -inf */ |

681 | n/a | else |

682 | n/a | return Py_NAN; /* log10(-ve) = nan */ |

683 | n/a | } |

684 | n/a | else if (Py_IS_NAN(x)) |

685 | n/a | return x; /* log10(nan) = nan */ |

686 | n/a | else if (x > 0.0) |

687 | n/a | return x; /* log10(inf) = inf */ |

688 | n/a | else { |

689 | n/a | errno = EDOM; |

690 | n/a | return Py_NAN; /* log10(-inf) = nan */ |

691 | n/a | } |

692 | n/a | } |

693 | n/a | |

694 | n/a | |

695 | n/a | /*[clinic input] |

696 | n/a | math.gcd |

697 | n/a | |

698 | n/a | x as a: object |

699 | n/a | y as b: object |

700 | n/a | / |

701 | n/a | |

702 | n/a | greatest common divisor of x and y |

703 | n/a | [clinic start generated code]*/ |

704 | n/a | |

705 | n/a | static PyObject * |

706 | n/a | math_gcd_impl(PyObject *module, PyObject *a, PyObject *b) |

707 | n/a | /*[clinic end generated code: output=7b2e0c151bd7a5d8 input=c2691e57fb2a98fa]*/ |

708 | n/a | { |

709 | n/a | PyObject *g; |

710 | n/a | |

711 | n/a | a = PyNumber_Index(a); |

712 | n/a | if (a == NULL) |

713 | n/a | return NULL; |

714 | n/a | b = PyNumber_Index(b); |

715 | n/a | if (b == NULL) { |

716 | n/a | Py_DECREF(a); |

717 | n/a | return NULL; |

718 | n/a | } |

719 | n/a | g = _PyLong_GCD(a, b); |

720 | n/a | Py_DECREF(a); |

721 | n/a | Py_DECREF(b); |

722 | n/a | return g; |

723 | n/a | } |

724 | n/a | |

725 | n/a | |

726 | n/a | /* Call is_error when errno != 0, and where x is the result libm |

727 | n/a | * returned. is_error will usually set up an exception and return |

728 | n/a | * true (1), but may return false (0) without setting up an exception. |

729 | n/a | */ |

730 | n/a | static int |

731 | n/a | is_error(double x) |

732 | n/a | { |

733 | n/a | int result = 1; /* presumption of guilt */ |

734 | n/a | assert(errno); /* non-zero errno is a precondition for calling */ |

735 | n/a | if (errno == EDOM) |

736 | n/a | PyErr_SetString(PyExc_ValueError, "math domain error"); |

737 | n/a | |

738 | n/a | else if (errno == ERANGE) { |

739 | n/a | /* ANSI C generally requires libm functions to set ERANGE |

740 | n/a | * on overflow, but also generally *allows* them to set |

741 | n/a | * ERANGE on underflow too. There's no consistency about |

742 | n/a | * the latter across platforms. |

743 | n/a | * Alas, C99 never requires that errno be set. |

744 | n/a | * Here we suppress the underflow errors (libm functions |

745 | n/a | * should return a zero on underflow, and +- HUGE_VAL on |

746 | n/a | * overflow, so testing the result for zero suffices to |

747 | n/a | * distinguish the cases). |

748 | n/a | * |

749 | n/a | * On some platforms (Ubuntu/ia64) it seems that errno can be |

750 | n/a | * set to ERANGE for subnormal results that do *not* underflow |

751 | n/a | * to zero. So to be safe, we'll ignore ERANGE whenever the |

752 | n/a | * function result is less than one in absolute value. |

753 | n/a | */ |

754 | n/a | if (fabs(x) < 1.0) |

755 | n/a | result = 0; |

756 | n/a | else |

757 | n/a | PyErr_SetString(PyExc_OverflowError, |

758 | n/a | "math range error"); |

759 | n/a | } |

760 | n/a | else |

761 | n/a | /* Unexpected math error */ |

762 | n/a | PyErr_SetFromErrno(PyExc_ValueError); |

763 | n/a | return result; |

764 | n/a | } |

765 | n/a | |

766 | n/a | /* |

767 | n/a | math_1 is used to wrap a libm function f that takes a double |

768 | n/a | argument and returns a double. |

769 | n/a | |

770 | n/a | The error reporting follows these rules, which are designed to do |

771 | n/a | the right thing on C89/C99 platforms and IEEE 754/non IEEE 754 |

772 | n/a | platforms. |

773 | n/a | |

774 | n/a | - a NaN result from non-NaN inputs causes ValueError to be raised |

775 | n/a | - an infinite result from finite inputs causes OverflowError to be |

776 | n/a | raised if can_overflow is 1, or raises ValueError if can_overflow |

777 | n/a | is 0. |

778 | n/a | - if the result is finite and errno == EDOM then ValueError is |

779 | n/a | raised |

780 | n/a | - if the result is finite and nonzero and errno == ERANGE then |

781 | n/a | OverflowError is raised |

782 | n/a | |

783 | n/a | The last rule is used to catch overflow on platforms which follow |

784 | n/a | C89 but for which HUGE_VAL is not an infinity. |

785 | n/a | |

786 | n/a | For the majority of one-argument functions these rules are enough |

787 | n/a | to ensure that Python's functions behave as specified in 'Annex F' |

788 | n/a | of the C99 standard, with the 'invalid' and 'divide-by-zero' |

789 | n/a | floating-point exceptions mapping to Python's ValueError and the |

790 | n/a | 'overflow' floating-point exception mapping to OverflowError. |

791 | n/a | math_1 only works for functions that don't have singularities *and* |

792 | n/a | the possibility of overflow; fortunately, that covers everything we |

793 | n/a | care about right now. |

794 | n/a | */ |

795 | n/a | |

796 | n/a | static PyObject * |

797 | n/a | math_1_to_whatever(PyObject *arg, double (*func) (double), |

798 | n/a | PyObject *(*from_double_func) (double), |

799 | n/a | int can_overflow) |

800 | n/a | { |

801 | n/a | double x, r; |

802 | n/a | x = PyFloat_AsDouble(arg); |

803 | n/a | if (x == -1.0 && PyErr_Occurred()) |

804 | n/a | return NULL; |

805 | n/a | errno = 0; |

806 | n/a | PyFPE_START_PROTECT("in math_1", return 0); |

807 | n/a | r = (*func)(x); |

808 | n/a | PyFPE_END_PROTECT(r); |

809 | n/a | if (Py_IS_NAN(r) && !Py_IS_NAN(x)) { |

810 | n/a | PyErr_SetString(PyExc_ValueError, |

811 | n/a | "math domain error"); /* invalid arg */ |

812 | n/a | return NULL; |

813 | n/a | } |

814 | n/a | if (Py_IS_INFINITY(r) && Py_IS_FINITE(x)) { |

815 | n/a | if (can_overflow) |

816 | n/a | PyErr_SetString(PyExc_OverflowError, |

817 | n/a | "math range error"); /* overflow */ |

818 | n/a | else |

819 | n/a | PyErr_SetString(PyExc_ValueError, |

820 | n/a | "math domain error"); /* singularity */ |

821 | n/a | return NULL; |

822 | n/a | } |

823 | n/a | if (Py_IS_FINITE(r) && errno && is_error(r)) |

824 | n/a | /* this branch unnecessary on most platforms */ |

825 | n/a | return NULL; |

826 | n/a | |

827 | n/a | return (*from_double_func)(r); |

828 | n/a | } |

829 | n/a | |

830 | n/a | /* variant of math_1, to be used when the function being wrapped is known to |

831 | n/a | set errno properly (that is, errno = EDOM for invalid or divide-by-zero, |

832 | n/a | errno = ERANGE for overflow). */ |

833 | n/a | |

834 | n/a | static PyObject * |

835 | n/a | math_1a(PyObject *arg, double (*func) (double)) |

836 | n/a | { |

837 | n/a | double x, r; |

838 | n/a | x = PyFloat_AsDouble(arg); |

839 | n/a | if (x == -1.0 && PyErr_Occurred()) |

840 | n/a | return NULL; |

841 | n/a | errno = 0; |

842 | n/a | PyFPE_START_PROTECT("in math_1a", return 0); |

843 | n/a | r = (*func)(x); |

844 | n/a | PyFPE_END_PROTECT(r); |

845 | n/a | if (errno && is_error(r)) |

846 | n/a | return NULL; |

847 | n/a | return PyFloat_FromDouble(r); |

848 | n/a | } |

849 | n/a | |

850 | n/a | /* |

851 | n/a | math_2 is used to wrap a libm function f that takes two double |

852 | n/a | arguments and returns a double. |

853 | n/a | |

854 | n/a | The error reporting follows these rules, which are designed to do |

855 | n/a | the right thing on C89/C99 platforms and IEEE 754/non IEEE 754 |

856 | n/a | platforms. |

857 | n/a | |

858 | n/a | - a NaN result from non-NaN inputs causes ValueError to be raised |

859 | n/a | - an infinite result from finite inputs causes OverflowError to be |

860 | n/a | raised. |

861 | n/a | - if the result is finite and errno == EDOM then ValueError is |

862 | n/a | raised |

863 | n/a | - if the result is finite and nonzero and errno == ERANGE then |

864 | n/a | OverflowError is raised |

865 | n/a | |

866 | n/a | The last rule is used to catch overflow on platforms which follow |

867 | n/a | C89 but for which HUGE_VAL is not an infinity. |

868 | n/a | |

869 | n/a | For most two-argument functions (copysign, fmod, hypot, atan2) |

870 | n/a | these rules are enough to ensure that Python's functions behave as |

871 | n/a | specified in 'Annex F' of the C99 standard, with the 'invalid' and |

872 | n/a | 'divide-by-zero' floating-point exceptions mapping to Python's |

873 | n/a | ValueError and the 'overflow' floating-point exception mapping to |

874 | n/a | OverflowError. |

875 | n/a | */ |

876 | n/a | |

877 | n/a | static PyObject * |

878 | n/a | math_1(PyObject *arg, double (*func) (double), int can_overflow) |

879 | n/a | { |

880 | n/a | return math_1_to_whatever(arg, func, PyFloat_FromDouble, can_overflow); |

881 | n/a | } |

882 | n/a | |

883 | n/a | static PyObject * |

884 | n/a | math_1_to_int(PyObject *arg, double (*func) (double), int can_overflow) |

885 | n/a | { |

886 | n/a | return math_1_to_whatever(arg, func, PyLong_FromDouble, can_overflow); |

887 | n/a | } |

888 | n/a | |

889 | n/a | static PyObject * |

890 | n/a | math_2(PyObject *args, double (*func) (double, double), const char *funcname) |

891 | n/a | { |

892 | n/a | PyObject *ox, *oy; |

893 | n/a | double x, y, r; |

894 | n/a | if (! PyArg_UnpackTuple(args, funcname, 2, 2, &ox, &oy)) |

895 | n/a | return NULL; |

896 | n/a | x = PyFloat_AsDouble(ox); |

897 | n/a | y = PyFloat_AsDouble(oy); |

898 | n/a | if ((x == -1.0 || y == -1.0) && PyErr_Occurred()) |

899 | n/a | return NULL; |

900 | n/a | errno = 0; |

901 | n/a | PyFPE_START_PROTECT("in math_2", return 0); |

902 | n/a | r = (*func)(x, y); |

903 | n/a | PyFPE_END_PROTECT(r); |

904 | n/a | if (Py_IS_NAN(r)) { |

905 | n/a | if (!Py_IS_NAN(x) && !Py_IS_NAN(y)) |

906 | n/a | errno = EDOM; |

907 | n/a | else |

908 | n/a | errno = 0; |

909 | n/a | } |

910 | n/a | else if (Py_IS_INFINITY(r)) { |

911 | n/a | if (Py_IS_FINITE(x) && Py_IS_FINITE(y)) |

912 | n/a | errno = ERANGE; |

913 | n/a | else |

914 | n/a | errno = 0; |

915 | n/a | } |

916 | n/a | if (errno && is_error(r)) |

917 | n/a | return NULL; |

918 | n/a | else |

919 | n/a | return PyFloat_FromDouble(r); |

920 | n/a | } |

921 | n/a | |

922 | n/a | #define FUNC1(funcname, func, can_overflow, docstring) \ |

923 | n/a | static PyObject * math_##funcname(PyObject *self, PyObject *args) { \ |

924 | n/a | return math_1(args, func, can_overflow); \ |

925 | n/a | }\ |

926 | n/a | PyDoc_STRVAR(math_##funcname##_doc, docstring); |

927 | n/a | |

928 | n/a | #define FUNC1A(funcname, func, docstring) \ |

929 | n/a | static PyObject * math_##funcname(PyObject *self, PyObject *args) { \ |

930 | n/a | return math_1a(args, func); \ |

931 | n/a | }\ |

932 | n/a | PyDoc_STRVAR(math_##funcname##_doc, docstring); |

933 | n/a | |

934 | n/a | #define FUNC2(funcname, func, docstring) \ |

935 | n/a | static PyObject * math_##funcname(PyObject *self, PyObject *args) { \ |

936 | n/a | return math_2(args, func, #funcname); \ |

937 | n/a | }\ |

938 | n/a | PyDoc_STRVAR(math_##funcname##_doc, docstring); |

939 | n/a | |

940 | n/a | FUNC1(acos, acos, 0, |

941 | n/a | "acos($module, x, /)\n--\n\n" |

942 | n/a | "Return the arc cosine (measured in radians) of x.") |

943 | n/a | FUNC1(acosh, m_acosh, 0, |

944 | n/a | "acosh($module, x, /)\n--\n\n" |

945 | n/a | "Return the inverse hyperbolic cosine of x.") |

946 | n/a | FUNC1(asin, asin, 0, |

947 | n/a | "asin($module, x, /)\n--\n\n" |

948 | n/a | "Return the arc sine (measured in radians) of x.") |

949 | n/a | FUNC1(asinh, m_asinh, 0, |

950 | n/a | "asinh($module, x, /)\n--\n\n" |

951 | n/a | "Return the inverse hyperbolic sine of x.") |

952 | n/a | FUNC1(atan, atan, 0, |

953 | n/a | "atan($module, x, /)\n--\n\n" |

954 | n/a | "Return the arc tangent (measured in radians) of x.") |

955 | n/a | FUNC2(atan2, m_atan2, |

956 | n/a | "atan2($module, y, x, /)\n--\n\n" |

957 | n/a | "Return the arc tangent (measured in radians) of y/x.\n\n" |

958 | n/a | "Unlike atan(y/x), the signs of both x and y are considered.") |

959 | n/a | FUNC1(atanh, m_atanh, 0, |

960 | n/a | "atanh($module, x, /)\n--\n\n" |

961 | n/a | "Return the inverse hyperbolic tangent of x.") |

962 | n/a | |

963 | n/a | /*[clinic input] |

964 | n/a | math.ceil |

965 | n/a | |

966 | n/a | x as number: object |

967 | n/a | / |

968 | n/a | |

969 | n/a | Return the ceiling of x as an Integral. |

970 | n/a | |

971 | n/a | This is the smallest integer >= x. |

972 | n/a | [clinic start generated code]*/ |

973 | n/a | |

974 | n/a | static PyObject * |

975 | n/a | math_ceil(PyObject *module, PyObject *number) |

976 | n/a | /*[clinic end generated code: output=6c3b8a78bc201c67 input=2725352806399cab]*/ |

977 | n/a | { |

978 | n/a | _Py_IDENTIFIER(__ceil__); |

979 | n/a | PyObject *method, *result; |

980 | n/a | |

981 | n/a | method = _PyObject_LookupSpecial(number, &PyId___ceil__); |

982 | n/a | if (method == NULL) { |

983 | n/a | if (PyErr_Occurred()) |

984 | n/a | return NULL; |

985 | n/a | return math_1_to_int(number, ceil, 0); |

986 | n/a | } |

987 | n/a | result = _PyObject_CallNoArg(method); |

988 | n/a | Py_DECREF(method); |

989 | n/a | return result; |

990 | n/a | } |

991 | n/a | |

992 | n/a | FUNC2(copysign, copysign, |

993 | n/a | "copysign($module, x, y, /)\n--\n\n" |

994 | n/a | "Return a float with the magnitude (absolute value) of x but the sign of y.\n\n" |

995 | n/a | "On platforms that support signed zeros, copysign(1.0, -0.0)\n" |

996 | n/a | "returns -1.0.\n") |

997 | n/a | FUNC1(cos, cos, 0, |

998 | n/a | "cos($module, x, /)\n--\n\n" |

999 | n/a | "Return the cosine of x (measured in radians).") |

1000 | n/a | FUNC1(cosh, cosh, 1, |

1001 | n/a | "cosh($module, x, /)\n--\n\n" |

1002 | n/a | "Return the hyperbolic cosine of x.") |

1003 | n/a | FUNC1A(erf, m_erf, |

1004 | n/a | "erf($module, x, /)\n--\n\n" |

1005 | n/a | "Error function at x.") |

1006 | n/a | FUNC1A(erfc, m_erfc, |

1007 | n/a | "erfc($module, x, /)\n--\n\n" |

1008 | n/a | "Complementary error function at x.") |

1009 | n/a | FUNC1(exp, exp, 1, |

1010 | n/a | "exp($module, x, /)\n--\n\n" |

1011 | n/a | "Return e raised to the power of x.") |

1012 | n/a | FUNC1(expm1, m_expm1, 1, |

1013 | n/a | "expm1($module, x, /)\n--\n\n" |

1014 | n/a | "Return exp(x)-1.\n\n" |

1015 | n/a | "This function avoids the loss of precision involved in the direct " |

1016 | n/a | "evaluation of exp(x)-1 for small x.") |

1017 | n/a | FUNC1(fabs, fabs, 0, |

1018 | n/a | "fabs($module, x, /)\n--\n\n" |

1019 | n/a | "Return the absolute value of the float x.") |

1020 | n/a | |

1021 | n/a | /*[clinic input] |

1022 | n/a | math.floor |

1023 | n/a | |

1024 | n/a | x as number: object |

1025 | n/a | / |

1026 | n/a | |

1027 | n/a | Return the floor of x as an Integral. |

1028 | n/a | |

1029 | n/a | This is the largest integer <= x. |

1030 | n/a | [clinic start generated code]*/ |

1031 | n/a | |

1032 | n/a | static PyObject * |

1033 | n/a | math_floor(PyObject *module, PyObject *number) |

1034 | n/a | /*[clinic end generated code: output=c6a65c4884884b8a input=63af6b5d7ebcc3d6]*/ |

1035 | n/a | { |

1036 | n/a | _Py_IDENTIFIER(__floor__); |

1037 | n/a | PyObject *method, *result; |

1038 | n/a | |

1039 | n/a | method = _PyObject_LookupSpecial(number, &PyId___floor__); |

1040 | n/a | if (method == NULL) { |

1041 | n/a | if (PyErr_Occurred()) |

1042 | n/a | return NULL; |

1043 | n/a | return math_1_to_int(number, floor, 0); |

1044 | n/a | } |

1045 | n/a | result = _PyObject_CallNoArg(method); |

1046 | n/a | Py_DECREF(method); |

1047 | n/a | return result; |

1048 | n/a | } |

1049 | n/a | |

1050 | n/a | FUNC1A(gamma, m_tgamma, |

1051 | n/a | "gamma($module, x, /)\n--\n\n" |

1052 | n/a | "Gamma function at x.") |

1053 | n/a | FUNC1A(lgamma, m_lgamma, |

1054 | n/a | "lgamma($module, x, /)\n--\n\n" |

1055 | n/a | "Natural logarithm of absolute value of Gamma function at x.") |

1056 | n/a | FUNC1(log1p, m_log1p, 0, |

1057 | n/a | "log1p($module, x, /)\n--\n\n" |

1058 | n/a | "Return the natural logarithm of 1+x (base e).\n\n" |

1059 | n/a | "The result is computed in a way which is accurate for x near zero.") |

1060 | n/a | FUNC1(sin, sin, 0, |

1061 | n/a | "sin($module, x, /)\n--\n\n" |

1062 | n/a | "Return the sine of x (measured in radians).") |

1063 | n/a | FUNC1(sinh, sinh, 1, |

1064 | n/a | "sinh($module, x, /)\n--\n\n" |

1065 | n/a | "Return the hyperbolic sine of x.") |

1066 | n/a | FUNC1(sqrt, sqrt, 0, |

1067 | n/a | "sqrt($module, x, /)\n--\n\n" |

1068 | n/a | "Return the square root of x.") |

1069 | n/a | FUNC1(tan, tan, 0, |

1070 | n/a | "tan($module, x, /)\n--\n\n" |

1071 | n/a | "Return the tangent of x (measured in radians).") |

1072 | n/a | FUNC1(tanh, tanh, 0, |

1073 | n/a | "tanh($module, x, /)\n--\n\n" |

1074 | n/a | "Return the hyperbolic tangent of x.") |

1075 | n/a | |

1076 | n/a | /* Precision summation function as msum() by Raymond Hettinger in |

1077 | n/a | <http://aspn.activestate.com/ASPN/Cookbook/Python/Recipe/393090>, |

1078 | n/a | enhanced with the exact partials sum and roundoff from Mark |

1079 | n/a | Dickinson's post at <http://bugs.python.org/file10357/msum4.py>. |

1080 | n/a | See those links for more details, proofs and other references. |

1081 | n/a | |

1082 | n/a | Note 1: IEEE 754R floating point semantics are assumed, |

1083 | n/a | but the current implementation does not re-establish special |

1084 | n/a | value semantics across iterations (i.e. handling -Inf + Inf). |

1085 | n/a | |

1086 | n/a | Note 2: No provision is made for intermediate overflow handling; |

1087 | n/a | therefore, sum([1e+308, 1e-308, 1e+308]) returns 1e+308 while |

1088 | n/a | sum([1e+308, 1e+308, 1e-308]) raises an OverflowError due to the |

1089 | n/a | overflow of the first partial sum. |

1090 | n/a | |

1091 | n/a | Note 3: The intermediate values lo, yr, and hi are declared volatile so |

1092 | n/a | aggressive compilers won't algebraically reduce lo to always be exactly 0.0. |

1093 | n/a | Also, the volatile declaration forces the values to be stored in memory as |

1094 | n/a | regular doubles instead of extended long precision (80-bit) values. This |

1095 | n/a | prevents double rounding because any addition or subtraction of two doubles |

1096 | n/a | can be resolved exactly into double-sized hi and lo values. As long as the |

1097 | n/a | hi value gets forced into a double before yr and lo are computed, the extra |

1098 | n/a | bits in downstream extended precision operations (x87 for example) will be |

1099 | n/a | exactly zero and therefore can be losslessly stored back into a double, |

1100 | n/a | thereby preventing double rounding. |

1101 | n/a | |

1102 | n/a | Note 4: A similar implementation is in Modules/cmathmodule.c. |

1103 | n/a | Be sure to update both when making changes. |

1104 | n/a | |

1105 | n/a | Note 5: The signature of math.fsum() differs from builtins.sum() |

1106 | n/a | because the start argument doesn't make sense in the context of |

1107 | n/a | accurate summation. Since the partials table is collapsed before |

1108 | n/a | returning a result, sum(seq2, start=sum(seq1)) may not equal the |

1109 | n/a | accurate result returned by sum(itertools.chain(seq1, seq2)). |

1110 | n/a | */ |

1111 | n/a | |

1112 | n/a | #define NUM_PARTIALS 32 /* initial partials array size, on stack */ |

1113 | n/a | |

1114 | n/a | /* Extend the partials array p[] by doubling its size. */ |

1115 | n/a | static int /* non-zero on error */ |

1116 | n/a | _fsum_realloc(double **p_ptr, Py_ssize_t n, |

1117 | n/a | double *ps, Py_ssize_t *m_ptr) |

1118 | n/a | { |

1119 | n/a | void *v = NULL; |

1120 | n/a | Py_ssize_t m = *m_ptr; |

1121 | n/a | |

1122 | n/a | m += m; /* double */ |

1123 | n/a | if (n < m && (size_t)m < ((size_t)PY_SSIZE_T_MAX / sizeof(double))) { |

1124 | n/a | double *p = *p_ptr; |

1125 | n/a | if (p == ps) { |

1126 | n/a | v = PyMem_Malloc(sizeof(double) * m); |

1127 | n/a | if (v != NULL) |

1128 | n/a | memcpy(v, ps, sizeof(double) * n); |

1129 | n/a | } |

1130 | n/a | else |

1131 | n/a | v = PyMem_Realloc(p, sizeof(double) * m); |

1132 | n/a | } |

1133 | n/a | if (v == NULL) { /* size overflow or no memory */ |

1134 | n/a | PyErr_SetString(PyExc_MemoryError, "math.fsum partials"); |

1135 | n/a | return 1; |

1136 | n/a | } |

1137 | n/a | *p_ptr = (double*) v; |

1138 | n/a | *m_ptr = m; |

1139 | n/a | return 0; |

1140 | n/a | } |

1141 | n/a | |

1142 | n/a | /* Full precision summation of a sequence of floats. |

1143 | n/a | |

1144 | n/a | def msum(iterable): |

1145 | n/a | partials = [] # sorted, non-overlapping partial sums |

1146 | n/a | for x in iterable: |

1147 | n/a | i = 0 |

1148 | n/a | for y in partials: |

1149 | n/a | if abs(x) < abs(y): |

1150 | n/a | x, y = y, x |

1151 | n/a | hi = x + y |

1152 | n/a | lo = y - (hi - x) |

1153 | n/a | if lo: |

1154 | n/a | partials[i] = lo |

1155 | n/a | i += 1 |

1156 | n/a | x = hi |

1157 | n/a | partials[i:] = [x] |

1158 | n/a | return sum_exact(partials) |

1159 | n/a | |

1160 | n/a | Rounded x+y stored in hi with the roundoff stored in lo. Together hi+lo |

1161 | n/a | are exactly equal to x+y. The inner loop applies hi/lo summation to each |

1162 | n/a | partial so that the list of partial sums remains exact. |

1163 | n/a | |

1164 | n/a | Sum_exact() adds the partial sums exactly and correctly rounds the final |

1165 | n/a | result (using the round-half-to-even rule). The items in partials remain |

1166 | n/a | non-zero, non-special, non-overlapping and strictly increasing in |

1167 | n/a | magnitude, but possibly not all having the same sign. |

1168 | n/a | |

1169 | n/a | Depends on IEEE 754 arithmetic guarantees and half-even rounding. |

1170 | n/a | */ |

1171 | n/a | |

1172 | n/a | /*[clinic input] |

1173 | n/a | math.fsum |

1174 | n/a | |

1175 | n/a | seq: object |

1176 | n/a | / |

1177 | n/a | |

1178 | n/a | Return an accurate floating point sum of values in the iterable seq. |

1179 | n/a | |

1180 | n/a | Assumes IEEE-754 floating point arithmetic. |

1181 | n/a | [clinic start generated code]*/ |

1182 | n/a | |

1183 | n/a | static PyObject * |

1184 | n/a | math_fsum(PyObject *module, PyObject *seq) |

1185 | n/a | /*[clinic end generated code: output=ba5c672b87fe34fc input=c51b7d8caf6f6e82]*/ |

1186 | n/a | { |

1187 | n/a | PyObject *item, *iter, *sum = NULL; |

1188 | n/a | Py_ssize_t i, j, n = 0, m = NUM_PARTIALS; |

1189 | n/a | double x, y, t, ps[NUM_PARTIALS], *p = ps; |

1190 | n/a | double xsave, special_sum = 0.0, inf_sum = 0.0; |

1191 | n/a | volatile double hi, yr, lo; |

1192 | n/a | |

1193 | n/a | iter = PyObject_GetIter(seq); |

1194 | n/a | if (iter == NULL) |

1195 | n/a | return NULL; |

1196 | n/a | |

1197 | n/a | PyFPE_START_PROTECT("fsum", Py_DECREF(iter); return NULL) |

1198 | n/a | |

1199 | n/a | for(;;) { /* for x in iterable */ |

1200 | n/a | assert(0 <= n && n <= m); |

1201 | n/a | assert((m == NUM_PARTIALS && p == ps) || |

1202 | n/a | (m > NUM_PARTIALS && p != NULL)); |

1203 | n/a | |

1204 | n/a | item = PyIter_Next(iter); |

1205 | n/a | if (item == NULL) { |

1206 | n/a | if (PyErr_Occurred()) |

1207 | n/a | goto _fsum_error; |

1208 | n/a | break; |

1209 | n/a | } |

1210 | n/a | x = PyFloat_AsDouble(item); |

1211 | n/a | Py_DECREF(item); |

1212 | n/a | if (PyErr_Occurred()) |

1213 | n/a | goto _fsum_error; |

1214 | n/a | |

1215 | n/a | xsave = x; |

1216 | n/a | for (i = j = 0; j < n; j++) { /* for y in partials */ |

1217 | n/a | y = p[j]; |

1218 | n/a | if (fabs(x) < fabs(y)) { |

1219 | n/a | t = x; x = y; y = t; |

1220 | n/a | } |

1221 | n/a | hi = x + y; |

1222 | n/a | yr = hi - x; |

1223 | n/a | lo = y - yr; |

1224 | n/a | if (lo != 0.0) |

1225 | n/a | p[i++] = lo; |

1226 | n/a | x = hi; |

1227 | n/a | } |

1228 | n/a | |

1229 | n/a | n = i; /* ps[i:] = [x] */ |

1230 | n/a | if (x != 0.0) { |

1231 | n/a | if (! Py_IS_FINITE(x)) { |

1232 | n/a | /* a nonfinite x could arise either as |

1233 | n/a | a result of intermediate overflow, or |

1234 | n/a | as a result of a nan or inf in the |

1235 | n/a | summands */ |

1236 | n/a | if (Py_IS_FINITE(xsave)) { |

1237 | n/a | PyErr_SetString(PyExc_OverflowError, |

1238 | n/a | "intermediate overflow in fsum"); |

1239 | n/a | goto _fsum_error; |

1240 | n/a | } |

1241 | n/a | if (Py_IS_INFINITY(xsave)) |

1242 | n/a | inf_sum += xsave; |

1243 | n/a | special_sum += xsave; |

1244 | n/a | /* reset partials */ |

1245 | n/a | n = 0; |

1246 | n/a | } |

1247 | n/a | else if (n >= m && _fsum_realloc(&p, n, ps, &m)) |

1248 | n/a | goto _fsum_error; |

1249 | n/a | else |

1250 | n/a | p[n++] = x; |

1251 | n/a | } |

1252 | n/a | } |

1253 | n/a | |

1254 | n/a | if (special_sum != 0.0) { |

1255 | n/a | if (Py_IS_NAN(inf_sum)) |

1256 | n/a | PyErr_SetString(PyExc_ValueError, |

1257 | n/a | "-inf + inf in fsum"); |

1258 | n/a | else |

1259 | n/a | sum = PyFloat_FromDouble(special_sum); |

1260 | n/a | goto _fsum_error; |

1261 | n/a | } |

1262 | n/a | |

1263 | n/a | hi = 0.0; |

1264 | n/a | if (n > 0) { |

1265 | n/a | hi = p[--n]; |

1266 | n/a | /* sum_exact(ps, hi) from the top, stop when the sum becomes |

1267 | n/a | inexact. */ |

1268 | n/a | while (n > 0) { |

1269 | n/a | x = hi; |

1270 | n/a | y = p[--n]; |

1271 | n/a | assert(fabs(y) < fabs(x)); |

1272 | n/a | hi = x + y; |

1273 | n/a | yr = hi - x; |

1274 | n/a | lo = y - yr; |

1275 | n/a | if (lo != 0.0) |

1276 | n/a | break; |

1277 | n/a | } |

1278 | n/a | /* Make half-even rounding work across multiple partials. |

1279 | n/a | Needed so that sum([1e-16, 1, 1e16]) will round-up the last |

1280 | n/a | digit to two instead of down to zero (the 1e-16 makes the 1 |

1281 | n/a | slightly closer to two). With a potential 1 ULP rounding |

1282 | n/a | error fixed-up, math.fsum() can guarantee commutativity. */ |

1283 | n/a | if (n > 0 && ((lo < 0.0 && p[n-1] < 0.0) || |

1284 | n/a | (lo > 0.0 && p[n-1] > 0.0))) { |

1285 | n/a | y = lo * 2.0; |

1286 | n/a | x = hi + y; |

1287 | n/a | yr = x - hi; |

1288 | n/a | if (y == yr) |

1289 | n/a | hi = x; |

1290 | n/a | } |

1291 | n/a | } |

1292 | n/a | sum = PyFloat_FromDouble(hi); |

1293 | n/a | |

1294 | n/a | _fsum_error: |

1295 | n/a | PyFPE_END_PROTECT(hi) |

1296 | n/a | Py_DECREF(iter); |

1297 | n/a | if (p != ps) |

1298 | n/a | PyMem_Free(p); |

1299 | n/a | return sum; |

1300 | n/a | } |

1301 | n/a | |

1302 | n/a | #undef NUM_PARTIALS |

1303 | n/a | |

1304 | n/a | |

1305 | n/a | /* Return the smallest integer k such that n < 2**k, or 0 if n == 0. |

1306 | n/a | * Equivalent to floor(lg(x))+1. Also equivalent to: bitwidth_of_type - |

1307 | n/a | * count_leading_zero_bits(x) |

1308 | n/a | */ |

1309 | n/a | |

1310 | n/a | /* XXX: This routine does more or less the same thing as |

1311 | n/a | * bits_in_digit() in Objects/longobject.c. Someday it would be nice to |

1312 | n/a | * consolidate them. On BSD, there's a library function called fls() |

1313 | n/a | * that we could use, and GCC provides __builtin_clz(). |

1314 | n/a | */ |

1315 | n/a | |

1316 | n/a | static unsigned long |

1317 | n/a | bit_length(unsigned long n) |

1318 | n/a | { |

1319 | n/a | unsigned long len = 0; |

1320 | n/a | while (n != 0) { |

1321 | n/a | ++len; |

1322 | n/a | n >>= 1; |

1323 | n/a | } |

1324 | n/a | return len; |

1325 | n/a | } |

1326 | n/a | |

1327 | n/a | static unsigned long |

1328 | n/a | count_set_bits(unsigned long n) |

1329 | n/a | { |

1330 | n/a | unsigned long count = 0; |

1331 | n/a | while (n != 0) { |

1332 | n/a | ++count; |

1333 | n/a | n &= n - 1; /* clear least significant bit */ |

1334 | n/a | } |

1335 | n/a | return count; |

1336 | n/a | } |

1337 | n/a | |

1338 | n/a | /* Divide-and-conquer factorial algorithm |

1339 | n/a | * |

1340 | n/a | * Based on the formula and pseudo-code provided at: |

1341 | n/a | * http://www.luschny.de/math/factorial/binarysplitfact.html |

1342 | n/a | * |

1343 | n/a | * Faster algorithms exist, but they're more complicated and depend on |

1344 | n/a | * a fast prime factorization algorithm. |

1345 | n/a | * |

1346 | n/a | * Notes on the algorithm |

1347 | n/a | * ---------------------- |

1348 | n/a | * |

1349 | n/a | * factorial(n) is written in the form 2**k * m, with m odd. k and m are |

1350 | n/a | * computed separately, and then combined using a left shift. |

1351 | n/a | * |

1352 | n/a | * The function factorial_odd_part computes the odd part m (i.e., the greatest |

1353 | n/a | * odd divisor) of factorial(n), using the formula: |

1354 | n/a | * |

1355 | n/a | * factorial_odd_part(n) = |

1356 | n/a | * |

1357 | n/a | * product_{i >= 0} product_{0 < j <= n / 2**i, j odd} j |

1358 | n/a | * |

1359 | n/a | * Example: factorial_odd_part(20) = |

1360 | n/a | * |

1361 | n/a | * (1) * |

1362 | n/a | * (1) * |

1363 | n/a | * (1 * 3 * 5) * |

1364 | n/a | * (1 * 3 * 5 * 7 * 9) |

1365 | n/a | * (1 * 3 * 5 * 7 * 9 * 11 * 13 * 15 * 17 * 19) |

1366 | n/a | * |

1367 | n/a | * Here i goes from large to small: the first term corresponds to i=4 (any |

1368 | n/a | * larger i gives an empty product), and the last term corresponds to i=0. |

1369 | n/a | * Each term can be computed from the last by multiplying by the extra odd |

1370 | n/a | * numbers required: e.g., to get from the penultimate term to the last one, |

1371 | n/a | * we multiply by (11 * 13 * 15 * 17 * 19). |

1372 | n/a | * |

1373 | n/a | * To see a hint of why this formula works, here are the same numbers as above |

1374 | n/a | * but with the even parts (i.e., the appropriate powers of 2) included. For |

1375 | n/a | * each subterm in the product for i, we multiply that subterm by 2**i: |

1376 | n/a | * |

1377 | n/a | * factorial(20) = |

1378 | n/a | * |

1379 | n/a | * (16) * |

1380 | n/a | * (8) * |

1381 | n/a | * (4 * 12 * 20) * |

1382 | n/a | * (2 * 6 * 10 * 14 * 18) * |

1383 | n/a | * (1 * 3 * 5 * 7 * 9 * 11 * 13 * 15 * 17 * 19) |

1384 | n/a | * |

1385 | n/a | * The factorial_partial_product function computes the product of all odd j in |

1386 | n/a | * range(start, stop) for given start and stop. It's used to compute the |

1387 | n/a | * partial products like (11 * 13 * 15 * 17 * 19) in the example above. It |

1388 | n/a | * operates recursively, repeatedly splitting the range into two roughly equal |

1389 | n/a | * pieces until the subranges are small enough to be computed using only C |

1390 | n/a | * integer arithmetic. |

1391 | n/a | * |

1392 | n/a | * The two-valuation k (i.e., the exponent of the largest power of 2 dividing |

1393 | n/a | * the factorial) is computed independently in the main math_factorial |

1394 | n/a | * function. By standard results, its value is: |

1395 | n/a | * |

1396 | n/a | * two_valuation = n//2 + n//4 + n//8 + .... |

1397 | n/a | * |

1398 | n/a | * It can be shown (e.g., by complete induction on n) that two_valuation is |

1399 | n/a | * equal to n - count_set_bits(n), where count_set_bits(n) gives the number of |

1400 | n/a | * '1'-bits in the binary expansion of n. |

1401 | n/a | */ |

1402 | n/a | |

1403 | n/a | /* factorial_partial_product: Compute product(range(start, stop, 2)) using |

1404 | n/a | * divide and conquer. Assumes start and stop are odd and stop > start. |

1405 | n/a | * max_bits must be >= bit_length(stop - 2). */ |

1406 | n/a | |

1407 | n/a | static PyObject * |

1408 | n/a | factorial_partial_product(unsigned long start, unsigned long stop, |

1409 | n/a | unsigned long max_bits) |

1410 | n/a | { |

1411 | n/a | unsigned long midpoint, num_operands; |

1412 | n/a | PyObject *left = NULL, *right = NULL, *result = NULL; |

1413 | n/a | |

1414 | n/a | /* If the return value will fit an unsigned long, then we can |

1415 | n/a | * multiply in a tight, fast loop where each multiply is O(1). |

1416 | n/a | * Compute an upper bound on the number of bits required to store |

1417 | n/a | * the answer. |

1418 | n/a | * |

1419 | n/a | * Storing some integer z requires floor(lg(z))+1 bits, which is |

1420 | n/a | * conveniently the value returned by bit_length(z). The |

1421 | n/a | * product x*y will require at most |

1422 | n/a | * bit_length(x) + bit_length(y) bits to store, based |

1423 | n/a | * on the idea that lg product = lg x + lg y. |

1424 | n/a | * |

1425 | n/a | * We know that stop - 2 is the largest number to be multiplied. From |

1426 | n/a | * there, we have: bit_length(answer) <= num_operands * |

1427 | n/a | * bit_length(stop - 2) |

1428 | n/a | */ |

1429 | n/a | |

1430 | n/a | num_operands = (stop - start) / 2; |

1431 | n/a | /* The "num_operands <= 8 * SIZEOF_LONG" check guards against the |

1432 | n/a | * unlikely case of an overflow in num_operands * max_bits. */ |

1433 | n/a | if (num_operands <= 8 * SIZEOF_LONG && |

1434 | n/a | num_operands * max_bits <= 8 * SIZEOF_LONG) { |

1435 | n/a | unsigned long j, total; |

1436 | n/a | for (total = start, j = start + 2; j < stop; j += 2) |

1437 | n/a | total *= j; |

1438 | n/a | return PyLong_FromUnsignedLong(total); |

1439 | n/a | } |

1440 | n/a | |

1441 | n/a | /* find midpoint of range(start, stop), rounded up to next odd number. */ |

1442 | n/a | midpoint = (start + num_operands) | 1; |

1443 | n/a | left = factorial_partial_product(start, midpoint, |

1444 | n/a | bit_length(midpoint - 2)); |

1445 | n/a | if (left == NULL) |

1446 | n/a | goto error; |

1447 | n/a | right = factorial_partial_product(midpoint, stop, max_bits); |

1448 | n/a | if (right == NULL) |

1449 | n/a | goto error; |

1450 | n/a | result = PyNumber_Multiply(left, right); |

1451 | n/a | |

1452 | n/a | error: |

1453 | n/a | Py_XDECREF(left); |

1454 | n/a | Py_XDECREF(right); |

1455 | n/a | return result; |

1456 | n/a | } |

1457 | n/a | |

1458 | n/a | /* factorial_odd_part: compute the odd part of factorial(n). */ |

1459 | n/a | |

1460 | n/a | static PyObject * |

1461 | n/a | factorial_odd_part(unsigned long n) |

1462 | n/a | { |

1463 | n/a | long i; |

1464 | n/a | unsigned long v, lower, upper; |

1465 | n/a | PyObject *partial, *tmp, *inner, *outer; |

1466 | n/a | |

1467 | n/a | inner = PyLong_FromLong(1); |

1468 | n/a | if (inner == NULL) |

1469 | n/a | return NULL; |

1470 | n/a | outer = inner; |

1471 | n/a | Py_INCREF(outer); |

1472 | n/a | |

1473 | n/a | upper = 3; |

1474 | n/a | for (i = bit_length(n) - 2; i >= 0; i--) { |

1475 | n/a | v = n >> i; |

1476 | n/a | if (v <= 2) |

1477 | n/a | continue; |

1478 | n/a | lower = upper; |

1479 | n/a | /* (v + 1) | 1 = least odd integer strictly larger than n / 2**i */ |

1480 | n/a | upper = (v + 1) | 1; |

1481 | n/a | /* Here inner is the product of all odd integers j in the range (0, |

1482 | n/a | n/2**(i+1)]. The factorial_partial_product call below gives the |

1483 | n/a | product of all odd integers j in the range (n/2**(i+1), n/2**i]. */ |

1484 | n/a | partial = factorial_partial_product(lower, upper, bit_length(upper-2)); |

1485 | n/a | /* inner *= partial */ |

1486 | n/a | if (partial == NULL) |

1487 | n/a | goto error; |

1488 | n/a | tmp = PyNumber_Multiply(inner, partial); |

1489 | n/a | Py_DECREF(partial); |

1490 | n/a | if (tmp == NULL) |

1491 | n/a | goto error; |

1492 | n/a | Py_DECREF(inner); |

1493 | n/a | inner = tmp; |

1494 | n/a | /* Now inner is the product of all odd integers j in the range (0, |

1495 | n/a | n/2**i], giving the inner product in the formula above. */ |

1496 | n/a | |

1497 | n/a | /* outer *= inner; */ |

1498 | n/a | tmp = PyNumber_Multiply(outer, inner); |

1499 | n/a | if (tmp == NULL) |

1500 | n/a | goto error; |

1501 | n/a | Py_DECREF(outer); |

1502 | n/a | outer = tmp; |

1503 | n/a | } |

1504 | n/a | Py_DECREF(inner); |

1505 | n/a | return outer; |

1506 | n/a | |

1507 | n/a | error: |

1508 | n/a | Py_DECREF(outer); |

1509 | n/a | Py_DECREF(inner); |

1510 | n/a | return NULL; |

1511 | n/a | } |

1512 | n/a | |

1513 | n/a | |

1514 | n/a | /* Lookup table for small factorial values */ |

1515 | n/a | |

1516 | n/a | static const unsigned long SmallFactorials[] = { |

1517 | n/a | 1, 1, 2, 6, 24, 120, 720, 5040, 40320, |

1518 | n/a | 362880, 3628800, 39916800, 479001600, |

1519 | n/a | #if SIZEOF_LONG >= 8 |

1520 | n/a | 6227020800, 87178291200, 1307674368000, |

1521 | n/a | 20922789888000, 355687428096000, 6402373705728000, |

1522 | n/a | 121645100408832000, 2432902008176640000 |

1523 | n/a | #endif |

1524 | n/a | }; |

1525 | n/a | |

1526 | n/a | /*[clinic input] |

1527 | n/a | math.factorial |

1528 | n/a | |

1529 | n/a | x as arg: object |

1530 | n/a | / |

1531 | n/a | |

1532 | n/a | Find x!. |

1533 | n/a | |

1534 | n/a | Raise a ValueError if x is negative or non-integral. |

1535 | n/a | [clinic start generated code]*/ |

1536 | n/a | |

1537 | n/a | static PyObject * |

1538 | n/a | math_factorial(PyObject *module, PyObject *arg) |

1539 | n/a | /*[clinic end generated code: output=6686f26fae00e9ca input=6d1c8105c0d91fb4]*/ |

1540 | n/a | { |

1541 | n/a | long x; |

1542 | n/a | int overflow; |

1543 | n/a | PyObject *result, *odd_part, *two_valuation; |

1544 | n/a | |

1545 | n/a | if (PyFloat_Check(arg)) { |

1546 | n/a | PyObject *lx; |

1547 | n/a | double dx = PyFloat_AS_DOUBLE((PyFloatObject *)arg); |

1548 | n/a | if (!(Py_IS_FINITE(dx) && dx == floor(dx))) { |

1549 | n/a | PyErr_SetString(PyExc_ValueError, |

1550 | n/a | "factorial() only accepts integral values"); |

1551 | n/a | return NULL; |

1552 | n/a | } |

1553 | n/a | lx = PyLong_FromDouble(dx); |

1554 | n/a | if (lx == NULL) |

1555 | n/a | return NULL; |

1556 | n/a | x = PyLong_AsLongAndOverflow(lx, &overflow); |

1557 | n/a | Py_DECREF(lx); |

1558 | n/a | } |

1559 | n/a | else |

1560 | n/a | x = PyLong_AsLongAndOverflow(arg, &overflow); |

1561 | n/a | |

1562 | n/a | if (x == -1 && PyErr_Occurred()) { |

1563 | n/a | return NULL; |

1564 | n/a | } |

1565 | n/a | else if (overflow == 1) { |

1566 | n/a | PyErr_Format(PyExc_OverflowError, |

1567 | n/a | "factorial() argument should not exceed %ld", |

1568 | n/a | LONG_MAX); |

1569 | n/a | return NULL; |

1570 | n/a | } |

1571 | n/a | else if (overflow == -1 || x < 0) { |

1572 | n/a | PyErr_SetString(PyExc_ValueError, |

1573 | n/a | "factorial() not defined for negative values"); |

1574 | n/a | return NULL; |

1575 | n/a | } |

1576 | n/a | |

1577 | n/a | /* use lookup table if x is small */ |

1578 | n/a | if (x < (long)Py_ARRAY_LENGTH(SmallFactorials)) |

1579 | n/a | return PyLong_FromUnsignedLong(SmallFactorials[x]); |

1580 | n/a | |

1581 | n/a | /* else express in the form odd_part * 2**two_valuation, and compute as |

1582 | n/a | odd_part << two_valuation. */ |

1583 | n/a | odd_part = factorial_odd_part(x); |

1584 | n/a | if (odd_part == NULL) |

1585 | n/a | return NULL; |

1586 | n/a | two_valuation = PyLong_FromLong(x - count_set_bits(x)); |

1587 | n/a | if (two_valuation == NULL) { |

1588 | n/a | Py_DECREF(odd_part); |

1589 | n/a | return NULL; |

1590 | n/a | } |

1591 | n/a | result = PyNumber_Lshift(odd_part, two_valuation); |

1592 | n/a | Py_DECREF(two_valuation); |

1593 | n/a | Py_DECREF(odd_part); |

1594 | n/a | return result; |

1595 | n/a | } |

1596 | n/a | |

1597 | n/a | |

1598 | n/a | /*[clinic input] |

1599 | n/a | math.trunc |

1600 | n/a | |

1601 | n/a | x: object |

1602 | n/a | / |

1603 | n/a | |

1604 | n/a | Truncates the Real x to the nearest Integral toward 0. |

1605 | n/a | |

1606 | n/a | Uses the __trunc__ magic method. |

1607 | n/a | [clinic start generated code]*/ |

1608 | n/a | |

1609 | n/a | static PyObject * |

1610 | n/a | math_trunc(PyObject *module, PyObject *x) |

1611 | n/a | /*[clinic end generated code: output=34b9697b707e1031 input=2168b34e0a09134d]*/ |

1612 | n/a | { |

1613 | n/a | _Py_IDENTIFIER(__trunc__); |

1614 | n/a | PyObject *trunc, *result; |

1615 | n/a | |

1616 | n/a | if (Py_TYPE(x)->tp_dict == NULL) { |

1617 | n/a | if (PyType_Ready(Py_TYPE(x)) < 0) |

1618 | n/a | return NULL; |

1619 | n/a | } |

1620 | n/a | |

1621 | n/a | trunc = _PyObject_LookupSpecial(x, &PyId___trunc__); |

1622 | n/a | if (trunc == NULL) { |

1623 | n/a | if (!PyErr_Occurred()) |

1624 | n/a | PyErr_Format(PyExc_TypeError, |

1625 | n/a | "type %.100s doesn't define __trunc__ method", |

1626 | n/a | Py_TYPE(x)->tp_name); |

1627 | n/a | return NULL; |

1628 | n/a | } |

1629 | n/a | result = _PyObject_CallNoArg(trunc); |

1630 | n/a | Py_DECREF(trunc); |

1631 | n/a | return result; |

1632 | n/a | } |

1633 | n/a | |

1634 | n/a | |

1635 | n/a | /*[clinic input] |

1636 | n/a | math.frexp |

1637 | n/a | |

1638 | n/a | x: double |

1639 | n/a | / |

1640 | n/a | |

1641 | n/a | Return the mantissa and exponent of x, as pair (m, e). |

1642 | n/a | |

1643 | n/a | m is a float and e is an int, such that x = m * 2.**e. |

1644 | n/a | If x is 0, m and e are both 0. Else 0.5 <= abs(m) < 1.0. |

1645 | n/a | [clinic start generated code]*/ |

1646 | n/a | |

1647 | n/a | static PyObject * |

1648 | n/a | math_frexp_impl(PyObject *module, double x) |

1649 | n/a | /*[clinic end generated code: output=03e30d252a15ad4a input=96251c9e208bc6e9]*/ |

1650 | n/a | { |

1651 | n/a | int i; |

1652 | n/a | /* deal with special cases directly, to sidestep platform |

1653 | n/a | differences */ |

1654 | n/a | if (Py_IS_NAN(x) || Py_IS_INFINITY(x) || !x) { |

1655 | n/a | i = 0; |

1656 | n/a | } |

1657 | n/a | else { |

1658 | n/a | PyFPE_START_PROTECT("in math_frexp", return 0); |

1659 | n/a | x = frexp(x, &i); |

1660 | n/a | PyFPE_END_PROTECT(x); |

1661 | n/a | } |

1662 | n/a | return Py_BuildValue("(di)", x, i); |

1663 | n/a | } |

1664 | n/a | |

1665 | n/a | |

1666 | n/a | /*[clinic input] |

1667 | n/a | math.ldexp |

1668 | n/a | |

1669 | n/a | x: double |

1670 | n/a | i: object |

1671 | n/a | / |

1672 | n/a | |

1673 | n/a | Return x * (2**i). |

1674 | n/a | |

1675 | n/a | This is essentially the inverse of frexp(). |

1676 | n/a | [clinic start generated code]*/ |

1677 | n/a | |

1678 | n/a | static PyObject * |

1679 | n/a | math_ldexp_impl(PyObject *module, double x, PyObject *i) |

1680 | n/a | /*[clinic end generated code: output=b6892f3c2df9cc6a input=17d5970c1a40a8c1]*/ |

1681 | n/a | { |

1682 | n/a | double r; |

1683 | n/a | long exp; |

1684 | n/a | int overflow; |

1685 | n/a | |

1686 | n/a | if (PyLong_Check(i)) { |

1687 | n/a | /* on overflow, replace exponent with either LONG_MAX |

1688 | n/a | or LONG_MIN, depending on the sign. */ |

1689 | n/a | exp = PyLong_AsLongAndOverflow(i, &overflow); |

1690 | n/a | if (exp == -1 && PyErr_Occurred()) |

1691 | n/a | return NULL; |

1692 | n/a | if (overflow) |

1693 | n/a | exp = overflow < 0 ? LONG_MIN : LONG_MAX; |

1694 | n/a | } |

1695 | n/a | else { |

1696 | n/a | PyErr_SetString(PyExc_TypeError, |

1697 | n/a | "Expected an int as second argument to ldexp."); |

1698 | n/a | return NULL; |

1699 | n/a | } |

1700 | n/a | |

1701 | n/a | if (x == 0. || !Py_IS_FINITE(x)) { |

1702 | n/a | /* NaNs, zeros and infinities are returned unchanged */ |

1703 | n/a | r = x; |

1704 | n/a | errno = 0; |

1705 | n/a | } else if (exp > INT_MAX) { |

1706 | n/a | /* overflow */ |

1707 | n/a | r = copysign(Py_HUGE_VAL, x); |

1708 | n/a | errno = ERANGE; |

1709 | n/a | } else if (exp < INT_MIN) { |

1710 | n/a | /* underflow to +-0 */ |

1711 | n/a | r = copysign(0., x); |

1712 | n/a | errno = 0; |

1713 | n/a | } else { |

1714 | n/a | errno = 0; |

1715 | n/a | PyFPE_START_PROTECT("in math_ldexp", return 0); |

1716 | n/a | r = ldexp(x, (int)exp); |

1717 | n/a | PyFPE_END_PROTECT(r); |

1718 | n/a | if (Py_IS_INFINITY(r)) |

1719 | n/a | errno = ERANGE; |

1720 | n/a | } |

1721 | n/a | |

1722 | n/a | if (errno && is_error(r)) |

1723 | n/a | return NULL; |

1724 | n/a | return PyFloat_FromDouble(r); |

1725 | n/a | } |

1726 | n/a | |

1727 | n/a | |

1728 | n/a | /*[clinic input] |

1729 | n/a | math.modf |

1730 | n/a | |

1731 | n/a | x: double |

1732 | n/a | / |

1733 | n/a | |

1734 | n/a | Return the fractional and integer parts of x. |

1735 | n/a | |

1736 | n/a | Both results carry the sign of x and are floats. |

1737 | n/a | [clinic start generated code]*/ |

1738 | n/a | |

1739 | n/a | static PyObject * |

1740 | n/a | math_modf_impl(PyObject *module, double x) |

1741 | n/a | /*[clinic end generated code: output=90cee0260014c3c0 input=b4cfb6786afd9035]*/ |

1742 | n/a | { |

1743 | n/a | double y; |

1744 | n/a | /* some platforms don't do the right thing for NaNs and |

1745 | n/a | infinities, so we take care of special cases directly. */ |

1746 | n/a | if (!Py_IS_FINITE(x)) { |

1747 | n/a | if (Py_IS_INFINITY(x)) |

1748 | n/a | return Py_BuildValue("(dd)", copysign(0., x), x); |

1749 | n/a | else if (Py_IS_NAN(x)) |

1750 | n/a | return Py_BuildValue("(dd)", x, x); |

1751 | n/a | } |

1752 | n/a | |

1753 | n/a | errno = 0; |

1754 | n/a | PyFPE_START_PROTECT("in math_modf", return 0); |

1755 | n/a | x = modf(x, &y); |

1756 | n/a | PyFPE_END_PROTECT(x); |

1757 | n/a | return Py_BuildValue("(dd)", x, y); |

1758 | n/a | } |

1759 | n/a | |

1760 | n/a | |

1761 | n/a | /* A decent logarithm is easy to compute even for huge ints, but libm can't |

1762 | n/a | do that by itself -- loghelper can. func is log or log10, and name is |

1763 | n/a | "log" or "log10". Note that overflow of the result isn't possible: an int |

1764 | n/a | can contain no more than INT_MAX * SHIFT bits, so has value certainly less |

1765 | n/a | than 2**(2**64 * 2**16) == 2**2**80, and log2 of that is 2**80, which is |

1766 | n/a | small enough to fit in an IEEE single. log and log10 are even smaller. |

1767 | n/a | However, intermediate overflow is possible for an int if the number of bits |

1768 | n/a | in that int is larger than PY_SSIZE_T_MAX. */ |

1769 | n/a | |

1770 | n/a | static PyObject* |

1771 | n/a | loghelper(PyObject* arg, double (*func)(double), const char *funcname) |

1772 | n/a | { |

1773 | n/a | /* If it is int, do it ourselves. */ |

1774 | n/a | if (PyLong_Check(arg)) { |

1775 | n/a | double x, result; |

1776 | n/a | Py_ssize_t e; |

1777 | n/a | |

1778 | n/a | /* Negative or zero inputs give a ValueError. */ |

1779 | n/a | if (Py_SIZE(arg) <= 0) { |

1780 | n/a | PyErr_SetString(PyExc_ValueError, |

1781 | n/a | "math domain error"); |

1782 | n/a | return NULL; |

1783 | n/a | } |

1784 | n/a | |

1785 | n/a | x = PyLong_AsDouble(arg); |

1786 | n/a | if (x == -1.0 && PyErr_Occurred()) { |

1787 | n/a | if (!PyErr_ExceptionMatches(PyExc_OverflowError)) |

1788 | n/a | return NULL; |

1789 | n/a | /* Here the conversion to double overflowed, but it's possible |

1790 | n/a | to compute the log anyway. Clear the exception and continue. */ |

1791 | n/a | PyErr_Clear(); |

1792 | n/a | x = _PyLong_Frexp((PyLongObject *)arg, &e); |

1793 | n/a | if (x == -1.0 && PyErr_Occurred()) |

1794 | n/a | return NULL; |

1795 | n/a | /* Value is ~= x * 2**e, so the log ~= log(x) + log(2) * e. */ |

1796 | n/a | result = func(x) + func(2.0) * e; |

1797 | n/a | } |

1798 | n/a | else |

1799 | n/a | /* Successfully converted x to a double. */ |

1800 | n/a | result = func(x); |

1801 | n/a | return PyFloat_FromDouble(result); |

1802 | n/a | } |

1803 | n/a | |

1804 | n/a | /* Else let libm handle it by itself. */ |

1805 | n/a | return math_1(arg, func, 0); |

1806 | n/a | } |

1807 | n/a | |

1808 | n/a | |

1809 | n/a | /*[clinic input] |

1810 | n/a | math.log |

1811 | n/a | |

1812 | n/a | x: object |

1813 | n/a | [ |

1814 | n/a | base: object(c_default="NULL") = math.e |

1815 | n/a | ] |

1816 | n/a | / |

1817 | n/a | |

1818 | n/a | Return the logarithm of x to the given base. |

1819 | n/a | |

1820 | n/a | If the base not specified, returns the natural logarithm (base e) of x. |

1821 | n/a | [clinic start generated code]*/ |

1822 | n/a | |

1823 | n/a | static PyObject * |

1824 | n/a | math_log_impl(PyObject *module, PyObject *x, int group_right_1, |

1825 | n/a | PyObject *base) |

1826 | n/a | /*[clinic end generated code: output=7b5a39e526b73fc9 input=0f62d5726cbfebbd]*/ |

1827 | n/a | { |

1828 | n/a | PyObject *num, *den; |

1829 | n/a | PyObject *ans; |

1830 | n/a | |

1831 | n/a | num = loghelper(x, m_log, "log"); |

1832 | n/a | if (num == NULL || base == NULL) |

1833 | n/a | return num; |

1834 | n/a | |

1835 | n/a | den = loghelper(base, m_log, "log"); |

1836 | n/a | if (den == NULL) { |

1837 | n/a | Py_DECREF(num); |

1838 | n/a | return NULL; |

1839 | n/a | } |

1840 | n/a | |

1841 | n/a | ans = PyNumber_TrueDivide(num, den); |

1842 | n/a | Py_DECREF(num); |

1843 | n/a | Py_DECREF(den); |

1844 | n/a | return ans; |

1845 | n/a | } |

1846 | n/a | |

1847 | n/a | |

1848 | n/a | /*[clinic input] |

1849 | n/a | math.log2 |

1850 | n/a | |

1851 | n/a | x: object |

1852 | n/a | / |

1853 | n/a | |

1854 | n/a | Return the base 2 logarithm of x. |

1855 | n/a | [clinic start generated code]*/ |

1856 | n/a | |

1857 | n/a | static PyObject * |

1858 | n/a | math_log2(PyObject *module, PyObject *x) |

1859 | n/a | /*[clinic end generated code: output=5425899a4d5d6acb input=08321262bae4f39b]*/ |

1860 | n/a | { |

1861 | n/a | return loghelper(x, m_log2, "log2"); |

1862 | n/a | } |

1863 | n/a | |

1864 | n/a | |

1865 | n/a | /*[clinic input] |

1866 | n/a | math.log10 |

1867 | n/a | |

1868 | n/a | x: object |

1869 | n/a | / |

1870 | n/a | |

1871 | n/a | Return the base 10 logarithm of x. |

1872 | n/a | [clinic start generated code]*/ |

1873 | n/a | |

1874 | n/a | static PyObject * |

1875 | n/a | math_log10(PyObject *module, PyObject *x) |

1876 | n/a | /*[clinic end generated code: output=be72a64617df9c6f input=b2469d02c6469e53]*/ |

1877 | n/a | { |

1878 | n/a | return loghelper(x, m_log10, "log10"); |

1879 | n/a | } |

1880 | n/a | |

1881 | n/a | |

1882 | n/a | /*[clinic input] |

1883 | n/a | math.fmod |

1884 | n/a | |

1885 | n/a | x: double |

1886 | n/a | y: double |

1887 | n/a | / |

1888 | n/a | |

1889 | n/a | Return fmod(x, y), according to platform C. |

1890 | n/a | |

1891 | n/a | x % y may differ. |

1892 | n/a | [clinic start generated code]*/ |

1893 | n/a | |

1894 | n/a | static PyObject * |

1895 | n/a | math_fmod_impl(PyObject *module, double x, double y) |

1896 | n/a | /*[clinic end generated code: output=7559d794343a27b5 input=4f84caa8cfc26a03]*/ |

1897 | n/a | { |

1898 | n/a | double r; |

1899 | n/a | /* fmod(x, +/-Inf) returns x for finite x. */ |

1900 | n/a | if (Py_IS_INFINITY(y) && Py_IS_FINITE(x)) |

1901 | n/a | return PyFloat_FromDouble(x); |

1902 | n/a | errno = 0; |

1903 | n/a | PyFPE_START_PROTECT("in math_fmod", return 0); |

1904 | n/a | r = fmod(x, y); |

1905 | n/a | PyFPE_END_PROTECT(r); |

1906 | n/a | if (Py_IS_NAN(r)) { |

1907 | n/a | if (!Py_IS_NAN(x) && !Py_IS_NAN(y)) |

1908 | n/a | errno = EDOM; |

1909 | n/a | else |

1910 | n/a | errno = 0; |

1911 | n/a | } |

1912 | n/a | if (errno && is_error(r)) |

1913 | n/a | return NULL; |

1914 | n/a | else |

1915 | n/a | return PyFloat_FromDouble(r); |

1916 | n/a | } |

1917 | n/a | |

1918 | n/a | |

1919 | n/a | /*[clinic input] |

1920 | n/a | math.hypot |

1921 | n/a | |

1922 | n/a | x: double |

1923 | n/a | y: double |

1924 | n/a | / |

1925 | n/a | |

1926 | n/a | Return the Euclidean distance, sqrt(x*x + y*y). |

1927 | n/a | [clinic start generated code]*/ |

1928 | n/a | |

1929 | n/a | static PyObject * |

1930 | n/a | math_hypot_impl(PyObject *module, double x, double y) |

1931 | n/a | /*[clinic end generated code: output=b7686e5be468ef87 input=7f8eea70406474aa]*/ |

1932 | n/a | { |

1933 | n/a | double r; |

1934 | n/a | /* hypot(x, +/-Inf) returns Inf, even if x is a NaN. */ |

1935 | n/a | if (Py_IS_INFINITY(x)) |

1936 | n/a | return PyFloat_FromDouble(fabs(x)); |

1937 | n/a | if (Py_IS_INFINITY(y)) |

1938 | n/a | return PyFloat_FromDouble(fabs(y)); |

1939 | n/a | errno = 0; |

1940 | n/a | PyFPE_START_PROTECT("in math_hypot", return 0); |

1941 | n/a | r = hypot(x, y); |

1942 | n/a | PyFPE_END_PROTECT(r); |

1943 | n/a | if (Py_IS_NAN(r)) { |

1944 | n/a | if (!Py_IS_NAN(x) && !Py_IS_NAN(y)) |

1945 | n/a | errno = EDOM; |

1946 | n/a | else |

1947 | n/a | errno = 0; |

1948 | n/a | } |

1949 | n/a | else if (Py_IS_INFINITY(r)) { |

1950 | n/a | if (Py_IS_FINITE(x) && Py_IS_FINITE(y)) |

1951 | n/a | errno = ERANGE; |

1952 | n/a | else |

1953 | n/a | errno = 0; |

1954 | n/a | } |

1955 | n/a | if (errno && is_error(r)) |

1956 | n/a | return NULL; |

1957 | n/a | else |

1958 | n/a | return PyFloat_FromDouble(r); |

1959 | n/a | } |

1960 | n/a | |

1961 | n/a | |

1962 | n/a | /* pow can't use math_2, but needs its own wrapper: the problem is |

1963 | n/a | that an infinite result can arise either as a result of overflow |

1964 | n/a | (in which case OverflowError should be raised) or as a result of |

1965 | n/a | e.g. 0.**-5. (for which ValueError needs to be raised.) |

1966 | n/a | */ |

1967 | n/a | |

1968 | n/a | /*[clinic input] |

1969 | n/a | math.pow |

1970 | n/a | |

1971 | n/a | x: double |

1972 | n/a | y: double |

1973 | n/a | / |

1974 | n/a | |

1975 | n/a | Return x**y (x to the power of y). |

1976 | n/a | [clinic start generated code]*/ |

1977 | n/a | |

1978 | n/a | static PyObject * |

1979 | n/a | math_pow_impl(PyObject *module, double x, double y) |

1980 | n/a | /*[clinic end generated code: output=fff93e65abccd6b0 input=c26f1f6075088bfd]*/ |

1981 | n/a | { |

1982 | n/a | double r; |

1983 | n/a | int odd_y; |

1984 | n/a | |

1985 | n/a | /* deal directly with IEEE specials, to cope with problems on various |

1986 | n/a | platforms whose semantics don't exactly match C99 */ |

1987 | n/a | r = 0.; /* silence compiler warning */ |

1988 | n/a | if (!Py_IS_FINITE(x) || !Py_IS_FINITE(y)) { |

1989 | n/a | errno = 0; |

1990 | n/a | if (Py_IS_NAN(x)) |

1991 | n/a | r = y == 0. ? 1. : x; /* NaN**0 = 1 */ |

1992 | n/a | else if (Py_IS_NAN(y)) |

1993 | n/a | r = x == 1. ? 1. : y; /* 1**NaN = 1 */ |

1994 | n/a | else if (Py_IS_INFINITY(x)) { |

1995 | n/a | odd_y = Py_IS_FINITE(y) && fmod(fabs(y), 2.0) == 1.0; |

1996 | n/a | if (y > 0.) |

1997 | n/a | r = odd_y ? x : fabs(x); |

1998 | n/a | else if (y == 0.) |

1999 | n/a | r = 1.; |

2000 | n/a | else /* y < 0. */ |

2001 | n/a | r = odd_y ? copysign(0., x) : 0.; |

2002 | n/a | } |

2003 | n/a | else if (Py_IS_INFINITY(y)) { |

2004 | n/a | if (fabs(x) == 1.0) |

2005 | n/a | r = 1.; |

2006 | n/a | else if (y > 0. && fabs(x) > 1.0) |

2007 | n/a | r = y; |

2008 | n/a | else if (y < 0. && fabs(x) < 1.0) { |

2009 | n/a | r = -y; /* result is +inf */ |

2010 | n/a | if (x == 0.) /* 0**-inf: divide-by-zero */ |

2011 | n/a | errno = EDOM; |

2012 | n/a | } |

2013 | n/a | else |

2014 | n/a | r = 0.; |

2015 | n/a | } |

2016 | n/a | } |

2017 | n/a | else { |

2018 | n/a | /* let libm handle finite**finite */ |

2019 | n/a | errno = 0; |

2020 | n/a | PyFPE_START_PROTECT("in math_pow", return 0); |

2021 | n/a | r = pow(x, y); |

2022 | n/a | PyFPE_END_PROTECT(r); |

2023 | n/a | /* a NaN result should arise only from (-ve)**(finite |

2024 | n/a | non-integer); in this case we want to raise ValueError. */ |

2025 | n/a | if (!Py_IS_FINITE(r)) { |

2026 | n/a | if (Py_IS_NAN(r)) { |

2027 | n/a | errno = EDOM; |

2028 | n/a | } |

2029 | n/a | /* |

2030 | n/a | an infinite result here arises either from: |

2031 | n/a | (A) (+/-0.)**negative (-> divide-by-zero) |

2032 | n/a | (B) overflow of x**y with x and y finite |

2033 | n/a | */ |

2034 | n/a | else if (Py_IS_INFINITY(r)) { |

2035 | n/a | if (x == 0.) |

2036 | n/a | errno = EDOM; |

2037 | n/a | else |

2038 | n/a | errno = ERANGE; |

2039 | n/a | } |

2040 | n/a | } |

2041 | n/a | } |

2042 | n/a | |

2043 | n/a | if (errno && is_error(r)) |

2044 | n/a | return NULL; |

2045 | n/a | else |

2046 | n/a | return PyFloat_FromDouble(r); |

2047 | n/a | } |

2048 | n/a | |

2049 | n/a | |

2050 | n/a | static const double degToRad = Py_MATH_PI / 180.0; |

2051 | n/a | static const double radToDeg = 180.0 / Py_MATH_PI; |

2052 | n/a | |

2053 | n/a | /*[clinic input] |

2054 | n/a | math.degrees |

2055 | n/a | |

2056 | n/a | x: double |

2057 | n/a | / |

2058 | n/a | |

2059 | n/a | Convert angle x from radians to degrees. |

2060 | n/a | [clinic start generated code]*/ |

2061 | n/a | |

2062 | n/a | static PyObject * |

2063 | n/a | math_degrees_impl(PyObject *module, double x) |

2064 | n/a | /*[clinic end generated code: output=7fea78b294acd12f input=81e016555d6e3660]*/ |

2065 | n/a | { |

2066 | n/a | return PyFloat_FromDouble(x * radToDeg); |

2067 | n/a | } |

2068 | n/a | |

2069 | n/a | |

2070 | n/a | /*[clinic input] |

2071 | n/a | math.radians |

2072 | n/a | |

2073 | n/a | x: double |

2074 | n/a | / |

2075 | n/a | |

2076 | n/a | Convert angle x from degrees to radians. |

2077 | n/a | [clinic start generated code]*/ |

2078 | n/a | |

2079 | n/a | static PyObject * |

2080 | n/a | math_radians_impl(PyObject *module, double x) |

2081 | n/a | /*[clinic end generated code: output=34daa47caf9b1590 input=91626fc489fe3d63]*/ |

2082 | n/a | { |

2083 | n/a | return PyFloat_FromDouble(x * degToRad); |

2084 | n/a | } |

2085 | n/a | |

2086 | n/a | |

2087 | n/a | /*[clinic input] |

2088 | n/a | math.isfinite |

2089 | n/a | |

2090 | n/a | x: double |

2091 | n/a | / |

2092 | n/a | |

2093 | n/a | Return True if x is neither an infinity nor a NaN, and False otherwise. |

2094 | n/a | [clinic start generated code]*/ |

2095 | n/a | |

2096 | n/a | static PyObject * |

2097 | n/a | math_isfinite_impl(PyObject *module, double x) |

2098 | n/a | /*[clinic end generated code: output=8ba1f396440c9901 input=46967d254812e54a]*/ |

2099 | n/a | { |

2100 | n/a | return PyBool_FromLong((long)Py_IS_FINITE(x)); |

2101 | n/a | } |

2102 | n/a | |

2103 | n/a | |

2104 | n/a | /*[clinic input] |

2105 | n/a | math.isnan |

2106 | n/a | |

2107 | n/a | x: double |

2108 | n/a | / |

2109 | n/a | |

2110 | n/a | Return True if x is a NaN (not a number), and False otherwise. |

2111 | n/a | [clinic start generated code]*/ |

2112 | n/a | |

2113 | n/a | static PyObject * |

2114 | n/a | math_isnan_impl(PyObject *module, double x) |

2115 | n/a | /*[clinic end generated code: output=f537b4d6df878c3e input=935891e66083f46a]*/ |

2116 | n/a | { |

2117 | n/a | return PyBool_FromLong((long)Py_IS_NAN(x)); |

2118 | n/a | } |

2119 | n/a | |

2120 | n/a | |

2121 | n/a | /*[clinic input] |

2122 | n/a | math.isinf |

2123 | n/a | |

2124 | n/a | x: double |

2125 | n/a | / |

2126 | n/a | |

2127 | n/a | Return True if x is a positive or negative infinity, and False otherwise. |

2128 | n/a | [clinic start generated code]*/ |

2129 | n/a | |

2130 | n/a | static PyObject * |

2131 | n/a | math_isinf_impl(PyObject *module, double x) |

2132 | n/a | /*[clinic end generated code: output=9f00cbec4de7b06b input=32630e4212cf961f]*/ |

2133 | n/a | { |

2134 | n/a | return PyBool_FromLong((long)Py_IS_INFINITY(x)); |

2135 | n/a | } |

2136 | n/a | |

2137 | n/a | |

2138 | n/a | /*[clinic input] |

2139 | n/a | math.isclose -> bool |

2140 | n/a | |

2141 | n/a | a: double |

2142 | n/a | b: double |

2143 | n/a | * |

2144 | n/a | rel_tol: double = 1e-09 |

2145 | n/a | maximum difference for being considered "close", relative to the |

2146 | n/a | magnitude of the input values |

2147 | n/a | abs_tol: double = 0.0 |

2148 | n/a | maximum difference for being considered "close", regardless of the |

2149 | n/a | magnitude of the input values |

2150 | n/a | |

2151 | n/a | Determine whether two floating point numbers are close in value. |

2152 | n/a | |

2153 | n/a | Return True if a is close in value to b, and False otherwise. |

2154 | n/a | |

2155 | n/a | For the values to be considered close, the difference between them |

2156 | n/a | must be smaller than at least one of the tolerances. |

2157 | n/a | |

2158 | n/a | -inf, inf and NaN behave similarly to the IEEE 754 Standard. That |

2159 | n/a | is, NaN is not close to anything, even itself. inf and -inf are |

2160 | n/a | only close to themselves. |

2161 | n/a | [clinic start generated code]*/ |

2162 | n/a | |

2163 | n/a | static int |

2164 | n/a | math_isclose_impl(PyObject *module, double a, double b, double rel_tol, |

2165 | n/a | double abs_tol) |

2166 | n/a | /*[clinic end generated code: output=b73070207511952d input=f28671871ea5bfba]*/ |

2167 | n/a | { |

2168 | n/a | double diff = 0.0; |

2169 | n/a | |

2170 | n/a | /* sanity check on the inputs */ |

2171 | n/a | if (rel_tol < 0.0 || abs_tol < 0.0 ) { |

2172 | n/a | PyErr_SetString(PyExc_ValueError, |

2173 | n/a | "tolerances must be non-negative"); |

2174 | n/a | return -1; |

2175 | n/a | } |

2176 | n/a | |

2177 | n/a | if ( a == b ) { |

2178 | n/a | /* short circuit exact equality -- needed to catch two infinities of |

2179 | n/a | the same sign. And perhaps speeds things up a bit sometimes. |

2180 | n/a | */ |

2181 | n/a | return 1; |

2182 | n/a | } |

2183 | n/a | |

2184 | n/a | /* This catches the case of two infinities of opposite sign, or |

2185 | n/a | one infinity and one finite number. Two infinities of opposite |

2186 | n/a | sign would otherwise have an infinite relative tolerance. |

2187 | n/a | Two infinities of the same sign are caught by the equality check |

2188 | n/a | above. |

2189 | n/a | */ |

2190 | n/a | |

2191 | n/a | if (Py_IS_INFINITY(a) || Py_IS_INFINITY(b)) { |

2192 | n/a | return 0; |

2193 | n/a | } |

2194 | n/a | |

2195 | n/a | /* now do the regular computation |

2196 | n/a | this is essentially the "weak" test from the Boost library |

2197 | n/a | */ |

2198 | n/a | |

2199 | n/a | diff = fabs(b - a); |

2200 | n/a | |

2201 | n/a | return (((diff <= fabs(rel_tol * b)) || |

2202 | n/a | (diff <= fabs(rel_tol * a))) || |

2203 | n/a | (diff <= abs_tol)); |

2204 | n/a | } |

2205 | n/a | |

2206 | n/a | |

2207 | n/a | static PyMethodDef math_methods[] = { |

2208 | n/a | {"acos", math_acos, METH_O, math_acos_doc}, |

2209 | n/a | {"acosh", math_acosh, METH_O, math_acosh_doc}, |

2210 | n/a | {"asin", math_asin, METH_O, math_asin_doc}, |

2211 | n/a | {"asinh", math_asinh, METH_O, math_asinh_doc}, |

2212 | n/a | {"atan", math_atan, METH_O, math_atan_doc}, |

2213 | n/a | {"atan2", math_atan2, METH_VARARGS, math_atan2_doc}, |

2214 | n/a | {"atanh", math_atanh, METH_O, math_atanh_doc}, |

2215 | n/a | MATH_CEIL_METHODDEF |

2216 | n/a | {"copysign", math_copysign, METH_VARARGS, math_copysign_doc}, |

2217 | n/a | {"cos", math_cos, METH_O, math_cos_doc}, |

2218 | n/a | {"cosh", math_cosh, METH_O, math_cosh_doc}, |

2219 | n/a | MATH_DEGREES_METHODDEF |

2220 | n/a | {"erf", math_erf, METH_O, math_erf_doc}, |

2221 | n/a | {"erfc", math_erfc, METH_O, math_erfc_doc}, |

2222 | n/a | {"exp", math_exp, METH_O, math_exp_doc}, |

2223 | n/a | {"expm1", math_expm1, METH_O, math_expm1_doc}, |

2224 | n/a | {"fabs", math_fabs, METH_O, math_fabs_doc}, |

2225 | n/a | MATH_FACTORIAL_METHODDEF |

2226 | n/a | MATH_FLOOR_METHODDEF |

2227 | n/a | MATH_FMOD_METHODDEF |

2228 | n/a | MATH_FREXP_METHODDEF |

2229 | n/a | MATH_FSUM_METHODDEF |

2230 | n/a | {"gamma", math_gamma, METH_O, math_gamma_doc}, |

2231 | n/a | MATH_GCD_METHODDEF |

2232 | n/a | MATH_HYPOT_METHODDEF |

2233 | n/a | MATH_ISCLOSE_METHODDEF |

2234 | n/a | MATH_ISFINITE_METHODDEF |

2235 | n/a | MATH_ISINF_METHODDEF |

2236 | n/a | MATH_ISNAN_METHODDEF |

2237 | n/a | MATH_LDEXP_METHODDEF |

2238 | n/a | {"lgamma", math_lgamma, METH_O, math_lgamma_doc}, |

2239 | n/a | MATH_LOG_METHODDEF |

2240 | n/a | {"log1p", math_log1p, METH_O, math_log1p_doc}, |

2241 | n/a | MATH_LOG10_METHODDEF |

2242 | n/a | MATH_LOG2_METHODDEF |

2243 | n/a | MATH_MODF_METHODDEF |

2244 | n/a | MATH_POW_METHODDEF |

2245 | n/a | MATH_RADIANS_METHODDEF |

2246 | n/a | {"sin", math_sin, METH_O, math_sin_doc}, |

2247 | n/a | {"sinh", math_sinh, METH_O, math_sinh_doc}, |

2248 | n/a | {"sqrt", math_sqrt, METH_O, math_sqrt_doc}, |

2249 | n/a | {"tan", math_tan, METH_O, math_tan_doc}, |

2250 | n/a | {"tanh", math_tanh, METH_O, math_tanh_doc}, |

2251 | n/a | MATH_TRUNC_METHODDEF |

2252 | n/a | {NULL, NULL} /* sentinel */ |

2253 | n/a | }; |

2254 | n/a | |

2255 | n/a | |

2256 | n/a | PyDoc_STRVAR(module_doc, |

2257 | n/a | "This module is always available. It provides access to the\n" |

2258 | n/a | "mathematical functions defined by the C standard."); |

2259 | n/a | |

2260 | n/a | |

2261 | n/a | static struct PyModuleDef mathmodule = { |

2262 | n/a | PyModuleDef_HEAD_INIT, |

2263 | n/a | "math", |

2264 | n/a | module_doc, |

2265 | n/a | -1, |

2266 | n/a | math_methods, |

2267 | n/a | NULL, |

2268 | n/a | NULL, |

2269 | n/a | NULL, |

2270 | n/a | NULL |

2271 | n/a | }; |

2272 | n/a | |

2273 | n/a | PyMODINIT_FUNC |

2274 | n/a | PyInit_math(void) |

2275 | n/a | { |

2276 | n/a | PyObject *m; |

2277 | n/a | |

2278 | n/a | m = PyModule_Create(&mathmodule); |

2279 | n/a | if (m == NULL) |

2280 | n/a | goto finally; |

2281 | n/a | |

2282 | n/a | PyModule_AddObject(m, "pi", PyFloat_FromDouble(Py_MATH_PI)); |

2283 | n/a | PyModule_AddObject(m, "e", PyFloat_FromDouble(Py_MATH_E)); |

2284 | n/a | PyModule_AddObject(m, "tau", PyFloat_FromDouble(Py_MATH_TAU)); /* 2pi */ |

2285 | n/a | PyModule_AddObject(m, "inf", PyFloat_FromDouble(m_inf())); |

2286 | n/a | #if !defined(PY_NO_SHORT_FLOAT_REPR) || defined(Py_NAN) |

2287 | n/a | PyModule_AddObject(m, "nan", PyFloat_FromDouble(m_nan())); |

2288 | n/a | #endif |

2289 | n/a | |

2290 | n/a | finally: |

2291 | n/a | return m; |

2292 | n/a | } |