1from mpmath.libmp import *
2from mpmath import *
3import random
4import time
5import math
6import cmath
7
8def mpc_ae(a, b, eps=eps):
9    res = True
10    res = res and a.real.ae(b.real, eps)
11    res = res and a.imag.ae(b.imag, eps)
12    return res
13
14#----------------------------------------------------------------------------
15# Constants and functions
16#
17
18tpi = "3.1415926535897932384626433832795028841971693993751058209749445923078\
191640628620899862803482534211706798"
20te = "2.71828182845904523536028747135266249775724709369995957496696762772407\
21663035354759457138217852516642743"
22tdegree = "0.017453292519943295769236907684886127134428718885417254560971914\
234017100911460344944368224156963450948221"
24teuler = "0.5772156649015328606065120900824024310421593359399235988057672348\
2584867726777664670936947063291746749516"
26tln2 = "0.693147180559945309417232121458176568075500134360255254120680009493\
27393621969694715605863326996418687542"
28tln10 = "2.30258509299404568401799145468436420760110148862877297603332790096\
29757260967735248023599720508959829834"
30tcatalan = "0.91596559417721901505460351493238411077414937428167213426649811\
319621763019776254769479356512926115106249"
32tkhinchin = "2.6854520010653064453097148354817956938203822939944629530511523\
334555721885953715200280114117493184769800"
34tglaisher = "1.2824271291006226368753425688697917277676889273250011920637400\
352174040630885882646112973649195820237439420646"
36tapery = "1.2020569031595942853997381615114499907649862923404988817922715553\
374183820578631309018645587360933525815"
38tphi = "1.618033988749894848204586834365638117720309179805762862135448622705\
3926046281890244970720720418939113748475"
40tmertens = "0.26149721284764278375542683860869585905156664826119920619206421\
413924924510897368209714142631434246651052"
42ttwinprime = "0.660161815846869573927812110014555778432623360284733413319448\
43423335405642304495277143760031413839867912"
44
45def test_constants():
46    for prec in [3, 7, 10, 15, 20, 37, 80, 100, 29]:
47        mp.dps = prec
48        assert pi == mpf(tpi)
49        assert e == mpf(te)
50        assert degree == mpf(tdegree)
51        assert euler == mpf(teuler)
52        assert ln2 == mpf(tln2)
53        assert ln10 == mpf(tln10)
54        assert catalan == mpf(tcatalan)
55        assert khinchin == mpf(tkhinchin)
56        assert glaisher == mpf(tglaisher)
57        assert phi == mpf(tphi)
58        if prec < 50:
59            assert mertens == mpf(tmertens)
60            assert twinprime == mpf(ttwinprime)
61    mp.dps = 15
62    assert pi >= -1
63    assert pi > 2
64    assert pi > 3
65    assert pi < 4
66
67def test_exact_sqrts():
68    for i in range(20000):
69        assert sqrt(mpf(i*i)) == i
70    random.seed(1)
71    for prec in [100, 300, 1000, 10000]:
72        mp.dps = prec
73        for i in range(20):
74            A = random.randint(10**(prec//2-2), 10**(prec//2-1))
75            assert sqrt(mpf(A*A)) == A
76    mp.dps = 15
77    for i in range(100):
78        for a in [1, 8, 25, 112307]:
79            assert sqrt(mpf((a*a, 2*i))) == mpf((a, i))
80            assert sqrt(mpf((a*a, -2*i))) == mpf((a, -i))
81
82def test_sqrt_rounding():
83    for i in [2, 3, 5, 6, 7, 8, 10, 11, 12, 13, 14, 15]:
84        i = from_int(i)
85        for dps in [7, 15, 83, 106, 2000]:
86            mp.dps = dps
87            a = mpf_pow_int(mpf_sqrt(i, mp.prec, round_down), 2, mp.prec, round_down)
88            b = mpf_pow_int(mpf_sqrt(i, mp.prec, round_up), 2, mp.prec, round_up)
89            assert mpf_lt(a, i)
90            assert mpf_gt(b, i)
91    random.seed(1234)
92    prec = 100
93    for rnd in [round_down, round_nearest, round_ceiling]:
94        for i in range(100):
95            a = mpf_rand(prec)
96            b = mpf_mul(a, a)
97            assert mpf_sqrt(b, prec, rnd) == a
98    # Test some extreme cases
99    mp.dps = 100
100    a = mpf(9) + 1e-90
101    b = mpf(9) - 1e-90
102    mp.dps = 15
103    assert sqrt(a, rounding='d') == 3
104    assert sqrt(a, rounding='n') == 3
105    assert sqrt(a, rounding='u') > 3
106    assert sqrt(b, rounding='d') < 3
107    assert sqrt(b, rounding='n') == 3
108    assert sqrt(b, rounding='u') == 3
109    # A worst case, from the MPFR test suite
110    assert sqrt(mpf('7.0503726185518891')) == mpf('2.655253776675949')
111
112def test_float_sqrt():
113    mp.dps = 15
114    # These should round identically
115    for x in [0, 1e-7, 0.1, 0.5, 1, 2, 3, 4, 5, 0.333, 76.19]:
116        assert sqrt(mpf(x)) == float(x)**0.5
117    assert sqrt(-1) == 1j
118    assert sqrt(-2).ae(cmath.sqrt(-2))
119    assert sqrt(-3).ae(cmath.sqrt(-3))
120    assert sqrt(-100).ae(cmath.sqrt(-100))
121    assert sqrt(1j).ae(cmath.sqrt(1j))
122    assert sqrt(-1j).ae(cmath.sqrt(-1j))
123    assert sqrt(math.pi + math.e*1j).ae(cmath.sqrt(math.pi + math.e*1j))
124    assert sqrt(math.pi - math.e*1j).ae(cmath.sqrt(math.pi - math.e*1j))
125
126def test_hypot():
127    assert hypot(0, 0) == 0
128    assert hypot(0, 0.33) == mpf(0.33)
129    assert hypot(0.33, 0) == mpf(0.33)
130    assert hypot(-0.33, 0) == mpf(0.33)
131    assert hypot(3, 4) == mpf(5)
132
133def test_exact_cbrt():
134    for i in range(0, 20000, 200):
135        assert cbrt(mpf(i*i*i)) == i
136    random.seed(1)
137    for prec in [100, 300, 1000, 10000]:
138        mp.dps = prec
139        A = random.randint(10**(prec//2-2), 10**(prec//2-1))
140        assert cbrt(mpf(A*A*A)) == A
141    mp.dps = 15
142
143def test_exp():
144    assert exp(0) == 1
145    assert exp(10000).ae(mpf('8.8068182256629215873e4342'))
146    assert exp(-10000).ae(mpf('1.1354838653147360985e-4343'))
147    a = exp(mpf((1, 8198646019315405, -53, 53)))
148    assert(a.bc == bitcount(a.man))
149    mp.prec = 67
150    a = exp(mpf((1, 1781864658064754565, -60, 61)))
151    assert(a.bc == bitcount(a.man))
152    mp.prec = 53
153    assert exp(ln2 * 10).ae(1024)
154    assert exp(2+2j).ae(cmath.exp(2+2j))
155
156def test_issue_73():
157    mp.dps = 512
158    a = exp(-1)
159    b = exp(1)
160    mp.dps = 15
161    assert (+a).ae(0.36787944117144233)
162    assert (+b).ae(2.7182818284590451)
163
164def test_log():
165    mp.dps = 15
166    assert log(1) == 0
167    for x in [0.5, 1.5, 2.0, 3.0, 100, 10**50, 1e-50]:
168        assert log(x).ae(math.log(x))
169        assert log(x, x) == 1
170    assert log(1024, 2) == 10
171    assert log(10**1234, 10) == 1234
172    assert log(2+2j).ae(cmath.log(2+2j))
173    # Accuracy near 1
174    assert (log(0.6+0.8j).real*10**17).ae(2.2204460492503131)
175    assert (log(0.6-0.8j).real*10**17).ae(2.2204460492503131)
176    assert (log(0.8-0.6j).real*10**17).ae(2.2204460492503131)
177    assert (log(1+1e-8j).real*10**16).ae(0.5)
178    assert (log(1-1e-8j).real*10**16).ae(0.5)
179    assert (log(-1+1e-8j).real*10**16).ae(0.5)
180    assert (log(-1-1e-8j).real*10**16).ae(0.5)
181    assert (log(1j+1e-8).real*10**16).ae(0.5)
182    assert (log(1j-1e-8).real*10**16).ae(0.5)
183    assert (log(-1j+1e-8).real*10**16).ae(0.5)
184    assert (log(-1j-1e-8).real*10**16).ae(0.5)
185    assert (log(1+1e-40j).real*10**80).ae(0.5)
186    assert (log(1j+1e-40).real*10**80).ae(0.5)
187    # Huge
188    assert log(ldexp(1.234,10**20)).ae(log(2)*1e20)
189    assert log(ldexp(1.234,10**200)).ae(log(2)*1e200)
190    # Some special values
191    assert log(mpc(0,0)) == mpc(-inf,0)
192    assert isnan(log(mpc(nan,0)).real)
193    assert isnan(log(mpc(nan,0)).imag)
194    assert isnan(log(mpc(0,nan)).real)
195    assert isnan(log(mpc(0,nan)).imag)
196    assert isnan(log(mpc(nan,1)).real)
197    assert isnan(log(mpc(nan,1)).imag)
198    assert isnan(log(mpc(1,nan)).real)
199    assert isnan(log(mpc(1,nan)).imag)
200
201def test_trig_hyperb_basic():
202    for x in (list(range(100)) + list(range(-100,0))):
203        t = x / 4.1
204        assert cos(mpf(t)).ae(math.cos(t))
205        assert sin(mpf(t)).ae(math.sin(t))
206        assert tan(mpf(t)).ae(math.tan(t))
207        assert cosh(mpf(t)).ae(math.cosh(t))
208        assert sinh(mpf(t)).ae(math.sinh(t))
209        assert tanh(mpf(t)).ae(math.tanh(t))
210    assert sin(1+1j).ae(cmath.sin(1+1j))
211    assert sin(-4-3.6j).ae(cmath.sin(-4-3.6j))
212    assert cos(1+1j).ae(cmath.cos(1+1j))
213    assert cos(-4-3.6j).ae(cmath.cos(-4-3.6j))
214
215def test_degrees():
216    assert cos(0*degree) == 1
217    assert cos(90*degree).ae(0)
218    assert cos(180*degree).ae(-1)
219    assert cos(270*degree).ae(0)
220    assert cos(360*degree).ae(1)
221    assert sin(0*degree) == 0
222    assert sin(90*degree).ae(1)
223    assert sin(180*degree).ae(0)
224    assert sin(270*degree).ae(-1)
225    assert sin(360*degree).ae(0)
226
227def random_complexes(N):
228    random.seed(1)
229    a = []
230    for i in range(N):
231        x1 = random.uniform(-10, 10)
232        y1 = random.uniform(-10, 10)
233        x2 = random.uniform(-10, 10)
234        y2 = random.uniform(-10, 10)
235        z1 = complex(x1, y1)
236        z2 = complex(x2, y2)
237        a.append((z1, z2))
238    return a
239
240def test_complex_powers():
241    for dps in [15, 30, 100]:
242        # Check accuracy for complex square root
243        mp.dps = dps
244        a = mpc(1j)**0.5
245        assert a.real == a.imag == mpf(2)**0.5 / 2
246    mp.dps = 15
247    random.seed(1)
248    for (z1, z2) in random_complexes(100):
249        assert (mpc(z1)**mpc(z2)).ae(z1**z2, 1e-12)
250    assert (e**(-pi*1j)).ae(-1)
251    mp.dps = 50
252    assert (e**(-pi*1j)).ae(-1)
253    mp.dps = 15
254
255def test_complex_sqrt_accuracy():
256    def test_mpc_sqrt(lst):
257        for a, b in lst:
258            z = mpc(a + j*b)
259            assert mpc_ae(sqrt(z*z), z)
260            z = mpc(-a + j*b)
261            assert mpc_ae(sqrt(z*z), -z)
262            z = mpc(a - j*b)
263            assert mpc_ae(sqrt(z*z), z)
264            z = mpc(-a - j*b)
265            assert mpc_ae(sqrt(z*z), -z)
266    random.seed(2)
267    N = 10
268    mp.dps = 30
269    dps = mp.dps
270    test_mpc_sqrt([(random.uniform(0, 10),random.uniform(0, 10)) for i in range(N)])
271    test_mpc_sqrt([(i + 0.1, (i + 0.2)*10**i) for i in range(N)])
272    mp.dps = 15
273
274def test_atan():
275    mp.dps = 15
276    assert atan(-2.3).ae(math.atan(-2.3))
277    assert atan(1e-50) == 1e-50
278    assert atan(1e50).ae(pi/2)
279    assert atan(-1e-50) == -1e-50
280    assert atan(-1e50).ae(-pi/2)
281    assert atan(10**1000).ae(pi/2)
282    for dps in [25, 70, 100, 300, 1000]:
283        mp.dps = dps
284        assert (4*atan(1)).ae(pi)
285    mp.dps = 15
286    pi2 = pi/2
287    assert atan(mpc(inf,-1)).ae(pi2)
288    assert atan(mpc(inf,0)).ae(pi2)
289    assert atan(mpc(inf,1)).ae(pi2)
290    assert atan(mpc(1,inf)).ae(pi2)
291    assert atan(mpc(0,inf)).ae(pi2)
292    assert atan(mpc(-1,inf)).ae(-pi2)
293    assert atan(mpc(-inf,1)).ae(-pi2)
294    assert atan(mpc(-inf,0)).ae(-pi2)
295    assert atan(mpc(-inf,-1)).ae(-pi2)
296    assert atan(mpc(-1,-inf)).ae(-pi2)
297    assert atan(mpc(0,-inf)).ae(-pi2)
298    assert atan(mpc(1,-inf)).ae(pi2)
299
300def test_atan2():
301    mp.dps = 15
302    assert atan2(1,1).ae(pi/4)
303    assert atan2(1,-1).ae(3*pi/4)
304    assert atan2(-1,-1).ae(-3*pi/4)
305    assert atan2(-1,1).ae(-pi/4)
306    assert atan2(-1,0).ae(-pi/2)
307    assert atan2(1,0).ae(pi/2)
308    assert atan2(0,0) == 0
309    assert atan2(inf,0).ae(pi/2)
310    assert atan2(-inf,0).ae(-pi/2)
311    assert isnan(atan2(inf,inf))
312    assert isnan(atan2(-inf,inf))
313    assert isnan(atan2(inf,-inf))
314    assert isnan(atan2(3,nan))
315    assert isnan(atan2(nan,3))
316    assert isnan(atan2(0,nan))
317    assert isnan(atan2(nan,0))
318    assert atan2(0,inf) == 0
319    assert atan2(0,-inf).ae(pi)
320    assert atan2(10,inf) == 0
321    assert atan2(-10,inf) == 0
322    assert atan2(-10,-inf).ae(-pi)
323    assert atan2(10,-inf).ae(pi)
324    assert atan2(inf,10).ae(pi/2)
325    assert atan2(inf,-10).ae(pi/2)
326    assert atan2(-inf,10).ae(-pi/2)
327    assert atan2(-inf,-10).ae(-pi/2)
328
329def test_areal_inverses():
330    assert asin(mpf(0)) == 0
331    assert asinh(mpf(0)) == 0
332    assert acosh(mpf(1)) == 0
333    assert isinstance(asin(mpf(0.5)), mpf)
334    assert isinstance(asin(mpf(2.0)), mpc)
335    assert isinstance(acos(mpf(0.5)), mpf)
336    assert isinstance(acos(mpf(2.0)), mpc)
337    assert isinstance(atanh(mpf(0.1)), mpf)
338    assert isinstance(atanh(mpf(1.1)), mpc)
339
340    random.seed(1)
341    for i in range(50):
342        x = random.uniform(0, 1)
343        assert asin(mpf(x)).ae(math.asin(x))
344        assert acos(mpf(x)).ae(math.acos(x))
345
346        x = random.uniform(-10, 10)
347        assert asinh(mpf(x)).ae(cmath.asinh(x).real)
348        assert isinstance(asinh(mpf(x)), mpf)
349        x = random.uniform(1, 10)
350        assert acosh(mpf(x)).ae(cmath.acosh(x).real)
351        assert isinstance(acosh(mpf(x)), mpf)
352        x = random.uniform(-10, 0.999)
353        assert isinstance(acosh(mpf(x)), mpc)
354
355        x = random.uniform(-1, 1)
356        assert atanh(mpf(x)).ae(cmath.atanh(x).real)
357        assert isinstance(atanh(mpf(x)), mpf)
358
359    dps = mp.dps
360    mp.dps = 300
361    assert isinstance(asin(0.5), mpf)
362    mp.dps = 1000
363    assert asin(1).ae(pi/2)
364    assert asin(-1).ae(-pi/2)
365    mp.dps = dps
366
367def test_invhyperb_inaccuracy():
368    mp.dps = 15
369    assert (asinh(1e-5)*10**5).ae(0.99999999998333333)
370    assert (asinh(1e-10)*10**10).ae(1)
371    assert (asinh(1e-50)*10**50).ae(1)
372    assert (asinh(-1e-5)*10**5).ae(-0.99999999998333333)
373    assert (asinh(-1e-10)*10**10).ae(-1)
374    assert (asinh(-1e-50)*10**50).ae(-1)
375    assert asinh(10**20).ae(46.744849040440862)
376    assert asinh(-10**20).ae(-46.744849040440862)
377    assert (tanh(1e-10)*10**10).ae(1)
378    assert (tanh(-1e-10)*10**10).ae(-1)
379    assert (atanh(1e-10)*10**10).ae(1)
380    assert (atanh(-1e-10)*10**10).ae(-1)
381
382def test_complex_functions():
383    for x in (list(range(10)) + list(range(-10,0))):
384        for y in (list(range(10)) + list(range(-10,0))):
385            z = complex(x, y)/4.3 + 0.01j
386            assert exp(mpc(z)).ae(cmath.exp(z))
387            assert log(mpc(z)).ae(cmath.log(z))
388            assert cos(mpc(z)).ae(cmath.cos(z))
389            assert sin(mpc(z)).ae(cmath.sin(z))
390            assert tan(mpc(z)).ae(cmath.tan(z))
391            assert sinh(mpc(z)).ae(cmath.sinh(z))
392            assert cosh(mpc(z)).ae(cmath.cosh(z))
393            assert tanh(mpc(z)).ae(cmath.tanh(z))
394
395def test_complex_inverse_functions():
396    mp.dps = 15
397    iv.dps = 15
398    for (z1, z2) in random_complexes(30):
399        # apparently cmath uses a different branch, so we
400        # can't use it for comparison
401        assert sinh(asinh(z1)).ae(z1)
402        #
403        assert acosh(z1).ae(cmath.acosh(z1))
404        assert atanh(z1).ae(cmath.atanh(z1))
405        assert atan(z1).ae(cmath.atan(z1))
406        # the reason we set a big eps here is that the cmath
407        # functions are inaccurate
408        assert asin(z1).ae(cmath.asin(z1), rel_eps=1e-12)
409        assert acos(z1).ae(cmath.acos(z1), rel_eps=1e-12)
410        one = mpf(1)
411    for i in range(-9, 10, 3):
412        for k in range(-9, 10, 3):
413            a = 0.9*j*10**k + 0.8*one*10**i
414            b = cos(acos(a))
415            assert b.ae(a)
416            b = sin(asin(a))
417            assert b.ae(a)
418    one = mpf(1)
419    err = 2*10**-15
420    for i in range(-9, 9, 3):
421        for k in range(-9, 9, 3):
422            a = -0.9*10**k + j*0.8*one*10**i
423            b = cosh(acosh(a))
424            assert b.ae(a, err)
425            b = sinh(asinh(a))
426            assert b.ae(a, err)
427
428def test_reciprocal_functions():
429    assert sec(3).ae(-1.01010866590799375)
430    assert csc(3).ae(7.08616739573718592)
431    assert cot(3).ae(-7.01525255143453347)
432    assert sech(3).ae(0.0993279274194332078)
433    assert csch(3).ae(0.0998215696688227329)
434    assert coth(3).ae(1.00496982331368917)
435    assert asec(3).ae(1.23095941734077468)
436    assert acsc(3).ae(0.339836909454121937)
437    assert acot(3).ae(0.321750554396642193)
438    assert asech(0.5).ae(1.31695789692481671)
439    assert acsch(3).ae(0.327450150237258443)
440    assert acoth(3).ae(0.346573590279972655)
441    assert acot(0).ae(1.5707963267948966192)
442    assert acoth(0).ae(1.5707963267948966192j)
443
444def test_ldexp():
445    mp.dps = 15
446    assert ldexp(mpf(2.5), 0) == 2.5
447    assert ldexp(mpf(2.5), -1) == 1.25
448    assert ldexp(mpf(2.5), 2) == 10
449    assert ldexp(mpf('inf'), 3) == mpf('inf')
450
451def test_frexp():
452    mp.dps = 15
453    assert frexp(0) == (0.0, 0)
454    assert frexp(9) == (0.5625, 4)
455    assert frexp(1) == (0.5, 1)
456    assert frexp(0.2) == (0.8, -2)
457    assert frexp(1000) == (0.9765625, 10)
458
459def test_aliases():
460    assert ln(7) == log(7)
461    assert log10(3.75) == log(3.75,10)
462    assert degrees(5.6) == 5.6 / degree
463    assert radians(5.6) == 5.6 * degree
464    assert power(-1,0.5) == j
465    assert fmod(25,7) == 4.0 and isinstance(fmod(25,7), mpf)
466
467def test_arg_sign():
468    assert arg(3) == 0
469    assert arg(-3).ae(pi)
470    assert arg(j).ae(pi/2)
471    assert arg(-j).ae(-pi/2)
472    assert arg(0) == 0
473    assert isnan(atan2(3,nan))
474    assert isnan(atan2(nan,3))
475    assert isnan(atan2(0,nan))
476    assert isnan(atan2(nan,0))
477    assert isnan(atan2(nan,nan))
478    assert arg(inf) == 0
479    assert arg(-inf).ae(pi)
480    assert isnan(arg(nan))
481    #assert arg(inf*j).ae(pi/2)
482    assert sign(0) == 0
483    assert sign(3) == 1
484    assert sign(-3) == -1
485    assert sign(inf) == 1
486    assert sign(-inf) == -1
487    assert isnan(sign(nan))
488    assert sign(j) == j
489    assert sign(-3*j) == -j
490    assert sign(1+j).ae((1+j)/sqrt(2))
491
492def test_misc_bugs():
493    # test that this doesn't raise an exception
494    mp.dps = 1000
495    log(1302)
496    mp.dps = 15
497
498def test_arange():
499    assert arange(10) == [mpf('0.0'), mpf('1.0'), mpf('2.0'), mpf('3.0'),
500                          mpf('4.0'), mpf('5.0'), mpf('6.0'), mpf('7.0'),
501                          mpf('8.0'), mpf('9.0')]
502    assert arange(-5, 5) == [mpf('-5.0'), mpf('-4.0'), mpf('-3.0'),
503                             mpf('-2.0'), mpf('-1.0'), mpf('0.0'),
504                             mpf('1.0'), mpf('2.0'), mpf('3.0'), mpf('4.0')]
505    assert arange(0, 1, 0.1) == [mpf('0.0'), mpf('0.10000000000000001'),
506                                 mpf('0.20000000000000001'),
507                                 mpf('0.30000000000000004'),
508                                 mpf('0.40000000000000002'),
509                                 mpf('0.5'), mpf('0.60000000000000009'),
510                                 mpf('0.70000000000000007'),
511                                 mpf('0.80000000000000004'),
512                                 mpf('0.90000000000000002')]
513    assert arange(17, -9, -3) == [mpf('17.0'), mpf('14.0'), mpf('11.0'),
514                                  mpf('8.0'), mpf('5.0'), mpf('2.0'),
515                                  mpf('-1.0'), mpf('-4.0'), mpf('-7.0')]
516    assert arange(0.2, 0.1, -0.1) == [mpf('0.20000000000000001')]
517    assert arange(0) == []
518    assert arange(1000, -1) == []
519    assert arange(-1.23, 3.21, -0.0000001) == []
520
521def test_linspace():
522    assert linspace(2, 9, 7) == [mpf('2.0'), mpf('3.166666666666667'),
523        mpf('4.3333333333333339'), mpf('5.5'), mpf('6.666666666666667'),
524        mpf('7.8333333333333339'), mpf('9.0')]
525    assert linspace(2, 9, 7, endpoint=0) == [mpf('2.0'), mpf('3.0'), mpf('4.0'),
526        mpf('5.0'), mpf('6.0'), mpf('7.0'), mpf('8.0')]
527    assert linspace(2, 7, 1) == [mpf(2)]
528
529def test_float_cbrt():
530    mp.dps = 30
531    for a in arange(0,10,0.1):
532        assert cbrt(a*a*a).ae(a, eps)
533    assert cbrt(-1).ae(0.5 + j*sqrt(3)/2)
534    one_third = mpf(1)/3
535    for a in arange(0,10,2.7) + [0.1 + 10**5]:
536        a = mpc(a + 1.1j)
537        r1 = cbrt(a)
538        mp.dps += 10
539        r2 = pow(a, one_third)
540        mp.dps -= 10
541        assert r1.ae(r2, eps)
542    mp.dps = 100
543    for n in range(100, 301, 100):
544        w = 10**n + j*10**-3
545        z = w*w*w
546        r = cbrt(z)
547        assert mpc_ae(r, w, eps)
548    mp.dps = 15
549
550def test_root():
551    mp.dps = 30
552    random.seed(1)
553    a = random.randint(0, 10000)
554    p = a*a*a
555    r = nthroot(mpf(p), 3)
556    assert r == a
557    for n in range(4, 10):
558        p = p*a
559        assert nthroot(mpf(p), n) == a
560    mp.dps = 40
561    for n in range(10, 5000, 100):
562        for a in [random.random()*10000, random.random()*10**100]:
563            r = nthroot(a, n)
564            r1 = pow(a, mpf(1)/n)
565            assert r.ae(r1)
566            r = nthroot(a, -n)
567            r1 = pow(a, -mpf(1)/n)
568            assert r.ae(r1)
569    # XXX: this is broken right now
570    # tests for nthroot rounding
571    for rnd in ['nearest', 'up', 'down']:
572        mp.rounding = rnd
573        for n in [-5, -3, 3, 5]:
574            prec = 50
575            for i in range(10):
576                mp.prec = prec
577                a = rand()
578                mp.prec = 2*prec
579                b = a**n
580                mp.prec = prec
581                r = nthroot(b, n)
582                assert r == a
583    mp.dps = 30
584    for n in range(3, 21):
585        a = (random.random() + j*random.random())
586        assert nthroot(a, n).ae(pow(a, mpf(1)/n))
587        assert mpc_ae(nthroot(a, n), pow(a, mpf(1)/n))
588        a = (random.random()*10**100 + j*random.random())
589        r = nthroot(a, n)
590        mp.dps += 4
591        r1 = pow(a, mpf(1)/n)
592        mp.dps -= 4
593        assert r.ae(r1)
594        assert mpc_ae(r, r1, eps)
595        r = nthroot(a, -n)
596        mp.dps += 4
597        r1 = pow(a, -mpf(1)/n)
598        mp.dps -= 4
599        assert r.ae(r1)
600        assert mpc_ae(r, r1, eps)
601    mp.dps = 15
602    assert nthroot(4, 1) == 4
603    assert nthroot(4, 0) == 1
604    assert nthroot(4, -1) == 0.25
605    assert nthroot(inf, 1) == inf
606    assert nthroot(inf, 2) == inf
607    assert nthroot(inf, 3) == inf
608    assert nthroot(inf, -1) == 0
609    assert nthroot(inf, -2) == 0
610    assert nthroot(inf, -3) == 0
611    assert nthroot(j, 1) == j
612    assert nthroot(j, 0) == 1
613    assert nthroot(j, -1) == -j
614    assert isnan(nthroot(nan, 1))
615    assert isnan(nthroot(nan, 0))
616    assert isnan(nthroot(nan, -1))
617    assert isnan(nthroot(inf, 0))
618    assert root(2,3) == nthroot(2,3)
619    assert root(16,4,0) == 2
620    assert root(16,4,1) == 2j
621    assert root(16,4,2) == -2
622    assert root(16,4,3) == -2j
623    assert root(16,4,4) == 2
624    assert root(-125,3,1) == -5
625
626def test_issue_136():
627    for dps in [20, 80]:
628        mp.dps = dps
629        r = nthroot(mpf('-1e-20'), 4)
630        assert r.ae(mpf(10)**(-5) * (1 + j) * mpf(2)**(-0.5))
631    mp.dps = 80
632    assert nthroot('-1e-3', 4).ae(mpf(10)**(-3./4) * (1 + j)/sqrt(2))
633    assert nthroot('-1e-6', 4).ae((1 + j)/(10 * sqrt(20)))
634    # Check that this doesn't take eternity to compute
635    mp.dps = 20
636    assert nthroot('-1e100000000', 4).ae((1+j)*mpf('1e25000000')/sqrt(2))
637    mp.dps = 15
638
639def test_mpcfun_real_imag():
640    mp.dps = 15
641    x = mpf(0.3)
642    y = mpf(0.4)
643    assert exp(mpc(x,0)) == exp(x)
644    assert exp(mpc(0,y)) == mpc(cos(y),sin(y))
645    assert cos(mpc(x,0)) == cos(x)
646    assert sin(mpc(x,0)) == sin(x)
647    assert cos(mpc(0,y)) == cosh(y)
648    assert sin(mpc(0,y)) == mpc(0,sinh(y))
649    assert cospi(mpc(x,0)) == cospi(x)
650    assert sinpi(mpc(x,0)) == sinpi(x)
651    assert cospi(mpc(0,y)).ae(cosh(pi*y))
652    assert sinpi(mpc(0,y)).ae(mpc(0,sinh(pi*y)))
653    c, s = cospi_sinpi(mpc(x,0))
654    assert c == cospi(x)
655    assert s == sinpi(x)
656    c, s = cospi_sinpi(mpc(0,y))
657    assert c.ae(cosh(pi*y))
658    assert s.ae(mpc(0,sinh(pi*y)))
659    c, s = cos_sin(mpc(x,0))
660    assert c == cos(x)
661    assert s == sin(x)
662    c, s = cos_sin(mpc(0,y))
663    assert c == cosh(y)
664    assert s == mpc(0,sinh(y))
665
666def test_perturbation_rounding():
667    mp.dps = 100
668    a = pi/10**50
669    b = -pi/10**50
670    c = 1 + a
671    d = 1 + b
672    mp.dps = 15
673    assert exp(a) == 1
674    assert exp(a, rounding='c') > 1
675    assert exp(b, rounding='c') == 1
676    assert exp(a, rounding='f') == 1
677    assert exp(b, rounding='f') < 1
678    assert cos(a) == 1
679    assert cos(a, rounding='c') == 1
680    assert cos(b, rounding='c') == 1
681    assert cos(a, rounding='f') < 1
682    assert cos(b, rounding='f') < 1
683    for f in [sin, atan, asinh, tanh]:
684        assert f(a) == +a
685        assert f(a, rounding='c') > a
686        assert f(a, rounding='f') < a
687        assert f(b) == +b
688        assert f(b, rounding='c') > b
689        assert f(b, rounding='f') < b
690    for f in [asin, tan, sinh, atanh]:
691        assert f(a) == +a
692        assert f(b) == +b
693        assert f(a, rounding='c') > a
694        assert f(b, rounding='c') > b
695        assert f(a, rounding='f') < a
696        assert f(b, rounding='f') < b
697    assert ln(c) == +a
698    assert ln(d) == +b
699    assert ln(c, rounding='c') > a
700    assert ln(c, rounding='f') < a
701    assert ln(d, rounding='c') > b
702    assert ln(d, rounding='f') < b
703    assert cosh(a) == 1
704    assert cosh(b) == 1
705    assert cosh(a, rounding='c') > 1
706    assert cosh(b, rounding='c') > 1
707    assert cosh(a, rounding='f') == 1
708    assert cosh(b, rounding='f') == 1
709
710def test_integer_parts():
711    assert floor(3.2) == 3
712    assert ceil(3.2) == 4
713    assert floor(3.2+5j) == 3+5j
714    assert ceil(3.2+5j) == 4+5j
715
716def test_complex_parts():
717    assert fabs('3') == 3
718    assert fabs(3+4j) == 5
719    assert re(3) == 3
720    assert re(1+4j) == 1
721    assert im(3) == 0
722    assert im(1+4j) == 4
723    assert conj(3) == 3
724    assert conj(3+4j) == 3-4j
725    assert mpf(3).conjugate() == 3
726
727def test_cospi_sinpi():
728    assert sinpi(0) == 0
729    assert sinpi(0.5) == 1
730    assert sinpi(1) == 0
731    assert sinpi(1.5) == -1
732    assert sinpi(2) == 0
733    assert sinpi(2.5) == 1
734    assert sinpi(-0.5) == -1
735    assert cospi(0) == 1
736    assert cospi(0.5) == 0
737    assert cospi(1) == -1
738    assert cospi(1.5) == 0
739    assert cospi(2) == 1
740    assert cospi(2.5) == 0
741    assert cospi(-0.5) == 0
742    assert cospi(100000000000.25).ae(sqrt(2)/2)
743    a = cospi(2+3j)
744    assert a.real.ae(cos((2+3j)*pi).real)
745    assert a.imag == 0
746    b = sinpi(2+3j)
747    assert b.imag.ae(sin((2+3j)*pi).imag)
748    assert b.real == 0
749    mp.dps = 35
750    x1 = mpf(10000) - mpf('1e-15')
751    x2 = mpf(10000) + mpf('1e-15')
752    x3 = mpf(10000.5) - mpf('1e-15')
753    x4 = mpf(10000.5) + mpf('1e-15')
754    x5 = mpf(10001) - mpf('1e-15')
755    x6 = mpf(10001) + mpf('1e-15')
756    x7 = mpf(10001.5) - mpf('1e-15')
757    x8 = mpf(10001.5) + mpf('1e-15')
758    mp.dps = 15
759    M = 10**15
760    assert (sinpi(x1)*M).ae(-pi)
761    assert (sinpi(x2)*M).ae(pi)
762    assert (cospi(x3)*M).ae(pi)
763    assert (cospi(x4)*M).ae(-pi)
764    assert (sinpi(x5)*M).ae(pi)
765    assert (sinpi(x6)*M).ae(-pi)
766    assert (cospi(x7)*M).ae(-pi)
767    assert (cospi(x8)*M).ae(pi)
768    assert 0.999 < cospi(x1, rounding='d') < 1
769    assert 0.999 < cospi(x2, rounding='d') < 1
770    assert 0.999 < sinpi(x3, rounding='d') < 1
771    assert 0.999 < sinpi(x4, rounding='d') < 1
772    assert -1 < cospi(x5, rounding='d') < -0.999
773    assert -1 < cospi(x6, rounding='d') < -0.999
774    assert -1 < sinpi(x7, rounding='d') < -0.999
775    assert -1 < sinpi(x8, rounding='d') < -0.999
776    assert (sinpi(1e-15)*M).ae(pi)
777    assert (sinpi(-1e-15)*M).ae(-pi)
778    assert cospi(1e-15) == 1
779    assert cospi(1e-15, rounding='d') < 1
780
781def test_expj():
782    assert expj(0) == 1
783    assert expj(1).ae(exp(j))
784    assert expj(j).ae(exp(-1))
785    assert expj(1+j).ae(exp(j*(1+j)))
786    assert expjpi(0) == 1
787    assert expjpi(1).ae(exp(j*pi))
788    assert expjpi(j).ae(exp(-pi))
789    assert expjpi(1+j).ae(exp(j*pi*(1+j)))
790    assert expjpi(-10**15 * j).ae('2.22579818340535731e+1364376353841841')
791
792def test_sinc():
793    assert sinc(0) == sincpi(0) == 1
794    assert sinc(inf) == sincpi(inf) == 0
795    assert sinc(-inf) == sincpi(-inf) == 0
796    assert sinc(2).ae(0.45464871341284084770)
797    assert sinc(2+3j).ae(0.4463290318402435457-2.7539470277436474940j)
798    assert sincpi(2) == 0
799    assert sincpi(1.5).ae(-0.212206590789193781)
800
801def test_fibonacci():
802    mp.dps = 15
803    assert [fibonacci(n) for n in range(-5, 10)] == \
804        [5, -3, 2, -1, 1, 0, 1, 1, 2, 3, 5, 8, 13, 21, 34]
805    assert fib(2.5).ae(1.4893065462657091)
806    assert fib(3+4j).ae(-5248.51130728372 - 14195.962288353j)
807    assert fib(1000).ae(4.3466557686937455e+208)
808    assert str(fib(10**100)) == '6.24499112864607e+2089876402499787337692720892375554168224592399182109535392875613974104853496745963277658556235103534'
809    mp.dps = 2100
810    a = fib(10000)
811    assert a % 10**10 == 9947366875
812    mp.dps = 15
813    assert fibonacci(inf) == inf
814    assert fib(3+0j) == 2
815
816def test_call_with_dps():
817    mp.dps = 15
818    assert abs(exp(1, dps=30)-e(dps=35)) < 1e-29
819
820def test_tanh():
821    mp.dps = 15
822    assert tanh(0) == 0
823    assert tanh(inf) == 1
824    assert tanh(-inf) == -1
825    assert isnan(tanh(nan))
826    assert tanh(mpc('inf', '0')) == 1
827
828def test_atanh():
829    mp.dps = 15
830    assert atanh(0) == 0
831    assert atanh(0.5).ae(0.54930614433405484570)
832    assert atanh(-0.5).ae(-0.54930614433405484570)
833    assert atanh(1) == inf
834    assert atanh(-1) == -inf
835    assert isnan(atanh(nan))
836    assert isinstance(atanh(1), mpf)
837    assert isinstance(atanh(-1), mpf)
838    # Limits at infinity
839    jpi2 = j*pi/2
840    assert atanh(inf).ae(-jpi2)
841    assert atanh(-inf).ae(jpi2)
842    assert atanh(mpc(inf,-1)).ae(-jpi2)
843    assert atanh(mpc(inf,0)).ae(-jpi2)
844    assert atanh(mpc(inf,1)).ae(jpi2)
845    assert atanh(mpc(1,inf)).ae(jpi2)
846    assert atanh(mpc(0,inf)).ae(jpi2)
847    assert atanh(mpc(-1,inf)).ae(jpi2)
848    assert atanh(mpc(-inf,1)).ae(jpi2)
849    assert atanh(mpc(-inf,0)).ae(jpi2)
850    assert atanh(mpc(-inf,-1)).ae(-jpi2)
851    assert atanh(mpc(-1,-inf)).ae(-jpi2)
852    assert atanh(mpc(0,-inf)).ae(-jpi2)
853    assert atanh(mpc(1,-inf)).ae(-jpi2)
854
855def test_expm1():
856    mp.dps = 15
857    assert expm1(0) == 0
858    assert expm1(3).ae(exp(3)-1)
859    assert expm1(inf) == inf
860    assert expm1(1e-50).ae(1e-50)
861    assert (expm1(1e-10)*1e10).ae(1.00000000005)
862
863def test_log1p():
864    mp.dps = 15
865    assert log1p(0) == 0
866    assert log1p(3).ae(log(1+3))
867    assert log1p(inf) == inf
868    assert log1p(1e-50).ae(1e-50)
869    assert (log1p(1e-10)*1e10).ae(0.99999999995)
870
871def test_powm1():
872    mp.dps = 15
873    assert powm1(2,3) == 7
874    assert powm1(-1,2) == 0
875    assert powm1(-1,0) == 0
876    assert powm1(-2,0) == 0
877    assert powm1(3+4j,0) == 0
878    assert powm1(0,1) == -1
879    assert powm1(0,0) == 0
880    assert powm1(1,0) == 0
881    assert powm1(1,2) == 0
882    assert powm1(1,3+4j) == 0
883    assert powm1(1,5) == 0
884    assert powm1(j,4) == 0
885    assert powm1(-j,4) == 0
886    assert (powm1(2,1e-100)*1e100).ae(ln2)
887    assert powm1(2,'1e-100000000000') != 0
888    assert (powm1(fadd(1,1e-100,exact=True), 5)*1e100).ae(5)
889
890def test_unitroots():
891    assert unitroots(1) == [1]
892    assert unitroots(2) == [1, -1]
893    a, b, c = unitroots(3)
894    assert a == 1
895    assert b.ae(-0.5 + 0.86602540378443864676j)
896    assert c.ae(-0.5 - 0.86602540378443864676j)
897    assert unitroots(1, primitive=True) == [1]
898    assert unitroots(2, primitive=True) == [-1]
899    assert unitroots(3, primitive=True) == unitroots(3)[1:]
900    assert unitroots(4, primitive=True) == [j, -j]
901    assert len(unitroots(17, primitive=True)) == 16
902    assert len(unitroots(16, primitive=True)) == 8
903
904def test_cyclotomic():
905    mp.dps = 15
906    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]
907    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]
908    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]
909    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]
910    assert cyclotomic(1624,j) == 1
911    assert cyclotomic(33600,j) == 1
912    u = sqrt(j, prec=500)
913    assert cyclotomic(8, u).ae(0)
914    assert cyclotomic(30, u).ae(5.8284271247461900976)
915    assert cyclotomic(2040, u).ae(1)
916    assert cyclotomic(0,2.5) == 1
917    assert cyclotomic(1,2.5) == 2.5-1
918    assert cyclotomic(2,2.5) == 2.5+1
919    assert cyclotomic(3,2.5) == 2.5**2 + 2.5 + 1
920    assert cyclotomic(7,2.5) == 406.234375
921