test_functions.py

from mpmath.libmp import *
from mpmath import *
import random
import time
import math
import cmath

def mpc_ae(a, b, eps=eps):
    res = True
    res = res and a.real.ae(b.real, eps)
    res = res and a.imag.ae(b.imag, eps)
    return res

#----------------------------------------------------------------------------
# Constants and functions
#

tpi = "3.1415926535897932384626433832795028841971693993751058209749445923078\
1640628620899862803482534211706798"
te = "2.71828182845904523536028747135266249775724709369995957496696762772407\
663035354759457138217852516642743"
tdegree = "0.017453292519943295769236907684886127134428718885417254560971914\
4017100911460344944368224156963450948221"
teuler = "0.5772156649015328606065120900824024310421593359399235988057672348\
84867726777664670936947063291746749516"
tln2 = "0.693147180559945309417232121458176568075500134360255254120680009493\
393621969694715605863326996418687542"
tln10 = "2.30258509299404568401799145468436420760110148862877297603332790096\
757260967735248023599720508959829834"
tcatalan = "0.91596559417721901505460351493238411077414937428167213426649811\
9621763019776254769479356512926115106249"
tkhinchin = "2.6854520010653064453097148354817956938203822939944629530511523\
4555721885953715200280114117493184769800"
tglaisher = "1.2824271291006226368753425688697917277676889273250011920637400\
2174040630885882646112973649195820237439420646"
tapery = "1.2020569031595942853997381615114499907649862923404988817922715553\
4183820578631309018645587360933525815"
tphi = "1.618033988749894848204586834365638117720309179805762862135448622705\
26046281890244970720720418939113748475"
tmertens = "0.26149721284764278375542683860869585905156664826119920619206421\
3924924510897368209714142631434246651052"
ttwinprime = "0.660161815846869573927812110014555778432623360284733413319448\
423335405642304495277143760031413839867912"

def test_constants():
    for prec in [3, 7, 10, 15, 20, 37, 80, 100, 29]:
        mp.dps = prec
        assert pi == mpf(tpi)
        assert e == mpf(te)
        assert degree == mpf(tdegree)
        assert euler == mpf(teuler)
        assert ln2 == mpf(tln2)
        assert ln10 == mpf(tln10)
        assert catalan == mpf(tcatalan)
        assert khinchin == mpf(tkhinchin)
        assert glaisher == mpf(tglaisher)
        assert phi == mpf(tphi)
        if prec < 50:
            assert mertens == mpf(tmertens)
            assert twinprime == mpf(ttwinprime)
    mp.dps = 15
    assert pi >= -1
    assert pi > 2
    assert pi > 3
    assert pi < 4

