text: 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
  