1""" 2This module implements computation of elementary transcendental 3functions (powers, logarithms, trigonometric and hyperbolic 4functions, inverse trigonometric and hyperbolic) for real 5floating-point numbers. 6 7For complex and interval implementations of the same functions, 8see libmpc and libmpi. 9 10""" 11 12import math 13from bisect import bisect 14 15from .backend import xrange 16from .backend import MPZ, MPZ_ZERO, MPZ_ONE, MPZ_TWO, MPZ_FIVE, BACKEND 17 18from .libmpf import ( 19 round_floor, round_ceiling, round_down, round_up, 20 round_nearest, round_fast, 21 ComplexResult, 22 bitcount, bctable, lshift, rshift, giant_steps, sqrt_fixed, 23 from_int, to_int, from_man_exp, to_fixed, to_float, from_float, 24 from_rational, normalize, 25 fzero, fone, fnone, fhalf, finf, fninf, fnan, 26 mpf_cmp, mpf_sign, mpf_abs, 27 mpf_pos, mpf_neg, mpf_add, mpf_sub, mpf_mul, mpf_div, mpf_shift, 28 mpf_rdiv_int, mpf_pow_int, mpf_sqrt, 29 reciprocal_rnd, negative_rnd, mpf_perturb, 30 isqrt_fast 31) 32 33from .libintmath import ifib 34 35 36#------------------------------------------------------------------------------- 37# Tuning parameters 38#------------------------------------------------------------------------------- 39 40# Cutoff for computing exp from cosh+sinh. This reduces the 41# number of terms by half, but also requires a square root which 42# is expensive with the pure-Python square root code. 43if BACKEND == 'python': 44 EXP_COSH_CUTOFF = 600 45else: 46 EXP_COSH_CUTOFF = 400 47# Cutoff for using more than 2 series 48EXP_SERIES_U_CUTOFF = 1500 49 50# Also basically determined by sqrt 51if BACKEND == 'python': 52 COS_SIN_CACHE_PREC = 400 53else: 54 COS_SIN_CACHE_PREC = 200 55COS_SIN_CACHE_STEP = 8 56cos_sin_cache = {} 57 58# Number of integer logarithms to cache (for zeta sums) 59MAX_LOG_INT_CACHE = 2000 60log_int_cache = {} 61 62LOG_TAYLOR_PREC = 2500 # Use Taylor series with caching up to this prec 63LOG_TAYLOR_SHIFT = 9 # Cache log values in steps of size 2^-N 64log_taylor_cache = {} 65# prec/size ratio of x for fastest convergence in AGM formula 66LOG_AGM_MAG_PREC_RATIO = 20 67 68ATAN_TAYLOR_PREC = 3000 # Same as for log 69ATAN_TAYLOR_SHIFT = 7 # steps of size 2^-N 70atan_taylor_cache = {} 71 72 73# ~= next power of two + 20 74cache_prec_steps = [22,22] 75for k in xrange(1, bitcount(LOG_TAYLOR_PREC)+1): 76 cache_prec_steps += [min(2**k,LOG_TAYLOR_PREC)+20] * 2**(k-1) 77 78 79#----------------------------------------------------------------------------# 80# # 81# Elementary mathematical constants # 82# # 83#----------------------------------------------------------------------------# 84 85def constant_memo(f): 86 """ 87 Decorator for caching computed values of mathematical 88 constants. This decorator should be applied to a 89 function taking a single argument prec as input and 90 returning a fixed-point value with the given precision. 91 """ 92 f.memo_prec = -1 93 f.memo_val = None 94 def g(prec, **kwargs): 95 memo_prec = f.memo_prec 96 if prec <= memo_prec: 97 return f.memo_val >> (memo_prec-prec) 98 newprec = int(prec*1.05+10) 99 f.memo_val = f(newprec, **kwargs) 100 f.memo_prec = newprec 101 return f.memo_val >> (newprec-prec) 102 g.__name__ = f.__name__ 103 g.__doc__ = f.__doc__ 104 return g 105 106def def_mpf_constant(fixed): 107 """ 108 Create a function that computes the mpf value for a mathematical 109 constant, given a function that computes the fixed-point value. 110 111 Assumptions: the constant is positive and has magnitude ~= 1; 112 the fixed-point function rounds to floor. 113 """ 114 def f(prec, rnd=round_fast): 115 wp = prec + 20 116 v = fixed(wp) 117 if rnd in (round_up, round_ceiling): 118 v += 1 119 return normalize(0, v, -wp, bitcount(v), prec, rnd) 120 f.__doc__ = fixed.__doc__ 121 return f 122 123def bsp_acot(q, a, b, hyperbolic): 124 if b - a == 1: 125 a1 = MPZ(2*a + 3) 126 if hyperbolic or a&1: 127 return MPZ_ONE, a1 * q**2, a1 128 else: 129 return -MPZ_ONE, a1 * q**2, a1 130 m = (a+b)//2 131 p1, q1, r1 = bsp_acot(q, a, m, hyperbolic) 132 p2, q2, r2 = bsp_acot(q, m, b, hyperbolic) 133 return q2*p1 + r1*p2, q1*q2, r1*r2 134 135# the acoth(x) series converges like the geometric series for x^2 136# N = ceil(p*log(2)/(2*log(x))) 137def acot_fixed(a, prec, hyperbolic): 138 """ 139 Compute acot(a) or acoth(a) for an integer a with binary splitting; see 140 http://numbers.computation.free.fr/Constants/Algorithms/splitting.html 141 """ 142 N = int(0.35 * prec/math.log(a) + 20) 143 p, q, r = bsp_acot(a, 0,N, hyperbolic) 144 return ((p+q)<<prec)//(q*a) 145 146def machin(coefs, prec, hyperbolic=False): 147 """ 148 Evaluate a Machin-like formula, i.e., a linear combination of 149 acot(n) or acoth(n) for specific integer values of n, using fixed- 150 point arithmetic. The input should be a list [(c, n), ...], giving 151 c*acot[h](n) + ... 152 """ 153 extraprec = 10 154 s = MPZ_ZERO 155 for a, b in coefs: 156 s += MPZ(a) * acot_fixed(MPZ(b), prec+extraprec, hyperbolic) 157 return (s >> extraprec) 158 159# Logarithms of integers are needed for various computations involving 160# logarithms, powers, radix conversion, etc 161 162@constant_memo 163def ln2_fixed(prec): 164 """ 165 Computes ln(2). This is done with a hyperbolic Machin-type formula, 166 with binary splitting at high precision. 167 """ 168 return machin([(18, 26), (-2, 4801), (8, 8749)], prec, True) 169 170@constant_memo 171def ln10_fixed(prec): 172 """ 173 Computes ln(10). This is done with a hyperbolic Machin-type formula. 174 """ 175 return machin([(46, 31), (34, 49), (20, 161)], prec, True) 176 177 178r""" 179For computation of pi, we use the Chudnovsky series: 180 181 oo 182 ___ k 183 1 \ (-1) (6 k)! (A + B k) 184 ----- = ) ----------------------- 185 12 pi /___ 3 3k+3/2 186 (3 k)! (k!) C 187 k = 0 188 189where A, B, and C are certain integer constants. This series adds roughly 19014 digits per term. Note that C^(3/2) can be extracted so that the 191series contains only rational terms. This makes binary splitting very 192efficient. 193 194The recurrence formulas for the binary splitting were taken from 195ftp://ftp.gmplib.org/pub/src/gmp-chudnovsky.c 196 197Previously, Machin's formula was used at low precision and the AGM iteration 198was used at high precision. However, the Chudnovsky series is essentially as 199fast as the Machin formula at low precision and in practice about 3x faster 200than the AGM at high precision (despite theoretically having a worse 201asymptotic complexity), so there is no reason not to use it in all cases. 202 203""" 204 205# Constants in Chudnovsky's series 206CHUD_A = MPZ(13591409) 207CHUD_B = MPZ(545140134) 208CHUD_C = MPZ(640320) 209CHUD_D = MPZ(12) 210 211def bs_chudnovsky(a, b, level, verbose): 212 """ 213 Computes the sum from a to b of the series in the Chudnovsky 214 formula. Returns g, p, q where p/q is the sum as an exact 215 fraction and g is a temporary value used to save work 216 for recursive calls. 217 """ 218 if b-a == 1: 219 g = MPZ((6*b-5)*(2*b-1)*(6*b-1)) 220 p = b**3 * CHUD_C**3 // 24 221 q = (-1)**b * g * (CHUD_A+CHUD_B*b) 222 else: 223 if verbose and level < 4: 224 print(" binary splitting", a, b) 225 mid = (a+b)//2 226 g1, p1, q1 = bs_chudnovsky(a, mid, level+1, verbose) 227 g2, p2, q2 = bs_chudnovsky(mid, b, level+1, verbose) 228 p = p1*p2 229 g = g1*g2 230 q = q1*p2 + q2*g1 231 return g, p, q 232 233@constant_memo 234def pi_fixed(prec, verbose=False, verbose_base=None): 235 """ 236 Compute floor(pi * 2**prec) as a big integer. 237 238 This is done using Chudnovsky's series (see comments in 239 libelefun.py for details). 240 """ 241 # The Chudnovsky series gives 14.18 digits per term 242 N = int(prec/3.3219280948/14.181647462 + 2) 243 if verbose: 244 print("binary splitting with N =", N) 245 g, p, q = bs_chudnovsky(0, N, 0, verbose) 246 sqrtC = isqrt_fast(CHUD_C<<(2*prec)) 247 v = p*CHUD_C*sqrtC//((q+CHUD_A*p)*CHUD_D) 248 return v 249 250def degree_fixed(prec): 251 return pi_fixed(prec)//180 252 253def bspe(a, b): 254 """ 255 Sum series for exp(1)-1 between a, b, returning the result 256 as an exact fraction (p, q). 257 """ 258 if b-a == 1: 259 return MPZ_ONE, MPZ(b) 260 m = (a+b)//2 261 p1, q1 = bspe(a, m) 262 p2, q2 = bspe(m, b) 263 return p1*q2+p2, q1*q2 264 265@constant_memo 266def e_fixed(prec): 267 """ 268 Computes exp(1). This is done using the ordinary Taylor series for 269 exp, with binary splitting. For a description of the algorithm, 270 see: 271 272 http://numbers.computation.free.fr/Constants/ 273 Algorithms/splitting.html 274 """ 275 # Slight overestimate of N needed for 1/N! < 2**(-prec) 276 # This could be tightened for large N. 277 N = int(1.1*prec/math.log(prec) + 20) 278 p, q = bspe(0,N) 279 return ((p+q)<<prec)//q 280 281@constant_memo 282def phi_fixed(prec): 283 """ 284 Computes the golden ratio, (1+sqrt(5))/2 285 """ 286 prec += 10 287 a = isqrt_fast(MPZ_FIVE<<(2*prec)) + (MPZ_ONE << prec) 288 return a >> 11 289 290mpf_phi = def_mpf_constant(phi_fixed) 291mpf_pi = def_mpf_constant(pi_fixed) 292mpf_e = def_mpf_constant(e_fixed) 293mpf_degree = def_mpf_constant(degree_fixed) 294mpf_ln2 = def_mpf_constant(ln2_fixed) 295mpf_ln10 = def_mpf_constant(ln10_fixed) 296 297 298@constant_memo 299def ln_sqrt2pi_fixed(prec): 300 wp = prec + 10 301 # ln(sqrt(2*pi)) = ln(2*pi)/2 302 return to_fixed(mpf_log(mpf_shift(mpf_pi(wp), 1), wp), prec-1) 303 304@constant_memo 305def sqrtpi_fixed(prec): 306 return sqrt_fixed(pi_fixed(prec), prec) 307 308mpf_sqrtpi = def_mpf_constant(sqrtpi_fixed) 309mpf_ln_sqrt2pi = def_mpf_constant(ln_sqrt2pi_fixed) 310 311 312#----------------------------------------------------------------------------# 313# # 314# Powers # 315# # 316#----------------------------------------------------------------------------# 317 318def mpf_pow(s, t, prec, rnd=round_fast): 319 """ 320 Compute s**t. Raises ComplexResult if s is negative and t is 321 fractional. 322 """ 323 ssign, sman, sexp, sbc = s 324 tsign, tman, texp, tbc = t 325 if ssign and texp < 0: 326 raise ComplexResult("negative number raised to a fractional power") 327 if texp >= 0: 328 return mpf_pow_int(s, (-1)**tsign * (tman<<texp), prec, rnd) 329 # s**(n/2) = sqrt(s)**n 330 if texp == -1: 331 if tman == 1: 332 if tsign: 333 return mpf_div(fone, mpf_sqrt(s, prec+10, 334 reciprocal_rnd[rnd]), prec, rnd) 335 return mpf_sqrt(s, prec, rnd) 336 else: 337 if tsign: 338 return mpf_pow_int(mpf_sqrt(s, prec+10, 339 reciprocal_rnd[rnd]), -tman, prec, rnd) 340 return mpf_pow_int(mpf_sqrt(s, prec+10, rnd), tman, prec, rnd) 341 # General formula: s**t = exp(t*log(s)) 342 # TODO: handle rnd direction of the logarithm carefully 343 c = mpf_log(s, prec+10, rnd) 344 return mpf_exp(mpf_mul(t, c), prec, rnd) 345 346def int_pow_fixed(y, n, prec): 347 """n-th power of a fixed point number with precision prec 348 349 Returns the power in the form man, exp, 350 man * 2**exp ~= y**n 351 """ 352 if n == 2: 353 return (y*y), 0 354 bc = bitcount(y) 355 exp = 0 356 workprec = 2 * (prec + 4*bitcount(n) + 4) 357 _, pm, pe, pbc = fone 358 while 1: 359 if n & 1: 360 pm = pm*y 361 pe = pe+exp 362 pbc += bc - 2 363 pbc = pbc + bctable[int(pm >> pbc)] 364 if pbc > workprec: 365 pm = pm >> (pbc-workprec) 366 pe += pbc - workprec 367 pbc = workprec 368 n -= 1 369 if not n: 370 break 371 y = y*y 372 exp = exp+exp 373 bc = bc + bc - 2 374 bc = bc + bctable[int(y >> bc)] 375 if bc > workprec: 376 y = y >> (bc-workprec) 377 exp += bc - workprec 378 bc = workprec 379 n = n // 2 380 return pm, pe 381 382# froot(s, n, prec, rnd) computes the real n-th root of a 383# positive mpf tuple s. 384# To compute the root we start from a 50-bit estimate for r 385# generated with ordinary floating-point arithmetic, and then refine 386# the value to full accuracy using the iteration 387 388# 1 / y \ 389# r = --- | (n-1) * r + ---------- | 390# n+1 n \ n r_n**(n-1) / 391 392# which is simply Newton's method applied to the equation r**n = y. 393# With giant_steps(start, prec+extra) = [p0,...,pm, prec+extra] 394# and y = man * 2**-shift one has 395# (man * 2**exp)**(1/n) = 396# y**(1/n) * 2**(start-prec/n) * 2**(p0-start) * ... * 2**(prec+extra-pm) * 397# 2**((exp+shift-(n-1)*prec)/n -extra)) 398# The last factor is accounted for in the last line of froot. 399 400def nthroot_fixed(y, n, prec, exp1): 401 start = 50 402 try: 403 y1 = rshift(y, prec - n*start) 404 r = MPZ(int(y1**(1.0/n))) 405 except OverflowError: 406 y1 = from_int(y1, start) 407 fn = from_int(n) 408 fn = mpf_rdiv_int(1, fn, start) 409 r = mpf_pow(y1, fn, start) 410 r = to_int(r) 411 extra = 10 412 extra1 = n 413 prevp = start 414 for p in giant_steps(start, prec+extra): 415 pm, pe = int_pow_fixed(r, n-1, prevp) 416 r2 = rshift(pm, (n-1)*prevp - p - pe - extra1) 417 B = lshift(y, 2*p-prec+extra1)//r2 418 r = (B + (n-1) * lshift(r, p-prevp))//n 419 prevp = p 420 return r 421 422def mpf_nthroot(s, n, prec, rnd=round_fast): 423 """nth-root of a positive number 424 425 Use the Newton method when faster, otherwise use x**(1/n) 426 """ 427 sign, man, exp, bc = s 428 if sign: 429 raise ComplexResult("nth root of a negative number") 430 if not man: 431 if s == fnan: 432 return fnan 433 if s == fzero: 434 if n > 0: 435 return fzero 436 if n == 0: 437 return fone 438 return finf 439 # Infinity 440 if not n: 441 return fnan 442 if n < 0: 443 return fzero 444 return finf 445 flag_inverse = False 446 if n < 2: 447 if n == 0: 448 return fone 449 if n == 1: 450 return mpf_pos(s, prec, rnd) 451 if n == -1: 452 return mpf_div(fone, s, prec, rnd) 453 # n < 0 454 rnd = reciprocal_rnd[rnd] 455 flag_inverse = True 456 extra_inverse = 5 457 prec += extra_inverse 458 n = -n 459 if n > 20 and (n >= 20000 or prec < int(233 + 28.3 * n**0.62)): 460 prec2 = prec + 10 461 fn = from_int(n) 462 nth = mpf_rdiv_int(1, fn, prec2) 463 r = mpf_pow(s, nth, prec2, rnd) 464 s = normalize(r[0], r[1], r[2], r[3], prec, rnd) 465 if flag_inverse: 466 return mpf_div(fone, s, prec-extra_inverse, rnd) 467 else: 468 return s 469 # Convert to a fixed-point number with prec2 bits. 470 prec2 = prec + 2*n - (prec%n) 471 # a few tests indicate that 472 # for 10 < n < 10**4 a bit more precision is needed 473 if n > 10: 474 prec2 += prec2//10 475 prec2 = prec2 - prec2%n 476 # Mantissa may have more bits than we need. Trim it down. 477 shift = bc - prec2 478 # Adjust exponents to make prec2 and exp+shift multiples of n. 479 sign1 = 0 480 es = exp+shift 481 if es < 0: 482 sign1 = 1 483 es = -es 484 if sign1: 485 shift += es%n 486 else: 487 shift -= es%n 488 man = rshift(man, shift) 489 extra = 10 490 exp1 = ((exp+shift-(n-1)*prec2)//n) - extra 491 rnd_shift = 0 492 if flag_inverse: 493 if rnd == 'u' or rnd == 'c': 494 rnd_shift = 1 495 else: 496 if rnd == 'd' or rnd == 'f': 497 rnd_shift = 1 498 man = nthroot_fixed(man+rnd_shift, n, prec2, exp1) 499 s = from_man_exp(man, exp1, prec, rnd) 500 if flag_inverse: 501 return mpf_div(fone, s, prec-extra_inverse, rnd) 502 else: 503 return s 504 505def mpf_cbrt(s, prec, rnd=round_fast): 506 """cubic root of a positive number""" 507 return mpf_nthroot(s, 3, prec, rnd) 508 509#----------------------------------------------------------------------------# 510# # 511# Logarithms # 512# # 513#----------------------------------------------------------------------------# 514 515 516def log_int_fixed(n, prec, ln2=None): 517 """ 518 Fast computation of log(n), caching the value for small n, 519 intended for zeta sums. 520 """ 521 if n in log_int_cache: 522 value, vprec = log_int_cache[n] 523 if vprec >= prec: 524 return value >> (vprec - prec) 525 wp = prec + 10 526 if wp <= LOG_TAYLOR_SHIFT: 527 if ln2 is None: 528 ln2 = ln2_fixed(wp) 529 r = bitcount(n) 530 x = n << (wp-r) 531 v = log_taylor_cached(x, wp) + r*ln2 532 else: 533 v = to_fixed(mpf_log(from_int(n), wp+5), wp) 534 if n < MAX_LOG_INT_CACHE: 535 log_int_cache[n] = (v, wp) 536 return v >> (wp-prec) 537 538def agm_fixed(a, b, prec): 539 """ 540 Fixed-point computation of agm(a,b), assuming 541 a, b both close to unit magnitude. 542 """ 543 i = 0 544 while 1: 545 anew = (a+b)>>1 546 if i > 4 and abs(a-anew) < 8: 547 return a 548 b = isqrt_fast(a*b) 549 a = anew 550 i += 1 551 return a 552 553def log_agm(x, prec): 554 """ 555 Fixed-point computation of -log(x) = log(1/x), suitable 556 for large precision. It is required that 0 < x < 1. The 557 algorithm used is the Sasaki-Kanada formula 558 559 -log(x) = pi/agm(theta2(x)^2,theta3(x)^2). [1] 560 561 For faster convergence in the theta functions, x should 562 be chosen closer to 0. 563 564 Guard bits must be added by the caller. 565 566 HYPOTHESIS: if x = 2^(-n), n bits need to be added to 567 account for the truncation to a fixed-point number, 568 and this is the only significant cancellation error. 569 570 The number of bits lost to roundoff is small and can be 571 considered constant. 572 573 [1] Richard P. Brent, "Fast Algorithms for High-Precision 574 Computation of Elementary Functions (extended abstract)", 575 http://wwwmaths.anu.edu.au/~brent/pd/RNC7-Brent.pdf 576 577 """ 578 x2 = (x*x) >> prec 579 # Compute jtheta2(x)**2 580 s = a = b = x2 581 while a: 582 b = (b*x2) >> prec 583 a = (a*b) >> prec 584 s += a 585 s += (MPZ_ONE<<prec) 586 s = (s*s)>>(prec-2) 587 s = (s*isqrt_fast(x<<prec))>>prec 588 # Compute jtheta3(x)**2 589 t = a = b = x 590 while a: 591 b = (b*x2) >> prec 592 a = (a*b) >> prec 593 t += a 594 t = (MPZ_ONE<<prec) + (t<<1) 595 t = (t*t)>>prec 596 # Final formula 597 p = agm_fixed(s, t, prec) 598 return (pi_fixed(prec) << prec) // p 599 600def log_taylor(x, prec, r=0): 601 """ 602 Fixed-point calculation of log(x). It is assumed that x is close 603 enough to 1 for the Taylor series to converge quickly. Convergence 604 can be improved by specifying r > 0 to compute 605 log(x^(1/2^r))*2^r, at the cost of performing r square roots. 606 607 The caller must provide sufficient guard bits. 608 """ 609 for i in xrange(r): 610 x = isqrt_fast(x<<prec) 611 one = MPZ_ONE << prec 612 v = ((x-one)<<prec)//(x+one) 613 sign = v < 0 614 if sign: 615 v = -v 616 v2 = (v*v) >> prec 617 v4 = (v2*v2) >> prec 618 s0 = v 619 s1 = v//3 620 v = (v*v4) >> prec 621 k = 5 622 while v: 623 s0 += v // k 624 k += 2 625 s1 += v // k 626 v = (v*v4) >> prec 627 k += 2 628 s1 = (s1*v2) >> prec 629 s = (s0+s1) << (1+r) 630 if sign: 631 return -s 632 return s 633 634def log_taylor_cached(x, prec): 635 """ 636 Fixed-point computation of log(x), assuming x in (0.5, 2) 637 and prec <= LOG_TAYLOR_PREC. 638 """ 639 n = x >> (prec-LOG_TAYLOR_SHIFT) 640 cached_prec = cache_prec_steps[prec] 641 dprec = cached_prec - prec 642 if (n, cached_prec) in log_taylor_cache: 643 a, log_a = log_taylor_cache[n, cached_prec] 644 else: 645 a = n << (cached_prec - LOG_TAYLOR_SHIFT) 646 log_a = log_taylor(a, cached_prec, 8) 647 log_taylor_cache[n, cached_prec] = (a, log_a) 648 a >>= dprec 649 log_a >>= dprec 650 u = ((x - a) << prec) // a 651 v = (u << prec) // ((MPZ_TWO << prec) + u) 652 v2 = (v*v) >> prec 653 v4 = (v2*v2) >> prec 654 s0 = v 655 s1 = v//3 656 v = (v*v4) >> prec 657 k = 5 658 while v: 659 s0 += v//k 660 k += 2 661 s1 += v//k 662 v = (v*v4) >> prec 663 k += 2 664 s1 = (s1*v2) >> prec 665 s = (s0+s1) << 1 666 return log_a + s 667 668def mpf_log(x, prec, rnd=round_fast): 669 """ 670 Compute the natural logarithm of the mpf value x. If x is negative, 671 ComplexResult is raised. 672 """ 673 sign, man, exp, bc = x 674 #------------------------------------------------------------------ 675 # Handle special values 676 if not man: 677 if x == fzero: return fninf 678 if x == finf: return finf 679 if x == fnan: return fnan 680 if sign: 681 raise ComplexResult("logarithm of a negative number") 682 wp = prec + 20 683 #------------------------------------------------------------------ 684 # Handle log(2^n) = log(n)*2. 685 # Here we catch the only possible exact value, log(1) = 0 686 if man == 1: 687 if not exp: 688 return fzero 689 return from_man_exp(exp*ln2_fixed(wp), -wp, prec, rnd) 690 mag = exp+bc 691 abs_mag = abs(mag) 692 #------------------------------------------------------------------ 693 # Handle x = 1+eps, where log(x) ~ x. We need to check for 694 # cancellation when moving to fixed-point math and compensate 695 # by increasing the precision. Note that abs_mag in (0, 1) <=> 696 # 0.5 < x < 2 and x != 1 697 if abs_mag <= 1: 698 # Calculate t = x-1 to measure distance from 1 in bits 699 tsign = 1-abs_mag 700 if tsign: 701 tman = (MPZ_ONE<<bc) - man 702 else: 703 tman = man - (MPZ_ONE<<(bc-1)) 704 tbc = bitcount(tman) 705 cancellation = bc - tbc 706 if cancellation > wp: 707 t = normalize(tsign, tman, abs_mag-bc, tbc, tbc, 'n') 708 return mpf_perturb(t, tsign, prec, rnd) 709 else: 710 wp += cancellation 711 # TODO: if close enough to 1, we could use Taylor series 712 # even in the AGM precision range, since the Taylor series 713 # converges rapidly 714 #------------------------------------------------------------------ 715 # Another special case: 716 # n*log(2) is a good enough approximation 717 if abs_mag > 10000: 718 if bitcount(abs_mag) > wp: 719 return from_man_exp(exp*ln2_fixed(wp), -wp, prec, rnd) 720 #------------------------------------------------------------------ 721 # General case. 722 # Perform argument reduction using log(x) = log(x*2^n) - n*log(2): 723 # If we are in the Taylor precision range, choose magnitude 0 or 1. 724 # If we are in the AGM precision range, choose magnitude -m for 725 # some large m; benchmarking on one machine showed m = prec/20 to be 726 # optimal between 1000 and 100,000 digits. 727 if wp <= LOG_TAYLOR_PREC: 728 m = log_taylor_cached(lshift(man, wp-bc), wp) 729 if mag: 730 m += mag*ln2_fixed(wp) 731 else: 732 optimal_mag = -wp//LOG_AGM_MAG_PREC_RATIO 733 n = optimal_mag - mag 734 x = mpf_shift(x, n) 735 wp += (-optimal_mag) 736 m = -log_agm(to_fixed(x, wp), wp) 737 m -= n*ln2_fixed(wp) 738 return from_man_exp(m, -wp, prec, rnd) 739 740def mpf_log_hypot(a, b, prec, rnd): 741 """ 742 Computes log(sqrt(a^2+b^2)) accurately. 743 """ 744 # If either a or b is inf/nan/0, assume it to be a 745 if not b[1]: 746 a, b = b, a 747 # a is inf/nan/0 748 if not a[1]: 749 # both are inf/nan/0 750 if not b[1]: 751 if a == b == fzero: 752 return fninf 753 if fnan in (a, b): 754 return fnan 755 # at least one term is (+/- inf)^2 756 return finf 757 # only a is inf/nan/0 758 if a == fzero: 759 # log(sqrt(0+b^2)) = log(|b|) 760 return mpf_log(mpf_abs(b), prec, rnd) 761 if a == fnan: 762 return fnan 763 return finf 764 # Exact 765 a2 = mpf_mul(a,a) 766 b2 = mpf_mul(b,b) 767 extra = 20 768 # Not exact 769 h2 = mpf_add(a2, b2, prec+extra) 770 cancelled = mpf_add(h2, fnone, 10) 771 mag_cancelled = cancelled[2]+cancelled[3] 772 # Just redo the sum exactly if necessary (could be smarter 773 # and avoid memory allocation when a or b is precisely 1 774 # and the other is tiny...) 775 if cancelled == fzero or mag_cancelled < -extra//2: 776 h2 = mpf_add(a2, b2, prec+extra-min(a2[2],b2[2])) 777 return mpf_shift(mpf_log(h2, prec, rnd), -1) 778 779 780#---------------------------------------------------------------------- 781# Inverse tangent 782# 783 784def atan_newton(x, prec): 785 if prec >= 100: 786 r = math.atan(int((x>>(prec-53)))/2.0**53) 787 else: 788 r = math.atan(int(x)/2.0**prec) 789 prevp = 50 790 r = MPZ(int(r * 2.0**53) >> (53-prevp)) 791 extra_p = 50 792 for wp in giant_steps(prevp, prec): 793 wp += extra_p 794 r = r << (wp-prevp) 795 cos, sin = cos_sin_fixed(r, wp) 796 tan = (sin << wp) // cos 797 a = ((tan-rshift(x, prec-wp)) << wp) // ((MPZ_ONE<<wp) + ((tan**2)>>wp)) 798 r = r - a 799 prevp = wp 800 return rshift(r, prevp-prec) 801 802def atan_taylor_get_cached(n, prec): 803 # Taylor series with caching wins up to huge precisions 804 # To avoid unnecessary precomputation at low precision, we 805 # do it in steps 806 # Round to next power of 2 807 prec2 = (1<<(bitcount(prec-1))) + 20 808 dprec = prec2 - prec 809 if (n, prec2) in atan_taylor_cache: 810 a, atan_a = atan_taylor_cache[n, prec2] 811 else: 812 a = n << (prec2 - ATAN_TAYLOR_SHIFT) 813 atan_a = atan_newton(a, prec2) 814 atan_taylor_cache[n, prec2] = (a, atan_a) 815 return (a >> dprec), (atan_a >> dprec) 816 817def atan_taylor(x, prec): 818 n = (x >> (prec-ATAN_TAYLOR_SHIFT)) 819 a, atan_a = atan_taylor_get_cached(n, prec) 820 d = x - a 821 s0 = v = (d << prec) // ((a**2 >> prec) + (a*d >> prec) + (MPZ_ONE << prec)) 822 v2 = (v**2 >> prec) 823 v4 = (v2 * v2) >> prec 824 s1 = v//3 825 v = (v * v4) >> prec 826 k = 5 827 while v: 828 s0 += v // k 829 k += 2 830 s1 += v // k 831 v = (v * v4) >> prec 832 k += 2 833 s1 = (s1 * v2) >> prec 834 s = s0 - s1 835 return atan_a + s 836 837def atan_inf(sign, prec, rnd): 838 if not sign: 839 return mpf_shift(mpf_pi(prec, rnd), -1) 840 return mpf_neg(mpf_shift(mpf_pi(prec, negative_rnd[rnd]), -1)) 841 842def mpf_atan(x, prec, rnd=round_fast): 843 sign, man, exp, bc = x 844 if not man: 845 if x == fzero: return fzero 846 if x == finf: return atan_inf(0, prec, rnd) 847 if x == fninf: return atan_inf(1, prec, rnd) 848 return fnan 849 mag = exp + bc 850 # Essentially infinity 851 if mag > prec+20: 852 return atan_inf(sign, prec, rnd) 853 # Essentially ~ x 854 if -mag > prec+20: 855 return mpf_perturb(x, 1-sign, prec, rnd) 856 wp = prec + 30 + abs(mag) 857 # For large x, use atan(x) = pi/2 - atan(1/x) 858 if mag >= 2: 859 x = mpf_rdiv_int(1, x, wp) 860 reciprocal = True 861 else: 862 reciprocal = False 863 t = to_fixed(x, wp) 864 if sign: 865 t = -t 866 if wp < ATAN_TAYLOR_PREC: 867 a = atan_taylor(t, wp) 868 else: 869 a = atan_newton(t, wp) 870 if reciprocal: 871 a = ((pi_fixed(wp)>>1)+1) - a 872 if sign: 873 a = -a 874 return from_man_exp(a, -wp, prec, rnd) 875 876# TODO: cleanup the special cases 877def mpf_atan2(y, x, prec, rnd=round_fast): 878 xsign, xman, xexp, xbc = x 879 ysign, yman, yexp, ybc = y 880 if not yman: 881 if y == fzero and x != fnan: 882 if mpf_sign(x) >= 0: 883 return fzero 884 return mpf_pi(prec, rnd) 885 if y in (finf, fninf): 886 if x in (finf, fninf): 887 return fnan 888 # pi/2 889 if y == finf: 890 return mpf_shift(mpf_pi(prec, rnd), -1) 891 # -pi/2 892 return mpf_neg(mpf_shift(mpf_pi(prec, negative_rnd[rnd]), -1)) 893 return fnan 894 if ysign: 895 return mpf_neg(mpf_atan2(mpf_neg(y), x, prec, negative_rnd[rnd])) 896 if not xman: 897 if x == fnan: 898 return fnan 899 if x == finf: 900 return fzero 901 if x == fninf: 902 return mpf_pi(prec, rnd) 903 if y == fzero: 904 return fzero 905 return mpf_shift(mpf_pi(prec, rnd), -1) 906 tquo = mpf_atan(mpf_div(y, x, prec+4), prec+4) 907 if xsign: 908 return mpf_add(mpf_pi(prec+4), tquo, prec, rnd) 909 else: 910 return mpf_pos(tquo, prec, rnd) 911 912def mpf_asin(x, prec, rnd=round_fast): 913 sign, man, exp, bc = x 914 if bc+exp > 0 and x not in (fone, fnone): 915 raise ComplexResult("asin(x) is real only for -1 <= x <= 1") 916 # asin(x) = 2*atan(x/(1+sqrt(1-x**2))) 917 wp = prec + 15 918 a = mpf_mul(x, x) 919 b = mpf_add(fone, mpf_sqrt(mpf_sub(fone, a, wp), wp), wp) 920 c = mpf_div(x, b, wp) 921 return mpf_shift(mpf_atan(c, prec, rnd), 1) 922 923def mpf_acos(x, prec, rnd=round_fast): 924 # acos(x) = 2*atan(sqrt(1-x**2)/(1+x)) 925 sign, man, exp, bc = x 926 if bc + exp > 0: 927 if x not in (fone, fnone): 928 raise ComplexResult("acos(x) is real only for -1 <= x <= 1") 929 if x == fnone: 930 return mpf_pi(prec, rnd) 931 wp = prec + 15 932 a = mpf_mul(x, x) 933 b = mpf_sqrt(mpf_sub(fone, a, wp), wp) 934 c = mpf_div(b, mpf_add(fone, x, wp), wp) 935 return mpf_shift(mpf_atan(c, prec, rnd), 1) 936 937def mpf_asinh(x, prec, rnd=round_fast): 938 wp = prec + 20 939 sign, man, exp, bc = x 940 mag = exp+bc 941 if mag < -8: 942 if mag < -wp: 943 return mpf_perturb(x, 1-sign, prec, rnd) 944 wp += (-mag) 945 # asinh(x) = log(x+sqrt(x**2+1)) 946 # use reflection symmetry to avoid cancellation 947 q = mpf_sqrt(mpf_add(mpf_mul(x, x), fone, wp), wp) 948 q = mpf_add(mpf_abs(x), q, wp) 949 if sign: 950 return mpf_neg(mpf_log(q, prec, negative_rnd[rnd])) 951 else: 952 return mpf_log(q, prec, rnd) 953 954def mpf_acosh(x, prec, rnd=round_fast): 955 # acosh(x) = log(x+sqrt(x**2-1)) 956 wp = prec + 15 957 if mpf_cmp(x, fone) == -1: 958 raise ComplexResult("acosh(x) is real only for x >= 1") 959 q = mpf_sqrt(mpf_add(mpf_mul(x,x), fnone, wp), wp) 960 return mpf_log(mpf_add(x, q, wp), prec, rnd) 961 962def mpf_atanh(x, prec, rnd=round_fast): 963 # atanh(x) = log((1+x)/(1-x))/2 964 sign, man, exp, bc = x 965 if (not man) and exp: 966 if x in (fzero, fnan): 967 return x 968 raise ComplexResult("atanh(x) is real only for -1 <= x <= 1") 969 mag = bc + exp 970 if mag > 0: 971 if mag == 1 and man == 1: 972 return [finf, fninf][sign] 973 raise ComplexResult("atanh(x) is real only for -1 <= x <= 1") 974 wp = prec + 15 975 if mag < -8: 976 if mag < -wp: 977 return mpf_perturb(x, sign, prec, rnd) 978 wp += (-mag) 979 a = mpf_add(x, fone, wp) 980 b = mpf_sub(fone, x, wp) 981 return mpf_shift(mpf_log(mpf_div(a, b, wp), prec, rnd), -1) 982 983def mpf_fibonacci(x, prec, rnd=round_fast): 984 sign, man, exp, bc = x 985 if not man: 986 if x == fninf: 987 return fnan 988 return x 989 # F(2^n) ~= 2^(2^n) 990 size = abs(exp+bc) 991 if exp >= 0: 992 # Exact 993 if size < 10 or size <= bitcount(prec): 994 return from_int(ifib(to_int(x)), prec, rnd) 995 # Use the modified Binet formula 996 wp = prec + size + 20 997 a = mpf_phi(wp) 998 b = mpf_add(mpf_shift(a, 1), fnone, wp) 999 u = mpf_pow(a, x, wp) 1000 v = mpf_cos_pi(x, wp) 1001 v = mpf_div(v, u, wp) 1002 u = mpf_sub(u, v, wp) 1003 u = mpf_div(u, b, prec, rnd) 1004 return u 1005 1006 1007#------------------------------------------------------------------------------- 1008# Exponential-type functions 1009#------------------------------------------------------------------------------- 1010 1011def exponential_series(x, prec, type=0): 1012 """ 1013 Taylor series for cosh/sinh or cos/sin. 1014 1015 type = 0 -- returns exp(x) (slightly faster than cosh+sinh) 1016 type = 1 -- returns (cosh(x), sinh(x)) 1017 type = 2 -- returns (cos(x), sin(x)) 1018 """ 1019 if x < 0: 1020 x = -x 1021 sign = 1 1022 else: 1023 sign = 0 1024 r = int(0.5*prec**0.5) 1025 xmag = bitcount(x) - prec 1026 r = max(0, xmag + r) 1027 extra = 10 + 2*max(r,-xmag) 1028 wp = prec + extra 1029 x <<= (extra - r) 1030 one = MPZ_ONE << wp 1031 alt = (type == 2) 1032 if prec < EXP_SERIES_U_CUTOFF: 1033 x2 = a = (x*x) >> wp 1034 x4 = (x2*x2) >> wp 1035 s0 = s1 = MPZ_ZERO 1036 k = 2 1037 while a: 1038 a //= (k-1)*k; s0 += a; k += 2 1039 a //= (k-1)*k; s1 += a; k += 2 1040 a = (a*x4) >> wp 1041 s1 = (x2*s1) >> wp 1042 if alt: 1043 c = s1 - s0 + one 1044 else: 1045 c = s1 + s0 + one 1046 else: 1047 u = int(0.3*prec**0.35) 1048 x2 = a = (x*x) >> wp 1049 xpowers = [one, x2] 1050 for i in xrange(1, u): 1051 xpowers.append((xpowers[-1]*x2)>>wp) 1052 sums = [MPZ_ZERO] * u 1053 k = 2 1054 while a: 1055 for i in xrange(u): 1056 a //= (k-1)*k 1057 if alt and k & 2: sums[i] -= a 1058 else: sums[i] += a 1059 k += 2 1060 a = (a*xpowers[-1]) >> wp 1061 for i in xrange(1, u): 1062 sums[i] = (sums[i]*xpowers[i]) >> wp 1063 c = sum(sums) + one 1064 if type == 0: 1065 s = isqrt_fast(c*c - (one<<wp)) 1066 if sign: 1067 v = c - s 1068 else: 1069 v = c + s 1070 for i in xrange(r): 1071 v = (v*v) >> wp 1072 return v >> extra 1073 else: 1074 # Repeatedly apply the double-angle formula 1075 # cosh(2*x) = 2*cosh(x)^2 - 1 1076 # cos(2*x) = 2*cos(x)^2 - 1 1077 pshift = wp-1 1078 for i in xrange(r): 1079 c = ((c*c) >> pshift) - one 1080 # With the abs, this is the same for sinh and sin 1081 s = isqrt_fast(abs((one<<wp) - c*c)) 1082 if sign: 1083 s = -s 1084 return (c>>extra), (s>>extra) 1085 1086def exp_basecase(x, prec): 1087 """ 1088 Compute exp(x) as a fixed-point number. Works for any x, 1089 but for speed should have |x| < 1. For an arbitrary number, 1090 use exp(x) = exp(x-m*log(2)) * 2^m where m = floor(x/log(2)). 1091 """ 1092 if prec > EXP_COSH_CUTOFF: 1093 return exponential_series(x, prec, 0) 1094 r = int(prec**0.5) 1095 prec += r 1096 s0 = s1 = (MPZ_ONE << prec) 1097 k = 2 1098 a = x2 = (x*x) >> prec 1099 while a: 1100 a //= k; s0 += a; k += 1 1101 a //= k; s1 += a; k += 1 1102 a = (a*x2) >> prec 1103 s1 = (s1*x) >> prec 1104 s = s0 + s1 1105 u = r 1106 while r: 1107 s = (s*s) >> prec 1108 r -= 1 1109 return s >> u 1110 1111def exp_expneg_basecase(x, prec): 1112 """ 1113 Computation of exp(x), exp(-x) 1114 """ 1115 if prec > EXP_COSH_CUTOFF: 1116 cosh, sinh = exponential_series(x, prec, 1) 1117 return cosh+sinh, cosh-sinh 1118 a = exp_basecase(x, prec) 1119 b = (MPZ_ONE << (prec+prec)) // a 1120 return a, b 1121 1122def cos_sin_basecase(x, prec): 1123 """ 1124 Compute cos(x), sin(x) as fixed-point numbers, assuming x 1125 in [0, pi/2). For an arbitrary number, use x' = x - m*(pi/2) 1126 where m = floor(x/(pi/2)) along with quarter-period symmetries. 1127 """ 1128 if prec > COS_SIN_CACHE_PREC: 1129 return exponential_series(x, prec, 2) 1130 precs = prec - COS_SIN_CACHE_STEP 1131 t = x >> precs 1132 n = int(t) 1133 if n not in cos_sin_cache: 1134 w = t<<(10+COS_SIN_CACHE_PREC-COS_SIN_CACHE_STEP) 1135 cos_t, sin_t = exponential_series(w, 10+COS_SIN_CACHE_PREC, 2) 1136 cos_sin_cache[n] = (cos_t>>10), (sin_t>>10) 1137 cos_t, sin_t = cos_sin_cache[n] 1138 offset = COS_SIN_CACHE_PREC - prec 1139 cos_t >>= offset 1140 sin_t >>= offset 1141 x -= t << precs 1142 cos = MPZ_ONE << prec 1143 sin = x 1144 k = 2 1145 a = -((x*x) >> prec) 1146 while a: 1147 a //= k; cos += a; k += 1; a = (a*x) >> prec 1148 a //= k; sin += a; k += 1; a = -((a*x) >> prec) 1149 return ((cos*cos_t-sin*sin_t) >> prec), ((sin*cos_t+cos*sin_t) >> prec) 1150 1151def mpf_exp(x, prec, rnd=round_fast): 1152 sign, man, exp, bc = x 1153 if man: 1154 mag = bc + exp 1155 wp = prec + 14 1156 if sign: 1157 man = -man 1158 # TODO: the best cutoff depends on both x and the precision. 1159 if prec > 600 and exp >= 0: 1160 # Need about log2(exp(n)) ~= 1.45*mag extra precision 1161 e = mpf_e(wp+int(1.45*mag)) 1162 return mpf_pow_int(e, man<<exp, prec, rnd) 1163 if mag < -wp: 1164 return mpf_perturb(fone, sign, prec, rnd) 1165 # |x| >= 2 1166 if mag > 1: 1167 # For large arguments: exp(2^mag*(1+eps)) = 1168 # exp(2^mag)*exp(2^mag*eps) = exp(2^mag)*(1 + 2^mag*eps + ...) 1169 # so about mag extra bits is required. 1170 wpmod = wp + mag 1171 offset = exp + wpmod 1172 if offset >= 0: 1173 t = man << offset 1174 else: 1175 t = man >> (-offset) 1176 lg2 = ln2_fixed(wpmod) 1177 n, t = divmod(t, lg2) 1178 n = int(n) 1179 t >>= mag 1180 else: 1181 offset = exp + wp 1182 if offset >= 0: 1183 t = man << offset 1184 else: 1185 t = man >> (-offset) 1186 n = 0 1187 man = exp_basecase(t, wp) 1188 return from_man_exp(man, n-wp, prec, rnd) 1189 if not exp: 1190 return fone 1191 if x == fninf: 1192 return fzero 1193 return x 1194 1195 1196def mpf_cosh_sinh(x, prec, rnd=round_fast, tanh=0): 1197 """Simultaneously compute (cosh(x), sinh(x)) for real x""" 1198 sign, man, exp, bc = x 1199 if (not man) and exp: 1200 if tanh: 1201 if x == finf: return fone 1202 if x == fninf: return fnone 1203 return fnan 1204 if x == finf: return (finf, finf) 1205 if x == fninf: return (finf, fninf) 1206 return fnan, fnan 1207 mag = exp+bc 1208 wp = prec+14 1209 if mag < -4: 1210 # Extremely close to 0, sinh(x) ~= x and cosh(x) ~= 1 1211 if mag < -wp: 1212 if tanh: 1213 return mpf_perturb(x, 1-sign, prec, rnd) 1214 cosh = mpf_perturb(fone, 0, prec, rnd) 1215 sinh = mpf_perturb(x, sign, prec, rnd) 1216 return cosh, sinh 1217 # Fix for cancellation when computing sinh 1218 wp += (-mag) 1219 # Does exp(-2*x) vanish? 1220 if mag > 10: 1221 if 3*(1<<(mag-1)) > wp: 1222 # XXX: rounding 1223 if tanh: 1224 return mpf_perturb([fone,fnone][sign], 1-sign, prec, rnd) 1225 c = s = mpf_shift(mpf_exp(mpf_abs(x), prec, rnd), -1) 1226 if sign: 1227 s = mpf_neg(s) 1228 return c, s 1229 # |x| > 1 1230 if mag > 1: 1231 wpmod = wp + mag 1232 offset = exp + wpmod 1233 if offset >= 0: 1234 t = man << offset 1235 else: 1236 t = man >> (-offset) 1237 lg2 = ln2_fixed(wpmod) 1238 n, t = divmod(t, lg2) 1239 n = int(n) 1240 t >>= mag 1241 else: 1242 offset = exp + wp 1243 if offset >= 0: 1244 t = man << offset 1245 else: 1246 t = man >> (-offset) 1247 n = 0 1248 a, b = exp_expneg_basecase(t, wp) 1249 # TODO: optimize division precision 1250 cosh = a + (b>>(2*n)) 1251 sinh = a - (b>>(2*n)) 1252 if sign: 1253 sinh = -sinh 1254 if tanh: 1255 man = (sinh << wp) // cosh 1256 return from_man_exp(man, -wp, prec, rnd) 1257 else: 1258 cosh = from_man_exp(cosh, n-wp-1, prec, rnd) 1259 sinh = from_man_exp(sinh, n-wp-1, prec, rnd) 1260 return cosh, sinh 1261 1262 1263def mod_pi2(man, exp, mag, wp): 1264 # Reduce to standard interval 1265 if mag > 0: 1266 i = 0 1267 while 1: 1268 cancellation_prec = 20 << i 1269 wpmod = wp + mag + cancellation_prec 1270 pi2 = pi_fixed(wpmod-1) 1271 pi4 = pi2 >> 1 1272 offset = wpmod + exp 1273 if offset >= 0: 1274 t = man << offset 1275 else: 1276 t = man >> (-offset) 1277 n, y = divmod(t, pi2) 1278 if y > pi4: 1279 small = pi2 - y 1280 else: 1281 small = y 1282 if small >> (wp+mag-10): 1283 n = int(n) 1284 t = y >> mag 1285 wp = wpmod - mag 1286 break 1287 i += 1 1288 else: 1289 wp += (-mag) 1290 offset = exp + wp 1291 if offset >= 0: 1292 t = man << offset 1293 else: 1294 t = man >> (-offset) 1295 n = 0 1296 return t, n, wp 1297 1298 1299def mpf_cos_sin(x, prec, rnd=round_fast, which=0, pi=False): 1300 """ 1301 which: 1302 0 -- return cos(x), sin(x) 1303 1 -- return cos(x) 1304 2 -- return sin(x) 1305 3 -- return tan(x) 1306 1307 if pi=True, compute for pi*x 1308 """ 1309 sign, man, exp, bc = x 1310 if not man: 1311 if exp: 1312 c, s = fnan, fnan 1313 else: 1314 c, s = fone, fzero 1315 if which == 0: return c, s 1316 if which == 1: return c 1317 if which == 2: return s 1318 if which == 3: return s 1319 1320 mag = bc + exp 1321 wp = prec + 10 1322 1323 # Extremely small? 1324 if mag < 0: 1325 if mag < -wp: 1326 if pi: 1327 x = mpf_mul(x, mpf_pi(wp)) 1328 c = mpf_perturb(fone, 1, prec, rnd) 1329 s = mpf_perturb(x, 1-sign, prec, rnd) 1330 if which == 0: return c, s 1331 if which == 1: return c 1332 if which == 2: return s 1333 if which == 3: return mpf_perturb(x, sign, prec, rnd) 1334 if pi: 1335 if exp >= -1: 1336 if exp == -1: 1337 c = fzero 1338 s = (fone, fnone)[bool(man & 2) ^ sign] 1339 elif exp == 0: 1340 c, s = (fnone, fzero) 1341 else: 1342 c, s = (fone, fzero) 1343 if which == 0: return c, s 1344 if which == 1: return c 1345 if which == 2: return s 1346 if which == 3: return mpf_div(s, c, prec, rnd) 1347 # Subtract nearest half-integer (= mod by pi/2) 1348 n = ((man >> (-exp-2)) + 1) >> 1 1349 man = man - (n << (-exp-1)) 1350 mag2 = bitcount(man) + exp 1351 wp = prec + 10 - mag2 1352 offset = exp + wp 1353 if offset >= 0: 1354 t = man << offset 1355 else: 1356 t = man >> (-offset) 1357 t = (t*pi_fixed(wp)) >> wp 1358 else: 1359 t, n, wp = mod_pi2(man, exp, mag, wp) 1360 c, s = cos_sin_basecase(t, wp) 1361 m = n & 3 1362 if m == 1: c, s = -s, c 1363 elif m == 2: c, s = -c, -s 1364 elif m == 3: c, s = s, -c 1365 if sign: 1366 s = -s 1367 if which == 0: 1368 c = from_man_exp(c, -wp, prec, rnd) 1369 s = from_man_exp(s, -wp, prec, rnd) 1370 return c, s 1371 if which == 1: 1372 return from_man_exp(c, -wp, prec, rnd) 1373 if which == 2: 1374 return from_man_exp(s, -wp, prec, rnd) 1375 if which == 3: 1376 return from_rational(s, c, prec, rnd) 1377 1378def mpf_cos(x, prec, rnd=round_fast): return mpf_cos_sin(x, prec, rnd, 1) 1379def mpf_sin(x, prec, rnd=round_fast): return mpf_cos_sin(x, prec, rnd, 2) 1380def mpf_tan(x, prec, rnd=round_fast): return mpf_cos_sin(x, prec, rnd, 3) 1381def mpf_cos_sin_pi(x, prec, rnd=round_fast): return mpf_cos_sin(x, prec, rnd, 0, 1) 1382def mpf_cos_pi(x, prec, rnd=round_fast): return mpf_cos_sin(x, prec, rnd, 1, 1) 1383def mpf_sin_pi(x, prec, rnd=round_fast): return mpf_cos_sin(x, prec, rnd, 2, 1) 1384def mpf_cosh(x, prec, rnd=round_fast): return mpf_cosh_sinh(x, prec, rnd)[0] 1385def mpf_sinh(x, prec, rnd=round_fast): return mpf_cosh_sinh(x, prec, rnd)[1] 1386def mpf_tanh(x, prec, rnd=round_fast): return mpf_cosh_sinh(x, prec, rnd, tanh=1) 1387 1388 1389# Low-overhead fixed-point versions 1390 1391def cos_sin_fixed(x, prec, pi2=None): 1392 if pi2 is None: 1393 pi2 = pi_fixed(prec-1) 1394 n, t = divmod(x, pi2) 1395 n = int(n) 1396 c, s = cos_sin_basecase(t, prec) 1397 m = n & 3 1398 if m == 0: return c, s 1399 if m == 1: return -s, c 1400 if m == 2: return -c, -s 1401 if m == 3: return s, -c 1402 1403def exp_fixed(x, prec, ln2=None): 1404 if ln2 is None: 1405 ln2 = ln2_fixed(prec) 1406 n, t = divmod(x, ln2) 1407 n = int(n) 1408 v = exp_basecase(t, prec) 1409 if n >= 0: 1410 return v << n 1411 else: 1412 return v >> (-n) 1413 1414 1415if BACKEND == 'sage': 1416 try: 1417 import sage.libs.mpmath.ext_libmp as _lbmp 1418 mpf_sqrt = _lbmp.mpf_sqrt 1419 mpf_exp = _lbmp.mpf_exp 1420 mpf_log = _lbmp.mpf_log 1421 mpf_cos = _lbmp.mpf_cos 1422 mpf_sin = _lbmp.mpf_sin 1423 mpf_pow = _lbmp.mpf_pow 1424 exp_fixed = _lbmp.exp_fixed 1425 cos_sin_fixed = _lbmp.cos_sin_fixed 1426 log_int_fixed = _lbmp.log_int_fixed 1427 except (ImportError, AttributeError): 1428 print("Warning: Sage imports in libelefun failed") 1429