test_functions.py 30 KB

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