sympycore.arithmetic.mpmath.lib.functions

Source Code for Module sympycore.arithmetic.mpmath.lib.functions

1 """Provides real valued functons. 2 3 Transcendental functions for real numbers: 4 * exp 5 * log 6 * sin/cos/tan 7 * sinh/cosh/tanh 8 9 """ 10 11 from util import * 12 from floatop import * 13 from squareroot import * 14 from constants import * 15 #from convert import * 16 17 18 """ 19 The exponential function has a rapidly convergent Maclaurin series: 20 21 exp(x) = 1 + x + x**2/2! + x**3/3! + x**4/4! + ... 22 23 The series can be summed very easily using fixed-point arithmetic. 24 The convergence can be improved further, using a trick due to 25 Richard P. Brent: instead of computing exp(x) directly, we choose a 26 small integer r (say, r=10) and compute exp(x/2**r)**(2**r). 27 28 The optimal value for r depends on the Python platform, the magnitude 29 of x and the target precision, and has to be estimated from 30 experimental timings. One test with x ~= 0.3 showed that 31 r = 2.2*prec**0.42 gave a good fit to the optimal values for r for 32 prec between 1 and 10000 bits, on one particular machine. 33 34 This optimization makes the summation about twice as fast at 35 low precision levels and much faster at high precision 36 (roughly five times faster at 1000 decimal digits). 37 38 If |x| is very large, we first rewrite it as t + n*log(2) with the 39 integer n chosen such that |t| <= log(2), and then calculate 40 exp(x) as exp(t)*(2**n), using the Maclaurin series for exp(t) 41 (the multiplication by 2**n just amounts to shifting the exponent). 42 """ 43

44 -def exp_series(x, prec):

45 r = int(2.2 * prec ** 0.42) 46 # XXX: more careful calculation of guard bits 47 guards = r + 3 48 if prec > 60: 49 guards += int(math.log(prec)) 50 prec2 = prec + guards 51 x = rshift_quick(x, r - guards) 52 s = (1 << prec2) + x 53 a = x 54 k = 2 55 # Sum exp(x/2**r) 56 while 1: 57 a = ((a*x) >> prec2) // k 58 if not a: break 59 s += a 60 k += 1 61 # Calculate s**(2**r) by repeated squaring 62 for j in range(r): 63 s = (s*s) >> prec2 64 return s >> guards

65

66 -def fexp(x, prec, rounding):

67 man, exp, bc = x 68 if bc == -1: 69 if x == fninf: 70 return fzero 71 return x 72 # extra precision needs to be similar in magnitude to log_2(|x|) 73 prec2 = prec + 6 + max(0, bc+exp) 74 t = make_fixed(x, prec2) 75 # abs(x) > 1? 76 if exp+bc > 1: #fcmp(fabs(x), fone) > 0: 77 lg2 = log2_fixed(prec2) 78 n, t = divmod(t, lg2) 79 else: 80 n = 0 81 return from_man_exp(exp_series(t, prec2), -prec2+n, prec, rounding)

82 83 84 """ 85 The basic strategy for computing log(x) is to set r = log(x) and use 86 Newton's method to solve the equation exp(r) = x. We set the initial 87 value r_0 to math.log(x) and then iterate r_{n+1} = r_n + exp(-r_n) - 1 88 until convergence. As with square roots, we increase the working 89 precision dynamically during the process so that only one full-precision 90 evaluation of exp is required. 91 92 log(x) is small for most inputs, so the r values can safely be 93 computed using fixed-point arithmetic. However, when x has a very 94 large or small exponent, we can improve performance through the 95 normalization log(t * 2**n) = log(t) + n*log(2), choosing n such 96 that 0.5 <= t <= 1 (for example). 97 98 There are some caveats: if x is extremely close to 1, the working 99 precision must be increased to maintain high relative precision in the 100 output (alternatively, the series approximation for log(1+x) could 101 be used in that case). 102 """ 103 104 # This function performs the Newton iteration using fixed-point 105 # arithmetic. x is assumed to have magnitude ~= 1

106 -def _log_newton(x, prec):

107 extra = 8 108 # 50-bit approximation 109 #r = int(_clog(Float((x, -prec), 64)) * 2.0**50) 110 fx = math.log(to_float((x, -prec, bitcount(x)))) 111 r = int(fx * 2.0**50) 112 prevp = 50 113 for p in giant_steps(50, prec+extra): 114 rb = lshift_quick(r, p-prevp) 115 e = exp_series(-rb, p) 116 r = rb + ((rshift_quick(x, prec-p)*e)>>p) - (1 << p) 117 prevp = p 118 return r >> extra

119

120 -def flog(x, prec, rounding):