def test_exact_sqrts():
    for i in range(20000):
        assert sqrt(mpf(i*i)) == i
    random.seed(1)
    for prec in [100, 300, 1000, 10000]:
        mp.dps = prec
        for i in range(20):
            A = random.randint(10**(prec//2-2), 10**(prec//2-1))
            assert sqrt(mpf(A*A)) == A
    mp.dps = 15
    for i in range(100):
        for a in [1, 8, 25, 112307]:
            assert sqrt(mpf((a*a, 2*i))) == mpf((a, i))
            assert sqrt(mpf((a*a, -2*i))) == mpf((a, -i))

def test_sqrt_rounding():
    for i in [2, 3, 5, 6, 7, 8, 10, 11, 12, 13, 14, 15]:
        i = from_int(i)
        for dps in [7, 15, 83, 106, 2000]:
            mp.dps = dps
            a = mpf_pow_int(mpf_sqrt(i, mp.prec, round_down), 2, mp.prec, round_down)
            b = mpf_pow_int(mpf_sqrt(i, mp.prec, round_up), 2, mp.prec, round_up)
            assert mpf_lt(a, i)
            assert mpf_gt(b, i)
    random.seed(1234)
    prec = 100
    for rnd in [round_down, round_nearest, round_ceiling]:
        for i in range(100):
            a = mpf_rand(prec)
            b = mpf_mul(a, a)
            assert mpf_sqrt(b, prec, rnd) == a
    # Test some extreme cases
    mp.dps = 100
    a = mpf(9) + 1e-90
    b = mpf(9) - 1e-90
    mp.dps = 15
    assert sqrt(a, rounding='d') == 3
    assert sqrt(a, rounding='n') == 3
    assert sqrt(a, rounding='u') > 3
    assert sqrt(b, rounding='d') < 3
    assert sqrt(b, rounding='n') == 3
    assert sqrt(b, rounding='u') == 3
    # A worst case, from the MPFR test suite
    assert sqrt(mpf('7.0503726185518891')) == mpf('2.655253776675949')

def test_float_sqrt():
    mp.dps = 15
    # These should round identically
    for x in [0, 1e-7, 0.1, 0.5, 1, 2, 3, 4, 5, 0.333, 76.19]:
        assert sqrt(mpf(x)) == float(x)**0.5
    assert sqrt(-1) == 1j
    assert sqrt(-2).ae(cmath.sqrt(-2))
    assert sqrt(-3).ae(cmath.sqrt(-3))
    assert sqrt(-100).ae(cmath.sqrt(-100))
    assert sqrt(1j).ae(cmath.sqrt(1j))
    assert sqrt(-1j).ae(cmath.sqrt(-1j))
    assert sqrt(math.pi + math.e*1j).ae(cmath.sqrt(math.pi + math.e*1j))
    assert sqrt(math.pi - math.e*1j).ae(cmath.sqrt(math.pi - math.e*1j))

def test_hypot():
    assert hypot(0, 0) == 0
    assert hypot(0, 0.33) == mpf(0.33)
    assert hypot(0.33, 0) == mpf(0.33)
    assert hypot(-0.33, 0) == mpf(0.33)
    assert hypot(3, 4) == mpf(5)

def test_exact_cbrt():
    for i in range(0, 20000, 200):
        assert cbrt(mpf(i*i*i)) == i
    random.seed(1)
    for prec in [100, 300, 1000, 10000]:
        mp.dps = prec
        A = random.randint(10**(prec//2-2), 10**(prec//2-1))
        assert cbrt(mpf(A*A*A)) == A
    mp.dps = 15

def test_exp():
    assert exp(0) == 1
    assert exp(10000).ae(mpf('8.8068182256629215873e4342'))
    assert exp(-10000).ae(mpf('1.1354838653147360985e-4343'))
    a = exp(mpf((1, 8198646019315405, -53, 53)))
    assert(a.bc == bitcount(a.man))
    mp.prec = 67
    a = exp(mpf((1, 1781864658064754565, -60, 61)))
    assert(a.bc == bitcount(a.man))
    mp.prec = 53
    assert exp(ln2 * 10).ae(1024)
    assert exp(2+2j).ae(cmath.exp(2+2j))

def test_issue_73():
    mp.dps = 512
    a = exp(-1)
    b = exp(1)
    mp.dps = 15
    assert (+a).ae(0.36787944117144233)
    assert (+b).ae(2.7182818284590451)

def test_log():
    mp.dps = 15
    assert log(1) == 0
    for x in [0.5, 1.5, 2.0, 3.0, 100, 10**50, 1e-50]:
        assert log(x).ae(math.log(x))
        assert log(x, x) == 1
    assert log(1024, 2) == 10
    assert log(10**1234, 10) == 1234
    assert log(2+2j).ae(cmath.log(2+2j))
    # Accuracy near 1
    assert (log(0.6+0.8j).real*10**17).ae(2.2204460492503131)
    assert (log(0.6-0.8j).real*10**17).ae(2.2204460492503131)
    assert (log(0.8-0.6j).real*10**17).ae(2.2204460492503131)
    assert (log(1+1e-8j).real*10**16).ae(0.5)
    assert (log(1-1e-8j).real*10**16).ae(0.5)
    assert (log(-1+1e-8j).real*10**16).ae(0.5)
    assert (log(-1-1e-8j).real*10**16).ae(0.5)
    assert (log(1j+1e-8).real*10**16).ae(0.5)
    assert (log(1j-1e-8).real*10**16).ae(0.5)
    assert (log(-1j+1e-8).real*10**16).ae(0.5)
    assert (log(-1j-1e-8).real*10**16).ae(0.5)
    assert (log(1+1e-40j).real*10**80).ae(0.5)
    assert (log(1j+1e-40).real*10**80).ae(0.5)
    # Huge
    assert log(ldexp(1.234,10**20)).ae(log(2)*1e20)
    assert log(ldexp(1.234,10**200)).ae(log(2)*1e200)
    # Some special values
    assert log(mpc(0,0)) == mpc(-inf,0)
    assert isnan(log(mpc(nan,0)).real)
    assert isnan(log(mpc(nan,0)).imag)
    assert isnan(log(mpc(0,nan)).real)
    assert isnan(log(mpc(0,nan)).imag)
    assert isnan(log(mpc(nan,1)).real)
    assert isnan(log(mpc(nan,1)).imag)
    assert isnan(log(mpc(1,nan)).real)
    assert isnan(log(mpc(1,nan)).imag)

def test_trig_hyperb_basic():
    for x in (list(range(100)) + list(range(-100,0))):
        t = x / 4.1
        assert cos(mpf(t)).ae(math.cos(t))
        assert sin(mpf(t)).ae(math.sin(t))
        assert tan(mpf(t)).ae(math.tan(t))
        assert cosh(mpf(t)).ae(math.cosh(t))
        assert sinh(mpf(t)).ae(math.sinh(t))
        assert tanh(mpf(t)).ae(math.tanh(t))
    assert sin(1+1j).ae(cmath.sin(1+1j))
    assert sin(-4-3.6j).ae(cmath.sin(-4-3.6j))
    assert cos(1+1j).ae(cmath.cos(1+1j))
    assert cos(-4-3.6j).ae(cmath.cos(-4-3.6j))

def test_degrees():
    assert cos(0*degree) == 1
    assert cos(90*degree).ae(0)
    assert cos(180*degree).ae(-1)
    assert cos(270*degree).ae(0)
    assert cos(360*degree).ae(1)
    assert sin(0*degree) == 0
    assert sin(90*degree).ae(1)
    assert sin(180*degree).ae(0)
    assert sin(270*degree).ae(-1)
    assert sin(360*degree).ae(0)

def random_complexes(N):
    random.seed(1)
    a = []
    for i in range(N):
        x1 = random.uniform(-10, 10)
        y1 = random.uniform(-10, 10)
        x2 = random.uniform(-10, 10)
        y2 = random.uniform(-10, 10)
        z1 = complex(x1, y1)
        z2 = complex(x2, y2)
        a.append((z1, z2))
    return a

def test_complex_powers():
    for dps in [15, 30, 100]:
        # Check accuracy for complex square root
        mp.dps = dps
        a = mpc(1j)**0.5
        assert a.real == a.imag == mpf(2)**0.5 / 2
    mp.dps = 15
    random.seed(1)
    for (z1, z2) in random_complexes(100):
        assert (mpc(z1)**mpc(z2)).ae(z1**z2, 1e-12)
    assert (e**(-pi*1j)).ae(-1)
    mp.dps = 50
    assert (e**(-pi*1j)).ae(-1)
    mp.dps = 15

def test_complex_sqrt_accuracy():
    def test_mpc_sqrt(lst):
        for a, b in lst:
            z = mpc(a + j*b)
            assert mpc_ae(sqrt(z*z), z)
            z = mpc(-a + j*b)
            assert mpc_ae(sqrt(z*z), -z)
            z = mpc(a - j*b)
            assert mpc_ae(sqrt(z*z), z)
            z = mpc(-a - j*b)
            assert mpc_ae(sqrt(z*z), -z)
    random.seed(2)
    N = 10
    mp.dps = 30
    dps = mp.dps
    test_mpc_sqrt([(random.uniform(0, 10),random.uniform(0, 10)) for i in range(N)])
    test_mpc_sqrt([(i + 0.1, (i + 0.2)*10**i) for i in range(N)])
    mp.dps = 15

def test_atan():
    mp.dps = 15
    assert atan(-2.3).ae(math.atan(-2.3))
    assert atan(1e-50) == 1e-50
    assert atan(1e50).ae(pi/2)
    assert atan(-1e-50) == -1e-50
    assert atan(-1e50).ae(-pi/2)
    assert atan(10**1000).ae(pi/2)
    for dps in [25, 70, 100, 300, 1000]:
        mp.dps = dps
        assert (4*atan(1)).ae(pi)
    mp.dps = 15
    pi2 = pi/2
    assert atan(mpc(inf,-1)).ae(pi2)
    assert atan(mpc(inf,0)).ae(pi2)
    assert atan(mpc(inf,1)).ae(pi2)
    assert atan(mpc(1,inf)).ae(pi2)
    assert atan(mpc(0,inf)).ae(pi2)
    assert atan(mpc(-1,inf)).ae(-pi2)
    assert atan(mpc(-inf,1)).ae(-pi2)
    assert atan(mpc(-inf,0)).ae(-pi2)
    assert atan(mpc(-inf,-1)).ae(-pi2)
    assert atan(mpc(-1,-inf)).ae(-pi2)
    assert atan(mpc(0,-inf)).ae(-pi2)
    assert atan(mpc(1,-inf)).ae(pi2)

def test_atan2():
    mp.dps = 15
    assert atan2(1,1).ae(pi/4)
    assert atan2(1,-1).ae(3*pi/4)
    assert atan2(-1,-1).ae(-3*pi/4)
    assert atan2(-1,1).ae(-pi/4)
    assert atan2(-1,0).ae(-pi/2)
    assert atan2(1,0).ae(pi/2)
    assert atan2(0,0) == 0
    assert atan2(inf,0).ae(pi/2)
    assert atan2(-inf,0).ae(-pi/2)
    assert isnan(atan2(inf,inf))
    assert isnan(atan2(-inf,inf))
    assert isnan(atan2(inf,-inf))
    assert isnan(atan2(3,nan))
    assert isnan(atan2(nan,3))
    assert isnan(atan2(0,nan))
    assert isnan(atan2(nan,0))
    assert atan2(0,inf) == 0
    assert atan2(0,-inf).ae(pi)
    assert atan2(10,inf) == 0
    assert atan2(-10,inf) == 0
    assert atan2(-10,-inf).ae(-pi)
    assert atan2(10,-inf).ae(pi)
    assert atan2(inf,10).ae(pi/2)
    assert atan2(inf,-10).ae(pi/2)
    assert atan2(-inf,10).ae(-pi/2)
    assert atan2(-inf,-10).ae(-pi/2)

def test_areal_inverses():
    assert asin(mpf(0)) == 0
    assert asinh(mpf(0)) == 0
    assert acosh(mpf(1)) == 0
    assert isinstance(asin(mpf(0.5)), mpf)
    assert isinstance(asin(mpf(2.0)), mpc)
    assert isinstance(acos(mpf(0.5)), mpf)
    assert isinstance(acos(mpf(2.0)), mpc)
    assert isinstance(atanh(mpf(0.1)), mpf)
    assert isinstance(atanh(mpf(1.1)), mpc)

    random.seed(1)
    for i in range(50):
        x = random.uniform(0, 1)
        assert asin(mpf(x)).ae(math.asin(x))
        assert acos(mpf(x)).ae(math.acos(x))

        x = random.uniform(-10, 10)
        assert asinh(mpf(x)).ae(cmath.asinh(x).real)
        assert isinstance(asinh(mpf(x)), mpf)
        x = random.uniform(1, 10)
        assert acosh(mpf(x)).ae(cmath.acosh(x).real)
        assert isinstance(acosh(mpf(x)), mpf)
        x = random.uniform(-10, 0.999)
        assert isinstance(acosh(mpf(x)), mpc)

        x = random.uniform(-1, 1)
        assert atanh(mpf(x)).ae(cmath.atanh(x).real)
        assert isinstance(atanh(mpf(x)), mpf)

    dps = mp.dps
    mp.dps = 300
    assert isinstance(asin(0.5), mpf)
    mp.dps = 1000
    assert asin(1).ae(pi/2)
    assert asin(-1).ae(-pi/2)
    mp.dps = dps

def test_invhyperb_inaccuracy():
    mp.dps = 15
    assert (asinh(1e-5)*10**5).ae(0.99999999998333333)
    assert (asinh(1e-10)*10**10).ae(1)
    assert (asinh(1e-50)*10**50).ae(1)
    assert (asinh(-1e-5)*10**5).ae(-0.99999999998333333)
    assert (asinh(-1e-10)*10**10).ae(-1)
    assert (asinh(-1e-50)*10**50).ae(-1)
    assert asinh(10**20).ae(46.744849040440862)
    assert asinh(-10**20).ae(-46.744849040440862)
    assert (tanh(1e-10)*10**10).ae(1)
    assert (tanh(-1e-10)*10**10).ae(-1)
    assert (atanh(1e-10)*10**10).ae(1)
    assert (atanh(-1e-10)*10**10).ae(-1)

def test_complex_functions():
    for x in (list(range(10)) + list(range(-10,0))):
        for y in (list(range(10)) + list(range(-10,0))):
            z = complex(x, y)/4.3 + 0.01j
            assert exp(mpc(z)).ae(cmath.exp(z))
            assert log(mpc(z)).ae(cmath.log(z))
            assert cos(mpc(z)).ae(cmath.cos(z))
            assert sin(mpc(z)).ae(cmath.sin(z))
            assert tan(mpc(z)).ae(cmath.tan(z))
            assert sinh(mpc(z)).ae(cmath.sinh(z))
            assert cosh(mpc(z)).ae(cmath.cosh(z))
            assert tanh(mpc(z)).ae(cmath.tanh(z))

def test_complex_inverse_functions():
    mp.dps = 15
    iv.dps = 15
    for (z1, z2) in random_complexes(30):
        # apparently cmath uses a different branch, so we
        # can't use it for comparison
        assert sinh(asinh(z1)).ae(z1)
        #
        assert acosh(z1).ae(cmath.acosh(z1))
        assert atanh(z1).ae(cmath.atanh(z1))
        assert atan(z1).ae(cmath.atan(z1))
        # the reason we set a big eps here is that the cmath
        # functions are inaccurate
        assert asin(z1).ae(cmath.asin(z1), rel_eps=1e-12)
        assert acos(z1).ae(cmath.acos(z1), rel_eps=1e-12)
        one = mpf(1)
    for i in range(-9, 10, 3):
        for k in range(-9, 10, 3):
            a = 0.9*j*10**k + 0.8*one*10**i
            b = cos(acos(a))
            assert b.ae(a)
            b = sin(asin(a))
            assert b.ae(a)
    one = mpf(1)
    err = 2*10**-15
    for i in range(-9, 9, 3):
        for k in range(-9, 9, 3):
            a = -0.9*10**k + j*0.8*one*10**i
            b = cosh(acosh(a))
            assert b.ae(a, err)
            b = sinh(asinh(a))
            assert b.ae(a, err)

def test_reciprocal_functions():
    assert sec(3).ae(-1.01010866590799375)
    assert csc(3).ae(7.08616739573718592)
    assert cot(3).ae(-7.01525255143453347)
    assert sech(3).ae(0.0993279274194332078)
    assert csch(3).ae(0.0998215696688227329)
    assert coth(3).ae(1.00496982331368917)
    assert asec(3).ae(1.23095941734077468)
    assert acsc(3).ae(0.339836909454121937)
    assert acot(3).ae(0.321750554396642193)
    assert asech(0.5).ae(1.31695789692481671)
    assert acsch(3).ae(0.327450150237258443)
    assert acoth(3).ae(0.346573590279972655)
    assert acot(0).ae(1.5707963267948966192)
    assert acoth(0).ae(1.5707963267948966192j)

def test_ldexp():
    mp.dps = 15
    assert ldexp(mpf(2.5), 0) == 2.5
    assert ldexp(mpf(2.5), -1) == 1.25
    assert ldexp(mpf(2.5), 2) == 10
    assert ldexp(mpf('inf'), 3) == mpf('inf')

def test_frexp():
    mp.dps = 15
    assert frexp(0) == (0.0, 0)
    assert frexp(9) == (0.5625, 4)
    assert frexp(1) == (0.5, 1)
    assert frexp(0.2) == (0.8, -2)
    assert frexp(1000) == (0.9765625, 10)

def test_aliases():
    assert ln(7) == log(7)
    assert log10(3.75) == log(3.75,10)
    assert degrees(5.6) == 5.6 / degree
    assert radians(5.6) == 5.6 * degree
    assert power(-1,0.5) == j
    assert fmod(25,7) == 4.0 and isinstance(fmod(25,7), mpf)

def test_arg_sign():
    assert arg(3) == 0
    assert arg(-3).ae(pi)
    assert arg(j).ae(pi/2)
    assert arg(-j).ae(-pi/2)
    assert arg(0) == 0
    assert isnan(atan2(3,nan))
    assert isnan(atan2(nan,3))
    assert isnan(atan2(0,nan))
    assert isnan(atan2(nan,0))
    assert isnan(atan2(nan,nan))
    assert arg(inf) == 0
    assert arg(-inf).ae(pi)
    assert isnan(arg(nan))
    #assert arg(inf*j).ae(pi/2)
    assert sign(0) == 0
    assert sign(3) == 1
    assert sign(-3) == -1
    assert sign(inf) == 1
    assert sign(-inf) == -1
    assert isnan(sign(nan))
    assert sign(j) == j
    assert sign(-3*j) == -j
    assert sign(1+j).ae((1+j)/sqrt(2))

def test_misc_bugs():
    # test that this doesn't raise an exception
    mp.dps = 1000
    log(1302)
    mp.dps = 15

def test_arange():
    assert arange(10) == [mpf('0.0'), mpf('1.0'), mpf('2.0'), mpf('3.0'),
                          mpf('4.0'), mpf('5.0'), mpf('6.0'), mpf('7.0'),
                          mpf('8.0'), mpf('9.0')]
    assert arange(-5, 5) == [mpf('-5.0'), mpf('-4.0'), mpf('-3.0'),
                             mpf('-2.0'), mpf('-1.0'), mpf('0.0'),
                             mpf('1.0'), mpf('2.0'), mpf('3.0'), mpf('4.0')]
    assert arange(0, 1, 0.1) == [mpf('0.0'), mpf('0.10000000000000001'),
                                 mpf('0.20000000000000001'),
                                 mpf('0.30000000000000004'),
                                 mpf('0.40000000000000002'),
                                 mpf('0.5'), mpf('0.60000000000000009'),
                                 mpf('0.70000000000000007'),
                                 mpf('0.80000000000000004'),
                                 mpf('0.90000000000000002')]
    assert arange(17, -9, -3) == [mpf('17.0'), mpf('14.0'), mpf('11.0'),
                                  mpf('8.0'), mpf('5.0'), mpf('2.0'),
                                  mpf('-1.0'), mpf('-4.0'), mpf('-7.0')]
    assert arange(0.2, 0.1, -0.1) == [mpf('0.20000000000000001')]
    assert arange(0) == []
    assert arange(1000, -1) == []
    assert arange(-1.23, 3.21, -0.0000001) == []

def test_linspace():
    assert linspace(2, 9, 7) == [mpf('2.0'), mpf('3.166666666666667'),
        mpf('4.3333333333333339'), mpf('5.5'), mpf('6.666666666666667'),
        mpf('7.8333333333333339'), mpf('9.0')]
    assert linspace(2, 9, 7, endpoint=0) == [mpf('2.0'), mpf('3.0'), mpf('4.0'),
        mpf('5.0'), mpf('6.0'), mpf('7.0'), mpf('8.0')]
    assert linspace(2, 7, 1) == [mpf(2)]

def test_float_cbrt():
    mp.dps = 30
    for a in arange(0,10,0.1):
        assert cbrt(a*a*a).ae(a, eps)
    assert cbrt(-1).ae(0.5 + j*sqrt(3)/2)
    one_third = mpf(1)/3
    for a in arange(0,10,2.7) + [0.1 + 10**5]:
        a = mpc(a + 1.1j)
        r1 = cbrt(a)
        mp.dps += 10
        r2 = pow(a, one_third)
        mp.dps -= 10
        assert r1.ae(r2, eps)
    mp.dps = 100
    for n in range(100, 301, 100):
        w = 10**n + j*10**-3
        z = w*w*w
        r = cbrt(z)
        assert mpc_ae(r, w, eps)
    mp.dps = 15

def test_root():
    mp.dps = 30
    random.seed(1)
    a = random.randint(0, 10000)
    p = a*a*a
    r = nthroot(mpf(p), 3)
    assert r == a
    for n in range(4, 10):
        p = p*a
        assert nthroot(mpf(p), n) == a
    mp.dps = 40
    for n in range(10, 5000, 100):
        for a in [random.random()*10000, random.random()*10**100]:
            r = nthroot(a, n)
            r1 = pow(a, mpf(1)/n)
            assert r.ae(r1)
            r = nthroot(a, -n)
            r1 = pow(a, -mpf(1)/n)
            assert r.ae(r1)
    # XXX: this is broken right now
    # tests for nthroot rounding
    for rnd in ['nearest', 'up', 'down']:
        mp.rounding = rnd
        for n in [-5, -3, 3, 5]:
            prec = 50
            for i in range(10):
                mp.prec = prec
                a = rand()
                mp.prec = 2*prec
                b = a**n
                mp.prec = prec
                r = nthroot(b, n)
                assert r == a
    mp.dps = 30
    for n in range(3, 21):
        a = (random.random() + j*random.random())
        assert nthroot(a, n).ae(pow(a, mpf(1)/n))
        assert mpc_ae(nthroot(a, n), pow(a, mpf(1)/n))
        a = (random.random()*10**100 + j*random.random())
        r = nthroot(a, n)
        mp.dps += 4
        r1 = pow(a, mpf(1)/n)
        mp.dps -= 4
        assert r.ae(r1)
        assert mpc_ae(r, r1, eps)
        r = nthroot(a, -n)
        mp.dps += 4
        r1 = pow(a, -mpf(1)/n)
        mp.dps -= 4
        assert r.ae(r1)
        assert mpc_ae(r, r1, eps)
    mp.dps = 15
    assert nthroot(4, 1) == 4
    assert nthroot(4, 0) == 1
    assert nthroot(4, -1) == 0.25
    assert nthroot(inf, 1) == inf
    assert nthroot(inf, 2) == inf
    assert nthroot(inf, 3) == inf
    assert nthroot(inf, -1) == 0
    assert nthroot(inf, -2) == 0
    assert nthroot(inf, -3) == 0
    assert nthroot(j, 1) == j
    assert nthroot(j, 0) == 1
    assert nthroot(j, -1) == -j
    assert isnan(nthroot(nan, 1))
    assert isnan(nthroot(nan, 0))
    assert isnan(nthroot(nan, -1))
    assert isnan(nthroot(inf, 0))
    assert root(2,3) == nthroot(2,3)
    assert root(16,4,0) == 2
    assert root(16,4,1) == 2j
    assert root(16,4,2) == -2
    assert root(16,4,3) == -2j
    assert root(16,4,4) == 2
    assert root(-125,3,1) == -5

def test_issue_136():
    for dps in [20, 80]:
        mp.dps = dps
        r = nthroot(mpf('-1e-20'), 4)
        assert r.ae(mpf(10)**(-5) * (1 + j) * mpf(2)**(-0.5))
    mp.dps = 80
    assert nthroot('-1e-3', 4).ae(mpf(10)**(-3./4) * (1 + j)/sqrt(2))
    assert nthroot('-1e-6', 4).ae((1 + j)/(10 * sqrt(20)))
    # Check that this doesn't take eternity to compute
    mp.dps = 20
    assert nthroot('-1e100000000', 4).ae((1+j)*mpf('1e25000000')/sqrt(2))
    mp.dps = 15

def test_mpcfun_real_imag():
    mp.dps = 15
    x = mpf(0.3)
    y = mpf(0.4)
    assert exp(mpc(x,0)) == exp(x)
    assert exp(mpc(0,y)) == mpc(cos(y),sin(y))
    assert cos(mpc(x,0)) == cos(x)
    assert sin(mpc(x,0)) == sin(x)
    assert cos(mpc(0,y)) == cosh(y)
    assert sin(mpc(0,y)) == mpc(0,sinh(y))
    assert cospi(mpc(x,0)) == cospi(x)
    assert sinpi(mpc(x,0)) == sinpi(x)
    assert cospi(mpc(0,y)).ae(cosh(pi*y))
    assert sinpi(mpc(0,y)).ae(mpc(0,sinh(pi*y)))
    c, s = cospi_sinpi(mpc(x,0))
    assert c == cospi(x)
    assert s == sinpi(x)
    c, s = cospi_sinpi(mpc(0,y))
    assert c.ae(cosh(pi*y))
    assert s.ae(mpc(0,sinh(pi*y)))
    c, s = cos_sin(mpc(x,0))
    assert c == cos(x)
    assert s == sin(x)
    c, s = cos_sin(mpc(0,y))
    assert c == cosh(y)
    assert s == mpc(0,sinh(y))

def test_perturbation_rounding():
    mp.dps = 100
    a = pi/10**50
    b = -pi/10**50
    c = 1 + a
    d = 1 + b
    mp.dps = 15
    assert exp(a) == 1
    assert exp(a, rounding='c') > 1
    assert exp(b, rounding='c') == 1
    assert exp(a, rounding='f') == 1
    assert exp(b, rounding='f') < 1
    assert cos(a) == 1
    assert cos(a, rounding='c') == 1
    assert cos(b, rounding='c') == 1
    assert cos(a, rounding='f') < 1
    assert cos(b, rounding='f') < 1
    for f in [sin, atan, asinh, tanh]:
        assert f(a) == +a
        assert f(a, rounding='c') > a
        assert f(a, rounding='f') < a
        assert f(b) == +b
        assert f(b, rounding='c') > b
        assert f(b, rounding='f') < b
    for f in [asin, tan, sinh, atanh]:
        assert f(a) == +a
        assert f(b) == +b
        assert f(a, rounding='c') > a
        assert f(b, rounding='c') > b
        assert f(a, rounding='f') < a
        assert f(b, rounding='f') < b
    assert ln(c) == +a
    assert ln(d) == +b
    assert ln(c, rounding='c') > a
    assert ln(c, rounding='f') < a
    assert ln(d, rounding='c') > b
    assert ln(d, rounding='f') < b
    assert cosh(a) == 1
    assert cosh(b) == 1
    assert cosh(a, rounding='c') > 1
    assert cosh(b, rounding='c') > 1
    assert cosh(a, rounding='f') == 1
    assert cosh(b, rounding='f') == 1

def test_integer_parts():
    assert floor(3.2) == 3
    assert ceil(3.2) == 4
    assert floor(3.2+5j) == 3+5j
    assert ceil(3.2+5j) == 4+5j

def test_complex_parts():
    assert fabs('3') == 3
    assert fabs(3+4j) == 5
    assert re(3) == 3
    assert re(1+4j) == 1
    assert im(3) == 0
    assert im(1+4j) == 4
    assert conj(3) == 3
    assert conj(3+4j) == 3-4j
    assert mpf(3).conjugate() == 3

def test_cospi_sinpi():
    assert sinpi(0) == 0
    assert sinpi(0.5) == 1
    assert sinpi(1) == 0
    assert sinpi(1.5) == -1
    assert sinpi(2) == 0
    assert sinpi(2.5) == 1
    assert sinpi(-0.5) == -1
    assert cospi(0) == 1
    assert cospi(0.5) == 0
    assert cospi(1) == -1
    assert cospi(1.5) == 0
    assert cospi(2) == 1
    assert cospi(2.5) == 0
    assert cospi(-0.5) == 0
    assert cospi(100000000000.25).ae(sqrt(2)/2)
    a = cospi(2+3j)
    assert a.real.ae(cos((2+3j)*pi).real)
    assert a.imag == 0
    b = sinpi(2+3j)
    assert b.imag.ae(sin((2+3j)*pi).imag)
    assert b.real == 0
    mp.dps = 35
    x1 = mpf(10000) - mpf('1e-15')
    x2 = mpf(10000) + mpf('1e-15')
    x3 = mpf(10000.5) - mpf('1e-15')
    x4 = mpf(10000.5) + mpf('1e-15')
    x5 = mpf(10001) - mpf('1e-15')
    x6 = mpf(10001) + mpf('1e-15')
    x7 = mpf(10001.5) - mpf('1e-15')
    x8 = mpf(10001.5) + mpf('1e-15')
    mp.dps = 15
    M = 10**15
    assert (sinpi(x1)*M).ae(-pi)
    assert (sinpi(x2)*M).ae(pi)
    assert (cospi(x3)*M).ae(pi)
    assert (cospi(x4)*M).ae(-pi)
    assert (sinpi(x5)*M).ae(pi)
    assert (sinpi(x6)*M).ae(-pi)
    assert (cospi(x7)*M).ae(-pi)
    assert (cospi(x8)*M).ae(pi)
    assert 0.999 < cospi(x1, rounding='d') < 1
    assert 0.999 < cospi(x2, rounding='d') < 1
    assert 0.999 < sinpi(x3, rounding='d') < 1
    assert 0.999 < sinpi(x4, rounding='d') < 1
    assert -1 < cospi(x5, rounding='d') < -0.999
    assert -1 < cospi(x6, rounding='d') < -0.999
    assert -1 < sinpi(x7, rounding='d') < -0.999
    assert -1 < sinpi(x8, rounding='d') < -0.999
    assert (sinpi(1e-15)*M).ae(pi)
    assert (sinpi(-1e-15)*M).ae(-pi)
    assert cospi(1e-15) == 1
    assert cospi(1e-15, rounding='d') < 1

def test_expj():
    assert expj(0) == 1
    assert expj(1).ae(exp(j))
    assert expj(j).ae(exp(-1))
    assert expj(1+j).ae(exp(j*(1+j)))
    assert expjpi(0) == 1
    assert expjpi(1).ae(exp(j*pi))
    assert expjpi(j).ae(exp(-pi))
    assert expjpi(1+j).ae(exp(j*pi*(1+j)))
    assert expjpi(-10**15 * j).ae('2.22579818340535731e+1364376353841841')

def test_sinc():
    assert sinc(0) == sincpi(0) == 1
    assert sinc(inf) == sincpi(inf) == 0
    assert sinc(-inf) == sincpi(-inf) == 0
    assert sinc(2).ae(0.45464871341284084770)
    assert sinc(2+3j).ae(0.4463290318402435457-2.7539470277436474940j)
    assert sincpi(2) == 0
    assert sincpi(1.5).ae(-0.212206590789193781)

def test_fibonacci():
    mp.dps = 15
    assert [fibonacci(n) for n in range(-5, 10)] == \
        [5, -3, 2, -1, 1, 0, 1, 1, 2, 3, 5, 8, 13, 21, 34]
    assert fib(2.5).ae(1.4893065462657091)
    assert fib(3+4j).ae(-5248.51130728372 - 14195.962288353j)
    assert fib(1000).ae(4.3466557686937455e+208)
    assert str(fib(10**100)) == '6.24499112864607e+2089876402499787337692720892375554168224592399182109535392875613974104853496745963277658556235103534'
    mp.dps = 2100
    a = fib(10000)
    assert a % 10**10 == 9947366875
    mp.dps = 15
    assert fibonacci(inf) == inf
    assert fib(3+0j) == 2

def test_call_with_dps():
    mp.dps = 15
    assert abs(exp(1, dps=30)-e(dps=35)) < 1e-29

def test_tanh():
    mp.dps = 15
    assert tanh(0) == 0
    assert tanh(inf) == 1
    assert tanh(-inf) == -1
    assert isnan(tanh(nan))
    assert tanh(mpc('inf', '0')) == 1

def test_atanh():
    mp.dps = 15
    assert atanh(0) == 0
    assert atanh(0.5).ae(0.54930614433405484570)
    assert atanh(-0.5).ae(-0.54930614433405484570)
    assert atanh(1) == inf
    assert atanh(-1) == -inf
    assert isnan(atanh(nan))
    assert isinstance(atanh(1), mpf)
    assert isinstance(atanh(-1), mpf)
    # Limits at infinity
    jpi2 = j*pi/2
    assert atanh(inf).ae(-jpi2)
    assert atanh(-inf).ae(jpi2)
    assert atanh(mpc(inf,-1)).ae(-jpi2)
    assert atanh(mpc(inf,0)).ae(-jpi2)
    assert atanh(mpc(inf,1)).ae(jpi2)
    assert atanh(mpc(1,inf)).ae(jpi2)
    assert atanh(mpc(0,inf)).ae(jpi2)
    assert atanh(mpc(-1,inf)).ae(jpi2)
    assert atanh(mpc(-inf,1)).ae(jpi2)
    assert atanh(mpc(-inf,0)).ae(jpi2)
    assert atanh(mpc(-inf,-1)).ae(-jpi2)
    assert atanh(mpc(-1,-inf)).ae(-jpi2)
    assert atanh(mpc(0,-inf)).ae(-jpi2)
    assert atanh(mpc(1,-inf)).ae(-jpi2)

def test_expm1():
    mp.dps = 15
    assert expm1(0) == 0
    assert expm1(3).ae(exp(3)-1)
    assert expm1(inf) == inf
    assert expm1(1e-50).ae(1e-50)
    assert (expm1(1e-10)*1e10).ae(1.00000000005)

def test_log1p():
    mp.dps = 15
    assert log1p(0) == 0
    assert log1p(3).ae(log(1+3))
    assert log1p(inf) == inf
    assert log1p(1e-50).ae(1e-50)
    assert (log1p(1e-10)*1e10).ae(0.99999999995)

def test_powm1():
    mp.dps = 15
    assert powm1(2,3) == 7
    assert powm1(-1,2) == 0
    assert powm1(-1,0) == 0
    assert powm1(-2,0) == 0
    assert powm1(3+4j,0) == 0
    assert powm1(0,1) == -1
    assert powm1(0,0) == 0
    assert powm1(1,0) == 0
    assert powm1(1,2) == 0
    assert powm1(1,3+4j) == 0
    assert powm1(1,5) == 0
    assert powm1(j,4) == 0
    assert powm1(-j,4) == 0
    assert (powm1(2,1e-100)*1e100).ae(ln2)
    assert powm1(2,'1e-100000000000') != 0
    assert (powm1(fadd(1,1e-100,exact=True), 5)*1e100).ae(5)

def test_unitroots():
    assert unitroots(1) == [1]
    assert unitroots(2) == [1, -1]
    a, b, c = unitroots(3)
    assert a == 1
    assert b.ae(-0.5 + 0.86602540378443864676j)
    assert c.ae(-0.5 - 0.86602540378443864676j)
    assert unitroots(1, primitive=True) == [1]
    assert unitroots(2, primitive=True) == [-1]
    assert unitroots(3, primitive=True) == unitroots(3)[1:]
    assert unitroots(4, primitive=True) == [j, -j]
    assert len(unitroots(17, primitive=True)) == 16
    assert len(unitroots(16, primitive=True)) == 8

def test_cyclotomic():
    mp.dps = 15
    assert [cyclotomic(n,1) for n in range(31)] == [1,0,2,3,2,5,1,7,2,3,1,11,1,13,1,1,2,17,1,19,1,1,1,23,1,5,1,3,1,29,1]
    assert [cyclotomic(n,-1) for n in range(31)] == [1,-2,0,1,2,1,3,1,2,1,5,1,1,1,7,1,2,1,3,1,1,1,11,1,1,1,13,1,1,1,1]
    assert [cyclotomic(n,j) for n in range(21)] == [1,-1+j,1+j,j,0,1,-j,j,2,-j,1,j,3,1,-j,1,2,1,j,j,5]
    assert [cyclotomic(n,-j) for n in range(21)] == [1,-1-j,1-j,-j,0,1,j,-j,2,j,1,-j,3,1,j,1,2,1,-j,-j,5]
    assert cyclotomic(1624,j) == 1
    assert cyclotomic(33600,j) == 1
    u = sqrt(j, prec=500)
    assert cyclotomic(8, u).ae(0)
    assert cyclotomic(30, u).ae(5.8284271247461900976)
    assert cyclotomic(2040, u).ae(1)
    assert cyclotomic(0,2.5) == 1
    assert cyclotomic(1,2.5) == 2.5-1
    assert cyclotomic(2,2.5) == 2.5+1
    assert cyclotomic(3,2.5) == 2.5**2 + 2.5 + 1
    assert cyclotomic(7,2.5) == 406.234375