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