121 if x == fzero: return fnan 122 if x == fone: return fzero 123 man, exp, bc = x 124 if bc == -1: 125 if x == finf: return finf 126 return fnan 127 if man < 0: return fnan 128 # Estimated precision needed for log(t) + n*log(2) 129 prec2 = prec + int(math.log(1+abs(bc+exp), 2)) + 10 130 # Watch out for the case when x is very close to 1 131 if -1 < bc + exp < 2: 132 near_one = fabs(fsub(x, fone, STANDARD_PREC, round_floor), STANDARD_PREC, round_floor) 133 if near_one == 0: 134 return fzero 135 # estimate how close 136 prec2 += -(near_one[1]) - bitcount(near_one[0]) 137 # Separate mantissa and exponent, calculate, join parts 138 t = rshift_quick(man, bc-prec2) 139 l = _log_newton(t, prec2) 140 a = (exp + bc) * log2_fixed(prec2) 141 return from_man_exp(l+a, -prec2, prec, rounding)

142 143 144 145 """ 146 We compute sin(x) around 0 from its Taylor series, and cos(x) around 0 147 from sqrt(1-sin(x)**2). This way we can simultaneously compute sin and 148 cos, which are often needed together (e.g. for the tangent function or 149 the complex exponential), with little extra cost compared to computing 150 just one of them. The main reason for computing sin first (and not sin 151 from cos) is to obtain high relative accuracy for x extremely close to 152 0, where the operation sqrt(1-cos(x)**2) can cause huge cancellations. 153 154 For any value of x, we can reduce it to the interval A = [-pi/4, pi/4] 155 (where the Taylor series converges quickly) by translations, changing 156 signs, and switching the roles of cos and sin: 157 158 A : sin(x) = sin(x) cos(x) = cos(x) 159 B : sin(x) = cos(x-pi/2) cos(x) = -sin(x-pi/2) 160 C : sin(x) = -sin(x-pi) cos(x) = -cos(x-pi) 161 D : sin(x) = -cos(x-3*pi/2) cos(x) = sin(x-3*pi/2) 162 163 | A | B | C | D | 164 v v v v v 165 166 1 | ____ .......... ____ 167 | _.. .. __ 168 | . __ . __ 169 | .. _ .. _ 170 | . __ . __ 171 -----| -.----------_-----------.-------------_----------- 172 | . _ .. _ . 173 | __ . __ . 174 | _ .. _ .. 175 | __ . __ . 176 | __ _.. .. 177 -1 | _________ .......... 178 0 pi 2*pi 179 180 181 TODO: could use cos series too when extremely close to 0 182 """ 183

184 -def _sin_series(x, prec):

185 x2 = (x*x) >> prec 186 s = a = x 187 k = 3 188 while a: 189 a = ((a * x2) >> prec) // (-k*(k-1)) 190 s += a 191 k += 2 192 return s

193

194 -def _trig_reduce(x, prec):

195 pi_ = pi_fixed(prec) 196 pi4 = pi_ >> 2 197 pi2 = pi_ >> 1 198 n, rem = divmod(x + pi4, pi2) 199 rem -= pi4 200 return n, rem

201

202 -def cos_sin(x, prec, rounding):

203 """Simultaneously compute (cos(x), sin(x)) for real x.""" 204 man, exp, bc = x 205 if bc == -1: 206 return (fnan, fnan) 207 bits_from_unit = abs(bc + exp) 208 prec2 = prec + bits_from_unit + 15 209 xf = make_fixed(x, prec2) 210 n, rx = _trig_reduce(xf, prec2) 211 case = n % 4 212 one = 1 << prec2 213 if case == 0: 214 s = _sin_series(rx, prec2) 215 c = sqrt_fixed(one - ((s*s)>>prec2), prec2) 216 elif case == 1: 217 c = -_sin_series(rx, prec2) 218 s = sqrt_fixed(one - ((c*c)>>prec2), prec2) 219 elif case == 2: 220 s = -_sin_series(rx, prec2) 221 c = -sqrt_fixed(one - ((s*s)>>prec2), prec2) 222 elif case == 3: 223 c = _sin_series(rx, prec2) 224 s = -sqrt_fixed(one - ((c*c)>>prec2), prec2) 225 c = from_man_exp(c, -prec2, prec, rounding) 226 s = from_man_exp(s, -prec2, prec, rounding) 227 return c, s

228

229 -def fcos(x, prec, rounding):

230 return cos_sin(x, prec, rounding)[0]

231

232 -def fsin(x, prec, rounding):

233 return cos_sin(x, prec, rounding)[1]

234

235 -def ftan(x, prec, rounding):

236 c, s = cos_sin(x, prec+6, round_floor) 237 return fdiv(s, c, prec, rounding)

238 239 240 #---------------------------------------------------------------------- 241 # Hyperbolic functions 242 # 243

244 -def _sinh_series(x, prec):

245 x2 = (x*x) >> prec 246 s = a = x 247 k = 3 248 while a: 249 a = ((a * x2) >> prec) // (k*(k-1)) 250 s += a 251 k += 2 252 return s

253

254 -def cosh_sinh(x, prec, rounding):

