Given integers x and M, M > 0, such that x/M is small in absolute value, compute an integer approximation to M*exp(x/M). For 0 <= x/M <= 2.4, the absolute error in the result is bounded by 60 (and is usually much smaller).
(x, M, L=8)
| 5899 | _log10_digits = _Log10Memoize().getdigits |
| 5900 | |
| 5901 | def _iexp(x, M, L=8): |
| 5902 | """Given integers x and M, M > 0, such that x/M is small in absolute |
| 5903 | value, compute an integer approximation to M*exp(x/M). For 0 <= |
| 5904 | x/M <= 2.4, the absolute error in the result is bounded by 60 (and |
| 5905 | is usually much smaller).""" |
| 5906 | |
| 5907 | # Algorithm: to compute exp(z) for a real number z, first divide z |
| 5908 | # by a suitable power R of 2 so that |z/2**R| < 2**-L. Then |
| 5909 | # compute expm1(z/2**R) = exp(z/2**R) - 1 using the usual Taylor |
| 5910 | # series |
| 5911 | # |
| 5912 | # expm1(x) = x + x**2/2! + x**3/3! + ... |
| 5913 | # |
| 5914 | # Now use the identity |
| 5915 | # |
| 5916 | # expm1(2x) = expm1(x)*(expm1(x)+2) |
| 5917 | # |
| 5918 | # R times to compute the sequence expm1(z/2**R), |
| 5919 | # expm1(z/2**(R-1)), ... , exp(z/2), exp(z). |
| 5920 | |
| 5921 | # Find R such that x/2**R/M <= 2**-L |
| 5922 | R = _nbits((x<<L)//M) |
| 5923 | |
| 5924 | # Taylor series. (2**L)**T > M |
| 5925 | T = -int(-10*len(str(M))//(3*L)) |
| 5926 | y = _div_nearest(x, T) |
| 5927 | Mshift = M<<R |
| 5928 | for i in range(T-1, 0, -1): |
| 5929 | y = _div_nearest(x*(Mshift + y), Mshift * i) |
| 5930 | |
| 5931 | # Expansion |
| 5932 | for k in range(R-1, -1, -1): |
| 5933 | Mshift = M<<(k+2) |
| 5934 | y = _div_nearest(y*(y+Mshift), Mshift) |
| 5935 | |
| 5936 | return M+y |
| 5937 | |
| 5938 | def _dexp(c, e, p): |
| 5939 | """Compute an approximation to exp(c*10**e), with p decimal places of |
no test coverage detected