add read me

venv/lib/python3.12/site-packages/mpmath/libmp/__init__.py

@@ -0,0 +1,77 @@
+from .libmpf import (prec_to_dps, dps_to_prec, repr_dps,
+  round_down, round_up, round_floor, round_ceiling, round_nearest,
+  to_pickable, from_pickable, ComplexResult,
+  fzero, fnzero, fone, fnone, ftwo, ften, fhalf, fnan, finf, fninf,
+  math_float_inf, round_int, normalize, normalize1,
+  from_man_exp, from_int, to_man_exp, to_int, mpf_ceil, mpf_floor,
+  mpf_nint, mpf_frac,
+  from_float, from_npfloat, from_Decimal, to_float, from_rational, to_rational, to_fixed,
+  mpf_rand, mpf_eq, mpf_hash, mpf_cmp, mpf_lt, mpf_le, mpf_gt, mpf_ge,
+  mpf_pos, mpf_neg, mpf_abs, mpf_sign, mpf_add, mpf_sub, mpf_sum,
+  mpf_mul, mpf_mul_int, mpf_shift, mpf_frexp,
+  mpf_div, mpf_rdiv_int, mpf_mod, mpf_pow_int,
+  mpf_perturb,
+  to_digits_exp, to_str, str_to_man_exp, from_str, from_bstr, to_bstr,
+  mpf_sqrt, mpf_hypot)
+
+from .libmpc import (mpc_one, mpc_zero, mpc_two, mpc_half,
+  mpc_is_inf, mpc_is_infnan, mpc_to_str, mpc_to_complex, mpc_hash,
+  mpc_conjugate, mpc_is_nonzero, mpc_add, mpc_add_mpf,
+  mpc_sub, mpc_sub_mpf, mpc_pos, mpc_neg, mpc_shift, mpc_abs,
+  mpc_arg, mpc_floor, mpc_ceil,  mpc_nint, mpc_frac, mpc_mul, mpc_square,
+  mpc_mul_mpf, mpc_mul_imag_mpf, mpc_mul_int,
+  mpc_div, mpc_div_mpf, mpc_reciprocal, mpc_mpf_div,
+  complex_int_pow, mpc_pow, mpc_pow_mpf, mpc_pow_int,
+  mpc_sqrt, mpc_nthroot, mpc_cbrt, mpc_exp, mpc_log, mpc_cos, mpc_sin,
+  mpc_tan, mpc_cos_pi, mpc_sin_pi, mpc_cosh, mpc_sinh, mpc_tanh,
+  mpc_atan, mpc_acos, mpc_asin, mpc_asinh, mpc_acosh, mpc_atanh,
+  mpc_fibonacci, mpf_expj, mpf_expjpi, mpc_expj, mpc_expjpi,
+  mpc_cos_sin, mpc_cos_sin_pi)
+
+from .libelefun import (ln2_fixed, mpf_ln2, ln10_fixed, mpf_ln10,
+  pi_fixed, mpf_pi, e_fixed, mpf_e, phi_fixed, mpf_phi,
+  degree_fixed, mpf_degree,
+  mpf_pow, mpf_nthroot, mpf_cbrt, log_int_fixed, agm_fixed,
+  mpf_log, mpf_log_hypot, mpf_exp, mpf_cos_sin, mpf_cos, mpf_sin, mpf_tan,
+  mpf_cos_sin_pi, mpf_cos_pi, mpf_sin_pi, mpf_cosh_sinh,
+  mpf_cosh, mpf_sinh, mpf_tanh, mpf_atan, mpf_atan2, mpf_asin,
+  mpf_acos, mpf_asinh, mpf_acosh, mpf_atanh, mpf_fibonacci)
+
+from .libhyper import (NoConvergence, make_hyp_summator,
+  mpf_erf, mpf_erfc, mpf_ei, mpc_ei, mpf_e1, mpc_e1, mpf_expint,
+  mpf_ci_si, mpf_ci, mpf_si, mpc_ci, mpc_si, mpf_besseljn,
+  mpc_besseljn, mpf_agm, mpf_agm1, mpc_agm, mpc_agm1,
+  mpf_ellipk, mpc_ellipk, mpf_ellipe, mpc_ellipe)
+
+from .gammazeta import (catalan_fixed, mpf_catalan,
+  khinchin_fixed, mpf_khinchin, glaisher_fixed, mpf_glaisher,
+  apery_fixed, mpf_apery, euler_fixed, mpf_euler, mertens_fixed,
+  mpf_mertens, twinprime_fixed, mpf_twinprime,
+  mpf_bernoulli, bernfrac, mpf_gamma_int,
+  mpf_factorial, mpc_factorial, mpf_gamma, mpc_gamma,
+  mpf_loggamma, mpc_loggamma, mpf_rgamma, mpc_rgamma,
+  mpf_harmonic, mpc_harmonic, mpf_psi0, mpc_psi0,
+  mpf_psi, mpc_psi, mpf_zeta_int, mpf_zeta, mpc_zeta,
+  mpf_altzeta, mpc_altzeta, mpf_zetasum, mpc_zetasum)
+
+from .libmpi import (mpi_str,
+  mpi_from_str, mpi_to_str,
+  mpi_eq, mpi_ne,
+  mpi_lt, mpi_le, mpi_gt, mpi_ge,
+  mpi_add, mpi_sub, mpi_delta, mpi_mid,
+  mpi_pos, mpi_neg, mpi_abs, mpi_mul, mpi_div, mpi_exp,
+  mpi_log, mpi_sqrt, mpi_pow_int, mpi_pow, mpi_cos_sin,
+  mpi_cos, mpi_sin, mpi_tan, mpi_cot,
+  mpi_atan, mpi_atan2,
+  mpci_pos, mpci_neg, mpci_add, mpci_sub, mpci_mul, mpci_div, mpci_pow,
+  mpci_abs, mpci_pow, mpci_exp, mpci_log, mpci_cos, mpci_sin,
+  mpi_gamma, mpci_gamma, mpi_loggamma, mpci_loggamma,
+  mpi_rgamma, mpci_rgamma, mpi_factorial, mpci_factorial)
+
+from .libintmath import (trailing, bitcount, numeral, bin_to_radix,
+  isqrt, isqrt_small, isqrt_fast, sqrt_fixed, sqrtrem, ifib, ifac,
+  list_primes, isprime, moebius, gcd, eulernum, stirling1, stirling2)
+
+from .backend import (gmpy, sage, BACKEND, STRICT, MPZ, MPZ_TYPE,
+  MPZ_ZERO, MPZ_ONE, MPZ_TWO, MPZ_THREE, MPZ_FIVE, int_types,
+  HASH_MODULUS, HASH_BITS)

venv/lib/python3.12/site-packages/mpmath/libmp/backend.py

@@ -0,0 +1,115 @@
+import os
+import sys
+
+#----------------------------------------------------------------------------#
+# Support GMPY for high-speed large integer arithmetic.                      #
+#                                                                            #
+# To allow an external module to handle arithmetic, we need to make sure     #
+# that all high-precision variables are declared of the correct type. MPZ    #
+# is the constructor for the high-precision type. It defaults to Python's    #
+# long type but can be assinged another type, typically gmpy.mpz.            #
+#                                                                            #
+# MPZ must be used for the mantissa component of an mpf and must be used     #
+# for internal fixed-point operations.                                       #
+#                                                                            #
+# Side-effects                                                               #
+# 1) "is" cannot be used to test for special values. Must use "==".          #
+# 2) There are bugs in GMPY prior to v1.02 so we must use v1.03 or later.    #
+#----------------------------------------------------------------------------#
+
+# So we can import it from this module
+gmpy = None
+sage = None
+sage_utils = None
+
+if sys.version_info[0] < 3:
+    python3 = False
+else:
+    python3 = True
+
+BACKEND = 'python'
+
+if not python3:
+    MPZ = long
+    xrange = xrange
+    basestring = basestring
+
+    def exec_(_code_, _globs_=None, _locs_=None):
+        """Execute code in a namespace."""
+        if _globs_ is None:
+            frame = sys._getframe(1)
+            _globs_ = frame.f_globals
+            if _locs_ is None:
+                _locs_ = frame.f_locals
+            del frame
+        elif _locs_ is None:
+            _locs_ = _globs_
+        exec("""exec _code_ in _globs_, _locs_""")
+else:
+    MPZ = int
+    xrange = range
+    basestring = str
+
+    import builtins
+    exec_ = getattr(builtins, "exec")
+
+# Define constants for calculating hash on Python 3.2.
+if sys.version_info >= (3, 2):
+    HASH_MODULUS = sys.hash_info.modulus
+    if sys.hash_info.width == 32:
+        HASH_BITS = 31
+    else:
+        HASH_BITS = 61
+else:
+    HASH_MODULUS = None
+    HASH_BITS = None
+
+if 'MPMATH_NOGMPY' not in os.environ:
+    try:
+        try:
+            import gmpy2 as gmpy
+        except ImportError:
+            try:
+                import gmpy
+            except ImportError:
+                raise ImportError
+        if gmpy.version() >= '1.03':
+            BACKEND = 'gmpy'
+            MPZ = gmpy.mpz
+    except:
+        pass
+
+if ('MPMATH_NOSAGE' not in os.environ and 'SAGE_ROOT' in os.environ or
+        'MPMATH_SAGE' in os.environ):
+    try:
+        import sage.all
+        import sage.libs.mpmath.utils as _sage_utils
+        sage = sage.all
+        sage_utils = _sage_utils
+        BACKEND = 'sage'
+        MPZ = sage.Integer
+    except:
+        pass
+
+if 'MPMATH_STRICT' in os.environ:
+    STRICT = True
+else:
+    STRICT = False
+
+MPZ_TYPE = type(MPZ(0))
+MPZ_ZERO = MPZ(0)
+MPZ_ONE = MPZ(1)
+MPZ_TWO = MPZ(2)
+MPZ_THREE = MPZ(3)
+MPZ_FIVE = MPZ(5)
+
+try:
+    if BACKEND == 'python':
+        int_types = (int, long)
+    else:
+        int_types = (int, long, MPZ_TYPE)
+except NameError:
+    if BACKEND == 'python':
+        int_types = (int,)
+    else:
+        int_types = (int, MPZ_TYPE)

venv/lib/python3.12/site-packages/mpmath/libmp/libintmath.py