255 """Simultaneously compute (cosh(x), sinh(x)) for real x""" 256 man, exp, bc = x 257 if bc == -1: 258 if x == finf: return (finf, finf) 259 if x == fninf: return (finf, fninf) 260 return fnan 261 262 high_bit = exp + bc 263 prec2 = prec + 6 264 265 if high_bit < -3: 266 # Extremely close to 0, sinh(x) ~= x and cosh(x) ~= 1 267 # TODO: support directed rounding 268 if high_bit < -prec-2: 269 return (fone, fpos(x, prec, rounding)) 270 271 # Avoid cancellation when computing sinh 272 # TODO: might be faster to use sinh series directly 273 prec2 += (-high_bit) + 4 274 275 # In the general case, we use 276 # cosh(x) = (exp(x) + exp(-x))/2 277 # sinh(x) = (exp(x) - exp(-x))/2 278 # and note that the exponential only needs to be computed once. 279 ep = fexp(x, prec2, round_floor) 280 em = fdiv(fone, ep, prec2, round_floor) 281 ch = fshift_exact(fadd(ep, em, prec, rounding), -1) 282 sh = fshift_exact(fsub(ep, em, prec, rounding), -1) 283 return ch, sh

284

285 -def fcosh(x, prec, rounding):

286 """Compute cosh(x) for a real argument x""" 287 return cosh_sinh(x, prec, rounding)[0]

288

289 -def fsinh(x, prec, rounding):

290 """Compute sinh(x) for a real argument x""" 291 return cosh_sinh(x, prec, rounding)[1]

292

293 -def ftanh(x, prec, rounding):

294 """Compute tanh(x) for a real argument x""" 295 ch, sh = cosh_sinh(x, prec+6, round_floor) 296 return fdiv(sh, ch, prec, rounding)

297 298 299 #---------------------------------------------------------------------- 300 # Inverse tangent 301 # 302 303 """ 304 Near x = 0, use atan(x) = x - x**3/3 + x**5/5 - ... 305 Near x = 1, use atan(x) = y/x * (1 + 2/3*y + 2*4/3/5*y**2 + ...) 306 where y = x**2/(1+x**2). 307 308 TODO: these series are not impressively fast. It is probably better 309 to calculate atan from tan, using Newton's method or even the 310 secant method. 311 """ 312

313 -def _atan_series_1(x, prec, rounding):

314 man, exp, bc = x 315 # Increase absolute precision when extremely close to 0 316 bc = bitcount(man) 317 diff = -(bc + exp) 318 prec2 = prec 319 if diff > 10: 320 if 3*diff - 4 > prec: # x**3 term vanishes; atan(x) ~x 321 return from_man_exp(man, exp, prec, rounding) 322 prec2 = prec + diff 323 prec2 += 15 # XXX: better estimate for number of guard bits 324 x = make_fixed(x, prec2) 325 x2 = (x*x)>>prec2; one = 1<<prec2; s=a=x 326 for n in xrange(1, 1000000): 327 a = (a*x2) >> prec2 328 s += a // ((-1)**n * (n+n+1)) 329 if -100 < a < 100: 330 break 331 return from_man_exp(s, -prec2, prec, rounding)

332

333 -def _atan_series_2(x, prec, rounding):

334 prec2 = prec + 15 335 x = make_fixed(x, prec2) 336 one = 1<<prec2; x2 = (x*x)>>prec2; y=(x2<<prec2)//(one+x2) 337 s = a = one 338 for n in xrange(1, 1000000): 339 a = ((a*y)>>prec2) * (2*n) // (2*n+1) 340 if a < 100: 341 break 342 s += a 343 return from_man_exp(y*s//x, -prec2, prec, rounding)

344 345 _cutoff_1 = (5, -3, 3) # ~0.6 346 _cutoff_2 = (3, -1, 2) # 1.5 347

348 -def fatan(x, prec, rounding):

349 man, exp, bc = x 350 if bc == -1: 351 if x == finf: return fshift_exact(fpi(prec, round_down), -1) 352 if x == fninf: return fneg(fshift_exact(fpi(prec, round_down), -1)) 353 return fnan 354 if man < 0: 355 t = fatan(fneg(x), prec+4, round_floor) 356 return from_man_exp(-t[0], t[1], prec, rounding) 357 if fcmp(x, _cutoff_1) < 0: 358 return _atan_series_1(x, prec, rounding) 359 if fcmp(x, _cutoff_2) < 0: 360 return _atan_series_2(x, prec, rounding) 361 # For large x, use atan(x) = pi/2 - atan(1/x) 362 if x[1] > 10*prec: 363 pi = fpi(prec, rounding) 364 pihalf = pi[0], pi[1]-1, pi[2] 365 else: 366 pi = fpi(prec+4, round_floor) 367 pihalf = pi[0], pi[1]-1, pi[2] 368 t = fatan(fdiv(fone, x, prec+4, round_floor), prec+4, round_floor) 369 return fsub(pihalf, t, prec, rounding)

370