123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905906907908909910911912913914915916917918919920 |
- 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
|