@@ -0,0 +1,584 @@
+"""
+Utility functions for integer math.
+
+TODO: rename, cleanup, perhaps move the gmpy wrapper code
+here from settings.py
+
+"""
+
+import math
+from bisect import bisect
+
+from .backend import xrange
+from .backend import BACKEND, gmpy, sage, sage_utils, MPZ, MPZ_ONE, MPZ_ZERO
+
+small_trailing = [0] * 256
+for j in range(1,8):
+    small_trailing[1<<j::1<<(j+1)] = [j] * (1<<(7-j))
+
+def giant_steps(start, target, n=2):
+    """
+    Return a list of integers ~=
+
+    [start, n*start, ..., target/n^2, target/n, target]
+
+    but conservatively rounded so that the quotient between two
+    successive elements is actually slightly less than n.
+
+    With n = 2, this describes suitable precision steps for a
+    quadratically convergent algorithm such as Newton's method;
+    with n = 3 steps for cubic convergence (Halley's method), etc.
+
+        >>> giant_steps(50,1000)
+        [66, 128, 253, 502, 1000]
+        >>> giant_steps(50,1000,4)
+        [65, 252, 1000]
+
+    """
+    L = [target]
+    while L[-1] > start*n:
+        L = L + [L[-1]//n + 2]
+    return L[::-1]
+
+def rshift(x, n):
+    """For an integer x, calculate x >> n with the fastest (floor)
+    rounding. Unlike the plain Python expression (x >> n), n is
+    allowed to be negative, in which case a left shift is performed."""
+    if n >= 0: return x >> n
+    else:      return x << (-n)
+
+def lshift(x, n):
+    """For an integer x, calculate x << n. Unlike the plain Python
+    expression (x << n), n is allowed to be negative, in which case a
+    right shift with default (floor) rounding is performed."""
+    if n >= 0: return x << n
+    else:      return x >> (-n)
+
+if BACKEND == 'sage':
+    import operator
+    rshift = operator.rshift
+    lshift = operator.lshift
+
+def python_trailing(n):
+    """Count the number of trailing zero bits in abs(n)."""
+    if not n:
+        return 0
+    low_byte = n & 0xff
+    if low_byte:
+        return small_trailing[low_byte]
+    t = 8
+    n >>= 8
+    while not n & 0xff:
+        n >>= 8
+        t += 8
+    return t + small_trailing[n & 0xff]
+
+if BACKEND == 'gmpy':
+    if gmpy.version() >= '2':
+        def gmpy_trailing(n):
+            """Count the number of trailing zero bits in abs(n) using gmpy."""
+            if n: return MPZ(n).bit_scan1()
+            else: return 0
+    else:
+        def gmpy_trailing(n):
+            """Count the number of trailing zero bits in abs(n) using gmpy."""
+            if n: return MPZ(n).scan1()
+            else: return 0
+
+# Small powers of 2
+powers = [1<<_ for _ in range(300)]
+
+def python_bitcount(n):
+    """Calculate bit size of the nonnegative integer n."""
+    bc = bisect(powers, n)
+    if bc != 300:
+        return bc
+    bc = int(math.log(n, 2)) - 4
+    return bc + bctable[n>>bc]
+
+def gmpy_bitcount(n):
+    """Calculate bit size of the nonnegative integer n."""
+    if n: return MPZ(n).numdigits(2)
+    else: return 0
+
+#def sage_bitcount(n):
+#    if n: return MPZ(n).nbits()
+#    else: return 0
+
+def sage_trailing(n):
+    return MPZ(n).trailing_zero_bits()
+
+if BACKEND == 'gmpy':
+    bitcount = gmpy_bitcount
+    trailing = gmpy_trailing
+elif BACKEND == 'sage':
+    sage_bitcount = sage_utils.bitcount
+    bitcount = sage_bitcount
+    trailing = sage_trailing
+else:
+    bitcount = python_bitcount
+    trailing = python_trailing
+
+if BACKEND == 'gmpy' and 'bit_length' in dir(gmpy):
+    bitcount = gmpy.bit_length
+
+# Used to avoid slow function calls as far as possible
+trailtable = [trailing(n) for n in range(256)]
+bctable = [bitcount(n) for n in range(1024)]
+
+# TODO: speed up for bases 2, 4, 8, 16, ...
+
+def bin_to_radix(x, xbits, base, bdigits):
+    """Changes radix of a fixed-point number; i.e., converts
+    x * 2**xbits to floor(x * 10**bdigits)."""
+    return x * (MPZ(base)**bdigits) >> xbits
+
+stddigits = '0123456789abcdefghijklmnopqrstuvwxyz'
+
+def small_numeral(n, base=10, digits=stddigits):
+    """Return the string numeral of a positive integer in an arbitrary
+    base. Most efficient for small input."""
+    if base == 10:
+        return str(n)
+    digs = []
+    while n:
+        n, digit = divmod(n, base)
+        digs.append(digits[digit])
+    return "".join(digs[::-1])
+
+def numeral_python(n, base=10, size=0, digits=stddigits):
+    """Represent the integer n as a string of digits in the given base.
+    Recursive division is used to make this function about 3x faster
+    than Python's str() for converting integers to decimal strings.
+
+    The 'size' parameters specifies the number of digits in n; this
+    number is only used to determine splitting points and need not be
+    exact."""
+    if n <= 0:
+        if not n:
+            return "0"
+        return "-" + numeral(-n, base, size, digits)
+    # Fast enough to do directly
+    if size < 250:
+        return small_numeral(n, base, digits)
+    # Divide in half
+    half = (size // 2) + (size & 1)
+    A, B = divmod(n, base**half)
+    ad = numeral(A, base, half, digits)
+    bd = numeral(B, base, half, digits).rjust(half, "0")
+    return ad + bd
+
+def numeral_gmpy(n, base=10, size=0, digits=stddigits):
+    """Represent the integer n as a string of digits in the given base.
+    Recursive division is used to make this function about 3x faster
+    than Python's str() for converting integers to decimal strings.
+
+    The 'size' parameters specifies the number of digits in n; this
+    number is only used to determine splitting points and need not be
+    exact."""
+    if n < 0:
+        return "-" + numeral(-n, base, size, digits)
+    # gmpy.digits() may cause a segmentation fault when trying to convert
+    # extremely large values to a string. The size limit may need to be
+    # adjusted on some platforms, but 1500000 works on Windows and Linux.
+    if size < 1500000:
+        return gmpy.digits(n, base)
+    # Divide in half
+    half = (size // 2) + (size & 1)
+    A, B = divmod(n, MPZ(base)**half)
+    ad = numeral(A, base, half, digits)
+    bd = numeral(B, base, half, digits).rjust(half, "0")
+    return ad + bd
+
+if BACKEND == "gmpy":
+    numeral = numeral_gmpy
+else:
+    numeral = numeral_python
+
+_1_800 = 1<<800
+_1_600 = 1<<600
+_1_400 = 1<<400
+_1_200 = 1<<200
+_1_100 = 1<<100
+_1_50 = 1<<50
+
+def isqrt_small_python(x):
+    """
+    Correctly (floor) rounded integer square root, using
+    division. Fast up to ~200 digits.
+    """
+    if not x:
+        return x
+    if x < _1_800:
+        # Exact with IEEE double precision arithmetic
+        if x < _1_50:
+            return int(x**0.5)
+        # Initial estimate can be any integer >= the true root; round up
+        r = int(x**0.5 * 1.00000000000001) + 1
+    else:
+        bc = bitcount(x)
+        n = bc//2
+        r = int((x>>(2*n-100))**0.5+2)<<(n-50)  # +2 is to round up
+    # The following iteration now precisely computes floor(sqrt(x))
+    # See e.g. Crandall & Pomerance, "Prime Numbers: A Computational
+    # Perspective"
+    while 1:
+        y = (r+x//r)>>1
+        if y >= r:
+            return r
+        r = y
+
+def isqrt_fast_python(x):
+    """
+    Fast approximate integer square root, computed using division-free
+    Newton iteration for large x. For random integers the result is almost
+    always correct (floor(sqrt(x))), but is 1 ulp too small with a roughly
+    0.1% probability. If x is very close to an exact square, the answer is
+    1 ulp wrong with high probability.
+
+    With 0 guard bits, the largest error over a set of 10^5 random
+    inputs of size 1-10^5 bits was 3 ulp. The use of 10 guard bits
+    almost certainly guarantees a max 1 ulp error.
+    """
+    # Use direct division-based iteration if sqrt(x) < 2^400
+    # Assume floating-point square root accurate to within 1 ulp, then:
+    # 0 Newton iterations good to 52 bits
+    # 1 Newton iterations good to 104 bits
+    # 2 Newton iterations good to 208 bits
+    # 3 Newton iterations good to 416 bits
+    if x < _1_800:
+        y = int(x**0.5)
+        if x >= _1_100:
+            y = (y + x//y) >> 1
+            if x >= _1_200:
+                y = (y + x//y) >> 1
+                if x >= _1_400:
+                    y = (y + x//y) >> 1
+        return y
+    bc = bitcount(x)
+    guard_bits = 10
+    x <<= 2*guard_bits
+    bc += 2*guard_bits
+    bc += (bc&1)
+    hbc = bc//2
+    startprec = min(50, hbc)
+    # Newton iteration for 1/sqrt(x), with floating-point starting value
+    r = int(2.0**(2*startprec) * (x >> (bc-2*startprec)) ** -0.5)
+    pp = startprec
+    for p in giant_steps(startprec, hbc):
+        # r**2, scaled from real size 2**(-bc) to 2**p
+        r2 = (r*r) >> (2*pp - p)
+        # x*r**2, scaled from real size ~1.0 to 2**p
+        xr2 = ((x >> (bc-p)) * r2) >> p
+        # New value of r, scaled from real size 2**(-bc/2) to 2**p
+        r = (r * ((3<<p) - xr2)) >> (pp+1)
+        pp = p
+    # (1/sqrt(x))*x = sqrt(x)
+    return (r*(x>>hbc)) >> (p+guard_bits)
+
+def sqrtrem_python(x):
+    """Correctly rounded integer (floor) square root with remainder."""
+    # to check cutoff:
+    # plot(lambda x: timing(isqrt, 2**int(x)), [0,2000])
+    if x < _1_600:
+        y = isqrt_small_python(x)
+        return y, x - y*y
+    y = isqrt_fast_python(x) + 1
+    rem = x - y*y
+    # Correct remainder
+    while rem < 0:
+        y -= 1
+        rem += (1+2*y)
+    else:
+        if rem:
+            while rem > 2*(1+y):
+                y += 1
+                rem -= (1+2*y)
+    return y, rem
+
+def isqrt_python(x):
+    """Integer square root with correct (floor) rounding."""
+    return sqrtrem_python(x)[0]
+
+def sqrt_fixed(x, prec):
+    return isqrt_fast(x<<prec)
+
+sqrt_fixed2 = sqrt_fixed
+
+if BACKEND == 'gmpy':
+    if gmpy.version() >= '2':
+        isqrt_small = isqrt_fast = isqrt = gmpy.isqrt
+        sqrtrem = gmpy.isqrt_rem
+    else:
+        isqrt_small = isqrt_fast = isqrt = gmpy.sqrt
+        sqrtrem = gmpy.sqrtrem
+elif BACKEND == 'sage':
+    isqrt_small = isqrt_fast = isqrt = \
+        getattr(sage_utils, "isqrt", lambda n: MPZ(n).isqrt())
+    sqrtrem = lambda n: MPZ(n).sqrtrem()
+else:
+    isqrt_small = isqrt_small_python
+    isqrt_fast = isqrt_fast_python
+    isqrt = isqrt_python
+    sqrtrem = sqrtrem_python
+
+
+def ifib(n, _cache={}):
+    """Computes the nth Fibonacci number as an integer, for
+    integer n."""
+    if n < 0:
+        return (-1)**(-n+1) * ifib(-n)
+    if n in _cache:
+        return _cache[n]
+    m = n
+    # Use Dijkstra's logarithmic algorithm
+    # The following implementation is basically equivalent to
+    # http://en.literateprograms.org/Fibonacci_numbers_(Scheme)
+    a, b, p, q = MPZ_ONE, MPZ_ZERO, MPZ_ZERO, MPZ_ONE
+    while n:
+        if n & 1:
+            aq = a*q
+            a, b = b*q+aq+a*p, b*p+aq
+            n -= 1
+        else:
+            qq = q*q
+            p, q = p*p+qq, qq+2*p*q
+            n >>= 1
+    if m < 250:
+        _cache[m] = b
+    return b
+
+MAX_FACTORIAL_CACHE = 1000
+
+def ifac(n, memo={0:1, 1:1}):
+    """Return n factorial (for integers n >= 0 only)."""
+    f = memo.get(n)
+    if f:
+        return f
+    k = len(memo)
+    p = memo[k-1]
+    MAX = MAX_FACTORIAL_CACHE
+    while k <= n:
+        p *= k
+        if k <= MAX:
+            memo[k] = p
+        k += 1
+    return p
+
+def ifac2(n, memo_pair=[{0:1}, {1:1}]):
+    """Return n!! (double factorial), integers n >= 0 only."""
+    memo = memo_pair[n&1]
+    f = memo.get(n)
+    if f:
+        return f
+    k = max(memo)
+    p = memo[k]
+    MAX = MAX_FACTORIAL_CACHE
+    while k < n:
+        k += 2
+        p *= k
+        if k <= MAX:
+            memo[k] = p
+    return p
+
+if BACKEND == 'gmpy':
+    ifac = gmpy.fac
+elif BACKEND == 'sage':
+    ifac = lambda n: int(sage.factorial(n))
+    ifib = sage.fibonacci
+
+def list_primes(n):
+    n = n + 1
+    sieve = list(xrange(n))
+    sieve[:2] = [0, 0]
+    for i in xrange(2, int(n**0.5)+1):
+        if sieve[i]:
+            for j in xrange(i**2, n, i):
+                sieve[j] = 0
+    return [p for p in sieve if p]
+
+if BACKEND == 'sage':
+    # Note: it is *VERY* important for performance that we convert
+    # the list to Python ints.
+    def list_primes(n):
+        return [int(_) for _ in sage.primes(n+1)]
+
+small_odd_primes = (3,5,7,11,13,17,19,23,29,31,37,41,43,47)
+small_odd_primes_set = set(small_odd_primes)
+
+def isprime(n):
+    """
+    Determines whether n is a prime number. A probabilistic test is
+    performed if n is very large. No special trick is used for detecting
+    perfect powers.
+
+        >>> sum(list_primes(100000))
+        454396537
+        >>> sum(n*isprime(n) for n in range(100000))
+        454396537
+
+    """
+    n = int(n)
+    if not n & 1:
+        return n == 2
+    if n < 50:
+        return n in small_odd_primes_set
+    for p in small_odd_primes:
+        if not n % p:
+            return False
+    m = n-1
+    s = trailing(m)
+    d = m >> s
+    def test(a):
+        x = pow(a,d,n)
+        if x == 1 or x == m:
+            return True
+        for r in xrange(1,s):
+            x = x**2 % n
+            if x == m:
+                return True
+        return False
+    # See http://primes.utm.edu/prove/prove2_3.html
+    if n < 1373653:
+        witnesses = [2,3]
+    elif n < 341550071728321:
+        witnesses = [2,3,5,7,11,13,17]
+    else:
+        witnesses = small_odd_primes
+    for a in witnesses:
+        if not test(a):
+            return False
+    return True
+
+def moebius(n):
+    """
+    Evaluates the Moebius function which is `mu(n) = (-1)^k` if `n`
+    is a product of `k` distinct primes and `mu(n) = 0` otherwise.
+
+    TODO: speed up using factorization
+    """
+    n = abs(int(n))
+    if n < 2:
+        return n
+    factors = []
+    for p in xrange(2, n+1):
+        if not (n % p):
+            if not (n % p**2):
+                return 0
+            if not sum(p % f for f in factors):
+                factors.append(p)
+    return (-1)**len(factors)
+
+def gcd(*args):
+    a = 0
+    for b in args:
+        if a:
+            while b:
+                a, b = b, a % b
+        else:
+            a = b
+    return a
+
+
+#  Comment by Juan Arias de Reyna:
+#
+#  I learn this method to compute EulerE[2n] from van de Lune.
+#
+#  We apply the formula   EulerE[2n] = (-1)^n 2**(-2n) sum_{j=0}^n a(2n,2j+1)
+#
+#  where the numbers a(n,j) vanish for  j > n+1 or j <= -1  and satisfies
+#
+#  a(0,-1) = a(0,0) = 0;  a(0,1)= 1; a(0,2) = a(0,3) = 0
+#
+#  a(n,j) = a(n-1,j)                              when n+j is even
+#  a(n,j) = (j-1) a(n-1,j-1) + (j+1) a(n-1,j+1)   when n+j is odd
+#
+#
+#  But we can use only one array unidimensional a(j) since to compute
+#  a(n,j) we only need to know a(n-1,k) where k and j are of different parity
+#  and we have not to conserve the used values.
+#
+#  We cached up the values of Euler numbers to sufficiently high order.
+#
+#  Important Observation: If we pretend to use the numbers
+#     EulerE[1], EulerE[2], ... , EulerE[n]
+#     it is convenient to compute first EulerE[n], since the algorithm
+#     computes first all
+#     the previous ones, and keeps them in the CACHE
+
+MAX_EULER_CACHE = 500
+
+def eulernum(m, _cache={0:MPZ_ONE}):
+    r"""
+    Computes the Euler numbers `E(n)`, which can be defined as
+    coefficients of the Taylor expansion of `1/cosh x`:
+
+    .. math ::
+
+        \frac{1}{\cosh x} = \sum_{n=0}^\infty \frac{E_n}{n!} x^n
+
+    Example::
+
+        >>> [int(eulernum(n)) for n in range(11)]
+        [1, 0, -1, 0, 5, 0, -61, 0, 1385, 0, -50521]
+        >>> [int(eulernum(n)) for n in range(11)]   # test cache
+        [1, 0, -1, 0, 5, 0, -61, 0, 1385, 0, -50521]
+
+    """
+    # for odd m > 1, the Euler numbers are zero
+    if m & 1:
+        return MPZ_ZERO
+    f = _cache.get(m)
+    if f:
+        return f
+    MAX = MAX_EULER_CACHE
+    n = m
+    a = [MPZ(_) for _ in [0,0,1,0,0,0]]
+    for  n in range(1, m+1):
+        for j in range(n+1, -1, -2):
+            a[j+1] = (j-1)*a[j] + (j+1)*a[j+2]
+        a.append(0)
+        suma = 0
+        for k in range(n+1, -1, -2):
+            suma += a[k+1]
+            if n <= MAX:
+                _cache[n] = ((-1)**(n//2))*(suma // 2**n)
+        if n == m:
+            return ((-1)**(n//2))*suma // 2**n
+
+def stirling1(n, k):
+    """
+    Stirling number of the first kind.
+    """
+    if n < 0 or k < 0:
+        raise ValueError
+    if k >= n:
+        return MPZ(n == k)
+    if k < 1:
+        return MPZ_ZERO
+    L = [MPZ_ZERO] * (k+1)
+    L[1] = MPZ_ONE
+    for m in xrange(2, n+1):
+        for j in xrange(min(k, m), 0, -1):
+            L[j] = (m-1) * L[j] + L[j-1]
+    return (-1)**(n+k) * L[k]
+
+def stirling2(n, k):
+    """
+    Stirling number of the second kind.
+    """
+    if n < 0 or k < 0:
+        raise ValueError
+    if k >= n:
+        return MPZ(n == k)
+    if k <= 1:
+        return MPZ(k == 1)
+    s = MPZ_ZERO
+    t = MPZ_ONE
+    for j in xrange(k+1):
+        if (k + j) & 1:
+            s -= t * MPZ(j)**n
+        else:
+            s += t * MPZ(j)**n
+        t = t * (k - j) // (j + 1)
+    return s // ifac(k)

venv/lib/python3.12/site-packages/mpmath/libmp/libmpc.py

@@ -0,0 +1,835 @@
+"""
+Low-level functions for complex arithmetic.
+"""
+
+import sys
+
+from .backend import MPZ, MPZ_ZERO, MPZ_ONE, MPZ_TWO, BACKEND
+
+from .libmpf import (\
+    round_floor, round_ceiling, round_down, round_up,
+    round_nearest, round_fast, bitcount,
+    bctable, normalize, normalize1, reciprocal_rnd, rshift, lshift, giant_steps,
+    negative_rnd,
+    to_str, to_fixed, from_man_exp, from_float, to_float, from_int, to_int,
+    fzero, fone, ftwo, fhalf, finf, fninf, fnan, fnone,
+    mpf_abs, mpf_pos, mpf_neg, mpf_add, mpf_sub, mpf_mul,
+    mpf_div, mpf_mul_int, mpf_shift, mpf_sqrt, mpf_hypot,
+    mpf_rdiv_int, mpf_floor, mpf_ceil, mpf_nint, mpf_frac,
+    mpf_sign, mpf_hash,
+    ComplexResult
+)
+
+from .libelefun import (\
+    mpf_pi, mpf_exp, mpf_log, mpf_cos_sin, mpf_cosh_sinh, mpf_tan, mpf_pow_int,
+    mpf_log_hypot,
+    mpf_cos_sin_pi, mpf_phi,
+    mpf_cos, mpf_sin, mpf_cos_pi, mpf_sin_pi,
+    mpf_atan, mpf_atan2, mpf_cosh, mpf_sinh, mpf_tanh,
+    mpf_asin, mpf_acos, mpf_acosh, mpf_nthroot, mpf_fibonacci
+)
+
+# An mpc value is a (real, imag) tuple
+mpc_one = fone, fzero
+mpc_zero = fzero, fzero
+mpc_two = ftwo, fzero
+mpc_half = (fhalf, fzero)
+
+_infs = (finf, fninf)
+_infs_nan = (finf, fninf, fnan)
+
+def mpc_is_inf(z):
+    """Check if either real or imaginary part is infinite"""
+    re, im = z
+    if re in _infs: return True
+    if im in _infs: return True
+    return False
+
+def mpc_is_infnan(z):
+    """Check if either real or imaginary part is infinite or nan"""
+    re, im = z
+    if re in _infs_nan: return True
+    if im in _infs_nan: return True
+    return False
+
+def mpc_to_str(z, dps, **kwargs):
+    re, im = z
+    rs = to_str(re, dps)
+    if im[0]:
+        return rs + " - " + to_str(mpf_neg(im), dps, **kwargs) + "j"
+    else:
+        return rs + " + " + to_str(im, dps, **kwargs) + "j"
+
+def mpc_to_complex(z, strict=False, rnd=round_fast):
+    re, im = z
+    return complex(to_float(re, strict, rnd), to_float(im, strict, rnd))
+
+def mpc_hash(z):
+    if sys.version_info >= (3, 2):
+        re, im = z
+        h = mpf_hash(re) + sys.hash_info.imag * mpf_hash(im)
+        # Need to reduce either module 2^32 or 2^64
+        h = h % (2**sys.hash_info.width)
+        return int(h)
+    else:
+        try:
+            return hash(mpc_to_complex(z, strict=True))
+        except OverflowError:
+            return hash(z)
+
+def mpc_conjugate(z, prec, rnd=round_fast):
+    re, im = z
+    return re, mpf_neg(im, prec, rnd)
+
+def mpc_is_nonzero(z):
+    return z != mpc_zero
+
+def mpc_add(z, w, prec, rnd=round_fast):
+    a, b = z
+    c, d = w
+    return mpf_add(a, c, prec, rnd), mpf_add(b, d, prec, rnd)
+
+def mpc_add_mpf(z, x, prec, rnd=round_fast):
+    a, b = z
+    return mpf_add(a, x, prec, rnd), b
+
+def mpc_sub(z, w, prec=0, rnd=round_fast):
+    a, b = z
+    c, d = w
+    return mpf_sub(a, c, prec, rnd), mpf_sub(b, d, prec, rnd)
+
+def mpc_sub_mpf(z, p, prec=0, rnd=round_fast):
+    a, b = z
+    return mpf_sub(a, p, prec, rnd), b
+
+def mpc_pos(z, prec, rnd=round_fast):
+    a, b = z
+    return mpf_pos(a, prec, rnd), mpf_pos(b, prec, rnd)
+
+def mpc_neg(z, prec=None, rnd=round_fast):
+    a, b = z
+    return mpf_neg(a, prec, rnd), mpf_neg(b, prec, rnd)
+
+def mpc_shift(z, n):
+    a, b = z
+    return mpf_shift(a, n), mpf_shift(b, n)
+
+def mpc_abs(z, prec, rnd=round_fast):
+    """Absolute value of a complex number, |a+bi|.
+    Returns an mpf value."""
+    a, b = z
+    return mpf_hypot(a, b, prec, rnd)
+
+def mpc_arg(z, prec, rnd=round_fast):
+    """Argument of a complex number. Returns an mpf value."""
+    a, b = z
+    return mpf_atan2(b, a, prec, rnd)
+
+def mpc_floor(z, prec, rnd=round_fast):
+    a, b = z
+    return mpf_floor(a, prec, rnd), mpf_floor(b, prec, rnd)
+
+def mpc_ceil(z, prec, rnd=round_fast):
+    a, b = z
+    return mpf_ceil(a, prec, rnd), mpf_ceil(b, prec, rnd)
+
+def mpc_nint(z, prec, rnd=round_fast):
+    a, b = z
+    return mpf_nint(a, prec, rnd), mpf_nint(b, prec, rnd)
+
+def mpc_frac(z, prec, rnd=round_fast):
+    a, b = z
+    return mpf_frac(a, prec, rnd), mpf_frac(b, prec, rnd)
+
+
+def mpc_mul(z, w, prec, rnd=round_fast):
+    """
+    Complex multiplication.
+
+    Returns the real and imaginary part of (a+bi)*(c+di), rounded to
+    the specified precision. The rounding mode applies to the real and
+    imaginary parts separately.
+    """
+    a, b = z
+    c, d = w
+    p = mpf_mul(a, c)
+    q = mpf_mul(b, d)
+    r = mpf_mul(a, d)
+    s = mpf_mul(b, c)
+    re = mpf_sub(p, q, prec, rnd)
+    im = mpf_add(r, s, prec, rnd)
+    return re, im
+
+def mpc_square(z, prec, rnd=round_fast):
+    # (a+b*I)**2 == a**2 - b**2 + 2*I*a*b
+    a, b = z
+    p = mpf_mul(a,a)
+    q = mpf_mul(b,b)
+    r = mpf_mul(a,b, prec, rnd)
+    re = mpf_sub(p, q, prec, rnd)
+    im = mpf_shift(r, 1)
+    return re, im
+
+def mpc_mul_mpf(z, p, prec, rnd=round_fast):
+    a, b = z
+    re = mpf_mul(a, p, prec, rnd)
+    im = mpf_mul(b, p, prec, rnd)
+    return re, im
+
+def mpc_mul_imag_mpf(z, x, prec, rnd=round_fast):
+    """
+    Multiply the mpc value z by I*x where x is an mpf value.
+    """
+    a, b = z
+    re = mpf_neg(mpf_mul(b, x, prec, rnd))
+    im = mpf_mul(a, x, prec, rnd)
+    return re, im
+
+def mpc_mul_int(z, n, prec, rnd=round_fast):
+    a, b = z
+    re = mpf_mul_int(a, n, prec, rnd)
+    im = mpf_mul_int(b, n, prec, rnd)
+    return re, im
+
+def mpc_div(z, w, prec, rnd=round_fast):
+    a, b = z
+    c, d = w
+    wp = prec + 10
+    # mag = c*c + d*d
+    mag = mpf_add(mpf_mul(c, c), mpf_mul(d, d), wp)
+    # (a*c+b*d)/mag, (b*c-a*d)/mag
+    t = mpf_add(mpf_mul(a,c), mpf_mul(b,d), wp)
+    u = mpf_sub(mpf_mul(b,c), mpf_mul(a,d), wp)
+    return mpf_div(t,mag,prec,rnd), mpf_div(u,mag,prec,rnd)
+
+def mpc_div_mpf(z, p, prec, rnd=round_fast):
+    """Calculate z/p where p is real"""
+    a, b = z
+    re = mpf_div(a, p, prec, rnd)
+    im = mpf_div(b, p, prec, rnd)
+    return re, im
+
+def mpc_reciprocal(z, prec, rnd=round_fast):
+    """Calculate 1/z efficiently"""
+    a, b = z
+    m = mpf_add(mpf_mul(a,a),mpf_mul(b,b),prec+10)
+    re = mpf_div(a, m, prec, rnd)
+    im = mpf_neg(mpf_div(b, m, prec, rnd))
+    return re, im
+
+def mpc_mpf_div(p, z, prec, rnd=round_fast):
+    """Calculate p/z where p is real efficiently"""
+    a, b = z
+    m = mpf_add(mpf_mul(a,a),mpf_mul(b,b), prec+10)
+    re = mpf_div(mpf_mul(a,p), m, prec, rnd)
+    im = mpf_div(mpf_neg(mpf_mul(b,p)), m, prec, rnd)
+    return re, im
+
+def complex_int_pow(a, b, n):
+    """Complex integer power: computes (a+b*I)**n exactly for
+    nonnegative n (a and b must be Python ints)."""
+    wre = 1
+    wim = 0
+    while n:
+        if n & 1:
+            wre, wim = wre*a - wim*b, wim*a + wre*b
+            n -= 1
+        a, b = a*a - b*b, 2*a*b
+        n //= 2
+    return wre, wim
+
+def mpc_pow(z, w, prec, rnd=round_fast):
+    if w[1] == fzero:
+        return mpc_pow_mpf(z, w[0], prec, rnd)
+    return mpc_exp(mpc_mul(mpc_log(z, prec+10), w, prec+10), prec, rnd)
+
+def mpc_pow_mpf(z, p, prec, rnd=round_fast):
+    psign, pman, pexp, pbc = p
+    if pexp >= 0:
+        return mpc_pow_int(z, (-1)**psign * (pman<<pexp), prec, rnd)
+    if pexp == -1:
+        sqrtz = mpc_sqrt(z, prec+10)
+        return mpc_pow_int(sqrtz, (-1)**psign * pman, prec, rnd)
+    return mpc_exp(mpc_mul_mpf(mpc_log(z, prec+10), p, prec+10), prec, rnd)
+
+def mpc_pow_int(z, n, prec, rnd=round_fast):
+    a, b = z
+    if b == fzero:
+        return mpf_pow_int(a, n, prec, rnd), fzero
+    if a == fzero:
+        v = mpf_pow_int(b, n, prec, rnd)
+        n %= 4
+        if n == 0:
+            return v, fzero
+        elif n == 1:
+            return fzero, v
+        elif n == 2:
+            return mpf_neg(v), fzero
+        elif n == 3:
+            return fzero, mpf_neg(v)
+    if n == 0: return mpc_one
+    if n == 1: return mpc_pos(z, prec, rnd)
+    if n == 2: return mpc_square(z, prec, rnd)
+    if n == -1: return mpc_reciprocal(z, prec, rnd)
+    if n < 0: return mpc_reciprocal(mpc_pow_int(z, -n, prec+4), prec, rnd)
+    asign, aman, aexp, abc = a
+    bsign, bman, bexp, bbc = b
+    if asign: aman = -aman
+    if bsign: bman = -bman
+    de = aexp - bexp
+    abs_de = abs(de)
+    exact_size = n*(abs_de + max(abc, bbc))
+    if exact_size < 10000:
+        if de > 0:
+            aman <<= de
+            aexp = bexp
+        else:
+            bman <<= (-de)
+            bexp = aexp
+        re, im = complex_int_pow(aman, bman, n)
+        re = from_man_exp(re, int(n*aexp), prec, rnd)
+        im = from_man_exp(im, int(n*bexp), prec, rnd)
+        return re, im
+    return mpc_exp(mpc_mul_int(mpc_log(z, prec+10), n, prec+10), prec, rnd)
+
+def mpc_sqrt(z, prec, rnd=round_fast):
+    """Complex square root (principal branch).
+
+    We have sqrt(a+bi) = sqrt((r+a)/2) + b/sqrt(2*(r+a))*i where
+    r = abs(a+bi), when a+bi is not a negative real number."""
+    a, b = z
+    if b == fzero:
+        if a == fzero:
+            return (a, b)
+        # When a+bi is a negative real number, we get a real sqrt times i
+        if a[0]:
+            im = mpf_sqrt(mpf_neg(a), prec, rnd)
+            return (fzero, im)
+        else:
+            re = mpf_sqrt(a, prec, rnd)
+            return (re, fzero)
+    wp = prec+20
+    if not a[0]:                               # case a positive
+        t  = mpf_add(mpc_abs((a, b), wp), a, wp)  # t = abs(a+bi) + a
+        u = mpf_shift(t, -1)                      # u = t/2
+        re = mpf_sqrt(u, prec, rnd)               # re = sqrt(u)
+        v = mpf_shift(t, 1)                       # v = 2*t
+        w  = mpf_sqrt(v, wp)                      # w = sqrt(v)
+        im = mpf_div(b, w, prec, rnd)             # im = b / w
+    else:                                      # case a negative
+        t = mpf_sub(mpc_abs((a, b), wp), a, wp)   # t = abs(a+bi) - a
+        u = mpf_shift(t, -1)                      # u = t/2
+        im = mpf_sqrt(u, prec, rnd)               # im = sqrt(u)
+        v = mpf_shift(t, 1)                       # v = 2*t
+        w  = mpf_sqrt(v, wp)                      # w = sqrt(v)
+        re = mpf_div(b, w, prec, rnd)             # re = b/w
+        if b[0]:
+            re = mpf_neg(re)
+            im = mpf_neg(im)
+    return re, im
+
+def mpc_nthroot_fixed(a, b, n, prec):
+    # a, b signed integers at fixed precision prec
+    start = 50
+    a1 = int(rshift(a, prec - n*start))
+    b1 = int(rshift(b, prec - n*start))
+    try:
+        r = (a1 + 1j * b1)**(1.0/n)
+        re = r.real
+        im = r.imag
+        re = MPZ(int(re))
+        im = MPZ(int(im))
+    except OverflowError:
+        a1 = from_int(a1, start)
+        b1 = from_int(b1, start)
+        fn = from_int(n)
+        nth = mpf_rdiv_int(1, fn, start)
+        re, im = mpc_pow((a1, b1), (nth, fzero), start)
+        re = to_int(re)
+        im = to_int(im)
+    extra = 10
+    prevp = start
+    extra1 = n
+    for p in giant_steps(start, prec+extra):
+        # this is slow for large n, unlike int_pow_fixed
+        re2, im2 = complex_int_pow(re, im, n-1)
+        re2 = rshift(re2, (n-1)*prevp - p - extra1)
+        im2 = rshift(im2, (n-1)*prevp - p - extra1)
+        r4 = (re2*re2 + im2*im2) >> (p + extra1)
+        ap = rshift(a, prec - p)
+        bp = rshift(b, prec - p)
+        rec = (ap * re2 + bp * im2) >> p
+        imc = (-ap * im2 + bp * re2) >> p
+        reb = (rec << p) // r4
+        imb = (imc << p) // r4
+        re = (reb + (n-1)*lshift(re, p-prevp))//n
+        im = (imb + (n-1)*lshift(im, p-prevp))//n
+        prevp = p
+    return re, im
+
+def mpc_nthroot(z, n, prec, rnd=round_fast):
+    """
+    Complex n-th root.
+
+    Use Newton method as in the real case when it is faster,
+    otherwise use z**(1/n)
+    """
+    a, b = z
+    if a[0] == 0 and b == fzero:
+        re = mpf_nthroot(a, n, prec, rnd)
+        return (re, fzero)
+    if n < 2:
+        if n == 0:
+            return mpc_one
+        if n == 1:
+            return mpc_pos((a, b), prec, rnd)
+        if n == -1:
+            return mpc_div(mpc_one, (a, b), prec, rnd)
+        inverse = mpc_nthroot((a, b), -n, prec+5, reciprocal_rnd[rnd])
+        return mpc_div(mpc_one, inverse, prec, rnd)
+    if n <= 20:
+        prec2 = int(1.2 * (prec + 10))
+        asign, aman, aexp, abc = a
+        bsign, bman, bexp, bbc = b
+        pf = mpc_abs((a,b), prec)
+        if pf[-2] + pf[-1] > -10  and pf[-2] + pf[-1] < prec:
+            af = to_fixed(a, prec2)
+            bf = to_fixed(b, prec2)
+            re, im = mpc_nthroot_fixed(af, bf, n, prec2)
+            extra = 10
+            re = from_man_exp(re, -prec2-extra, prec2, rnd)
+            im = from_man_exp(im, -prec2-extra, prec2, rnd)
+            return re, im
+    fn = from_int(n)
+    prec2 = prec+10 + 10
+    nth = mpf_rdiv_int(1, fn, prec2)
+    re, im = mpc_pow((a, b), (nth, fzero), prec2, rnd)
+    re = normalize(re[0], re[1], re[2], re[3], prec, rnd)
+    im = normalize(im[0], im[1], im[2], im[3], prec, rnd)
+    return re, im
+
+def mpc_cbrt(z, prec, rnd=round_fast):
+    """
+    Complex cubic root.
+    """
+    return mpc_nthroot(z, 3, prec, rnd)
+
+def mpc_exp(z, prec, rnd=round_fast):
+    """
+    Complex exponential function.
+
+    We use the direct formula exp(a+bi) = exp(a) * (cos(b) + sin(b)*i)
+    for the computation. This formula is very nice because it is
+    pefectly stable; since we just do real multiplications, the only
+    numerical errors that can creep in are single-ulp rounding errors.
+
+    The formula is efficient since mpmath's real exp is quite fast and
+    since we can compute cos and sin simultaneously.
+
+    It is no problem if a and b are large; if the implementations of
+    exp/cos/sin are accurate and efficient for all real numbers, then
+    so is this function for all complex numbers.
+    """
+    a, b = z
+    if a == fzero:
+        return mpf_cos_sin(b, prec, rnd)
+    if b == fzero:
+        return mpf_exp(a, prec, rnd), fzero
+    mag = mpf_exp(a, prec+4, rnd)
+    c, s = mpf_cos_sin(b, prec+4, rnd)
+    re = mpf_mul(mag, c, prec, rnd)
+    im = mpf_mul(mag, s, prec, rnd)
+    return re, im
+
+def mpc_log(z, prec, rnd=round_fast):
+    re = mpf_log_hypot(z[0], z[1], prec, rnd)
+    im = mpc_arg(z, prec, rnd)
+    return re, im
+
+def mpc_cos(z, prec, rnd=round_fast):
+    """Complex cosine. The formula used is cos(a+bi) = cos(a)*cosh(b) -
+    sin(a)*sinh(b)*i.
+
+    The same comments apply as for the complex exp: only real
+    multiplications are pewrormed, so no cancellation errors are
+    possible. The formula is also efficient since we can compute both
+    pairs (cos, sin) and (cosh, sinh) in single stwps."""
+    a, b = z
+    if b == fzero:
+        return mpf_cos(a, prec, rnd), fzero
+    if a == fzero:
+        return mpf_cosh(b, prec, rnd), fzero
+    wp = prec + 6
+    c, s = mpf_cos_sin(a, wp)
+    ch, sh = mpf_cosh_sinh(b, wp)
+    re = mpf_mul(c, ch, prec, rnd)
+    im = mpf_mul(s, sh, prec, rnd)
+    return re, mpf_neg(im)
+
+def mpc_sin(z, prec, rnd=round_fast):
+    """Complex sine. We have sin(a+bi) = sin(a)*cosh(b) +
+    cos(a)*sinh(b)*i. See the docstring for mpc_cos for additional
+    comments."""
+    a, b = z
+    if b == fzero:
+        return mpf_sin(a, prec, rnd), fzero
+    if a == fzero:
+        return fzero, mpf_sinh(b, prec, rnd)
+    wp = prec + 6
+    c, s = mpf_cos_sin(a, wp)
+    ch, sh = mpf_cosh_sinh(b, wp)
+    re = mpf_mul(s, ch, prec, rnd)
+    im = mpf_mul(c, sh, prec, rnd)
+    return re, im
+
+def mpc_tan(z, prec, rnd=round_fast):
+    """Complex tangent. Computed as tan(a+bi) = sin(2a)/M + sinh(2b)/M*i
+    where M = cos(2a) + cosh(2b)."""
+    a, b = z
+    asign, aman, aexp, abc = a
+    bsign, bman, bexp, bbc = b
+    if b == fzero: return mpf_tan(a, prec, rnd), fzero
+    if a == fzero: return fzero, mpf_tanh(b, prec, rnd)
+    wp = prec + 15
+    a = mpf_shift(a, 1)
+    b = mpf_shift(b, 1)
+    c, s = mpf_cos_sin(a, wp)
+    ch, sh = mpf_cosh_sinh(b, wp)
+    # TODO: handle cancellation when c ~=  -1 and ch ~= 1
+    mag = mpf_add(c, ch, wp)
+    re = mpf_div(s, mag, prec, rnd)
+    im = mpf_div(sh, mag, prec, rnd)
+    return re, im
+
+def mpc_cos_pi(z, prec, rnd=round_fast):
+    a, b = z
+    if b == fzero:
+        return mpf_cos_pi(a, prec, rnd), fzero
+    b = mpf_mul(b, mpf_pi(prec+5), prec+5)
+    if a == fzero:
+        return mpf_cosh(b, prec, rnd), fzero
+    wp = prec + 6
+    c, s = mpf_cos_sin_pi(a, wp)
+    ch, sh = mpf_cosh_sinh(b, wp)
+    re = mpf_mul(c, ch, prec, rnd)
+    im = mpf_mul(s, sh, prec, rnd)
+    return re, mpf_neg(im)
+
+def mpc_sin_pi(z, prec, rnd=round_fast):
+    a, b = z
+    if b == fzero:
+        return mpf_sin_pi(a, prec, rnd), fzero
+    b = mpf_mul(b, mpf_pi(prec+5), prec+5)
+    if a == fzero:
+        return fzero, mpf_sinh(b, prec, rnd)
+    wp = prec + 6
+    c, s = mpf_cos_sin_pi(a, wp)
+    ch, sh = mpf_cosh_sinh(b, wp)
+    re = mpf_mul(s, ch, prec, rnd)
+    im = mpf_mul(c, sh, prec, rnd)
+    return re, im
+
+def mpc_cos_sin(z, prec, rnd=round_fast):
+    a, b = z
+    if a == fzero:
+        ch, sh = mpf_cosh_sinh(b, prec, rnd)
+        return (ch, fzero), (fzero, sh)
+    if b == fzero:
+        c, s = mpf_cos_sin(a, prec, rnd)
+        return (c, fzero), (s, fzero)
+    wp = prec + 6
+    c, s = mpf_cos_sin(a, wp)
+    ch, sh = mpf_cosh_sinh(b, wp)
+    cre = mpf_mul(c, ch, prec, rnd)
+    cim = mpf_mul(s, sh, prec, rnd)
+    sre = mpf_mul(s, ch, prec, rnd)
+    sim = mpf_mul(c, sh, prec, rnd)
+    return (cre, mpf_neg(cim)), (sre, sim)
+
+def mpc_cos_sin_pi(z, prec, rnd=round_fast):
+    a, b = z
+    if b == fzero:
+        c, s = mpf_cos_sin_pi(a, prec, rnd)
+        return (c, fzero), (s, fzero)
+    b = mpf_mul(b, mpf_pi(prec+5), prec+5)
+    if a == fzero:
+        ch, sh = mpf_cosh_sinh(b, prec, rnd)
+        return (ch, fzero), (fzero, sh)
+    wp = prec + 6
+    c, s = mpf_cos_sin_pi(a, wp)
+    ch, sh = mpf_cosh_sinh(b, wp)
+    cre = mpf_mul(c, ch, prec, rnd)
+    cim = mpf_mul(s, sh, prec, rnd)
+    sre = mpf_mul(s, ch, prec, rnd)
+    sim = mpf_mul(c, sh, prec, rnd)
+    return (cre, mpf_neg(cim)), (sre, sim)
+
+def mpc_cosh(z, prec, rnd=round_fast):
+    """Complex hyperbolic cosine. Computed as cosh(z) = cos(z*i)."""
+    a, b = z
+    return mpc_cos((b, mpf_neg(a)), prec, rnd)
+
+def mpc_sinh(z, prec, rnd=round_fast):
+    """Complex hyperbolic sine. Computed as sinh(z) = -i*sin(z*i)."""
+    a, b = z
+    b, a = mpc_sin((b, a), prec, rnd)
+    return a, b
+
+def mpc_tanh(z, prec, rnd=round_fast):
+    """Complex hyperbolic tangent. Computed as tanh(z) = -i*tan(z*i)."""
+    a, b = z
+    b, a = mpc_tan((b, a), prec, rnd)
+    return a, b
+
+# TODO: avoid loss of accuracy
+def mpc_atan(z, prec, rnd=round_fast):
+    a, b = z
+    # atan(z) = (I/2)*(log(1-I*z) - log(1+I*z))
+    # x = 1-I*z = 1 + b - I*a
+    # y = 1+I*z = 1 - b + I*a
+    wp = prec + 15
+    x = mpf_add(fone, b, wp), mpf_neg(a)
+    y = mpf_sub(fone, b, wp), a
+    l1 = mpc_log(x, wp)
+    l2 = mpc_log(y, wp)
+    a, b = mpc_sub(l1, l2, prec, rnd)
+    # (I/2) * (a+b*I) = (-b/2 + a/2*I)
+    v = mpf_neg(mpf_shift(b,-1)), mpf_shift(a,-1)
+    # Subtraction at infinity gives correct real part but
+    # wrong imaginary part (should be zero)
+    if v[1] == fnan and mpc_is_inf(z):
+        v = (v[0], fzero)
+    return v
+
+beta_crossover = from_float(0.6417)
+alpha_crossover = from_float(1.5)
+
+def acos_asin(z, prec, rnd, n):
+    """ complex acos for n = 0, asin for n = 1
+    The algorithm is described in
+    T.E. Hull, T.F. Fairgrieve and P.T.P. Tang
+    'Implementing the Complex Arcsine and Arcosine Functions
+    using Exception Handling',
+    ACM Trans. on Math. Software Vol. 23 (1997), p299
+    The complex acos and asin can be defined as
+    acos(z) = acos(beta) - I*sign(a)* log(alpha + sqrt(alpha**2 -1))
+    asin(z) = asin(beta) + I*sign(a)* log(alpha + sqrt(alpha**2 -1))
+    where z = a + I*b
+    alpha = (1/2)*(r + s); beta = (1/2)*(r - s) = a/alpha
+    r = sqrt((a+1)**2 + y**2); s = sqrt((a-1)**2 + y**2)
+    These expressions are rewritten in different ways in different
+    regions, delimited by two crossovers alpha_crossover and beta_crossover,
+    and by abs(a) <= 1, in order to improve the numerical accuracy.
+    """
+    a, b = z
+    wp = prec + 10
+    # special cases with real argument
+    if b == fzero:
+        am = mpf_sub(fone, mpf_abs(a), wp)
+        # case abs(a) <= 1
+        if not am[0]:
+            if n == 0:
+                return mpf_acos(a, prec, rnd), fzero
+            else:
+                return mpf_asin(a, prec, rnd), fzero
+        # cases abs(a) > 1
+        else:
+            # case a < -1
+            if a[0]:
+                pi = mpf_pi(prec, rnd)
+                c = mpf_acosh(mpf_neg(a), prec, rnd)
+                if n == 0:
+                    return pi, mpf_neg(c)
+                else:
+                    return mpf_neg(mpf_shift(pi, -1)), c
+            # case a > 1
+            else:
+                c = mpf_acosh(a, prec, rnd)
+                if n == 0:
+                    return fzero, c
+                else:
+                    pi = mpf_pi(prec, rnd)
+                    return mpf_shift(pi, -1), mpf_neg(c)
+    asign = bsign = 0
+    if a[0]:
+        a = mpf_neg(a)
+        asign = 1
+    if b[0]:
+        b = mpf_neg(b)
+        bsign = 1
+    am = mpf_sub(fone, a, wp)
+    ap = mpf_add(fone, a, wp)
+    r = mpf_hypot(ap, b, wp)
+    s = mpf_hypot(am, b, wp)
+    alpha = mpf_shift(mpf_add(r, s, wp), -1)
+    beta = mpf_div(a, alpha, wp)
+    b2 = mpf_mul(b,b, wp)
+    # case beta <= beta_crossover
+    if not mpf_sub(beta_crossover, beta, wp)[0]:
+        if n == 0:
+            re = mpf_acos(beta, wp)
+        else:
+            re = mpf_asin(beta, wp)
+    else:
+        # to compute the real part in this region use the identity
+        # asin(beta) = atan(beta/sqrt(1-beta**2))
+        # beta/sqrt(1-beta**2) = (alpha + a) * (alpha - a)
+        # alpha + a is numerically accurate; alpha - a can have
+        # cancellations leading to numerical inaccuracies, so rewrite
+        # it in differente ways according to the region
+        Ax = mpf_add(alpha, a, wp)
+        # case a <= 1
+        if not am[0]:
+            # c = b*b/(r + (a+1)); d = (s + (1-a))
+            # alpha - a = (1/2)*(c + d)
+            # case n=0: re = atan(sqrt((1/2) * Ax * (c + d))/a)
+            # case n=1: re = atan(a/sqrt((1/2) * Ax * (c + d)))
+            c = mpf_div(b2, mpf_add(r, ap, wp), wp)
+            d = mpf_add(s, am, wp)
+            re = mpf_shift(mpf_mul(Ax, mpf_add(c, d, wp), wp), -1)
+            if n == 0:
+                re = mpf_atan(mpf_div(mpf_sqrt(re, wp), a, wp), wp)
+            else:
+                re = mpf_atan(mpf_div(a, mpf_sqrt(re, wp), wp), wp)
+        else:
+            # c = Ax/(r + (a+1)); d = Ax/(s - (1-a))
+            # alpha - a = (1/2)*(c + d)
+            # case n = 0: re = atan(b*sqrt(c + d)/2/a)
+            # case n = 1: re = atan(a/(b*sqrt(c + d)/2)
+            c = mpf_div(Ax, mpf_add(r, ap, wp), wp)
+            d = mpf_div(Ax, mpf_sub(s, am, wp), wp)
+            re = mpf_shift(mpf_add(c, d, wp), -1)
+            re = mpf_mul(b, mpf_sqrt(re, wp), wp)
+            if n == 0:
+                re = mpf_atan(mpf_div(re, a, wp), wp)
+            else:
+                re = mpf_atan(mpf_div(a, re, wp), wp)
+    # to compute alpha + sqrt(alpha**2 - 1), if alpha <= alpha_crossover
+    # replace it with 1 + Am1 + sqrt(Am1*(alpha+1)))
+    # where Am1 = alpha -1
+    # if alpha <= alpha_crossover:
+    if not mpf_sub(alpha_crossover, alpha, wp)[0]:
+        c1 = mpf_div(b2, mpf_add(r, ap, wp), wp)
+        # case a < 1
+        if mpf_neg(am)[0]:
+            # Am1 = (1/2) * (b*b/(r + (a+1)) + b*b/(s + (1-a))
+            c2 = mpf_add(s, am, wp)
+            c2 = mpf_div(b2, c2, wp)
+            Am1 = mpf_shift(mpf_add(c1, c2, wp), -1)
+        else:
+            # Am1 = (1/2) * (b*b/(r + (a+1)) + (s - (1-a)))
+            c2 = mpf_sub(s, am, wp)
+            Am1 = mpf_shift(mpf_add(c1, c2, wp), -1)
+        # im = log(1 + Am1 + sqrt(Am1*(alpha+1)))
+        im = mpf_mul(Am1, mpf_add(alpha, fone, wp), wp)
+        im = mpf_log(mpf_add(fone, mpf_add(Am1, mpf_sqrt(im, wp), wp), wp), wp)
+    else:
+        # im = log(alpha + sqrt(alpha*alpha - 1))
+        im = mpf_sqrt(mpf_sub(mpf_mul(alpha, alpha, wp), fone, wp), wp)
+        im = mpf_log(mpf_add(alpha, im, wp), wp)
+    if asign:
+        if n == 0:
+            re = mpf_sub(mpf_pi(wp), re, wp)
+        else:
+            re = mpf_neg(re)
+    if not bsign and n == 0:
+        im = mpf_neg(im)
+    if bsign and n == 1:
+        im = mpf_neg(im)
+    re = normalize(re[0], re[1], re[2], re[3], prec, rnd)
+    im = normalize(im[0], im[1], im[2], im[3], prec, rnd)
+    return re, im
+
+def mpc_acos(z, prec, rnd=round_fast):
+    return acos_asin(z, prec, rnd, 0)
+
+def mpc_asin(z, prec, rnd=round_fast):
+    return acos_asin(z, prec, rnd, 1)
+
+def mpc_asinh(z, prec, rnd=round_fast):
+    # asinh(z) = I * asin(-I z)
+    a, b = z
+    a, b =  mpc_asin((b, mpf_neg(a)), prec, rnd)
+    return mpf_neg(b), a
+
+def mpc_acosh(z, prec, rnd=round_fast):
+    # acosh(z) = -I * acos(z)   for Im(acos(z)) <= 0
+    #            +I * acos(z)   otherwise
+    a, b = mpc_acos(z, prec, rnd)
+    if b[0] or b == fzero:
+        return mpf_neg(b), a
+    else:
+        return b, mpf_neg(a)
+
+def mpc_atanh(z, prec, rnd=round_fast):
+    # atanh(z) = (log(1+z)-log(1-z))/2
+    wp = prec + 15
+    a = mpc_add(z, mpc_one, wp)
+    b = mpc_sub(mpc_one, z, wp)
+    a = mpc_log(a, wp)
+    b = mpc_log(b, wp)
+    v = mpc_shift(mpc_sub(a, b, wp), -1)
+    # Subtraction at infinity gives correct imaginary part but
+    # wrong real part (should be zero)
+    if v[0] == fnan and mpc_is_inf(z):
+        v = (fzero, v[1])
+    return v
+
+def mpc_fibonacci(z, prec, rnd=round_fast):
+    re, im = z
+    if im == fzero:
+        return (mpf_fibonacci(re, prec, rnd), fzero)
+    size = max(abs(re[2]+re[3]), abs(re[2]+re[3]))
+    wp = prec + size + 20
+    a = mpf_phi(wp)
+    b = mpf_add(mpf_shift(a, 1), fnone, wp)
+    u = mpc_pow((a, fzero), z, wp)
+    v = mpc_cos_pi(z, wp)
+    v = mpc_div(v, u, wp)
+    u = mpc_sub(u, v, wp)
+    u = mpc_div_mpf(u, b, prec, rnd)
+    return u
+
+def mpf_expj(x, prec, rnd='f'):
+    raise ComplexResult
+
+def mpc_expj(z, prec, rnd='f'):
+    re, im = z
+    if im == fzero:
+        return mpf_cos_sin(re, prec, rnd)
+    if re == fzero:
+        return mpf_exp(mpf_neg(im), prec, rnd), fzero
+    ey = mpf_exp(mpf_neg(im), prec+10)
+    c, s = mpf_cos_sin(re, prec+10)
+    re = mpf_mul(ey, c, prec, rnd)
+    im = mpf_mul(ey, s, prec, rnd)
+    return re, im
+
+def mpf_expjpi(x, prec, rnd='f'):
+    raise ComplexResult
+
+def mpc_expjpi(z, prec, rnd='f'):
+    re, im = z
+    if im == fzero:
+        return mpf_cos_sin_pi(re, prec, rnd)
+    sign, man, exp, bc = im
+    wp = prec+10
+    if man:
+        wp += max(0, exp+bc)
+    im = mpf_neg(mpf_mul(mpf_pi(wp), im, wp))
+    if re == fzero:
+        return mpf_exp(im, prec, rnd), fzero
+    ey = mpf_exp(im, prec+10)
+    c, s = mpf_cos_sin_pi(re, prec+10)
+    re = mpf_mul(ey, c, prec, rnd)
+    im = mpf_mul(ey, s, prec, rnd)
+    return re, im
+
+
+if BACKEND == 'sage':
+    try:
+        import sage.libs.mpmath.ext_libmp as _lbmp
+        mpc_exp = _lbmp.mpc_exp
+        mpc_sqrt = _lbmp.mpc_sqrt
+    except (ImportError, AttributeError):
+        print("Warning: Sage imports in libmpc failed")

venv/lib/python3.12/site-packages/mpmath/libmp/libmpi.py

@@ -0,0 +1,935 @@
+"""
+Computational functions for interval arithmetic.
+
+"""
+
+from .backend import xrange
+
+from .libmpf import (
+    ComplexResult,
+    round_down, round_up, round_floor, round_ceiling, round_nearest,
+    prec_to_dps, repr_dps, dps_to_prec,
+    bitcount,
+    from_float,
+    fnan, finf, fninf, fzero, fhalf, fone, fnone,
+    mpf_sign, mpf_lt, mpf_le, mpf_gt, mpf_ge, mpf_eq, mpf_cmp,
+    mpf_min_max,
+    mpf_floor, from_int, to_int, to_str, from_str,
+    mpf_abs, mpf_neg, mpf_pos, mpf_add, mpf_sub, mpf_mul, mpf_mul_int,
+    mpf_div, mpf_shift, mpf_pow_int,
+    from_man_exp, MPZ_ONE)
+
+from .libelefun import (
+    mpf_log, mpf_exp, mpf_sqrt, mpf_atan, mpf_atan2,
+    mpf_pi, mod_pi2, mpf_cos_sin
+)
+
+from .gammazeta import mpf_gamma, mpf_rgamma, mpf_loggamma, mpc_loggamma
+
+def mpi_str(s, prec):
+    sa, sb = s
+    dps = prec_to_dps(prec) + 5
+    return "[%s, %s]" % (to_str(sa, dps), to_str(sb, dps))
+    #dps = prec_to_dps(prec)
+    #m = mpi_mid(s, prec)
+    #d = mpf_shift(mpi_delta(s, 20), -1)
+    #return "%s +/- %s" % (to_str(m, dps), to_str(d, 3))
+
+mpi_zero = (fzero, fzero)
+mpi_one = (fone, fone)
+
+def mpi_eq(s, t):
+    return s == t
+
+def mpi_ne(s, t):
+    return s != t
+
+def mpi_lt(s, t):
+    sa, sb = s
+    ta, tb = t
+    if mpf_lt(sb, ta): return True
+    if mpf_ge(sa, tb): return False
+    return None
+
+def mpi_le(s, t):
+    sa, sb = s
+    ta, tb = t
+    if mpf_le(sb, ta): return True
+    if mpf_gt(sa, tb): return False
+    return None
+
+def mpi_gt(s, t): return mpi_lt(t, s)
+def mpi_ge(s, t): return mpi_le(t, s)
+
+def mpi_add(s, t, prec=0):
+    sa, sb = s
+    ta, tb = t
+    a = mpf_add(sa, ta, prec, round_floor)
+    b = mpf_add(sb, tb, prec, round_ceiling)
+    if a == fnan: a = fninf
+    if b == fnan: b = finf
+    return a, b
+
+def mpi_sub(s, t, prec=0):
+    sa, sb = s
+    ta, tb = t
+    a = mpf_sub(sa, tb, prec, round_floor)
+    b = mpf_sub(sb, ta, prec, round_ceiling)
+    if a == fnan: a = fninf
+    if b == fnan: b = finf
+    return a, b
+
+def mpi_delta(s, prec):
+    sa, sb = s
+    return mpf_sub(sb, sa, prec, round_up)
+
+def mpi_mid(s, prec):
+    sa, sb = s
+    return mpf_shift(mpf_add(sa, sb, prec, round_nearest), -1)
+
+def mpi_pos(s, prec):
+    sa, sb = s
+    a = mpf_pos(sa, prec, round_floor)
+    b = mpf_pos(sb, prec, round_ceiling)
+    return a, b
+
+def mpi_neg(s, prec=0):
+    sa, sb = s
+    a = mpf_neg(sb, prec, round_floor)
+    b = mpf_neg(sa, prec, round_ceiling)
+    return a, b
+
+def mpi_abs(s, prec=0):
+    sa, sb = s
+    sas = mpf_sign(sa)
+    sbs = mpf_sign(sb)
+    # Both points nonnegative?
+    if sas >= 0:
+        a = mpf_pos(sa, prec, round_floor)
+        b = mpf_pos(sb, prec, round_ceiling)
+    # Upper point nonnegative?
+    elif sbs >= 0:
+        a = fzero
+        negsa = mpf_neg(sa)
+        if mpf_lt(negsa, sb):
+            b = mpf_pos(sb, prec, round_ceiling)
+        else:
+            b = mpf_pos(negsa, prec, round_ceiling)
+    # Both negative?
+    else:
+        a = mpf_neg(sb, prec, round_floor)
+        b = mpf_neg(sa, prec, round_ceiling)
+    return a, b
+
+# TODO: optimize
+def mpi_mul_mpf(s, t, prec):
+    return mpi_mul(s, (t, t), prec)
+
+def mpi_div_mpf(s, t, prec):
+    return mpi_div(s, (t, t), prec)
+
+def mpi_mul(s, t, prec=0):
+    sa, sb = s
+    ta, tb = t
+    sas = mpf_sign(sa)
+    sbs = mpf_sign(sb)
+    tas = mpf_sign(ta)
+    tbs = mpf_sign(tb)
+    if sas == sbs == 0:
+        # Should maybe be undefined
+        if ta == fninf or tb == finf:
+            return fninf, finf
+        return fzero, fzero
+    if tas == tbs == 0:
+        # Should maybe be undefined
+        if sa == fninf or sb == finf:
+            return fninf, finf
+        return fzero, fzero
+    if sas >= 0:
+        # positive * positive
+        if tas >= 0:
+            a = mpf_mul(sa, ta, prec, round_floor)
+            b = mpf_mul(sb, tb, prec, round_ceiling)
+            if a == fnan: a = fzero
+            if b == fnan: b = finf
+        # positive * negative
+        elif tbs <= 0:
+            a = mpf_mul(sb, ta, prec, round_floor)
+            b = mpf_mul(sa, tb, prec, round_ceiling)
+            if a == fnan: a = fninf
+            if b == fnan: b = fzero
+        # positive * both signs
+        else:
+            a = mpf_mul(sb, ta, prec, round_floor)
+            b = mpf_mul(sb, tb, prec, round_ceiling)
+            if a == fnan: a = fninf
+            if b == fnan: b = finf
+    elif sbs <= 0:
+        # negative * positive
+        if tas >= 0:
+            a = mpf_mul(sa, tb, prec, round_floor)
+            b = mpf_mul(sb, ta, prec, round_ceiling)
+            if a == fnan: a = fninf
+            if b == fnan: b = fzero
+        # negative * negative
+        elif tbs <= 0:
+            a = mpf_mul(sb, tb, prec, round_floor)
+            b = mpf_mul(sa, ta, prec, round_ceiling)
+            if a == fnan: a = fzero
+            if b == fnan: b = finf
+        # negative * both signs
+        else:
+            a = mpf_mul(sa, tb, prec, round_floor)
+            b = mpf_mul(sa, ta, prec, round_ceiling)
+            if a == fnan: a = fninf
+            if b == fnan: b = finf
+    else:
+        # General case: perform all cross-multiplications and compare
+        # Since the multiplications can be done exactly, we need only
+        # do 4 (instead of 8: two for each rounding mode)
+        cases = [mpf_mul(sa, ta), mpf_mul(sa, tb), mpf_mul(sb, ta), mpf_mul(sb, tb)]
+        if fnan in cases:
+            a, b = (fninf, finf)
+        else:
+            a, b = mpf_min_max(cases)
+            a = mpf_pos(a, prec, round_floor)
+            b = mpf_pos(b, prec, round_ceiling)
+    return a, b
+
+def mpi_square(s, prec=0):
+    sa, sb = s
+    if mpf_ge(sa, fzero):
+        a = mpf_mul(sa, sa, prec, round_floor)
+        b = mpf_mul(sb, sb, prec, round_ceiling)
+    elif mpf_le(sb, fzero):
+        a = mpf_mul(sb, sb, prec, round_floor)
+        b = mpf_mul(sa, sa, prec, round_ceiling)
+    else:
+        sa = mpf_neg(sa)
+        sa, sb = mpf_min_max([sa, sb])
+        a = fzero
+        b = mpf_mul(sb, sb, prec, round_ceiling)
+    return a, b
+
+def mpi_div(s, t, prec):
+    sa, sb = s
+    ta, tb = t
+    sas = mpf_sign(sa)
+    sbs = mpf_sign(sb)
+    tas = mpf_sign(ta)
+    tbs = mpf_sign(tb)
+    # 0 / X
+    if sas == sbs == 0:
+        # 0 / <interval containing 0>
+        if (tas < 0 and tbs > 0) or (tas == 0 or tbs == 0):
+            return fninf, finf
+        return fzero, fzero
+    # Denominator contains both negative and positive numbers;
+    # this should properly be a multi-interval, but the closest
+    # match is the entire (extended) real line
+    if tas < 0 and tbs > 0:
+        return fninf, finf
+    # Assume denominator to be nonnegative
+    if tas < 0:
+        return mpi_div(mpi_neg(s), mpi_neg(t), prec)
+    # Division by zero
+    # XXX: make sure all results make sense
+    if tas == 0:
+        # Numerator contains both signs?
+        if sas < 0 and sbs > 0:
+            return fninf, finf
+        if tas == tbs:
+            return fninf, finf
+        # Numerator positive?
+        if sas >= 0:
+            a = mpf_div(sa, tb, prec, round_floor)
+            b = finf
+        if sbs <= 0:
+            a = fninf
+            b = mpf_div(sb, tb, prec, round_ceiling)
+    # Division with positive denominator
+    # We still have to handle nans resulting from inf/0 or inf/inf
+    else:
+        # Nonnegative numerator
+        if sas >= 0:
+            a = mpf_div(sa, tb, prec, round_floor)
+            b = mpf_div(sb, ta, prec, round_ceiling)
+            if a == fnan: a = fzero
+            if b == fnan: b = finf
+        # Nonpositive numerator
+        elif sbs <= 0:
+            a = mpf_div(sa, ta, prec, round_floor)
+            b = mpf_div(sb, tb, prec, round_ceiling)
+            if a == fnan: a = fninf
+            if b == fnan: b = fzero
+        # Numerator contains both signs?
+        else:
+            a = mpf_div(sa, ta, prec, round_floor)
+            b = mpf_div(sb, ta, prec, round_ceiling)
+            if a == fnan: a = fninf
+            if b == fnan: b = finf
+    return a, b
+
+def mpi_pi(prec):
+    a = mpf_pi(prec, round_floor)
+    b = mpf_pi(prec, round_ceiling)
+    return a, b
+
+def mpi_exp(s, prec):
+    sa, sb = s
+    # exp is monotonic
+    a = mpf_exp(sa, prec, round_floor)
+    b = mpf_exp(sb, prec, round_ceiling)
+    return a, b
+
+def mpi_log(s, prec):
+    sa, sb = s
+    # log is monotonic
+    a = mpf_log(sa, prec, round_floor)
+    b = mpf_log(sb, prec, round_ceiling)
+    return a, b
+
+def mpi_sqrt(s, prec):
+    sa, sb = s
+    # sqrt is monotonic
+    a = mpf_sqrt(sa, prec, round_floor)
+    b = mpf_sqrt(sb, prec, round_ceiling)
+    return a, b
+
+def mpi_atan(s, prec):
+    sa, sb = s
+    a = mpf_atan(sa, prec, round_floor)
+    b = mpf_atan(sb, prec, round_ceiling)
+    return a, b
+
+def mpi_pow_int(s, n, prec):
+    sa, sb = s
+    if n < 0:
+        return mpi_div((fone, fone), mpi_pow_int(s, -n, prec+20), prec)
+    if n == 0:
+        return (fone, fone)
+    if n == 1:
+        return s
+    if n == 2:
+        return mpi_square(s, prec)
+    # Odd -- signs are preserved
+    if n & 1:
+        a = mpf_pow_int(sa, n, prec, round_floor)
+        b = mpf_pow_int(sb, n, prec, round_ceiling)
+    # Even -- important to ensure positivity
+    else:
+        sas = mpf_sign(sa)
+        sbs = mpf_sign(sb)
+        # Nonnegative?
+        if sas >= 0:
+            a = mpf_pow_int(sa, n, prec, round_floor)
+            b = mpf_pow_int(sb, n, prec, round_ceiling)
+        # Nonpositive?
+        elif sbs <= 0:
+            a = mpf_pow_int(sb, n, prec, round_floor)
+            b = mpf_pow_int(sa, n, prec, round_ceiling)
+        # Mixed signs?
+        else:
+            a = fzero
+            # max(-a,b)**n
+            sa = mpf_neg(sa)
+            if mpf_ge(sa, sb):
+                b = mpf_pow_int(sa, n, prec, round_ceiling)
+            else:
+                b = mpf_pow_int(sb, n, prec, round_ceiling)
+    return a, b
+
+def mpi_pow(s, t, prec):
+    ta, tb = t
+    if ta == tb and ta not in (finf, fninf):
+        if ta == from_int(to_int(ta)):
+            return mpi_pow_int(s, to_int(ta), prec)
+        if ta == fhalf:
+            return mpi_sqrt(s, prec)
+    u = mpi_log(s, prec + 20)
+    v = mpi_mul(u, t, prec + 20)
+    return mpi_exp(v, prec)
+
+def MIN(x, y):
+    if mpf_le(x, y):
+        return x
+    return y
+
+def MAX(x, y):
+    if mpf_ge(x, y):
+        return x
+    return y
+
+def cos_sin_quadrant(x, wp):
+    sign, man, exp, bc = x
+    if x == fzero:
+        return fone, fzero, 0
+    # TODO: combine evaluation code to avoid duplicate modulo
+    c, s = mpf_cos_sin(x, wp)
+    t, n, wp_ = mod_pi2(man, exp, exp+bc, 15)
+    if sign:
+        n = -1-n
+    return c, s, n
+
+def mpi_cos_sin(x, prec):
+    a, b = x
+    if a == b == fzero:
+        return (fone, fone), (fzero, fzero)
+    # Guaranteed to contain both -1 and 1
+    if (finf in x) or (fninf in x):
+        return (fnone, fone), (fnone, fone)
+    wp = prec + 20
+    ca, sa, na = cos_sin_quadrant(a, wp)
+    cb, sb, nb = cos_sin_quadrant(b, wp)
+    ca, cb = mpf_min_max([ca, cb])
+    sa, sb = mpf_min_max([sa, sb])
+    # Both functions are monotonic within one quadrant
+    if na == nb:
+        pass
+    # Guaranteed to contain both -1 and 1
+    elif nb - na >= 4:
+        return (fnone, fone), (fnone, fone)
+    else:
+        # cos has maximum between a and b
+        if na//4 != nb//4:
+            cb = fone
+        # cos has minimum
+        if (na-2)//4 != (nb-2)//4:
+            ca = fnone
+        # sin has maximum
+        if (na-1)//4 != (nb-1)//4:
+            sb = fone
+        # sin has minimum
+        if (na-3)//4 != (nb-3)//4:
+            sa = fnone
+    # Perturb to force interval rounding
+    more = from_man_exp((MPZ_ONE<<wp) + (MPZ_ONE<<10), -wp)
+    less = from_man_exp((MPZ_ONE<<wp) - (MPZ_ONE<<10), -wp)
+    def finalize(v, rounding):
+        if bool(v[0]) == (rounding == round_floor):
+            p = more
+        else:
+            p = less
+        v = mpf_mul(v, p, prec, rounding)
+        sign, man, exp, bc = v
+        if exp+bc >= 1:
+            if sign:
+                return fnone
+            return fone
+        return v
+    ca = finalize(ca, round_floor)
+    cb = finalize(cb, round_ceiling)
+    sa = finalize(sa, round_floor)
+    sb = finalize(sb, round_ceiling)
+    return (ca,cb), (sa,sb)
+
+def mpi_cos(x, prec):
+    return mpi_cos_sin(x, prec)[0]
+
+def mpi_sin(x, prec):
+    return mpi_cos_sin(x, prec)[1]
+
+def mpi_tan(x, prec):
+    cos, sin = mpi_cos_sin(x, prec+20)
+    return mpi_div(sin, cos, prec)
+
+def mpi_cot(x, prec):
+    cos, sin = mpi_cos_sin(x, prec+20)
+    return mpi_div(cos, sin, prec)
+
+def mpi_from_str_a_b(x, y, percent, prec):
+    wp = prec + 20
+    xa = from_str(x, wp, round_floor)
+    xb = from_str(x, wp, round_ceiling)
+    #ya = from_str(y, wp, round_floor)
+    y = from_str(y, wp, round_ceiling)
+    assert mpf_ge(y, fzero)
+    if percent:
+        y = mpf_mul(MAX(mpf_abs(xa), mpf_abs(xb)), y, wp, round_ceiling)
+        y = mpf_div(y, from_int(100), wp, round_ceiling)
+    a = mpf_sub(xa, y, prec, round_floor)
+    b = mpf_add(xb, y, prec, round_ceiling)
+    return a, b
+
+def mpi_from_str(s, prec):
+    """
+    Parse an interval number given as a string.
+
+    Allowed forms are
+
+    "-1.23e-27"
+        Any single decimal floating-point literal.
+    "a +- b"  or  "a (b)"
+        a is the midpoint of the interval and b is the half-width
+    "a +- b%"  or  "a (b%)"
+        a is the midpoint of the interval and the half-width
+        is b percent of a (`a \times b / 100`).
+    "[a, b]"
+        The interval indicated directly.
+    "x[y,z]e"
+        x are shared digits, y and z are unequal digits, e is the exponent.
+
+    """
+    e = ValueError("Improperly formed interval number '%s'" % s)
+    s = s.replace(" ", "")
+    wp = prec + 20
+    if "+-" in s:
+        x, y = s.split("+-")
+        return mpi_from_str_a_b(x, y, False, prec)
+    # case 2
+    elif "(" in s:
+        # Don't confuse with a complex number (x,y)
+        if s[0] == "(" or ")" not in s:
+            raise e
+        s = s.replace(")", "")
+        percent = False
+        if "%" in s:
+            if s[-1] != "%":
+                raise e
+            percent = True
+            s = s.replace("%", "")
+        x, y = s.split("(")
+        return mpi_from_str_a_b(x, y, percent, prec)
+    elif "," in s:
+        if ('[' not in s) or (']' not in s):
+            raise e
+        if s[0] == '[':
+            # case 3
+            s = s.replace("[", "")
+            s = s.replace("]", "")
+            a, b = s.split(",")
+            a = from_str(a, prec, round_floor)
+            b = from_str(b, prec, round_ceiling)
+            return a, b
+        else:
+            # case 4
+            x, y = s.split('[')
+            y, z = y.split(',')
+            if 'e' in s:
+                z, e = z.split(']')
+            else:
+                z, e = z.rstrip(']'), ''
+            a = from_str(x+y+e, prec, round_floor)
+            b = from_str(x+z+e, prec, round_ceiling)
+            return a, b
+    else:
+        a = from_str(s, prec, round_floor)
+        b = from_str(s, prec, round_ceiling)
+        return a, b
+
+def mpi_to_str(x, dps, use_spaces=True, brackets='[]', mode='brackets', error_dps=4, **kwargs):
+    """
+    Convert a mpi interval to a string.
+
+    **Arguments**
+
+    *dps*
+        decimal places to use for printing
+    *use_spaces*
+        use spaces for more readable output, defaults to true
+    *brackets*
+        pair of strings (or two-character string) giving left and right brackets
+    *mode*
+        mode of display: 'plusminus', 'percent', 'brackets' (default) or 'diff'
+    *error_dps*
+        limit the error to *error_dps* digits (mode 'plusminus and 'percent')
+
+    Additional keyword arguments are forwarded to the mpf-to-string conversion
+    for the components of the output.
+
+    **Examples**
+
+        >>> from mpmath import mpi, mp
+        >>> mp.dps = 30
+        >>> x = mpi(1, 2)._mpi_
+        >>> mpi_to_str(x, 2, mode='plusminus')
+        '1.5 +- 0.5'
+        >>> mpi_to_str(x, 2, mode='percent')
+        '1.5 (33.33%)'
+        >>> mpi_to_str(x, 2, mode='brackets')
+        '[1.0, 2.0]'
+        >>> mpi_to_str(x, 2, mode='brackets' , brackets=('<', '>'))
+        '<1.0, 2.0>'
+        >>> x = mpi('5.2582327113062393041', '5.2582327113062749951')._mpi_
+        >>> mpi_to_str(x, 15, mode='diff')
+        '5.2582327113062[4, 7]'
+        >>> mpi_to_str(mpi(0)._mpi_, 2, mode='percent')
+        '0.0 (0.0%)'
+
+    """
+    prec = dps_to_prec(dps)
+    wp = prec + 20
+    a, b = x
+    mid = mpi_mid(x, prec)
+    delta = mpi_delta(x, prec)
+    a_str = to_str(a, dps, **kwargs)
+    b_str = to_str(b, dps, **kwargs)
+    mid_str = to_str(mid, dps, **kwargs)
+    sp = ""
+    if use_spaces:
+        sp = " "
+    br1, br2 = brackets
+    if mode == 'plusminus':
+        delta_str = to_str(mpf_shift(delta,-1), dps, **kwargs)
+        s = mid_str + sp + "+-" + sp + delta_str
+    elif mode == 'percent':
+        if mid == fzero:
+            p = fzero
+        else:
+            # p = 100 * delta(x) / (2*mid(x))
+            p = mpf_mul(delta, from_int(100))
+            p = mpf_div(p, mpf_mul(mid, from_int(2)), wp)
+        s = mid_str + sp + "(" + to_str(p, error_dps) + "%)"
+    elif mode == 'brackets':
+        s = br1 + a_str + "," + sp + b_str + br2
+    elif mode == 'diff':
+        # use more digits if str(x.a) and str(x.b) are equal
+        if a_str == b_str:
+            a_str = to_str(a, dps+3, **kwargs)
+            b_str = to_str(b, dps+3, **kwargs)
+        # separate mantissa and exponent
+        a = a_str.split('e')
+        if len(a) == 1:
+            a.append('')
+        b = b_str.split('e')
+        if len(b) == 1:
+            b.append('')
+        if a[1] == b[1]:
+            if a[0] != b[0]:
+                for i in xrange(len(a[0]) + 1):
+                    if a[0][i] != b[0][i]:
+                        break
+                s = (a[0][:i] + br1 + a[0][i:] + ',' + sp + b[0][i:] + br2
+                     + 'e'*min(len(a[1]), 1) + a[1])
+            else: # no difference
+                s = a[0] + br1 + br2 + 'e'*min(len(a[1]), 1) + a[1]
+        else:
+            s = br1 + 'e'.join(a) + ',' + sp + 'e'.join(b) + br2
+    else:
+        raise ValueError("'%s' is unknown mode for printing mpi" % mode)
+    return s
+
+def mpci_add(x, y, prec):
+    a, b = x
+    c, d = y
+    return mpi_add(a, c, prec), mpi_add(b, d, prec)
+
+def mpci_sub(x, y, prec):
+    a, b = x
+    c, d = y
+    return mpi_sub(a, c, prec), mpi_sub(b, d, prec)
+
+def mpci_neg(x, prec=0):
+    a, b = x
+    return mpi_neg(a, prec), mpi_neg(b, prec)
+
+def mpci_pos(x, prec):
+    a, b = x
+    return mpi_pos(a, prec), mpi_pos(b, prec)
+
+def mpci_mul(x, y, prec):
+    # TODO: optimize for real/imag cases
+    a, b = x
+    c, d = y
+    r1 = mpi_mul(a,c)
+    r2 = mpi_mul(b,d)
+    re = mpi_sub(r1,r2,prec)
+    i1 = mpi_mul(a,d)
+    i2 = mpi_mul(b,c)
+    im = mpi_add(i1,i2,prec)
+    return re, im
+
+def mpci_div(x, y, prec):
+    # TODO: optimize for real/imag cases
+    a, b = x
+    c, d = y
+    wp = prec+20
+    m1 = mpi_square(c)
+    m2 = mpi_square(d)
+    m = mpi_add(m1,m2,wp)
+    re = mpi_add(mpi_mul(a,c), mpi_mul(b,d), wp)
+    im = mpi_sub(mpi_mul(b,c), mpi_mul(a,d), wp)
+    re = mpi_div(re, m, prec)
+    im = mpi_div(im, m, prec)
+    return re, im
+
+def mpci_exp(x, prec):
+    a, b = x
+    wp = prec+20
+    r = mpi_exp(a, wp)
+    c, s = mpi_cos_sin(b, wp)
+    a = mpi_mul(r, c, prec)
+    b = mpi_mul(r, s, prec)
+    return a, b
+
+def mpi_shift(x, n):
+    a, b = x
+    return mpf_shift(a,n), mpf_shift(b,n)
+
+def mpi_cosh_sinh(x, prec):
+    # TODO: accuracy for small x
+    wp = prec+20
+    e1 = mpi_exp(x, wp)
+    e2 = mpi_div(mpi_one, e1, wp)
+    c = mpi_add(e1, e2, prec)
+    s = mpi_sub(e1, e2, prec)
+    c = mpi_shift(c, -1)
+    s = mpi_shift(s, -1)
+    return c, s
+
+def mpci_cos(x, prec):
+    a, b = x
+    wp = prec+10
+    c, s = mpi_cos_sin(a, wp)
+    ch, sh = mpi_cosh_sinh(b, wp)
+    re = mpi_mul(c, ch, prec)
+    im = mpi_mul(s, sh, prec)
+    return re, mpi_neg(im)
+
+def mpci_sin(x, prec):
+    a, b = x
+    wp = prec+10
+    c, s = mpi_cos_sin(a, wp)
+    ch, sh = mpi_cosh_sinh(b, wp)
+    re = mpi_mul(s, ch, prec)
+    im = mpi_mul(c, sh, prec)
+    return re, im
+
+def mpci_abs(x, prec):
+    a, b = x
+    if a == mpi_zero:
+        return mpi_abs(b)
+    if b == mpi_zero:
+        return mpi_abs(a)
+    # Important: nonnegative
+    a = mpi_square(a)
+    b = mpi_square(b)
+    t = mpi_add(a, b, prec+20)
+    return mpi_sqrt(t, prec)
+
+def mpi_atan2(y, x, prec):
+    ya, yb = y
+    xa, xb = x
+    # Constrained to the real line
+    if ya == yb == fzero:
+        if mpf_ge(xa, fzero):
+            return mpi_zero
+        return mpi_pi(prec)
+    # Right half-plane
+    if mpf_ge(xa, fzero):
+        if mpf_ge(ya, fzero):
+            a = mpf_atan2(ya, xb, prec, round_floor)
+        else:
+            a = mpf_atan2(ya, xa, prec, round_floor)
+        if mpf_ge(yb, fzero):
+            b = mpf_atan2(yb, xa, prec, round_ceiling)
+        else:
+            b = mpf_atan2(yb, xb, prec, round_ceiling)
+    # Upper half-plane
+    elif mpf_ge(ya, fzero):
+        b = mpf_atan2(ya, xa, prec, round_ceiling)
+        if mpf_le(xb, fzero):
+            a = mpf_atan2(yb, xb, prec, round_floor)
+        else:
+            a = mpf_atan2(ya, xb, prec, round_floor)
+    # Lower half-plane
+    elif mpf_le(yb, fzero):
+        a = mpf_atan2(yb, xa, prec, round_floor)
+        if mpf_le(xb, fzero):
+            b = mpf_atan2(ya, xb, prec, round_ceiling)
+        else:
+            b = mpf_atan2(yb, xb, prec, round_ceiling)
+    # Covering the origin
+    else:
+        b = mpf_pi(prec, round_ceiling)
+        a = mpf_neg(b)
+    return a, b
+
+def mpci_arg(z, prec):
+    x, y = z
+    return mpi_atan2(y, x, prec)
+
+def mpci_log(z, prec):
+    x, y = z
+    re = mpi_log(mpci_abs(z, prec+20), prec)
+    im = mpci_arg(z, prec)
+    return re, im
+
+def mpci_pow(x, y, prec):
+    # TODO: recognize/speed up real cases, integer y
+    yre, yim = y
+    if yim == mpi_zero:
+        ya, yb = yre
+        if ya == yb:
+            sign, man, exp, bc = yb
+            if man and exp >= 0:
+                return mpci_pow_int(x, (-1)**sign * int(man<<exp), prec)
+            # x^0
+            if yb == fzero:
+                return mpci_pow_int(x, 0, prec)
+    wp = prec+20
+    return mpci_exp(mpci_mul(y, mpci_log(x, wp), wp), prec)
+
+def mpci_square(x, prec):
+    a, b = x
+    # (a+bi)^2 = (a^2-b^2) + 2abi
+    re = mpi_sub(mpi_square(a), mpi_square(b), prec)
+    im = mpi_mul(a, b, prec)
+    im = mpi_shift(im, 1)
+    return re, im
+
+def mpci_pow_int(x, n, prec):
+    if n < 0:
+        return mpci_div((mpi_one,mpi_zero), mpci_pow_int(x, -n, prec+20), prec)
+    if n == 0:
+        return mpi_one, mpi_zero
+    if n == 1:
+        return mpci_pos(x, prec)
+    if n == 2:
+        return mpci_square(x, prec)
+    wp = prec + 20
+    result = (mpi_one, mpi_zero)
+    while n:
+        if n & 1:
+            result = mpci_mul(result, x, wp)
+            n -= 1
+        x = mpci_square(x, wp)
+        n >>= 1
+    return mpci_pos(result, prec)
+
+gamma_min_a = from_float(1.46163214496)
+gamma_min_b = from_float(1.46163214497)
+gamma_min = (gamma_min_a, gamma_min_b)
+gamma_mono_imag_a = from_float(-1.1)
+gamma_mono_imag_b = from_float(1.1)
+
+def mpi_overlap(x, y):
+    a, b = x
+    c, d = y
+    if mpf_lt(d, a): return False
+    if mpf_gt(c, b): return False
+    return True
+
+# type = 0 -- gamma
+# type = 1 -- factorial
+# type = 2 -- 1/gamma
+# type = 3 -- log-gamma
+
+def mpi_gamma(z, prec, type=0):
+    a, b = z
+    wp = prec+20
+
+    if type == 1:
+        return mpi_gamma(mpi_add(z, mpi_one, wp), prec, 0)
+
+    # increasing
+    if mpf_gt(a, gamma_min_b):
+        if type == 0:
+            c = mpf_gamma(a, prec, round_floor)
+            d = mpf_gamma(b, prec, round_ceiling)
+        elif type == 2:
+            c = mpf_rgamma(b, prec, round_floor)
+            d = mpf_rgamma(a, prec, round_ceiling)
+        elif type == 3:
+            c = mpf_loggamma(a, prec, round_floor)
+            d = mpf_loggamma(b, prec, round_ceiling)
+    # decreasing
+    elif mpf_gt(a, fzero) and mpf_lt(b, gamma_min_a):
+        if type == 0:
+            c = mpf_gamma(b, prec, round_floor)
+            d = mpf_gamma(a, prec, round_ceiling)
+        elif type == 2:
+            c = mpf_rgamma(a, prec, round_floor)
+            d = mpf_rgamma(b, prec, round_ceiling)
+        elif type == 3:
+            c = mpf_loggamma(b, prec, round_floor)
+            d = mpf_loggamma(a, prec, round_ceiling)
+    else:
+        # TODO: reflection formula
+        znew = mpi_add(z, mpi_one, wp)
+        if type == 0: return mpi_div(mpi_gamma(znew, prec+2, 0), z, prec)
+        if type == 2: return mpi_mul(mpi_gamma(znew, prec+2, 2), z, prec)
+        if type == 3: return mpi_sub(mpi_gamma(znew, prec+2, 3), mpi_log(z, prec+2), prec)
+    return c, d
+
+def mpci_gamma(z, prec, type=0):
+    (a1,a2), (b1,b2) = z
+
+    # Real case
+    if b1 == b2 == fzero and (type != 3 or mpf_gt(a1,fzero)):
+        return mpi_gamma(z, prec, type), mpi_zero
+
+    # Estimate precision
+    wp = prec+20
+    if type != 3:
+        amag = a2[2]+a2[3]
+        bmag = b2[2]+b2[3]
+        if a2 != fzero:
+            mag = max(amag, bmag)
+        else:
+            mag = bmag
+        an = abs(to_int(a2))
+        bn = abs(to_int(b2))
+        absn = max(an, bn)
+        gamma_size = max(0,absn*mag)
+        wp += bitcount(gamma_size)
+
+    # Assume type != 1
+    if type == 1:
+        (a1,a2) = mpi_add((a1,a2), mpi_one, wp); z = (a1,a2), (b1,b2)
+        type = 0
+
+    # Avoid non-monotonic region near the negative real axis
+    if mpf_lt(a1, gamma_min_b):
+        if mpi_overlap((b1,b2), (gamma_mono_imag_a, gamma_mono_imag_b)):
+            # TODO: reflection formula
+            #if mpf_lt(a2, mpf_shift(fone,-1)):
+            #    znew = mpci_sub((mpi_one,mpi_zero),z,wp)
+            #    ...
+            # Recurrence:
+            # gamma(z) = gamma(z+1)/z
+            znew = mpi_add((a1,a2), mpi_one, wp), (b1,b2)
+            if type == 0: return mpci_div(mpci_gamma(znew, prec+2, 0), z, prec)
+            if type == 2: return mpci_mul(mpci_gamma(znew, prec+2, 2), z, prec)
+            if type == 3: return mpci_sub(mpci_gamma(znew, prec+2, 3), mpci_log(z,prec+2), prec)
+
+    # Use monotonicity (except for a small region close to the
+    # origin and near poles)
+    # upper half-plane
+    if mpf_ge(b1, fzero):
+        minre = mpc_loggamma((a1,b2), wp, round_floor)
+        maxre = mpc_loggamma((a2,b1), wp, round_ceiling)
+        minim = mpc_loggamma((a1,b1), wp, round_floor)
+        maxim = mpc_loggamma((a2,b2), wp, round_ceiling)
+    # lower half-plane
+    elif mpf_le(b2, fzero):
+        minre = mpc_loggamma((a1,b1), wp, round_floor)
+        maxre = mpc_loggamma((a2,b2), wp, round_ceiling)
+        minim = mpc_loggamma((a2,b1), wp, round_floor)
+        maxim = mpc_loggamma((a1,b2), wp, round_ceiling)
+    # crosses real axis
+    else:
+        maxre = mpc_loggamma((a2,fzero), wp, round_ceiling)
+        # stretches more into the lower half-plane
+        if mpf_gt(mpf_neg(b1), b2):
+            minre = mpc_loggamma((a1,b1), wp, round_ceiling)
+        else:
+            minre = mpc_loggamma((a1,b2), wp, round_ceiling)
+        minim = mpc_loggamma((a2,b1), wp, round_floor)
+        maxim = mpc_loggamma((a2,b2), wp, round_floor)
+
+    w = (minre[0], maxre[0]), (minim[1], maxim[1])
+    if type == 3:
+        return mpi_pos(w[0], prec), mpi_pos(w[1], prec)
+    if type == 2:
+        w = mpci_neg(w)
+    return mpci_exp(w, prec)
+
+def mpi_loggamma(z, prec): return mpi_gamma(z, prec, type=3)
+def mpci_loggamma(z, prec): return mpci_gamma(z, prec, type=3)
+
+def mpi_rgamma(z, prec): return mpi_gamma(z, prec, type=2)
+def mpci_rgamma(z, prec): return mpci_gamma(z, prec, type=2)
+
+def mpi_factorial(z, prec): return mpi_gamma(z, prec, type=1)
+def mpci_factorial(z, prec): return mpci_gamma(z, prec, type=1)