test_functions2.py 95 KB

12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879808182838485868788899091929394959697989910010110210310410510610710810911011111211311411511611711811912012112212312412512612712812913013113213313413513613713813914014114214314414514614714814915015115215315415515615715815916016116216316416516616716816917017117217317417517617717817918018118218318418518618718818919019119219319419519619719819920020120220320420520620720820921021121221321421521621721821922022122222322422522622722822923023123223323423523623723823924024124224324424524624724824925025125225325425525625725825926026126226326426526626726826927027127227327427527627727827928028128228328428528628728828929029129229329429529629729829930030130230330430530630730830931031131231331431531631731831932032132232332432532632732832933033133233333433533633733833934034134234334434534634734834935035135235335435535635735835936036136236336436536636736836937037137237337437537637737837938038138238338438538638738838939039139239339439539639739839940040140240340440540640740840941041141241341441541641741841942042142242342442542642742842943043143243343443543643743843944044144244344444544644744844945045145245345445545645745845946046146246346446546646746846947047147247347447547647747847948048148248348448548648748848949049149249349449549649749849950050150250350450550650750850951051151251351451551651751851952052152252352452552652752852953053153253353453553653753853954054154254354454554654754854955055155255355455555655755855956056156256356456556656756856957057157257357457557657757857958058158258358458558658758858959059159259359459559659759859960060160260360460560660760860961061161261361461561661761861962062162262362462562662762862963063163263363463563663763863964064164264364464564664764864965065165265365465565665765865966066166266366466566666766866967067167267367467567667767867968068168268368468568668768868969069169269369469569669769869970070170270370470570670770870971071171271371471571671771871972072172272372472572672772872973073173273373473573673773873974074174274374474574674774874975075175275375475575675775875976076176276376476576676776876977077177277377477577677777877978078178278378478578678778878979079179279379479579679779879980080180280380480580680780880981081181281381481581681781881982082182282382482582682782882983083183283383483583683783883984084184284384484584684784884985085185285385485585685785885986086186286386486586686786886987087187287387487587687787887988088188288388488588688788888989089189289389489589689789889990090190290390490590690790890991091191291391491591691791891992092192292392492592692792892993093193293393493593693793893994094194294394494594694794894995095195295395495595695795895996096196296396496596696796896997097197297397497597697797897998098198298398498598698798898999099199299399499599699799899910001001100210031004100510061007100810091010101110121013101410151016101710181019102010211022102310241025102610271028102910301031103210331034103510361037103810391040104110421043104410451046104710481049105010511052105310541055105610571058105910601061106210631064106510661067106810691070107110721073107410751076107710781079108010811082108310841085108610871088108910901091109210931094109510961097109810991100110111021103110411051106110711081109111011111112111311141115111611171118111911201121112211231124112511261127112811291130113111321133113411351136113711381139114011411142114311441145114611471148114911501151115211531154115511561157115811591160116111621163116411651166116711681169117011711172117311741175117611771178117911801181118211831184118511861187118811891190119111921193119411951196119711981199120012011202120312041205120612071208120912101211121212131214121512161217121812191220122112221223122412251226122712281229123012311232123312341235123612371238123912401241124212431244124512461247124812491250125112521253125412551256125712581259126012611262126312641265126612671268126912701271127212731274127512761277127812791280128112821283128412851286128712881289129012911292129312941295129612971298129913001301130213031304130513061307130813091310131113121313131413151316131713181319132013211322132313241325132613271328132913301331133213331334133513361337133813391340134113421343134413451346134713481349135013511352135313541355135613571358135913601361136213631364136513661367136813691370137113721373137413751376137713781379138013811382138313841385138613871388138913901391139213931394139513961397139813991400140114021403140414051406140714081409141014111412141314141415141614171418141914201421142214231424142514261427142814291430143114321433143414351436143714381439144014411442144314441445144614471448144914501451145214531454145514561457145814591460146114621463146414651466146714681469147014711472147314741475147614771478147914801481148214831484148514861487148814891490149114921493149414951496149714981499150015011502150315041505150615071508150915101511151215131514151515161517151815191520152115221523152415251526152715281529153015311532153315341535153615371538153915401541154215431544154515461547154815491550155115521553155415551556155715581559156015611562156315641565156615671568156915701571157215731574157515761577157815791580158115821583158415851586158715881589159015911592159315941595159615971598159916001601160216031604160516061607160816091610161116121613161416151616161716181619162016211622162316241625162616271628162916301631163216331634163516361637163816391640164116421643164416451646164716481649165016511652165316541655165616571658165916601661166216631664166516661667166816691670167116721673167416751676167716781679168016811682168316841685168616871688168916901691169216931694169516961697169816991700170117021703170417051706170717081709171017111712171317141715171617171718171917201721172217231724172517261727172817291730173117321733173417351736173717381739174017411742174317441745174617471748174917501751175217531754175517561757175817591760176117621763176417651766176717681769177017711772177317741775177617771778177917801781178217831784178517861787178817891790179117921793179417951796179717981799180018011802180318041805180618071808180918101811181218131814181518161817181818191820182118221823182418251826182718281829183018311832183318341835183618371838183918401841184218431844184518461847184818491850185118521853185418551856185718581859186018611862186318641865186618671868186918701871187218731874187518761877187818791880188118821883188418851886188718881889189018911892189318941895189618971898189919001901190219031904190519061907190819091910191119121913191419151916191719181919192019211922192319241925192619271928192919301931193219331934193519361937193819391940194119421943194419451946194719481949195019511952195319541955195619571958195919601961196219631964196519661967196819691970197119721973197419751976197719781979198019811982198319841985198619871988198919901991199219931994199519961997199819992000200120022003200420052006200720082009201020112012201320142015201620172018201920202021202220232024202520262027202820292030203120322033203420352036203720382039204020412042204320442045204620472048204920502051205220532054205520562057205820592060206120622063206420652066206720682069207020712072207320742075207620772078207920802081208220832084208520862087208820892090209120922093209420952096209720982099210021012102210321042105210621072108210921102111211221132114211521162117211821192120212121222123212421252126212721282129213021312132213321342135213621372138213921402141214221432144214521462147214821492150215121522153215421552156215721582159216021612162216321642165216621672168216921702171217221732174217521762177217821792180218121822183218421852186218721882189219021912192219321942195219621972198219922002201220222032204220522062207220822092210221122122213221422152216221722182219222022212222222322242225222622272228222922302231223222332234223522362237223822392240224122422243224422452246224722482249225022512252225322542255225622572258225922602261226222632264226522662267226822692270227122722273227422752276227722782279228022812282228322842285228622872288228922902291229222932294229522962297229822992300230123022303230423052306230723082309231023112312231323142315231623172318231923202321232223232324232523262327232823292330233123322333233423352336233723382339234023412342234323442345234623472348234923502351235223532354235523562357235823592360236123622363236423652366236723682369237023712372237323742375237623772378237923802381238223832384
  1. import math
  2. import pytest
  3. from mpmath import *
  4. def test_bessel():
  5. mp.dps = 15
  6. assert j0(1).ae(0.765197686557966551)
  7. assert j0(pi).ae(-0.304242177644093864)
  8. assert j0(1000).ae(0.0247866861524201746)
  9. assert j0(-25).ae(0.0962667832759581162)
  10. assert j1(1).ae(0.440050585744933516)
  11. assert j1(pi).ae(0.284615343179752757)
  12. assert j1(1000).ae(0.00472831190708952392)
  13. assert j1(-25).ae(0.125350249580289905)
  14. assert besselj(5,1).ae(0.000249757730211234431)
  15. assert besselj(5+0j,1).ae(0.000249757730211234431)
  16. assert besselj(5,pi).ae(0.0521411843671184747)
  17. assert besselj(5,1000).ae(0.00502540694523318607)
  18. assert besselj(5,-25).ae(0.0660079953984229934)
  19. assert besselj(-3,2).ae(-0.128943249474402051)
  20. assert besselj(-4,2).ae(0.0339957198075684341)
  21. assert besselj(3,3+2j).ae(0.424718794929639595942 + 0.625665327745785804812j)
  22. assert besselj(0.25,4).ae(-0.374760630804249715)
  23. assert besselj(1+2j,3+4j).ae(0.319247428741872131 - 0.669557748880365678j)
  24. assert (besselj(3, 10**10) * 10**5).ae(0.76765081748139204023)
  25. assert bessely(-0.5, 0) == 0
  26. assert bessely(0.5, 0) == -inf
  27. assert bessely(1.5, 0) == -inf
  28. assert bessely(0,0) == -inf
  29. assert bessely(-0.4, 0) == -inf
  30. assert bessely(-0.6, 0) == inf
  31. assert bessely(-1, 0) == inf
  32. assert bessely(-1.4, 0) == inf
  33. assert bessely(-1.6, 0) == -inf
  34. assert bessely(-1, 0) == inf
  35. assert bessely(-2, 0) == -inf
  36. assert bessely(-3, 0) == inf
  37. assert bessely(0.5, 0) == -inf
  38. assert bessely(1, 0) == -inf
  39. assert bessely(1.5, 0) == -inf
  40. assert bessely(2, 0) == -inf
  41. assert bessely(2.5, 0) == -inf
  42. assert bessely(3, 0) == -inf
  43. assert bessely(0,0.5).ae(-0.44451873350670655715)
  44. assert bessely(1,0.5).ae(-1.4714723926702430692)
  45. assert bessely(-1,0.5).ae(1.4714723926702430692)
  46. assert bessely(3.5,0.5).ae(-138.86400867242488443)
  47. assert bessely(0,3+4j).ae(4.6047596915010138655-8.8110771408232264208j)
  48. assert bessely(0,j).ae(-0.26803248203398854876+1.26606587775200833560j)
  49. assert (bessely(3, 10**10) * 10**5).ae(0.21755917537013204058)
  50. assert besseli(0,0) == 1
  51. assert besseli(1,0) == 0
  52. assert besseli(2,0) == 0
  53. assert besseli(-1,0) == 0
  54. assert besseli(-2,0) == 0
  55. assert besseli(0,0.5).ae(1.0634833707413235193)
  56. assert besseli(1,0.5).ae(0.25789430539089631636)
  57. assert besseli(-1,0.5).ae(0.25789430539089631636)
  58. assert besseli(3.5,0.5).ae(0.00068103597085793815863)
  59. assert besseli(0,3+4j).ae(-3.3924877882755196097-1.3239458916287264815j)
  60. assert besseli(0,j).ae(besselj(0,1))
  61. assert (besseli(3, 10**10) * mpf(10)**(-4342944813)).ae(4.2996028505491271875)
  62. assert besselk(0,0) == inf
  63. assert besselk(1,0) == inf
  64. assert besselk(2,0) == inf
  65. assert besselk(-1,0) == inf
  66. assert besselk(-2,0) == inf
  67. assert besselk(0,0.5).ae(0.92441907122766586178)
  68. assert besselk(1,0.5).ae(1.6564411200033008937)
  69. assert besselk(-1,0.5).ae(1.6564411200033008937)
  70. assert besselk(3.5,0.5).ae(207.48418747548460607)
  71. assert besselk(0,3+4j).ae(-0.007239051213570155013+0.026510418350267677215j)
  72. assert besselk(0,j).ae(-0.13863371520405399968-1.20196971531720649914j)
  73. assert (besselk(3, 10**10) * mpf(10)**4342944824).ae(1.1628981033356187851)
  74. # test for issue 331, bug reported by Michael Hartmann
  75. for n in range(10,100,10):
  76. mp.dps = n
  77. assert besseli(91.5,24.7708).ae("4.00830632138673963619656140653537080438462342928377020695738635559218797348548092636896796324190271316137982810144874264e-41")
  78. def test_bessel_zeros():
  79. mp.dps = 15
  80. assert besseljzero(0,1).ae(2.40482555769577276869)
  81. assert besseljzero(2,1).ae(5.1356223018406825563)
  82. assert besseljzero(1,50).ae(157.86265540193029781)
  83. assert besseljzero(10,1).ae(14.475500686554541220)
  84. assert besseljzero(0.5,3).ae(9.4247779607693797153)
  85. assert besseljzero(2,1,1).ae(3.0542369282271403228)
  86. assert besselyzero(0,1).ae(0.89357696627916752158)
  87. assert besselyzero(2,1).ae(3.3842417671495934727)
  88. assert besselyzero(1,50).ae(156.29183520147840108)
  89. assert besselyzero(10,1).ae(12.128927704415439387)
  90. assert besselyzero(0.5,3).ae(7.8539816339744830962)
  91. assert besselyzero(2,1,1).ae(5.0025829314460639452)
  92. def test_hankel():
  93. mp.dps = 15
  94. assert hankel1(0,0.5).ae(0.93846980724081290423-0.44451873350670655715j)
  95. assert hankel1(1,0.5).ae(0.2422684576748738864-1.4714723926702430692j)
  96. assert hankel1(-1,0.5).ae(-0.2422684576748738864+1.4714723926702430692j)
  97. assert hankel1(1.5,0.5).ae(0.0917016996256513026-2.5214655504213378514j)
  98. assert hankel1(1.5,3+4j).ae(0.0066806866476728165382-0.0036684231610839127106j)
  99. assert hankel2(0,0.5).ae(0.93846980724081290423+0.44451873350670655715j)
  100. assert hankel2(1,0.5).ae(0.2422684576748738864+1.4714723926702430692j)
  101. assert hankel2(-1,0.5).ae(-0.2422684576748738864-1.4714723926702430692j)
  102. assert hankel2(1.5,0.5).ae(0.0917016996256513026+2.5214655504213378514j)
  103. assert hankel2(1.5,3+4j).ae(14.783528526098567526-7.397390270853446512j)
  104. def test_struve():
  105. mp.dps = 15
  106. assert struveh(2,3).ae(0.74238666967748318564)
  107. assert struveh(-2.5,3).ae(0.41271003220971599344)
  108. assert struvel(2,3).ae(1.7476573277362782744)
  109. assert struvel(-2.5,3).ae(1.5153394466819651377)
  110. def test_whittaker():
  111. mp.dps = 15
  112. assert whitm(2,3,4).ae(49.753745589025246591)
  113. assert whitw(2,3,4).ae(14.111656223052932215)
  114. def test_kelvin():
  115. mp.dps = 15
  116. assert ber(2,3).ae(0.80836846563726819091)
  117. assert ber(3,4).ae(-0.28262680167242600233)
  118. assert ber(-3,2).ae(-0.085611448496796363669)
  119. assert bei(2,3).ae(-0.89102236377977331571)
  120. assert bei(-3,2).ae(-0.14420994155731828415)
  121. assert ker(2,3).ae(0.12839126695733458928)
  122. assert ker(-3,2).ae(-0.29802153400559142783)
  123. assert ker(0.5,3).ae(-0.085662378535217097524)
  124. assert kei(2,3).ae(0.036804426134164634000)
  125. assert kei(-3,2).ae(0.88682069845786731114)
  126. assert kei(0.5,3).ae(0.013633041571314302948)
  127. def test_hyper_misc():
  128. mp.dps = 15
  129. assert hyp0f1(1,0) == 1
  130. assert hyp1f1(1,2,0) == 1
  131. assert hyp1f2(1,2,3,0) == 1
  132. assert hyp2f1(1,2,3,0) == 1
  133. assert hyp2f2(1,2,3,4,0) == 1
  134. assert hyp2f3(1,2,3,4,5,0) == 1
  135. # Degenerate case: 0F0
  136. assert hyper([],[],0) == 1
  137. assert hyper([],[],-2).ae(exp(-2))
  138. # Degenerate case: 1F0
  139. assert hyper([2],[],1.5) == 4
  140. #
  141. assert hyp2f1((1,3),(2,3),(5,6),mpf(27)/32).ae(1.6)
  142. assert hyp2f1((1,4),(1,2),(3,4),mpf(80)/81).ae(1.8)
  143. assert hyp2f1((2,3),(1,1),(3,2),(2+j)/3).ae(1.327531603558679093+0.439585080092769253j)
  144. mp.dps = 25
  145. v = mpc('1.2282306665029814734863026', '-0.1225033830118305184672133')
  146. assert hyper([(3,4),2+j,1],[1,5,j/3],mpf(1)/5+j/8).ae(v)
  147. mp.dps = 15
  148. def test_elliptic_integrals():
  149. mp.dps = 15
  150. assert ellipk(0).ae(pi/2)
  151. assert ellipk(0.5).ae(gamma(0.25)**2/(4*sqrt(pi)))
  152. assert ellipk(1) == inf
  153. assert ellipk(1+0j) == inf
  154. assert ellipk(-1).ae('1.3110287771460599052')
  155. assert ellipk(-2).ae('1.1714200841467698589')
  156. assert isinstance(ellipk(-2), mpf)
  157. assert isinstance(ellipe(-2), mpf)
  158. assert ellipk(-50).ae('0.47103424540873331679')
  159. mp.dps = 30
  160. n1 = +fraction(99999,100000)
  161. n2 = +fraction(100001,100000)
  162. mp.dps = 15
  163. assert ellipk(n1).ae('7.1427724505817781901')
  164. assert ellipk(n2).ae(mpc('7.1427417367963090109', '-1.5707923998261688019'))
  165. assert ellipe(n1).ae('1.0000332138990829170')
  166. v = ellipe(n2)
  167. assert v.real.ae('0.999966786328145474069137')
  168. assert (v.imag*10**6).ae('7.853952181727432')
  169. assert ellipk(2).ae(mpc('1.3110287771460599052', '-1.3110287771460599052'))
  170. assert ellipk(50).ae(mpc('0.22326753950210985451', '-0.47434723226254522087'))
  171. assert ellipk(3+4j).ae(mpc('0.91119556380496500866', '0.63133428324134524388'))
  172. assert ellipk(3-4j).ae(mpc('0.91119556380496500866', '-0.63133428324134524388'))
  173. assert ellipk(-3+4j).ae(mpc('0.95357894880405122483', '0.23093044503746114444'))
  174. assert ellipk(-3-4j).ae(mpc('0.95357894880405122483', '-0.23093044503746114444'))
  175. assert isnan(ellipk(nan))
  176. assert isnan(ellipe(nan))
  177. assert ellipk(inf) == 0
  178. assert isinstance(ellipk(inf), mpc)
  179. assert ellipk(-inf) == 0
  180. assert ellipk(1+0j) == inf
  181. assert ellipe(0).ae(pi/2)
  182. assert ellipe(0.5).ae(pi**(mpf(3)/2)/gamma(0.25)**2 +gamma(0.25)**2/(8*sqrt(pi)))
  183. assert ellipe(1) == 1
  184. assert ellipe(1+0j) == 1
  185. assert ellipe(inf) == mpc(0,inf)
  186. assert ellipe(-inf) == inf
  187. assert ellipe(3+4j).ae(1.4995535209333469543-1.5778790079127582745j)
  188. assert ellipe(3-4j).ae(1.4995535209333469543+1.5778790079127582745j)
  189. assert ellipe(-3+4j).ae(2.5804237855343377803-0.8306096791000413778j)
  190. assert ellipe(-3-4j).ae(2.5804237855343377803+0.8306096791000413778j)
  191. assert ellipe(2).ae(0.59907011736779610372+0.59907011736779610372j)
  192. assert ellipe('1e-1000000000').ae(pi/2)
  193. assert ellipk('1e-1000000000').ae(pi/2)
  194. assert ellipe(-pi).ae(2.4535865983838923)
  195. mp.dps = 50
  196. assert ellipk(1/pi).ae('1.724756270009501831744438120951614673874904182624739673')
  197. assert ellipe(1/pi).ae('1.437129808135123030101542922290970050337425479058225712')
  198. assert ellipk(-10*pi).ae('0.5519067523886233967683646782286965823151896970015484512')
  199. assert ellipe(-10*pi).ae('5.926192483740483797854383268707108012328213431657645509')
  200. v = ellipk(pi)
  201. assert v.real.ae('0.973089521698042334840454592642137667227167622330325225')
  202. assert v.imag.ae('-1.156151296372835303836814390793087600271609993858798016')
  203. v = ellipe(pi)
  204. assert v.real.ae('0.4632848917264710404078033487934663562998345622611263332')
  205. assert v.imag.ae('1.0637961621753130852473300451583414489944099504180510966')
  206. mp.dps = 15
  207. def test_exp_integrals():
  208. mp.dps = 15
  209. x = +e
  210. z = e + sqrt(3)*j
  211. assert ei(x).ae(8.21168165538361560)
  212. assert li(x).ae(1.89511781635593676)
  213. assert si(x).ae(1.82104026914756705)
  214. assert ci(x).ae(0.213958001340379779)
  215. assert shi(x).ae(4.11520706247846193)
  216. assert chi(x).ae(4.09647459290515367)
  217. assert fresnels(x).ae(0.437189718149787643)
  218. assert fresnelc(x).ae(0.401777759590243012)
  219. assert airyai(x).ae(0.0108502401568586681)
  220. assert airybi(x).ae(8.98245748585468627)
  221. assert ei(z).ae(3.72597969491314951 + 7.34213212314224421j)
  222. assert li(z).ae(2.28662658112562502 + 1.50427225297269364j)
  223. assert si(z).ae(2.48122029237669054 + 0.12684703275254834j)
  224. assert ci(z).ae(0.169255590269456633 - 0.892020751420780353j)
  225. assert shi(z).ae(1.85810366559344468 + 3.66435842914920263j)
  226. assert chi(z).ae(1.86787602931970484 + 3.67777369399304159j)
  227. assert fresnels(z/3).ae(0.034534397197008182 + 0.754859844188218737j)
  228. assert fresnelc(z/3).ae(1.261581645990027372 + 0.417949198775061893j)
  229. assert airyai(z).ae(-0.0162552579839056062 - 0.0018045715700210556j)
  230. assert airybi(z).ae(-4.98856113282883371 + 2.08558537872180623j)
  231. assert li(0) == 0.0
  232. assert li(1) == -inf
  233. assert li(inf) == inf
  234. assert isinstance(li(0.7), mpf)
  235. assert si(inf).ae(pi/2)
  236. assert si(-inf).ae(-pi/2)
  237. assert ci(inf) == 0
  238. assert ci(0) == -inf
  239. assert isinstance(ei(-0.7), mpf)
  240. assert airyai(inf) == 0
  241. assert airybi(inf) == inf
  242. assert airyai(-inf) == 0
  243. assert airybi(-inf) == 0
  244. assert fresnels(inf) == 0.5
  245. assert fresnelc(inf) == 0.5
  246. assert fresnels(-inf) == -0.5
  247. assert fresnelc(-inf) == -0.5
  248. assert shi(0) == 0
  249. assert shi(inf) == inf
  250. assert shi(-inf) == -inf
  251. assert chi(0) == -inf
  252. assert chi(inf) == inf
  253. def test_ei():
  254. mp.dps = 15
  255. assert ei(0) == -inf
  256. assert ei(inf) == inf
  257. assert ei(-inf) == -0.0
  258. assert ei(20+70j).ae(6.1041351911152984397e6 - 2.7324109310519928872e6j)
  259. # tests for the asymptotic expansion
  260. # values checked with Mathematica ExpIntegralEi
  261. mp.dps = 50
  262. r = ei(20000)
  263. s = '3.8781962825045010930273870085501819470698476975019e+8681'
  264. assert str(r) == s
  265. r = ei(-200)
  266. s = '-6.8852261063076355977108174824557929738368086933303e-90'
  267. assert str(r) == s
  268. r =ei(20000 + 10*j)
  269. sre = '-3.255138234032069402493850638874410725961401274106e+8681'
  270. sim = '-2.1081929993474403520785942429469187647767369645423e+8681'
  271. assert str(r.real) == sre and str(r.imag) == sim
  272. mp.dps = 15
  273. # More asymptotic expansions
  274. assert chi(-10**6+100j).ae('1.3077239389562548386e+434288 + 7.6808956999707408158e+434287j')
  275. assert shi(-10**6+100j).ae('-1.3077239389562548386e+434288 - 7.6808956999707408158e+434287j')
  276. mp.dps = 15
  277. assert ei(10j).ae(-0.0454564330044553726+3.2291439210137706686j)
  278. assert ei(100j).ae(-0.0051488251426104921+3.1330217936839529126j)
  279. u = ei(fmul(10**20, j, exact=True))
  280. assert u.real.ae(-6.4525128526578084421345e-21, abs_eps=0, rel_eps=8*eps)
  281. assert u.imag.ae(pi)
  282. assert ei(-10j).ae(-0.0454564330044553726-3.2291439210137706686j)
  283. assert ei(-100j).ae(-0.0051488251426104921-3.1330217936839529126j)
  284. u = ei(fmul(-10**20, j, exact=True))
  285. assert u.real.ae(-6.4525128526578084421345e-21, abs_eps=0, rel_eps=8*eps)
  286. assert u.imag.ae(-pi)
  287. assert ei(10+10j).ae(-1576.1504265768517448+436.9192317011328140j)
  288. u = ei(-10+10j)
  289. assert u.real.ae(7.6698978415553488362543e-7, abs_eps=0, rel_eps=8*eps)
  290. assert u.imag.ae(3.141595611735621062025)
  291. def test_e1():
  292. mp.dps = 15
  293. assert e1(0) == inf
  294. assert e1(inf) == 0
  295. assert e1(-inf) == mpc(-inf, -pi)
  296. assert e1(10j).ae(0.045456433004455372635 + 0.087551267423977430100j)
  297. assert e1(100j).ae(0.0051488251426104921444 - 0.0085708599058403258790j)
  298. assert e1(fmul(10**20, j, exact=True)).ae(6.4525128526578084421e-21 - 7.6397040444172830039e-21j, abs_eps=0, rel_eps=8*eps)
  299. assert e1(-10j).ae(0.045456433004455372635 - 0.087551267423977430100j)
  300. assert e1(-100j).ae(0.0051488251426104921444 + 0.0085708599058403258790j)
  301. assert e1(fmul(-10**20, j, exact=True)).ae(6.4525128526578084421e-21 + 7.6397040444172830039e-21j, abs_eps=0, rel_eps=8*eps)
  302. def test_expint():
  303. mp.dps = 15
  304. assert expint(0,0) == inf
  305. assert expint(0,1).ae(1/e)
  306. assert expint(0,1.5).ae(2/exp(1.5)/3)
  307. assert expint(1,1).ae(-ei(-1))
  308. assert expint(2,0).ae(1)
  309. assert expint(3,0).ae(1/2.)
  310. assert expint(4,0).ae(1/3.)
  311. assert expint(-2, 0.5).ae(26/sqrt(e))
  312. assert expint(-1,-1) == 0
  313. assert expint(-2,-1).ae(-e)
  314. assert expint(5.5, 0).ae(2/9.)
  315. assert expint(2.00000001,0).ae(100000000./100000001)
  316. assert expint(2+3j,4-j).ae(0.0023461179581675065414+0.0020395540604713669262j)
  317. assert expint('1.01', '1e-1000').ae(99.9999999899412802)
  318. assert expint('1.000000000001', 3.5).ae(0.00697013985754701819446)
  319. assert expint(2,3).ae(3*ei(-3)+exp(-3))
  320. assert (expint(10,20)*10**10).ae(0.694439055541231353)
  321. assert expint(3,inf) == 0
  322. assert expint(3.2,inf) == 0
  323. assert expint(3.2+2j,inf) == 0
  324. assert expint(1,3j).ae(-0.11962978600800032763 + 0.27785620120457163717j)
  325. assert expint(1,3).ae(0.013048381094197037413)
  326. assert expint(1,-3).ae(-ei(3)-pi*j)
  327. #assert expint(3) == expint(1,3)
  328. assert expint(1,-20).ae(-25615652.66405658882 - 3.1415926535897932385j)
  329. assert expint(1000000,0).ae(1./999999)
  330. assert expint(0,2+3j).ae(-0.025019798357114678171 + 0.027980439405104419040j)
  331. assert expint(-1,2+3j).ae(-0.022411973626262070419 + 0.038058922011377716932j)
  332. assert expint(-1.5,0) == inf
  333. def test_trig_integrals():
  334. mp.dps = 30
  335. assert si(mpf(1)/1000000).ae('0.000000999999999999944444444444446111')
  336. assert ci(mpf(1)/1000000).ae('-13.2382948930629912435014366276')
  337. assert si(10**10).ae('1.5707963267075846569685111517747537')
  338. assert ci(10**10).ae('-4.87506025174822653785729773959e-11')
  339. assert si(10**100).ae(pi/2)
  340. assert (ci(10**100)*10**100).ae('-0.372376123661276688262086695553')
  341. assert si(-3) == -si(3)
  342. assert ci(-3).ae(ci(3) + pi*j)
  343. # Test complex structure
  344. mp.dps = 15
  345. assert mp.ci(50).ae(-0.0056283863241163054402)
  346. assert mp.ci(50+2j).ae(-0.018378282946133067149+0.070352808023688336193j)
  347. assert mp.ci(20j).ae(1.28078263320282943611e7+1.5707963267949j)
  348. assert mp.ci(-2+20j).ae(-4.050116856873293505e6+1.207476188206989909e7j)
  349. assert mp.ci(-50+2j).ae(-0.0183782829461330671+3.0712398455661049023j)
  350. assert mp.ci(-50).ae(-0.0056283863241163054+3.1415926535897932385j)
  351. assert mp.ci(-50-2j).ae(-0.0183782829461330671-3.0712398455661049023j)
  352. assert mp.ci(-2-20j).ae(-4.050116856873293505e6-1.207476188206989909e7j)
  353. assert mp.ci(-20j).ae(1.28078263320282943611e7-1.5707963267949j)
  354. assert mp.ci(50-2j).ae(-0.018378282946133067149-0.070352808023688336193j)
  355. assert mp.si(50).ae(1.5516170724859358947)
  356. assert mp.si(50+2j).ae(1.497884414277228461-0.017515007378437448j)
  357. assert mp.si(20j).ae(1.2807826332028294459e7j)
  358. assert mp.si(-2+20j).ae(-1.20747603112735722103e7-4.050116856873293554e6j)
  359. assert mp.si(-50+2j).ae(-1.497884414277228461-0.017515007378437448j)
  360. assert mp.si(-50).ae(-1.5516170724859358947)
  361. assert mp.si(-50-2j).ae(-1.497884414277228461+0.017515007378437448j)
  362. assert mp.si(-2-20j).ae(-1.20747603112735722103e7+4.050116856873293554e6j)
  363. assert mp.si(-20j).ae(-1.2807826332028294459e7j)
  364. assert mp.si(50-2j).ae(1.497884414277228461+0.017515007378437448j)
  365. assert mp.chi(50j).ae(-0.0056283863241163054+1.5707963267948966192j)
  366. assert mp.chi(-2+50j).ae(-0.0183782829461330671+1.6411491348185849554j)
  367. assert mp.chi(-20).ae(1.28078263320282943611e7+3.1415926535898j)
  368. assert mp.chi(-20-2j).ae(-4.050116856873293505e6+1.20747571696809187053e7j)
  369. assert mp.chi(-2-50j).ae(-0.0183782829461330671-1.6411491348185849554j)
  370. assert mp.chi(-50j).ae(-0.0056283863241163054-1.5707963267948966192j)
  371. assert mp.chi(2-50j).ae(-0.0183782829461330671-1.500443518771208283j)
  372. assert mp.chi(20-2j).ae(-4.050116856873293505e6-1.20747603112735722951e7j)
  373. assert mp.chi(20).ae(1.2807826332028294361e7)
  374. assert mp.chi(2+50j).ae(-0.0183782829461330671+1.500443518771208283j)
  375. assert mp.shi(50j).ae(1.5516170724859358947j)
  376. assert mp.shi(-2+50j).ae(0.017515007378437448+1.497884414277228461j)
  377. assert mp.shi(-20).ae(-1.2807826332028294459e7)
  378. assert mp.shi(-20-2j).ae(4.050116856873293554e6-1.20747603112735722103e7j)
  379. assert mp.shi(-2-50j).ae(0.017515007378437448-1.497884414277228461j)
  380. assert mp.shi(-50j).ae(-1.5516170724859358947j)
  381. assert mp.shi(2-50j).ae(-0.017515007378437448-1.497884414277228461j)
  382. assert mp.shi(20-2j).ae(-4.050116856873293554e6-1.20747603112735722103e7j)
  383. assert mp.shi(20).ae(1.2807826332028294459e7)
  384. assert mp.shi(2+50j).ae(-0.017515007378437448+1.497884414277228461j)
  385. def ae(x,y,tol=1e-12):
  386. return abs(x-y) <= abs(y)*tol
  387. assert fp.ci(fp.inf) == 0
  388. assert ae(fp.ci(fp.ninf), fp.pi*1j)
  389. assert ae(fp.si(fp.inf), fp.pi/2)
  390. assert ae(fp.si(fp.ninf), -fp.pi/2)
  391. assert fp.si(0) == 0
  392. assert ae(fp.ci(50), -0.0056283863241163054402)
  393. assert ae(fp.ci(50+2j), -0.018378282946133067149+0.070352808023688336193j)
  394. assert ae(fp.ci(20j), 1.28078263320282943611e7+1.5707963267949j)
  395. assert ae(fp.ci(-2+20j), -4.050116856873293505e6+1.207476188206989909e7j)
  396. assert ae(fp.ci(-50+2j), -0.0183782829461330671+3.0712398455661049023j)
  397. assert ae(fp.ci(-50), -0.0056283863241163054+3.1415926535897932385j)
  398. assert ae(fp.ci(-50-2j), -0.0183782829461330671-3.0712398455661049023j)
  399. assert ae(fp.ci(-2-20j), -4.050116856873293505e6-1.207476188206989909e7j)
  400. assert ae(fp.ci(-20j), 1.28078263320282943611e7-1.5707963267949j)
  401. assert ae(fp.ci(50-2j), -0.018378282946133067149-0.070352808023688336193j)
  402. assert ae(fp.si(50), 1.5516170724859358947)
  403. assert ae(fp.si(50+2j), 1.497884414277228461-0.017515007378437448j)
  404. assert ae(fp.si(20j), 1.2807826332028294459e7j)
  405. assert ae(fp.si(-2+20j), -1.20747603112735722103e7-4.050116856873293554e6j)
  406. assert ae(fp.si(-50+2j), -1.497884414277228461-0.017515007378437448j)
  407. assert ae(fp.si(-50), -1.5516170724859358947)
  408. assert ae(fp.si(-50-2j), -1.497884414277228461+0.017515007378437448j)
  409. assert ae(fp.si(-2-20j), -1.20747603112735722103e7+4.050116856873293554e6j)
  410. assert ae(fp.si(-20j), -1.2807826332028294459e7j)
  411. assert ae(fp.si(50-2j), 1.497884414277228461+0.017515007378437448j)
  412. assert ae(fp.chi(50j), -0.0056283863241163054+1.5707963267948966192j)
  413. assert ae(fp.chi(-2+50j), -0.0183782829461330671+1.6411491348185849554j)
  414. assert ae(fp.chi(-20), 1.28078263320282943611e7+3.1415926535898j)
  415. assert ae(fp.chi(-20-2j), -4.050116856873293505e6+1.20747571696809187053e7j)
  416. assert ae(fp.chi(-2-50j), -0.0183782829461330671-1.6411491348185849554j)
  417. assert ae(fp.chi(-50j), -0.0056283863241163054-1.5707963267948966192j)
  418. assert ae(fp.chi(2-50j), -0.0183782829461330671-1.500443518771208283j)
  419. assert ae(fp.chi(20-2j), -4.050116856873293505e6-1.20747603112735722951e7j)
  420. assert ae(fp.chi(20), 1.2807826332028294361e7)
  421. assert ae(fp.chi(2+50j), -0.0183782829461330671+1.500443518771208283j)
  422. assert ae(fp.shi(50j), 1.5516170724859358947j)
  423. assert ae(fp.shi(-2+50j), 0.017515007378437448+1.497884414277228461j)
  424. assert ae(fp.shi(-20), -1.2807826332028294459e7)
  425. assert ae(fp.shi(-20-2j), 4.050116856873293554e6-1.20747603112735722103e7j)
  426. assert ae(fp.shi(-2-50j), 0.017515007378437448-1.497884414277228461j)
  427. assert ae(fp.shi(-50j), -1.5516170724859358947j)
  428. assert ae(fp.shi(2-50j), -0.017515007378437448-1.497884414277228461j)
  429. assert ae(fp.shi(20-2j), -4.050116856873293554e6-1.20747603112735722103e7j)
  430. assert ae(fp.shi(20), 1.2807826332028294459e7)
  431. assert ae(fp.shi(2+50j), -0.017515007378437448+1.497884414277228461j)
  432. def test_airy():
  433. mp.dps = 15
  434. assert (airyai(10)*10**10).ae(1.1047532552898687)
  435. assert (airybi(10)/10**9).ae(0.45564115354822515)
  436. assert (airyai(1000)*10**9158).ae(9.306933063179556004)
  437. assert (airybi(1000)/10**9154).ae(5.4077118391949465477)
  438. assert airyai(-1000).ae(0.055971895773019918842)
  439. assert airybi(-1000).ae(-0.083264574117080633012)
  440. assert (airyai(100+100j)*10**188).ae(2.9099582462207032076 + 2.353013591706178756j)
  441. assert (airybi(100+100j)/10**185).ae(1.7086751714463652039 - 3.1416590020830804578j)
  442. def test_hyper_0f1():
  443. mp.dps = 15
  444. v = 8.63911136507950465
  445. assert hyper([],[(1,3)],1.5).ae(v)
  446. assert hyper([],[1/3.],1.5).ae(v)
  447. assert hyp0f1(1/3.,1.5).ae(v)
  448. assert hyp0f1((1,3),1.5).ae(v)
  449. # Asymptotic expansion
  450. assert hyp0f1(3,1e9).ae('4.9679055380347771271e+27455')
  451. assert hyp0f1(3,1e9j).ae('-2.1222788784457702157e+19410 + 5.0840597555401854116e+19410j')
  452. def test_hyper_1f1():
  453. mp.dps = 15
  454. v = 1.2917526488617656673
  455. assert hyper([(1,2)],[(3,2)],0.7).ae(v)
  456. assert hyper([(1,2)],[(3,2)],0.7+0j).ae(v)
  457. assert hyper([0.5],[(3,2)],0.7).ae(v)
  458. assert hyper([0.5],[1.5],0.7).ae(v)
  459. assert hyper([0.5],[(3,2)],0.7+0j).ae(v)
  460. assert hyper([0.5],[1.5],0.7+0j).ae(v)
  461. assert hyper([(1,2)],[1.5+0j],0.7).ae(v)
  462. assert hyper([0.5+0j],[1.5],0.7).ae(v)
  463. assert hyper([0.5+0j],[1.5+0j],0.7+0j).ae(v)
  464. assert hyp1f1(0.5,1.5,0.7).ae(v)
  465. assert hyp1f1((1,2),1.5,0.7).ae(v)
  466. # Asymptotic expansion
  467. assert hyp1f1(2,3,1e10).ae('2.1555012157015796988e+4342944809')
  468. assert (hyp1f1(2,3,1e10j)*10**10).ae(-0.97501205020039745852 - 1.7462392454512132074j)
  469. # Shouldn't use asymptotic expansion
  470. assert hyp1f1(-2, 1, 10000).ae(49980001)
  471. # Bug
  472. assert hyp1f1(1j,fraction(1,3),0.415-69.739j).ae(25.857588206024346592 + 15.738060264515292063j)
  473. def test_hyper_2f1():
  474. mp.dps = 15
  475. v = 1.0652207633823291032
  476. assert hyper([(1,2), (3,4)], [2], 0.3).ae(v)
  477. assert hyper([(1,2), 0.75], [2], 0.3).ae(v)
  478. assert hyper([0.5, 0.75], [2.0], 0.3).ae(v)
  479. assert hyper([0.5, 0.75], [2.0], 0.3+0j).ae(v)
  480. assert hyper([0.5+0j, (3,4)], [2.0], 0.3+0j).ae(v)
  481. assert hyper([0.5+0j, (3,4)], [2.0], 0.3).ae(v)
  482. assert hyper([0.5, (3,4)], [2.0+0j], 0.3).ae(v)
  483. assert hyper([0.5+0j, 0.75+0j], [2.0+0j], 0.3+0j).ae(v)
  484. v = 1.09234681096223231717 + 0.18104859169479360380j
  485. assert hyper([(1,2),0.75+j], [2], 0.5).ae(v)
  486. assert hyper([0.5,0.75+j], [2.0], 0.5).ae(v)
  487. assert hyper([0.5,0.75+j], [2.0], 0.5+0j).ae(v)
  488. assert hyper([0.5,0.75+j], [2.0+0j], 0.5+0j).ae(v)
  489. v = 0.9625 - 0.125j
  490. assert hyper([(3,2),-1],[4], 0.1+j/3).ae(v)
  491. assert hyper([1.5,-1.0],[4], 0.1+j/3).ae(v)
  492. assert hyper([1.5,-1.0],[4+0j], 0.1+j/3).ae(v)
  493. assert hyper([1.5+0j,-1.0+0j],[4+0j], 0.1+j/3).ae(v)
  494. v = 1.02111069501693445001 - 0.50402252613466859521j
  495. assert hyper([(2,10),(3,10)],[(4,10)],1.5).ae(v)
  496. assert hyper([0.2,(3,10)],[0.4+0j],1.5).ae(v)
  497. assert hyper([0.2,(3,10)],[0.4+0j],1.5+0j).ae(v)
  498. v = 0.76922501362865848528 + 0.32640579593235886194j
  499. assert hyper([(2,10),(3,10)],[(4,10)],4+2j).ae(v)
  500. assert hyper([0.2,(3,10)],[0.4+0j],4+2j).ae(v)
  501. assert hyper([0.2,(3,10)],[(4,10)],4+2j).ae(v)
  502. def test_hyper_2f1_hard():
  503. mp.dps = 15
  504. # Singular cases
  505. assert hyp2f1(2,-1,-1,3).ae(7)
  506. assert hyp2f1(2,-1,-1,3,eliminate_all=True).ae(0.25)
  507. assert hyp2f1(2,-2,-2,3).ae(34)
  508. assert hyp2f1(2,-2,-2,3,eliminate_all=True).ae(0.25)
  509. assert hyp2f1(2,-2,-3,3) == 14
  510. assert hyp2f1(2,-3,-2,3) == inf
  511. assert hyp2f1(2,-1.5,-1.5,3) == 0.25
  512. assert hyp2f1(1,2,3,0) == 1
  513. assert hyp2f1(0,1,0,0) == 1
  514. assert hyp2f1(0,0,0,0) == 1
  515. assert isnan(hyp2f1(1,1,0,0))
  516. assert hyp2f1(2,-1,-5, 0.25+0.25j).ae(1.1+0.1j)
  517. assert hyp2f1(2,-5,-5, 0.25+0.25j, eliminate=False).ae(163./128 + 125./128*j)
  518. assert hyp2f1(0.7235, -1, -5, 0.3).ae(1.04341)
  519. assert hyp2f1(0.7235, -5, -5, 0.3, eliminate=False).ae(1.2939225017815903812)
  520. assert hyp2f1(-1,-2,4,1) == 1.5
  521. assert hyp2f1(1,2,-3,1) == inf
  522. assert hyp2f1(-2,-2,1,1) == 6
  523. assert hyp2f1(1,-2,-4,1).ae(5./3)
  524. assert hyp2f1(0,-6,-4,1) == 1
  525. assert hyp2f1(0,-3,-4,1) == 1
  526. assert hyp2f1(0,0,0,1) == 1
  527. assert hyp2f1(1,0,0,1,eliminate=False) == 1
  528. assert hyp2f1(1,1,0,1) == inf
  529. assert hyp2f1(1,-6,-4,1) == inf
  530. assert hyp2f1(-7.2,-0.5,-4.5,1) == 0
  531. assert hyp2f1(-7.2,-1,-2,1).ae(-2.6)
  532. assert hyp2f1(1,-0.5,-4.5, 1) == inf
  533. assert hyp2f1(1,0.5,-4.5, 1) == -inf
  534. # Check evaluation on / close to unit circle
  535. z = exp(j*pi/3)
  536. w = (nthroot(2,3)+1)*exp(j*pi/12)/nthroot(3,4)**3
  537. assert hyp2f1('1/2','1/6','1/3', z).ae(w)
  538. assert hyp2f1('1/2','1/6','1/3', z.conjugate()).ae(w.conjugate())
  539. assert hyp2f1(0.25, (1,3), 2, '0.999').ae(1.06826449496030635)
  540. assert hyp2f1(0.25, (1,3), 2, '1.001').ae(1.06867299254830309446-0.00001446586793975874j)
  541. assert hyp2f1(0.25, (1,3), 2, -1).ae(0.96656584492524351673)
  542. assert hyp2f1(0.25, (1,3), 2, j).ae(0.99041766248982072266+0.03777135604180735522j)
  543. assert hyp2f1(2,3,5,'0.99').ae(27.699347904322690602)
  544. assert hyp2f1((3,2),-0.5,3,'0.99').ae(0.68403036843911661388)
  545. assert hyp2f1(2,3,5,1j).ae(0.37290667145974386127+0.59210004902748285917j)
  546. assert fsum([hyp2f1((7,10),(2,3),(-1,2), 0.95*exp(j*k)) for k in range(1,15)]).ae(52.851400204289452922+6.244285013912953225j)
  547. assert fsum([hyp2f1((7,10),(2,3),(-1,2), 1.05*exp(j*k)) for k in range(1,15)]).ae(54.506013786220655330-3.000118813413217097j)
  548. assert fsum([hyp2f1((7,10),(2,3),(-1,2), exp(j*k)) for k in range(1,15)]).ae(55.792077935955314887+1.731986485778500241j)
  549. assert hyp2f1(2,2.5,-3.25,0.999).ae(218373932801217082543180041.33)
  550. # Branches
  551. assert hyp2f1(1,1,2,1.01).ae(4.5595744415723676911-3.1104877758314784539j)
  552. assert hyp2f1(1,1,2,1.01+0.1j).ae(2.4149427480552782484+1.4148224796836938829j)
  553. assert hyp2f1(1,1,2,3+4j).ae(0.14576709331407297807+0.48379185417980360773j)
  554. assert hyp2f1(1,1,2,4).ae(-0.27465307216702742285 - 0.78539816339744830962j)
  555. assert hyp2f1(1,1,2,-4).ae(0.40235947810852509365)
  556. # Other:
  557. # Cancellation with a large parameter involved (bug reported on sage-devel)
  558. assert hyp2f1(112, (51,10), (-9,10), -0.99999).ae(-1.6241361047970862961e-24, abs_eps=0, rel_eps=eps*16)
  559. def test_hyper_3f2_etc():
  560. assert hyper([1,2,3],[1.5,8],-1).ae(0.67108992351533333030)
  561. assert hyper([1,2,3,4],[5,6,7], -1).ae(0.90232988035425506008)
  562. assert hyper([1,2,3],[1.25,5], 1).ae(28.924181329701905701)
  563. assert hyper([1,2,3,4],[5,6,7],5).ae(1.5192307344006649499-1.1529845225075537461j)
  564. assert hyper([1,2,3,4,5],[6,7,8,9],-1).ae(0.96288759462882357253)
  565. assert hyper([1,2,3,4,5],[6,7,8,9],1).ae(1.0428697385885855841)
  566. assert hyper([1,2,3,4,5],[6,7,8,9],5).ae(1.33980653631074769423-0.07143405251029226699j)
  567. assert hyper([1,2.79,3.08,4.37],[5.2,6.1,7.3],5).ae(1.0996321464692607231-1.7748052293979985001j)
  568. assert hyper([1,1,1],[1,2],1) == inf
  569. assert hyper([1,1,1],[2,(101,100)],1).ae(100.01621213528313220)
  570. # slow -- covered by doctests
  571. #assert hyper([1,1,1],[2,3],0.9999).ae(1.2897972005319693905)
  572. def test_hyper_u():
  573. mp.dps = 15
  574. assert hyperu(2,-3,0).ae(0.05)
  575. assert hyperu(2,-3.5,0).ae(4./99)
  576. assert hyperu(2,0,0) == 0.5
  577. assert hyperu(-5,1,0) == -120
  578. assert hyperu(-5,2,0) == inf
  579. assert hyperu(-5,-2,0) == 0
  580. assert hyperu(7,7,3).ae(0.00014681269365593503986) #exp(3)*gammainc(-6,3)
  581. assert hyperu(2,-3,4).ae(0.011836478100271995559)
  582. assert hyperu(3,4,5).ae(1./125)
  583. assert hyperu(2,3,0.0625) == 256
  584. assert hyperu(-1,2,0.25+0.5j) == -1.75+0.5j
  585. assert hyperu(0.5,1.5,7.25).ae(2/sqrt(29))
  586. assert hyperu(2,6,pi).ae(0.55804439825913399130)
  587. assert (hyperu((3,2),8,100+201j)*10**4).ae(-0.3797318333856738798 - 2.9974928453561707782j)
  588. assert (hyperu((5,2),(-1,2),-5000)*10**10).ae(-5.6681877926881664678j)
  589. # XXX: fails because of undetected cancellation in low level series code
  590. # Alternatively: could use asymptotic series here, if convergence test
  591. # tweaked back to recognize this one
  592. #assert (hyperu((5,2),(-1,2),-500)*10**7).ae(-1.82526906001593252847j)
  593. def test_hyper_2f0():
  594. mp.dps = 15
  595. assert hyper([1,2],[],3) == hyp2f0(1,2,3)
  596. assert hyp2f0(2,3,7).ae(0.0116108068639728714668 - 0.0073727413865865802130j)
  597. assert hyp2f0(2,3,0) == 1
  598. assert hyp2f0(0,0,0) == 1
  599. assert hyp2f0(-1,-1,1).ae(2)
  600. assert hyp2f0(-4,1,1.5).ae(62.5)
  601. assert hyp2f0(-4,1,50).ae(147029801)
  602. assert hyp2f0(-4,1,0.0001).ae(0.99960011997600240000)
  603. assert hyp2f0(0.5,0.25,0.001).ae(1.0001251174078538115)
  604. assert hyp2f0(0.5,0.25,3+4j).ae(0.85548875824755163518 + 0.21636041283392292973j)
  605. # Important: cancellation check
  606. assert hyp2f0((1,6),(5,6),-0.02371708245126284498).ae(0.996785723120804309)
  607. # Should be exact; polynomial case
  608. assert hyp2f0(-2,1,0.5+0.5j,zeroprec=200) == 0
  609. assert hyp2f0(1,-2,0.5+0.5j,zeroprec=200) == 0
  610. # There used to be a bug in thresholds that made one of the following hang
  611. for d in [15, 50, 80]:
  612. mp.dps = d
  613. assert hyp2f0(1.5, 0.5, 0.009).ae('1.006867007239309717945323585695344927904000945829843527398772456281301440034218290443367270629519483 + 1.238277162240704919639384945859073461954721356062919829456053965502443570466701567100438048602352623e-46j')
  614. def test_hyper_1f2():
  615. mp.dps = 15
  616. assert hyper([1],[2,3],4) == hyp1f2(1,2,3,4)
  617. a1,b1,b2 = (1,10),(2,3),1./16
  618. assert hyp1f2(a1,b1,b2,10).ae(298.7482725554557568)
  619. assert hyp1f2(a1,b1,b2,100).ae(224128961.48602947604)
  620. assert hyp1f2(a1,b1,b2,1000).ae(1.1669528298622675109e+27)
  621. assert hyp1f2(a1,b1,b2,10000).ae(2.4780514622487212192e+86)
  622. assert hyp1f2(a1,b1,b2,100000).ae(1.3885391458871523997e+274)
  623. assert hyp1f2(a1,b1,b2,1000000).ae('9.8851796978960318255e+867')
  624. assert hyp1f2(a1,b1,b2,10**7).ae('1.1505659189516303646e+2746')
  625. assert hyp1f2(a1,b1,b2,10**8).ae('1.4672005404314334081e+8685')
  626. assert hyp1f2(a1,b1,b2,10**20).ae('3.6888217332150976493e+8685889636')
  627. assert hyp1f2(a1,b1,b2,10*j).ae(-16.163252524618572878 - 44.321567896480184312j)
  628. assert hyp1f2(a1,b1,b2,100*j).ae(61938.155294517848171 + 637349.45215942348739j)
  629. assert hyp1f2(a1,b1,b2,1000*j).ae(8455057657257695958.7 + 6261969266997571510.6j)
  630. assert hyp1f2(a1,b1,b2,10000*j).ae(-8.9771211184008593089e+60 + 4.6550528111731631456e+59j)
  631. assert hyp1f2(a1,b1,b2,100000*j).ae(2.6398091437239324225e+193 + 4.1658080666870618332e+193j)
  632. assert hyp1f2(a1,b1,b2,1000000*j).ae('3.5999042951925965458e+613 + 1.5026014707128947992e+613j')
  633. assert hyp1f2(a1,b1,b2,10**7*j).ae('-8.3208715051623234801e+1939 - 3.6752883490851869429e+1941j')
  634. assert hyp1f2(a1,b1,b2,10**8*j).ae('2.0724195707891484454e+6140 - 1.3276619482724266387e+6141j')
  635. assert hyp1f2(a1,b1,b2,10**20*j).ae('-1.1734497974795488504e+6141851462 + 1.1498106965385471542e+6141851462j')
  636. def test_hyper_2f3():
  637. mp.dps = 15
  638. assert hyper([1,2],[3,4,5],6) == hyp2f3(1,2,3,4,5,6)
  639. a1,a2,b1,b2,b3 = (1,10),(2,3),(3,10), 2, 1./16
  640. # Check asymptotic expansion
  641. assert hyp2f3(a1,a2,b1,b2,b3,10).ae(128.98207160698659976)
  642. assert hyp2f3(a1,a2,b1,b2,b3,1000).ae(6.6309632883131273141e25)
  643. assert hyp2f3(a1,a2,b1,b2,b3,10000).ae(4.6863639362713340539e84)
  644. assert hyp2f3(a1,a2,b1,b2,b3,100000).ae(8.6632451236103084119e271)
  645. assert hyp2f3(a1,a2,b1,b2,b3,10**6).ae('2.0291718386574980641e865')
  646. assert hyp2f3(a1,a2,b1,b2,b3,10**7).ae('7.7639836665710030977e2742')
  647. assert hyp2f3(a1,a2,b1,b2,b3,10**8).ae('3.2537462584071268759e8681')
  648. assert hyp2f3(a1,a2,b1,b2,b3,10**20).ae('1.2966030542911614163e+8685889627')
  649. assert hyp2f3(a1,a2,b1,b2,b3,10*j).ae(-18.551602185587547854 - 13.348031097874113552j)
  650. assert hyp2f3(a1,a2,b1,b2,b3,100*j).ae(78634.359124504488695 + 74459.535945281973996j)
  651. assert hyp2f3(a1,a2,b1,b2,b3,1000*j).ae(597682550276527901.59 - 65136194809352613.078j)
  652. assert hyp2f3(a1,a2,b1,b2,b3,10000*j).ae(-1.1779696326238582496e+59 + 1.2297607505213133872e+59j)
  653. assert hyp2f3(a1,a2,b1,b2,b3,100000*j).ae(2.9844228969804380301e+191 + 7.5587163231490273296e+190j)
  654. assert hyp2f3(a1,a2,b1,b2,b3,1000000*j).ae('7.4859161049322370311e+610 - 2.8467477015940090189e+610j')
  655. assert hyp2f3(a1,a2,b1,b2,b3,10**7*j).ae('-1.7477645579418800826e+1938 - 1.7606522995808116405e+1938j')
  656. assert hyp2f3(a1,a2,b1,b2,b3,10**8*j).ae('-1.6932731942958401784e+6137 - 2.4521909113114629368e+6137j')
  657. assert hyp2f3(a1,a2,b1,b2,b3,10**20*j).ae('-2.0988815677627225449e+6141851451 + 5.7708223542739208681e+6141851452j')
  658. def test_hyper_2f2():
  659. mp.dps = 15
  660. assert hyper([1,2],[3,4],5) == hyp2f2(1,2,3,4,5)
  661. a1,a2,b1,b2 = (3,10),4,(1,2),1./16
  662. assert hyp2f2(a1,a2,b1,b2,10).ae(448225936.3377556696)
  663. assert hyp2f2(a1,a2,b1,b2,10000).ae('1.2012553712966636711e+4358')
  664. assert hyp2f2(a1,a2,b1,b2,-20000).ae(-0.04182343755661214626)
  665. assert hyp2f2(a1,a2,b1,b2,10**20).ae('1.1148680024303263661e+43429448190325182840')
  666. def test_orthpoly():
  667. mp.dps = 15
  668. assert jacobi(-4,2,3,0.7).ae(22800./4913)
  669. assert jacobi(3,2,4,5.5) == 4133.125
  670. assert jacobi(1.5,5/6.,4,0).ae(-1.0851951434075508417)
  671. assert jacobi(-2, 1, 2, 4).ae(-0.16)
  672. assert jacobi(2, -1, 2.5, 4).ae(34.59375)
  673. #assert jacobi(2, -1, 2, 4) == 28.5
  674. assert legendre(5, 7) == 129367
  675. assert legendre(0.5,0).ae(0.53935260118837935667)
  676. assert legendre(-1,-1) == 1
  677. assert legendre(0,-1) == 1
  678. assert legendre(0, 1) == 1
  679. assert legendre(1, -1) == -1
  680. assert legendre(7, 1) == 1
  681. assert legendre(7, -1) == -1
  682. assert legendre(8,1.5).ae(15457523./32768)
  683. assert legendre(j,-j).ae(2.4448182735671431011 + 0.6928881737669934843j)
  684. assert chebyu(5,1) == 6
  685. assert chebyt(3,2) == 26
  686. assert legendre(3.5,-1) == inf
  687. assert legendre(4.5,-1) == -inf
  688. assert legendre(3.5+1j,-1) == mpc(inf,inf)
  689. assert legendre(4.5+1j,-1) == mpc(-inf,-inf)
  690. assert laguerre(4, -2, 3).ae(-1.125)
  691. assert laguerre(3, 1+j, 0.5).ae(0.2291666666666666667 + 2.5416666666666666667j)
  692. def test_hermite():
  693. mp.dps = 15
  694. assert hermite(-2, 0).ae(0.5)
  695. assert hermite(-1, 0).ae(0.88622692545275801365)
  696. assert hermite(0, 0).ae(1)
  697. assert hermite(1, 0) == 0
  698. assert hermite(2, 0).ae(-2)
  699. assert hermite(0, 2).ae(1)
  700. assert hermite(1, 2).ae(4)
  701. assert hermite(1, -2).ae(-4)
  702. assert hermite(2, -2).ae(14)
  703. assert hermite(0.5, 0).ae(0.69136733903629335053)
  704. assert hermite(9, 0) == 0
  705. assert hermite(4,4).ae(3340)
  706. assert hermite(3,4).ae(464)
  707. assert hermite(-4,4).ae(0.00018623860287512396181)
  708. assert hermite(-3,4).ae(0.0016540169879668766270)
  709. assert hermite(9, 2.5j).ae(13638725j)
  710. assert hermite(9, -2.5j).ae(-13638725j)
  711. assert hermite(9, 100).ae(511078883759363024000)
  712. assert hermite(9, -100).ae(-511078883759363024000)
  713. assert hermite(9, 100j).ae(512922083920643024000j)
  714. assert hermite(9, -100j).ae(-512922083920643024000j)
  715. assert hermite(-9.5, 2.5j).ae(-2.9004951258126778174e-6 + 1.7601372934039951100e-6j)
  716. assert hermite(-9.5, -2.5j).ae(-2.9004951258126778174e-6 - 1.7601372934039951100e-6j)
  717. assert hermite(-9.5, 100).ae(1.3776300722767084162e-22, abs_eps=0, rel_eps=eps)
  718. assert hermite(-9.5, -100).ae('1.3106082028470671626e4355')
  719. assert hermite(-9.5, 100j).ae(-9.7900218581864768430e-23 - 9.7900218581864768430e-23j, abs_eps=0, rel_eps=eps)
  720. assert hermite(-9.5, -100j).ae(-9.7900218581864768430e-23 + 9.7900218581864768430e-23j, abs_eps=0, rel_eps=eps)
  721. assert hermite(2+3j, -1-j).ae(851.3677063883687676 - 1496.4373467871007997j)
  722. def test_gegenbauer():
  723. mp.dps = 15
  724. assert gegenbauer(1,2,3).ae(12)
  725. assert gegenbauer(2,3,4).ae(381)
  726. assert gegenbauer(0,0,0) == 0
  727. assert gegenbauer(2,-1,3) == 0
  728. assert gegenbauer(-7, 0.5, 3).ae(8989)
  729. assert gegenbauer(1, -0.5, 3).ae(-3)
  730. assert gegenbauer(1, -1.5, 3).ae(-9)
  731. assert gegenbauer(1, -0.5, 3).ae(-3)
  732. assert gegenbauer(-0.5, -0.5, 3).ae(-2.6383553159023906245)
  733. assert gegenbauer(2+3j, 1-j, 3+4j).ae(14.880536623203696780 + 20.022029711598032898j)
  734. #assert gegenbauer(-2, -0.5, 3).ae(-12)
  735. def test_legenp():
  736. mp.dps = 15
  737. assert legenp(2,0,4) == legendre(2,4)
  738. assert legenp(-2, -1, 0.5).ae(0.43301270189221932338)
  739. assert legenp(-2, -1, 0.5, type=3).ae(0.43301270189221932338j)
  740. assert legenp(-2, 1, 0.5).ae(-0.86602540378443864676)
  741. assert legenp(2+j, 3+4j, -j).ae(134742.98773236786148 + 429782.72924463851745j)
  742. assert legenp(2+j, 3+4j, -j, type=3).ae(802.59463394152268507 - 251.62481308942906447j)
  743. assert legenp(2,4,3).ae(0)
  744. assert legenp(2,4,3,type=3).ae(0)
  745. assert legenp(2,1,0.5).ae(-1.2990381056766579701)
  746. assert legenp(2,1,0.5,type=3).ae(1.2990381056766579701j)
  747. assert legenp(3,2,3).ae(-360)
  748. assert legenp(3,3,3).ae(240j*2**0.5)
  749. assert legenp(3,4,3).ae(0)
  750. assert legenp(0,0.5,2).ae(0.52503756790433198939 - 0.52503756790433198939j)
  751. assert legenp(-1,-0.5,2).ae(0.60626116232846498110 + 0.60626116232846498110j)
  752. assert legenp(-2,0.5,2).ae(1.5751127037129959682 - 1.5751127037129959682j)
  753. assert legenp(-2,0.5,-0.5).ae(-0.85738275810499171286)
  754. def test_legenq():
  755. mp.dps = 15
  756. f = legenq
  757. # Evaluation at poles
  758. assert isnan(f(3,2,1))
  759. assert isnan(f(3,2,-1))
  760. assert isnan(f(3,2,1,type=3))
  761. assert isnan(f(3,2,-1,type=3))
  762. # Evaluation at 0
  763. assert f(0,1,0,type=2).ae(-1)
  764. assert f(-2,2,0,type=2,zeroprec=200).ae(0)
  765. assert f(1.5,3,0,type=2).ae(-2.2239343475841951023)
  766. assert f(0,1,0,type=3).ae(j)
  767. assert f(-2,2,0,type=3,zeroprec=200).ae(0)
  768. assert f(1.5,3,0,type=3).ae(2.2239343475841951022*(1-1j))
  769. # Standard case, degree 0
  770. assert f(0,0,-1.5).ae(-0.8047189562170501873 + 1.5707963267948966192j)
  771. assert f(0,0,-0.5).ae(-0.54930614433405484570)
  772. assert f(0,0,0,zeroprec=200).ae(0)
  773. assert f(0,0,0.5).ae(0.54930614433405484570)
  774. assert f(0,0,1.5).ae(0.8047189562170501873 - 1.5707963267948966192j)
  775. assert f(0,0,-1.5,type=3).ae(-0.80471895621705018730)
  776. assert f(0,0,-0.5,type=3).ae(-0.5493061443340548457 - 1.5707963267948966192j)
  777. assert f(0,0,0,type=3).ae(-1.5707963267948966192j)
  778. assert f(0,0,0.5,type=3).ae(0.5493061443340548457 - 1.5707963267948966192j)
  779. assert f(0,0,1.5,type=3).ae(0.80471895621705018730)
  780. # Standard case, degree 1
  781. assert f(1,0,-1.5).ae(0.2070784343255752810 - 2.3561944901923449288j)
  782. assert f(1,0,-0.5).ae(-0.72534692783297257715)
  783. assert f(1,0,0).ae(-1)
  784. assert f(1,0,0.5).ae(-0.72534692783297257715)
  785. assert f(1,0,1.5).ae(0.2070784343255752810 - 2.3561944901923449288j)
  786. # Standard case, degree 2
  787. assert f(2,0,-1.5).ae(-0.0635669991240192885 + 4.5160394395353277803j)
  788. assert f(2,0,-0.5).ae(0.81866326804175685571)
  789. assert f(2,0,0,zeroprec=200).ae(0)
  790. assert f(2,0,0.5).ae(-0.81866326804175685571)
  791. assert f(2,0,1.5).ae(0.0635669991240192885 - 4.5160394395353277803j)
  792. # Misc orders and degrees
  793. assert f(2,3,1.5,type=2).ae(-5.7243340223994616228j)
  794. assert f(2,3,1.5,type=3).ae(-5.7243340223994616228)
  795. assert f(2,3,0.5,type=2).ae(-12.316805742712016310)
  796. assert f(2,3,0.5,type=3).ae(-12.316805742712016310j)
  797. assert f(2,3,-1.5,type=2).ae(-5.7243340223994616228j)
  798. assert f(2,3,-1.5,type=3).ae(5.7243340223994616228)
  799. assert f(2,3,-0.5,type=2).ae(-12.316805742712016310)
  800. assert f(2,3,-0.5,type=3).ae(-12.316805742712016310j)
  801. assert f(2+3j, 3+4j, 0.5, type=3).ae(0.0016119404873235186807 - 0.0005885900510718119836j)
  802. assert f(2+3j, 3+4j, -1.5, type=3).ae(0.008451400254138808670 + 0.020645193304593235298j)
  803. assert f(-2.5,1,-1.5).ae(3.9553395527435335749j)
  804. assert f(-2.5,1,-0.5).ae(1.9290561746445456908)
  805. assert f(-2.5,1,0).ae(1.2708196271909686299)
  806. assert f(-2.5,1,0.5).ae(-0.31584812990742202869)
  807. assert f(-2.5,1,1.5).ae(-3.9553395527435335742 + 0.2993235655044701706j)
  808. assert f(-2.5,1,-1.5,type=3).ae(0.29932356550447017254j)
  809. assert f(-2.5,1,-0.5,type=3).ae(-0.3158481299074220287 - 1.9290561746445456908j)
  810. assert f(-2.5,1,0,type=3).ae(1.2708196271909686292 - 1.2708196271909686299j)
  811. assert f(-2.5,1,0.5,type=3).ae(1.9290561746445456907 + 0.3158481299074220287j)
  812. assert f(-2.5,1,1.5,type=3).ae(-0.29932356550447017254)
  813. def test_agm():
  814. mp.dps = 15
  815. assert agm(0,0) == 0
  816. assert agm(0,1) == 0
  817. assert agm(1,1) == 1
  818. assert agm(7,7) == 7
  819. assert agm(j,j) == j
  820. assert (1/agm(1,sqrt(2))).ae(0.834626841674073186)
  821. assert agm(1,2).ae(1.4567910310469068692)
  822. assert agm(1,3).ae(1.8636167832448965424)
  823. assert agm(1,j).ae(0.599070117367796104+0.599070117367796104j)
  824. assert agm(2) == agm(1,2)
  825. assert agm(-3,4).ae(0.63468509766550907+1.3443087080896272j)
  826. def test_gammainc():
  827. mp.dps = 15
  828. assert gammainc(2,5).ae(6*exp(-5))
  829. assert gammainc(2,0,5).ae(1-6*exp(-5))
  830. assert gammainc(2,3,5).ae(-6*exp(-5)+4*exp(-3))
  831. assert gammainc(-2.5,-0.5).ae(-0.9453087204829418812-5.3164237738936178621j)
  832. assert gammainc(0,2,4).ae(0.045121158298212213088)
  833. assert gammainc(0,3).ae(0.013048381094197037413)
  834. assert gammainc(0,2+j,1-j).ae(0.00910653685850304839-0.22378752918074432574j)
  835. assert gammainc(0,1-j).ae(0.00028162445198141833+0.17932453503935894015j)
  836. assert gammainc(3,4,5,True).ae(0.11345128607046320253)
  837. assert gammainc(3.5,0,inf).ae(gamma(3.5))
  838. assert gammainc(-150.5,500).ae('6.9825435345798951153e-627')
  839. assert gammainc(-150.5,800).ae('4.6885137549474089431e-788')
  840. assert gammainc(-3.5, -20.5).ae(0.27008820585226911 - 1310.31447140574997636j)
  841. assert gammainc(-3.5, -200.5).ae(0.27008820585226911 - 5.3264597096208368435e76j) # XXX real part
  842. assert gammainc(0,0,2) == inf
  843. assert gammainc(1,b=1).ae(0.6321205588285576784)
  844. assert gammainc(3,2,2) == 0
  845. assert gammainc(2,3+j,3-j).ae(-0.28135485191849314194j)
  846. assert gammainc(4+0j,1).ae(5.8860710587430771455)
  847. # GH issue #301
  848. assert gammainc(-1,-1).ae(-0.8231640121031084799 + 3.1415926535897932385j)
  849. assert gammainc(-2,-1).ae(1.7707229202810768576 - 1.5707963267948966192j)
  850. assert gammainc(-3,-1).ae(-1.4963349162467073643 + 0.5235987755982988731j)
  851. assert gammainc(-4,-1).ae(1.05365418617643814992 - 0.13089969389957471827j)
  852. # Regularized upper gamma
  853. assert isnan(gammainc(0, 0, regularized=True))
  854. assert gammainc(-1, 0, regularized=True) == inf
  855. assert gammainc(1, 0, regularized=True) == 1
  856. assert gammainc(0, 5, regularized=True) == 0
  857. assert gammainc(0, 2+3j, regularized=True) == 0
  858. assert gammainc(0, 5000, regularized=True) == 0
  859. assert gammainc(0, 10**30, regularized=True) == 0
  860. assert gammainc(-1, 5, regularized=True) == 0
  861. assert gammainc(-1, 5000, regularized=True) == 0
  862. assert gammainc(-1, 10**30, regularized=True) == 0
  863. assert gammainc(-1, -5, regularized=True) == 0
  864. assert gammainc(-1, -5000, regularized=True) == 0
  865. assert gammainc(-1, -10**30, regularized=True) == 0
  866. assert gammainc(-1, 3+4j, regularized=True) == 0
  867. assert gammainc(1, 5, regularized=True).ae(exp(-5))
  868. assert gammainc(1, 5000, regularized=True).ae(exp(-5000))
  869. assert gammainc(1, 10**30, regularized=True).ae(exp(-10**30))
  870. assert gammainc(1, 3+4j, regularized=True).ae(exp(-3-4j))
  871. assert gammainc(-1000000,2).ae('1.3669297209397347754e-301037', abs_eps=0, rel_eps=8*eps)
  872. assert gammainc(-1000000,2,regularized=True) == 0
  873. assert gammainc(-1000000,3+4j).ae('-1.322575609404222361e-698979 - 4.9274570591854533273e-698978j', abs_eps=0, rel_eps=8*eps)
  874. assert gammainc(-1000000,3+4j,regularized=True) == 0
  875. assert gammainc(2+3j, 4+5j, regularized=True).ae(0.085422013530993285774-0.052595379150390078503j)
  876. assert gammainc(1000j, 1000j, regularized=True).ae(0.49702647628921131761 + 0.00297355675013575341j)
  877. # Generalized
  878. assert gammainc(3,4,2) == -gammainc(3,2,4)
  879. assert gammainc(4, 2, 3).ae(1.2593494302978947396)
  880. assert gammainc(4, 2, 3, regularized=True).ae(0.20989157171631578993)
  881. assert gammainc(0, 2, 3).ae(0.035852129613864082155)
  882. assert gammainc(0, 2, 3, regularized=True) == 0
  883. assert gammainc(-1, 2, 3).ae(0.015219822548487616132)
  884. assert gammainc(-1, 2, 3, regularized=True) == 0
  885. assert gammainc(0, 2, 3).ae(0.035852129613864082155)
  886. assert gammainc(0, 2, 3, regularized=True) == 0
  887. # Should use upper gammas
  888. assert gammainc(5, 10000, 12000).ae('1.1359381951461801687e-4327', abs_eps=0, rel_eps=8*eps)
  889. # Should use lower gammas
  890. assert gammainc(10000, 2, 3).ae('8.1244514125995785934e4765')
  891. # GH issue 306
  892. assert gammainc(3,-1-1j) == 0
  893. assert gammainc(3,-1+1j) == 0
  894. assert gammainc(2,-1) == 0
  895. assert gammainc(2,-1+0j) == 0
  896. assert gammainc(2+0j,-1) == 0
  897. def test_gammainc_expint_n():
  898. # These tests are intended to check all cases of the low-level code
  899. # for upper gamma and expint with small integer index.
  900. # Need to cover positive/negative arguments; small/large/huge arguments
  901. # for both positive and negative indices, as well as indices 0 and 1
  902. # which may be special-cased
  903. mp.dps = 15
  904. assert expint(-3,3.5).ae(0.021456366563296693987)
  905. assert expint(-2,3.5).ae(0.014966633183073309405)
  906. assert expint(-1,3.5).ae(0.011092916359219041088)
  907. assert expint(0,3.5).ae(0.0086278238349481430685)
  908. assert expint(1,3.5).ae(0.0069701398575483929193)
  909. assert expint(2,3.5).ae(0.0058018939208991255223)
  910. assert expint(3,3.5).ae(0.0049453773495857807058)
  911. assert expint(-3,-3.5).ae(-4.6618170604073311319)
  912. assert expint(-2,-3.5).ae(-5.5996974157555515963)
  913. assert expint(-1,-3.5).ae(-6.7582555017739415818)
  914. assert expint(0,-3.5).ae(-9.4615577024835182145)
  915. assert expint(1,-3.5).ae(-13.925353995152335292 - 3.1415926535897932385j)
  916. assert expint(2,-3.5).ae(-15.62328702434085977 - 10.995574287564276335j)
  917. assert expint(3,-3.5).ae(-10.783026313250347722 - 19.242255003237483586j)
  918. assert expint(-3,350).ae(2.8614825451252838069e-155, abs_eps=0, rel_eps=8*eps)
  919. assert expint(-2,350).ae(2.8532837224504675901e-155, abs_eps=0, rel_eps=8*eps)
  920. assert expint(-1,350).ae(2.8451316155828634555e-155, abs_eps=0, rel_eps=8*eps)
  921. assert expint(0,350).ae(2.8370258275042797989e-155, abs_eps=0, rel_eps=8*eps)
  922. assert expint(1,350).ae(2.8289659656701459404e-155, abs_eps=0, rel_eps=8*eps)
  923. assert expint(2,350).ae(2.8209516419468505006e-155, abs_eps=0, rel_eps=8*eps)
  924. assert expint(3,350).ae(2.8129824725501272171e-155, abs_eps=0, rel_eps=8*eps)
  925. assert expint(-3,-350).ae(-2.8528796154044839443e+149)
  926. assert expint(-2,-350).ae(-2.8610072121701264351e+149)
  927. assert expint(-1,-350).ae(-2.8691813842677537647e+149)
  928. assert expint(0,-350).ae(-2.8774025343659421709e+149)
  929. u = expint(1,-350)
  930. assert u.ae(-2.8856710698020863568e+149)
  931. assert u.imag.ae(-3.1415926535897932385)
  932. u = expint(2,-350)
  933. assert u.ae(-2.8939874026504650534e+149)
  934. assert u.imag.ae(-1099.5574287564276335)
  935. u = expint(3,-350)
  936. assert u.ae(-2.9023519497915044349e+149)
  937. assert u.imag.ae(-192422.55003237483586)
  938. assert expint(-3,350000000000000000000000).ae('2.1592908471792544286e-152003068666138139677919', abs_eps=0, rel_eps=8*eps)
  939. assert expint(-2,350000000000000000000000).ae('2.1592908471792544286e-152003068666138139677919', abs_eps=0, rel_eps=8*eps)
  940. assert expint(-1,350000000000000000000000).ae('2.1592908471792544286e-152003068666138139677919', abs_eps=0, rel_eps=8*eps)
  941. assert expint(0,350000000000000000000000).ae('2.1592908471792544286e-152003068666138139677919', abs_eps=0, rel_eps=8*eps)
  942. assert expint(1,350000000000000000000000).ae('2.1592908471792544286e-152003068666138139677919', abs_eps=0, rel_eps=8*eps)
  943. assert expint(2,350000000000000000000000).ae('2.1592908471792544286e-152003068666138139677919', abs_eps=0, rel_eps=8*eps)
  944. assert expint(3,350000000000000000000000).ae('2.1592908471792544286e-152003068666138139677919', abs_eps=0, rel_eps=8*eps)
  945. assert expint(-3,-350000000000000000000000).ae('-3.7805306852415755699e+152003068666138139677871')
  946. assert expint(-2,-350000000000000000000000).ae('-3.7805306852415755699e+152003068666138139677871')
  947. assert expint(-1,-350000000000000000000000).ae('-3.7805306852415755699e+152003068666138139677871')
  948. assert expint(0,-350000000000000000000000).ae('-3.7805306852415755699e+152003068666138139677871')
  949. u = expint(1,-350000000000000000000000)
  950. assert u.ae('-3.7805306852415755699e+152003068666138139677871')
  951. assert u.imag.ae(-3.1415926535897932385)
  952. u = expint(2,-350000000000000000000000)
  953. assert u.imag.ae(-1.0995574287564276335e+24)
  954. assert u.ae('-3.7805306852415755699e+152003068666138139677871')
  955. u = expint(3,-350000000000000000000000)
  956. assert u.imag.ae(-1.9242255003237483586e+47)
  957. assert u.ae('-3.7805306852415755699e+152003068666138139677871')
  958. # Small case; no branch cut
  959. assert gammainc(-3,3.5).ae(0.00010020262545203707109)
  960. assert gammainc(-2,3.5).ae(0.00040370427343557393517)
  961. assert gammainc(-1,3.5).ae(0.0016576839773997501492)
  962. assert gammainc(0,3.5).ae(0.0069701398575483929193)
  963. assert gammainc(1,3.5).ae(0.03019738342231850074)
  964. assert gammainc(2,3.5).ae(0.13588822540043325333)
  965. assert gammainc(3,3.5).ae(0.64169439772426814072)
  966. # Small case; with branch cut
  967. assert gammainc(-3,-3.5).ae(0.03595832954467563286 + 0.52359877559829887308j)
  968. assert gammainc(-2,-3.5).ae(-0.88024704597962022221 - 1.5707963267948966192j)
  969. assert gammainc(-1,-3.5).ae(4.4637962926688170771 + 3.1415926535897932385j)
  970. assert gammainc(0,-3.5).ae(-13.925353995152335292 - 3.1415926535897932385j)
  971. assert gammainc(1,-3.5).ae(33.115451958692313751)
  972. assert gammainc(2,-3.5).ae(-82.788629896730784377)
  973. assert gammainc(3,-3.5).ae(240.08702670051927469)
  974. # Asymptotic case; no branch cut
  975. assert gammainc(-3,350).ae(6.5424095113340358813e-163, abs_eps=0, rel_eps=8*eps)
  976. assert gammainc(-2,350).ae(2.296312222489899769e-160, abs_eps=0, rel_eps=8*eps)
  977. assert gammainc(-1,350).ae(8.059861834133858573e-158, abs_eps=0, rel_eps=8*eps)
  978. assert gammainc(0,350).ae(2.8289659656701459404e-155, abs_eps=0, rel_eps=8*eps)
  979. assert gammainc(1,350).ae(9.9295903962649792963e-153, abs_eps=0, rel_eps=8*eps)
  980. assert gammainc(2,350).ae(3.485286229089007733e-150, abs_eps=0, rel_eps=8*eps)
  981. assert gammainc(3,350).ae(1.2233453960006379793e-147, abs_eps=0, rel_eps=8*eps)
  982. # Asymptotic case; branch cut
  983. u = gammainc(-3,-350)
  984. assert u.ae(6.7889565783842895085e+141)
  985. assert u.imag.ae(0.52359877559829887308)
  986. u = gammainc(-2,-350)
  987. assert u.ae(-2.3692668977889832121e+144)
  988. assert u.imag.ae(-1.5707963267948966192)
  989. u = gammainc(-1,-350)
  990. assert u.ae(8.2685354361441858669e+146)
  991. assert u.imag.ae(3.1415926535897932385)
  992. u = gammainc(0,-350)
  993. assert u.ae(-2.8856710698020863568e+149)
  994. assert u.imag.ae(-3.1415926535897932385)
  995. u = gammainc(1,-350)
  996. assert u.ae(1.0070908870280797598e+152)
  997. assert u.imag == 0
  998. u = gammainc(2,-350)
  999. assert u.ae(-3.5147471957279983618e+154)
  1000. assert u.imag == 0
  1001. u = gammainc(3,-350)
  1002. assert u.ae(1.2266568422179417091e+157)
  1003. assert u.imag == 0
  1004. # Extreme asymptotic case
  1005. assert gammainc(-3,350000000000000000000000).ae('5.0362468738874738859e-152003068666138139677990', abs_eps=0, rel_eps=8*eps)
  1006. assert gammainc(-2,350000000000000000000000).ae('1.7626864058606158601e-152003068666138139677966', abs_eps=0, rel_eps=8*eps)
  1007. assert gammainc(-1,350000000000000000000000).ae('6.1694024205121555102e-152003068666138139677943', abs_eps=0, rel_eps=8*eps)
  1008. assert gammainc(0,350000000000000000000000).ae('2.1592908471792544286e-152003068666138139677919', abs_eps=0, rel_eps=8*eps)
  1009. assert gammainc(1,350000000000000000000000).ae('7.5575179651273905e-152003068666138139677896', abs_eps=0, rel_eps=8*eps)
  1010. assert gammainc(2,350000000000000000000000).ae('2.645131287794586675e-152003068666138139677872', abs_eps=0, rel_eps=8*eps)
  1011. assert gammainc(3,350000000000000000000000).ae('9.2579595072810533625e-152003068666138139677849', abs_eps=0, rel_eps=8*eps)
  1012. u = gammainc(-3,-350000000000000000000000)
  1013. assert u.ae('8.8175642804468234866e+152003068666138139677800')
  1014. assert u.imag.ae(0.52359877559829887308)
  1015. u = gammainc(-2,-350000000000000000000000)
  1016. assert u.ae('-3.0861474981563882203e+152003068666138139677824')
  1017. assert u.imag.ae(-1.5707963267948966192)
  1018. u = gammainc(-1,-350000000000000000000000)
  1019. assert u.ae('1.0801516243547358771e+152003068666138139677848')
  1020. assert u.imag.ae(3.1415926535897932385)
  1021. u = gammainc(0,-350000000000000000000000)
  1022. assert u.ae('-3.7805306852415755699e+152003068666138139677871')
  1023. assert u.imag.ae(-3.1415926535897932385)
  1024. assert gammainc(1,-350000000000000000000000).ae('1.3231857398345514495e+152003068666138139677895')
  1025. assert gammainc(2,-350000000000000000000000).ae('-4.6311500894209300731e+152003068666138139677918')
  1026. assert gammainc(3,-350000000000000000000000).ae('1.6209025312973255256e+152003068666138139677942')
  1027. def test_incomplete_beta():
  1028. mp.dps = 15
  1029. assert betainc(-2,-3,0.5,0.75).ae(63.4305673311255413583969)
  1030. assert betainc(4.5,0.5+2j,2.5,6).ae(0.2628801146130621387903065 + 0.5162565234467020592855378j)
  1031. assert betainc(4,5,0,6).ae(90747.77142857142857142857)
  1032. def test_erf():
  1033. mp.dps = 15
  1034. assert erf(0) == 0
  1035. assert erf(1).ae(0.84270079294971486934)
  1036. assert erf(3+4j).ae(-120.186991395079444098 - 27.750337293623902498j)
  1037. assert erf(-4-3j).ae(-0.99991066178539168236 + 0.00004972026054496604j)
  1038. assert erf(pi).ae(0.99999112385363235839)
  1039. assert erf(1j).ae(1.6504257587975428760j)
  1040. assert erf(-1j).ae(-1.6504257587975428760j)
  1041. assert isinstance(erf(1), mpf)
  1042. assert isinstance(erf(-1), mpf)
  1043. assert isinstance(erf(0), mpf)
  1044. assert isinstance(erf(0j), mpc)
  1045. assert erf(inf) == 1
  1046. assert erf(-inf) == -1
  1047. assert erfi(0) == 0
  1048. assert erfi(1/pi).ae(0.371682698493894314)
  1049. assert erfi(inf) == inf
  1050. assert erfi(-inf) == -inf
  1051. assert erf(1+0j) == erf(1)
  1052. assert erfc(1+0j) == erfc(1)
  1053. assert erf(0.2+0.5j).ae(1 - erfc(0.2+0.5j))
  1054. assert erfc(0) == 1
  1055. assert erfc(1).ae(1-erf(1))
  1056. assert erfc(-1).ae(1-erf(-1))
  1057. assert erfc(1/pi).ae(1-erf(1/pi))
  1058. assert erfc(-10) == 2
  1059. assert erfc(-1000000) == 2
  1060. assert erfc(-inf) == 2
  1061. assert erfc(inf) == 0
  1062. assert isnan(erfc(nan))
  1063. assert (erfc(10**4)*mpf(10)**43429453).ae('3.63998738656420')
  1064. assert erf(8+9j).ae(-1072004.2525062051158 + 364149.91954310255423j)
  1065. assert erfc(8+9j).ae(1072005.2525062051158 - 364149.91954310255423j)
  1066. assert erfc(-8-9j).ae(-1072003.2525062051158 + 364149.91954310255423j)
  1067. mp.dps = 50
  1068. # This one does not use the asymptotic series
  1069. assert (erfc(10)*10**45).ae('2.0884875837625447570007862949577886115608181193212')
  1070. # This one does
  1071. assert (erfc(50)*10**1088).ae('2.0709207788416560484484478751657887929322509209954')
  1072. mp.dps = 15
  1073. assert str(erfc(10**50)) == '3.66744826532555e-4342944819032518276511289189166050822943970058036665661144537831658646492088707747292249493384317534'
  1074. assert erfinv(0) == 0
  1075. assert erfinv(0.5).ae(0.47693627620446987338)
  1076. assert erfinv(-0.5).ae(-0.47693627620446987338)
  1077. assert erfinv(1) == inf
  1078. assert erfinv(-1) == -inf
  1079. assert erf(erfinv(0.95)).ae(0.95)
  1080. assert erf(erfinv(0.999999999995)).ae(0.999999999995)
  1081. assert erf(erfinv(-0.999999999995)).ae(-0.999999999995)
  1082. mp.dps = 50
  1083. assert erf(erfinv('0.99999999999999999999999999999995')).ae('0.99999999999999999999999999999995')
  1084. assert erf(erfinv('0.999999999999999999999999999999995')).ae('0.999999999999999999999999999999995')
  1085. assert erf(erfinv('-0.999999999999999999999999999999995')).ae('-0.999999999999999999999999999999995')
  1086. mp.dps = 15
  1087. # Complex asymptotic expansions
  1088. v = erfc(50j)
  1089. assert v.real == 1
  1090. assert v.imag.ae('-6.1481820666053078736e+1083')
  1091. assert erfc(-100+5j).ae(2)
  1092. assert (erfc(100+5j)*10**4335).ae(2.3973567853824133572 - 3.9339259530609420597j)
  1093. assert erfc(100+100j).ae(0.00065234366376857698698 - 0.0039357263629214118437j)
  1094. def test_pdf():
  1095. mp.dps = 15
  1096. assert npdf(-inf) == 0
  1097. assert npdf(inf) == 0
  1098. assert npdf(5,0,2).ae(npdf(5+4,4,2))
  1099. assert quadts(lambda x: npdf(x,-0.5,0.8), [-inf, inf]) == 1
  1100. assert ncdf(0) == 0.5
  1101. assert ncdf(3,3) == 0.5
  1102. assert ncdf(-inf) == 0
  1103. assert ncdf(inf) == 1
  1104. assert ncdf(10) == 1
  1105. # Verify that this is computed accurately
  1106. assert (ncdf(-10)*10**24).ae(7.619853024160526)
  1107. def test_lambertw():
  1108. mp.dps = 15
  1109. assert lambertw(0) == 0
  1110. assert lambertw(0+0j) == 0
  1111. assert lambertw(inf) == inf
  1112. assert isnan(lambertw(nan))
  1113. assert lambertw(inf,1).real == inf
  1114. assert lambertw(inf,1).imag.ae(2*pi)
  1115. assert lambertw(-inf,1).real == inf
  1116. assert lambertw(-inf,1).imag.ae(3*pi)
  1117. assert lambertw(0,-1) == -inf
  1118. assert lambertw(0,1) == -inf
  1119. assert lambertw(0,3) == -inf
  1120. assert lambertw(e).ae(1)
  1121. assert lambertw(1).ae(0.567143290409783873)
  1122. assert lambertw(-pi/2).ae(j*pi/2)
  1123. assert lambertw(-log(2)/2).ae(-log(2))
  1124. assert lambertw(0.25).ae(0.203888354702240164)
  1125. assert lambertw(-0.25).ae(-0.357402956181388903)
  1126. assert lambertw(-1./10000,0).ae(-0.000100010001500266719)
  1127. assert lambertw(-0.25,-1).ae(-2.15329236411034965)
  1128. assert lambertw(0.25,-1).ae(-3.00899800997004620-4.07652978899159763j)
  1129. assert lambertw(-0.25,-1).ae(-2.15329236411034965)
  1130. assert lambertw(0.25,1).ae(-3.00899800997004620+4.07652978899159763j)
  1131. assert lambertw(-0.25,1).ae(-3.48973228422959210+7.41405453009603664j)
  1132. assert lambertw(-4).ae(0.67881197132094523+1.91195078174339937j)
  1133. assert lambertw(-4,1).ae(-0.66743107129800988+7.76827456802783084j)
  1134. assert lambertw(-4,-1).ae(0.67881197132094523-1.91195078174339937j)
  1135. assert lambertw(1000).ae(5.24960285240159623)
  1136. assert lambertw(1000,1).ae(4.91492239981054535+5.44652615979447070j)
  1137. assert lambertw(1000,-1).ae(4.91492239981054535-5.44652615979447070j)
  1138. assert lambertw(1000,5).ae(3.5010625305312892+29.9614548941181328j)
  1139. assert lambertw(3+4j).ae(1.281561806123775878+0.533095222020971071j)
  1140. assert lambertw(-0.4+0.4j).ae(-0.10396515323290657+0.61899273315171632j)
  1141. assert lambertw(3+4j,1).ae(-0.11691092896595324+5.61888039871282334j)
  1142. assert lambertw(3+4j,-1).ae(0.25856740686699742-3.85211668616143559j)
  1143. assert lambertw(-0.5,-1).ae(-0.794023632344689368-0.770111750510379110j)
  1144. assert lambertw(-1./10000,1).ae(-11.82350837248724344+6.80546081842002101j)
  1145. assert lambertw(-1./10000,-1).ae(-11.6671145325663544)
  1146. assert lambertw(-1./10000,-2).ae(-11.82350837248724344-6.80546081842002101j)
  1147. assert lambertw(-1./100000,4).ae(-14.9186890769540539+26.1856750178782046j)
  1148. assert lambertw(-1./100000,5).ae(-15.0931437726379218666+32.5525721210262290086j)
  1149. assert lambertw((2+j)/10).ae(0.173704503762911669+0.071781336752835511j)
  1150. assert lambertw((2+j)/10,1).ae(-3.21746028349820063+4.56175438896292539j)
  1151. assert lambertw((2+j)/10,-1).ae(-3.03781405002993088-3.53946629633505737j)
  1152. assert lambertw((2+j)/10,4).ae(-4.6878509692773249+23.8313630697683291j)
  1153. assert lambertw(-(2+j)/10).ae(-0.226933772515757933-0.164986470020154580j)
  1154. assert lambertw(-(2+j)/10,1).ae(-2.43569517046110001+0.76974067544756289j)
  1155. assert lambertw(-(2+j)/10,-1).ae(-3.54858738151989450-6.91627921869943589j)
  1156. assert lambertw(-(2+j)/10,4).ae(-4.5500846928118151+20.6672982215434637j)
  1157. mp.dps = 50
  1158. assert lambertw(pi).ae('1.073658194796149172092178407024821347547745350410314531')
  1159. mp.dps = 15
  1160. # Former bug in generated branch
  1161. assert lambertw(-0.5+0.002j).ae(-0.78917138132659918344 + 0.76743539379990327749j)
  1162. assert lambertw(-0.5-0.002j).ae(-0.78917138132659918344 - 0.76743539379990327749j)
  1163. assert lambertw(-0.448+0.4j).ae(-0.11855133765652382241 + 0.66570534313583423116j)
  1164. assert lambertw(-0.448-0.4j).ae(-0.11855133765652382241 - 0.66570534313583423116j)
  1165. assert lambertw(-0.65475+0.0001j).ae(-0.61053421111385310898+1.0396534993944097723803j)
  1166. # Huge branch index
  1167. w = lambertw(1,10**20)
  1168. assert w.real.ae(-47.889578926290259164)
  1169. assert w.imag.ae(6.2831853071795864769e+20)
  1170. def test_lambertw_hard():
  1171. def check(x,y):
  1172. y = convert(y)
  1173. type_ok = True
  1174. if isinstance(y, mpf):
  1175. type_ok = isinstance(x, mpf)
  1176. real_ok = abs(x.real-y.real) <= abs(y.real)*8*eps
  1177. imag_ok = abs(x.imag-y.imag) <= abs(y.imag)*8*eps
  1178. #print x, y, abs(x.real-y.real), abs(x.imag-y.imag)
  1179. return real_ok and imag_ok
  1180. # Evaluation near 0
  1181. mp.dps = 15
  1182. assert check(lambertw(1e-10), 9.999999999000000000e-11)
  1183. assert check(lambertw(-1e-10), -1.000000000100000000e-10)
  1184. assert check(lambertw(1e-10j), 9.999999999999999999733e-21 + 9.99999999999999999985e-11j)
  1185. assert check(lambertw(-1e-10j), 9.999999999999999999733e-21 - 9.99999999999999999985e-11j)
  1186. assert check(lambertw(1e-10,1), -26.303186778379041559 + 3.265093911703828397j)
  1187. assert check(lambertw(-1e-10,1), -26.326236166739163892 + 6.526183280686333315j)
  1188. assert check(lambertw(1e-10j,1), -26.312931726911421551 + 4.896366881798013421j)
  1189. assert check(lambertw(-1e-10j,1), -26.297238779529035066 + 1.632807161345576513j)
  1190. assert check(lambertw(1e-10,-1), -26.303186778379041559 - 3.265093911703828397j)
  1191. assert check(lambertw(-1e-10,-1), -26.295238819246925694)
  1192. assert check(lambertw(1e-10j,-1), -26.297238779529035028 - 1.6328071613455765135j)
  1193. assert check(lambertw(-1e-10j,-1), -26.312931726911421551 - 4.896366881798013421j)
  1194. # Test evaluation very close to the branch point -1/e
  1195. # on the -1, 0, and 1 branches
  1196. add = lambda x, y: fadd(x,y,exact=True)
  1197. sub = lambda x, y: fsub(x,y,exact=True)
  1198. addj = lambda x, y: fadd(x,fmul(y,1j,exact=True),exact=True)
  1199. subj = lambda x, y: fadd(x,fmul(y,-1j,exact=True),exact=True)
  1200. mp.dps = 1500
  1201. a = -1/e + 10*eps
  1202. d3 = mpf('1e-3')
  1203. d10 = mpf('1e-10')
  1204. d20 = mpf('1e-20')
  1205. d40 = mpf('1e-40')
  1206. d80 = mpf('1e-80')
  1207. d300 = mpf('1e-300')
  1208. d1000 = mpf('1e-1000')
  1209. mp.dps = 15
  1210. # ---- Branch 0 ----
  1211. # -1/e + eps
  1212. assert check(lambertw(add(a,d3)), -0.92802015005456704876)
  1213. assert check(lambertw(add(a,d10)), -0.99997668374140088071)
  1214. assert check(lambertw(add(a,d20)), -0.99999999976683560186)
  1215. assert lambertw(add(a,d40)) == -1
  1216. assert lambertw(add(a,d80)) == -1
  1217. assert lambertw(add(a,d300)) == -1
  1218. assert lambertw(add(a,d1000)) == -1
  1219. # -1/e - eps
  1220. assert check(lambertw(sub(a,d3)), -0.99819016149860989001+0.07367191188934638577j)
  1221. assert check(lambertw(sub(a,d10)), -0.9999999998187812114595992+0.0000233164398140346109194j)
  1222. assert check(lambertw(sub(a,d20)), -0.99999999999999999998187+2.331643981597124203344e-10j)
  1223. assert check(lambertw(sub(a,d40)), -1.0+2.33164398159712420336e-20j)
  1224. assert check(lambertw(sub(a,d80)), -1.0+2.33164398159712420336e-40j)
  1225. assert check(lambertw(sub(a,d300)), -1.0+2.33164398159712420336e-150j)
  1226. assert check(lambertw(sub(a,d1000)), mpc(-1,'2.33164398159712420336e-500'))
  1227. # -1/e + eps*j
  1228. assert check(lambertw(addj(a,d3)), -0.94790387486938526634+0.05036819639190132490j)
  1229. assert check(lambertw(addj(a,d10)), -0.9999835127872943680999899+0.0000164870314895821225256j)
  1230. assert check(lambertw(addj(a,d20)), -0.999999999835127872929987+1.64872127051890935830e-10j)
  1231. assert check(lambertw(addj(a,d40)), -0.9999999999999999999835+1.6487212707001281468305e-20j)
  1232. assert check(lambertw(addj(a,d80)), -1.0 + 1.64872127070012814684865e-40j)
  1233. assert check(lambertw(addj(a,d300)), -1.0 + 1.64872127070012814684865e-150j)
  1234. assert check(lambertw(addj(a,d1000)), mpc(-1.0,'1.64872127070012814684865e-500'))
  1235. # -1/e - eps*j
  1236. assert check(lambertw(subj(a,d3)), -0.94790387486938526634-0.05036819639190132490j)
  1237. assert check(lambertw(subj(a,d10)), -0.9999835127872943680999899-0.0000164870314895821225256j)
  1238. assert check(lambertw(subj(a,d20)), -0.999999999835127872929987-1.64872127051890935830e-10j)
  1239. assert check(lambertw(subj(a,d40)), -0.9999999999999999999835-1.6487212707001281468305e-20j)
  1240. assert check(lambertw(subj(a,d80)), -1.0 - 1.64872127070012814684865e-40j)
  1241. assert check(lambertw(subj(a,d300)), -1.0 - 1.64872127070012814684865e-150j)
  1242. assert check(lambertw(subj(a,d1000)), mpc(-1.0,'-1.64872127070012814684865e-500'))
  1243. # ---- Branch 1 ----
  1244. assert check(lambertw(addj(a,d3),1), -3.088501303219933378005990 + 7.458676867597474813950098j)
  1245. assert check(lambertw(addj(a,d80),1), -3.088843015613043855957087 + 7.461489285654254556906117j)
  1246. assert check(lambertw(addj(a,d300),1), -3.088843015613043855957087 + 7.461489285654254556906117j)
  1247. assert check(lambertw(addj(a,d1000),1), -3.088843015613043855957087 + 7.461489285654254556906117j)
  1248. assert check(lambertw(subj(a,d3),1), -1.0520914180450129534365906 + 0.0539925638125450525673175j)
  1249. assert check(lambertw(subj(a,d10),1), -1.0000164872127056318529390 + 0.000016487393927159250398333077j)
  1250. assert check(lambertw(subj(a,d20),1), -1.0000000001648721270700128 + 1.64872127088134693542628e-10j)
  1251. assert check(lambertw(subj(a,d40),1), -1.000000000000000000016487 + 1.64872127070012814686677e-20j)
  1252. assert check(lambertw(subj(a,d80),1), -1.0 + 1.64872127070012814684865e-40j)
  1253. assert check(lambertw(subj(a,d300),1), -1.0 + 1.64872127070012814684865e-150j)
  1254. assert check(lambertw(subj(a,d1000),1), mpc(-1.0, '1.64872127070012814684865e-500'))
  1255. # ---- Branch -1 ----
  1256. # -1/e + eps
  1257. assert check(lambertw(add(a,d3),-1), -1.075608941186624989414945)
  1258. assert check(lambertw(add(a,d10),-1), -1.000023316621036696460620)
  1259. assert check(lambertw(add(a,d20),-1), -1.000000000233164398177834)
  1260. assert lambertw(add(a,d40),-1) == -1
  1261. assert lambertw(add(a,d80),-1) == -1
  1262. assert lambertw(add(a,d300),-1) == -1
  1263. assert lambertw(add(a,d1000),-1) == -1
  1264. # -1/e - eps
  1265. assert check(lambertw(sub(a,d3),-1), -0.99819016149860989001-0.07367191188934638577j)
  1266. assert check(lambertw(sub(a,d10),-1), -0.9999999998187812114595992-0.0000233164398140346109194j)
  1267. assert check(lambertw(sub(a,d20),-1), -0.99999999999999999998187-2.331643981597124203344e-10j)
  1268. assert check(lambertw(sub(a,d40),-1), -1.0-2.33164398159712420336e-20j)
  1269. assert check(lambertw(sub(a,d80),-1), -1.0-2.33164398159712420336e-40j)
  1270. assert check(lambertw(sub(a,d300),-1), -1.0-2.33164398159712420336e-150j)
  1271. assert check(lambertw(sub(a,d1000),-1), mpc(-1,'-2.33164398159712420336e-500'))
  1272. # -1/e + eps*j
  1273. assert check(lambertw(addj(a,d3),-1), -1.0520914180450129534365906 - 0.0539925638125450525673175j)
  1274. assert check(lambertw(addj(a,d10),-1), -1.0000164872127056318529390 - 0.0000164873939271592503983j)
  1275. assert check(lambertw(addj(a,d20),-1), -1.0000000001648721270700 - 1.64872127088134693542628e-10j)
  1276. assert check(lambertw(addj(a,d40),-1), -1.00000000000000000001648 - 1.6487212707001281468667726e-20j)
  1277. assert check(lambertw(addj(a,d80),-1), -1.0 - 1.64872127070012814684865e-40j)
  1278. assert check(lambertw(addj(a,d300),-1), -1.0 - 1.64872127070012814684865e-150j)
  1279. assert check(lambertw(addj(a,d1000),-1), mpc(-1.0,'-1.64872127070012814684865e-500'))
  1280. # -1/e - eps*j
  1281. assert check(lambertw(subj(a,d3),-1), -3.088501303219933378005990-7.458676867597474813950098j)
  1282. assert check(lambertw(subj(a,d10),-1), -3.088843015579260686911033-7.461489285372968780020716j)
  1283. assert check(lambertw(subj(a,d20),-1), -3.088843015613043855953708-7.461489285654254556877988j)
  1284. assert check(lambertw(subj(a,d40),-1), -3.088843015613043855957087-7.461489285654254556906117j)
  1285. assert check(lambertw(subj(a,d80),-1), -3.088843015613043855957087 - 7.461489285654254556906117j)
  1286. assert check(lambertw(subj(a,d300),-1), -3.088843015613043855957087 - 7.461489285654254556906117j)
  1287. assert check(lambertw(subj(a,d1000),-1), -3.088843015613043855957087 - 7.461489285654254556906117j)
  1288. # One more case, testing higher precision
  1289. mp.dps = 500
  1290. x = -1/e + mpf('1e-13')
  1291. ans = "-0.99999926266961377166355784455394913638782494543377383"\
  1292. "744978844374498153493943725364881490261187530235150668593869563"\
  1293. "168276697689459394902153960200361935311512317183678882"
  1294. mp.dps = 15
  1295. assert lambertw(x).ae(ans)
  1296. mp.dps = 50
  1297. assert lambertw(x).ae(ans)
  1298. mp.dps = 150
  1299. assert lambertw(x).ae(ans)
  1300. def test_meijerg():
  1301. mp.dps = 15
  1302. assert meijerg([[2,3],[1]],[[0.5,2],[3,4]], 2.5).ae(4.2181028074787439386)
  1303. assert meijerg([[],[1+j]],[[1],[1]], 3+4j).ae(271.46290321152464592 - 703.03330399954820169j)
  1304. assert meijerg([[0.25],[1]],[[0.5],[2]],0) == 0
  1305. assert meijerg([[0],[]],[[0,0,'1/3','2/3'], []], '2/27').ae(2.2019391389653314120)
  1306. # Verify 1/z series being used
  1307. assert meijerg([[-3],[-0.5]], [[-1],[-2.5]], -0.5).ae(-1.338096165935754898687431)
  1308. assert meijerg([[1-(-1)],[1-(-2.5)]], [[1-(-3)],[1-(-0.5)]], -2.0).ae(-1.338096165935754898687431)
  1309. assert meijerg([[-3],[-0.5]], [[-1],[-2.5]], -1).ae(-(pi+4)/(4*pi))
  1310. a = 2.5
  1311. b = 1.25
  1312. for z in [mpf(0.25), mpf(2)]:
  1313. x1 = hyp1f1(a,b,z)
  1314. x2 = gamma(b)/gamma(a)*meijerg([[1-a],[]],[[0],[1-b]],-z)
  1315. x3 = gamma(b)/gamma(a)*meijerg([[1-0],[1-(1-b)]],[[1-(1-a)],[]],-1/z)
  1316. assert x1.ae(x2)
  1317. assert x1.ae(x3)
  1318. def test_appellf1():
  1319. mp.dps = 15
  1320. assert appellf1(2,-2,1,1,2,3).ae(-1.75)
  1321. assert appellf1(2,1,-2,1,2,3).ae(-8)
  1322. assert appellf1(2,1,-2,1,0.5,0.25).ae(1.5)
  1323. assert appellf1(-2,1,3,2,3,3).ae(19)
  1324. assert appellf1(1,2,3,4,0.5,0.125).ae( 1.53843285792549786518)
  1325. def test_coulomb():
  1326. # Note: most tests are doctests
  1327. # Test for a bug:
  1328. mp.dps = 15
  1329. assert coulombg(mpc(-5,0),2,3).ae(20.087729487721430394)
  1330. def test_hyper_param_accuracy():
  1331. mp.dps = 15
  1332. As = [n+1e-10 for n in range(-5,-1)]
  1333. Bs = [n+1e-10 for n in range(-12,-5)]
  1334. assert hyper(As,Bs,10).ae(-381757055858.652671927)
  1335. assert legenp(0.5, 100, 0.25).ae(-2.4124576567211311755e+144)
  1336. assert (hyp1f1(1000,1,-100)*10**24).ae(5.2589445437370169113)
  1337. assert (hyp2f1(10, -900, 10.5, 0.99)*10**24).ae(1.9185370579660768203)
  1338. assert (hyp2f1(1000,1.5,-3.5,-1.5)*10**385).ae(-2.7367529051334000764)
  1339. assert hyp2f1(-5, 10, 3, 0.5, zeroprec=500) == 0
  1340. assert (hyp1f1(-10000, 1000, 100)*10**424).ae(-3.1046080515824859974)
  1341. assert (hyp2f1(1000,1.5,-3.5,-0.75,maxterms=100000)*10**231).ae(-4.0534790813913998643)
  1342. assert legenp(2, 3, 0.25) == 0
  1343. pytest.raises(ValueError, lambda: hypercomb(lambda a: [([],[],[],[],[a],[-a],0.5)], [3]))
  1344. assert hypercomb(lambda a: [([],[],[],[],[a],[-a],0.5)], [3], infprec=200) == inf
  1345. assert meijerg([[],[]],[[0,0,0,0],[]],0.1).ae(1.5680822343832351418)
  1346. assert (besselk(400,400)*10**94).ae(1.4387057277018550583)
  1347. mp.dps = 5
  1348. (hyp1f1(-5000.5, 1500, 100)*10**185).ae(8.5185229673381935522)
  1349. (hyp1f1(-5000, 1500, 100)*10**185).ae(9.1501213424563944311)
  1350. mp.dps = 15
  1351. (hyp1f1(-5000.5, 1500, 100)*10**185).ae(8.5185229673381935522)
  1352. (hyp1f1(-5000, 1500, 100)*10**185).ae(9.1501213424563944311)
  1353. assert hyp0f1(fadd(-20,'1e-100',exact=True), 0.25).ae(1.85014429040102783e+49)
  1354. assert hyp0f1((-20*10**100+1, 10**100), 0.25).ae(1.85014429040102783e+49)
  1355. def test_hypercomb_zero_pow():
  1356. # check that 0^0 = 1
  1357. assert hypercomb(lambda a: (([0],[a],[],[],[],[],0),), [0]) == 1
  1358. assert meijerg([[-1.5],[]],[[0],[-0.75]],0).ae(1.4464090846320771425)
  1359. def test_spherharm():
  1360. mp.dps = 15
  1361. t = 0.5; r = 0.25
  1362. assert spherharm(0,0,t,r).ae(0.28209479177387814347)
  1363. assert spherharm(1,-1,t,r).ae(0.16048941205971996369 - 0.04097967481096344271j)
  1364. assert spherharm(1,0,t,r).ae(0.42878904414183579379)
  1365. assert spherharm(1,1,t,r).ae(-0.16048941205971996369 - 0.04097967481096344271j)
  1366. assert spherharm(2,-2,t,r).ae(0.077915886919031181734 - 0.042565643022253962264j)
  1367. assert spherharm(2,-1,t,r).ae(0.31493387233497459884 - 0.08041582001959297689j)
  1368. assert spherharm(2,0,t,r).ae(0.41330596756220761898)
  1369. assert spherharm(2,1,t,r).ae(-0.31493387233497459884 - 0.08041582001959297689j)
  1370. assert spherharm(2,2,t,r).ae(0.077915886919031181734 + 0.042565643022253962264j)
  1371. assert spherharm(3,-3,t,r).ae(0.033640236589690881646 - 0.031339125318637082197j)
  1372. assert spherharm(3,-2,t,r).ae(0.18091018743101461963 - 0.09883168583167010241j)
  1373. assert spherharm(3,-1,t,r).ae(0.42796713930907320351 - 0.10927795157064962317j)
  1374. assert spherharm(3,0,t,r).ae(0.27861659336351639787)
  1375. assert spherharm(3,1,t,r).ae(-0.42796713930907320351 - 0.10927795157064962317j)
  1376. assert spherharm(3,2,t,r).ae(0.18091018743101461963 + 0.09883168583167010241j)
  1377. assert spherharm(3,3,t,r).ae(-0.033640236589690881646 - 0.031339125318637082197j)
  1378. assert spherharm(0,-1,t,r) == 0
  1379. assert spherharm(0,-2,t,r) == 0
  1380. assert spherharm(0,1,t,r) == 0
  1381. assert spherharm(0,2,t,r) == 0
  1382. assert spherharm(1,2,t,r) == 0
  1383. assert spherharm(1,3,t,r) == 0
  1384. assert spherharm(1,-2,t,r) == 0
  1385. assert spherharm(1,-3,t,r) == 0
  1386. assert spherharm(2,3,t,r) == 0
  1387. assert spherharm(2,4,t,r) == 0
  1388. assert spherharm(2,-3,t,r) == 0
  1389. assert spherharm(2,-4,t,r) == 0
  1390. assert spherharm(3,4.5,0.5,0.25).ae(-22.831053442240790148 + 10.910526059510013757j)
  1391. assert spherharm(2+3j, 1-j, 1+j, 3+4j).ae(-2.6582752037810116935 - 1.0909214905642160211j)
  1392. assert spherharm(-6,2.5,t,r).ae(0.39383644983851448178 + 0.28414687085358299021j)
  1393. assert spherharm(-3.5, 3, 0.5, 0.25).ae(0.014516852987544698924 - 0.015582769591477628495j)
  1394. assert spherharm(-3, 3, 0.5, 0.25) == 0
  1395. assert spherharm(-6, 3, 0.5, 0.25).ae(-0.16544349818782275459 - 0.15412657723253924562j)
  1396. assert spherharm(-6, 1.5, 0.5, 0.25).ae(0.032208193499767402477 + 0.012678000924063664921j)
  1397. assert spherharm(3,0,0,1).ae(0.74635266518023078283)
  1398. assert spherharm(3,-2,0,1) == 0
  1399. assert spherharm(3,-2,1,1).ae(-0.16270707338254028971 - 0.35552144137546777097j)
  1400. def test_qfunctions():
  1401. mp.dps = 15
  1402. assert qp(2,3,100).ae('2.7291482267247332183e2391')
  1403. def test_issue_239():
  1404. mp.prec = 150
  1405. x = ldexp(2476979795053773,-52)
  1406. assert betainc(206, 385, 0, 0.55, 1).ae('0.99999999999999999999996570910644857895771110649954')
  1407. mp.dps = 15
  1408. pytest.raises(ValueError, lambda: hyp2f1(-5,5,0.5,0.5))
  1409. # Extra stress testing for Bessel functions
  1410. # Reference zeros generated with the aid of scipy.special
  1411. # jn_zero, jnp_zero, yn_zero, ynp_zero
  1412. V = 15
  1413. M = 15
  1414. jn_small_zeros = \
  1415. [[2.4048255576957728,
  1416. 5.5200781102863106,
  1417. 8.6537279129110122,
  1418. 11.791534439014282,
  1419. 14.930917708487786,
  1420. 18.071063967910923,
  1421. 21.211636629879259,
  1422. 24.352471530749303,
  1423. 27.493479132040255,
  1424. 30.634606468431975,
  1425. 33.775820213573569,
  1426. 36.917098353664044,
  1427. 40.058425764628239,
  1428. 43.19979171317673,
  1429. 46.341188371661814],
  1430. [3.8317059702075123,
  1431. 7.0155866698156188,
  1432. 10.173468135062722,
  1433. 13.323691936314223,
  1434. 16.470630050877633,
  1435. 19.615858510468242,
  1436. 22.760084380592772,
  1437. 25.903672087618383,
  1438. 29.046828534916855,
  1439. 32.189679910974404,
  1440. 35.332307550083865,
  1441. 38.474766234771615,
  1442. 41.617094212814451,
  1443. 44.759318997652822,
  1444. 47.901460887185447],
  1445. [5.1356223018406826,
  1446. 8.4172441403998649,
  1447. 11.619841172149059,
  1448. 14.795951782351261,
  1449. 17.959819494987826,
  1450. 21.116997053021846,
  1451. 24.270112313573103,
  1452. 27.420573549984557,
  1453. 30.569204495516397,
  1454. 33.7165195092227,
  1455. 36.86285651128381,
  1456. 40.008446733478192,
  1457. 43.153453778371463,
  1458. 46.297996677236919,
  1459. 49.442164110416873],
  1460. [6.3801618959239835,
  1461. 9.7610231299816697,
  1462. 13.015200721698434,
  1463. 16.223466160318768,
  1464. 19.409415226435012,
  1465. 22.582729593104442,
  1466. 25.748166699294978,
  1467. 28.908350780921758,
  1468. 32.064852407097709,
  1469. 35.218670738610115,
  1470. 38.370472434756944,
  1471. 41.520719670406776,
  1472. 44.669743116617253,
  1473. 47.817785691533302,
  1474. 50.965029906205183],
  1475. [7.5883424345038044,
  1476. 11.064709488501185,
  1477. 14.37253667161759,
  1478. 17.615966049804833,
  1479. 20.826932956962388,
  1480. 24.01901952477111,
  1481. 27.199087765981251,
  1482. 30.371007667117247,
  1483. 33.537137711819223,
  1484. 36.699001128744649,
  1485. 39.857627302180889,
  1486. 43.01373772335443,
  1487. 46.167853512924375,
  1488. 49.320360686390272,
  1489. 52.471551398458023],
  1490. [8.771483815959954,
  1491. 12.338604197466944,
  1492. 15.700174079711671,
  1493. 18.980133875179921,
  1494. 22.217799896561268,
  1495. 25.430341154222704,
  1496. 28.626618307291138,
  1497. 31.811716724047763,
  1498. 34.988781294559295,
  1499. 38.159868561967132,
  1500. 41.326383254047406,
  1501. 44.489319123219673,
  1502. 47.649399806697054,
  1503. 50.80716520300633,
  1504. 53.963026558378149],
  1505. [9.9361095242176849,
  1506. 13.589290170541217,
  1507. 17.003819667816014,
  1508. 20.320789213566506,
  1509. 23.58608443558139,
  1510. 26.820151983411405,
  1511. 30.033722386570469,
  1512. 33.233041762847123,
  1513. 36.422019668258457,
  1514. 39.603239416075404,
  1515. 42.778481613199507,
  1516. 45.949015998042603,
  1517. 49.11577372476426,
  1518. 52.279453903601052,
  1519. 55.440592068853149],
  1520. [11.086370019245084,
  1521. 14.821268727013171,
  1522. 18.287582832481726,
  1523. 21.641541019848401,
  1524. 24.934927887673022,
  1525. 28.191188459483199,
  1526. 31.42279419226558,
  1527. 34.637089352069324,
  1528. 37.838717382853611,
  1529. 41.030773691585537,
  1530. 44.21540850526126,
  1531. 47.394165755570512,
  1532. 50.568184679795566,
  1533. 53.738325371963291,
  1534. 56.905249991978781],
  1535. [12.225092264004655,
  1536. 16.037774190887709,
  1537. 19.554536430997055,
  1538. 22.94517313187462,
  1539. 26.266814641176644,
  1540. 29.54565967099855,
  1541. 32.795800037341462,
  1542. 36.025615063869571,
  1543. 39.240447995178135,
  1544. 42.443887743273558,
  1545. 45.638444182199141,
  1546. 48.825930381553857,
  1547. 52.007691456686903,
  1548. 55.184747939289049,
  1549. 58.357889025269694],
  1550. [13.354300477435331,
  1551. 17.241220382489128,
  1552. 20.807047789264107,
  1553. 24.233885257750552,
  1554. 27.583748963573006,
  1555. 30.885378967696675,
  1556. 34.154377923855096,
  1557. 37.400099977156589,
  1558. 40.628553718964528,
  1559. 43.843801420337347,
  1560. 47.048700737654032,
  1561. 50.245326955305383,
  1562. 53.435227157042058,
  1563. 56.619580266508436,
  1564. 59.799301630960228],
  1565. [14.475500686554541,
  1566. 18.433463666966583,
  1567. 22.046985364697802,
  1568. 25.509450554182826,
  1569. 28.887375063530457,
  1570. 32.211856199712731,
  1571. 35.499909205373851,
  1572. 38.761807017881651,
  1573. 42.004190236671805,
  1574. 45.231574103535045,
  1575. 48.447151387269394,
  1576. 51.653251668165858,
  1577. 54.851619075963349,
  1578. 58.043587928232478,
  1579. 61.230197977292681],
  1580. [15.589847884455485,
  1581. 19.61596690396692,
  1582. 23.275853726263409,
  1583. 26.773322545509539,
  1584. 30.17906117878486,
  1585. 33.526364075588624,
  1586. 36.833571341894905,
  1587. 40.111823270954241,
  1588. 43.368360947521711,
  1589. 46.608132676274944,
  1590. 49.834653510396724,
  1591. 53.050498959135054,
  1592. 56.257604715114484,
  1593. 59.457456908388002,
  1594. 62.651217388202912],
  1595. [16.698249933848246,
  1596. 20.789906360078443,
  1597. 24.494885043881354,
  1598. 28.026709949973129,
  1599. 31.45996003531804,
  1600. 34.829986990290238,
  1601. 38.156377504681354,
  1602. 41.451092307939681,
  1603. 44.721943543191147,
  1604. 47.974293531269048,
  1605. 51.211967004101068,
  1606. 54.437776928325074,
  1607. 57.653844811906946,
  1608. 60.8618046824805,
  1609. 64.062937824850136],
  1610. [17.801435153282442,
  1611. 21.95624406783631,
  1612. 25.705103053924724,
  1613. 29.270630441874802,
  1614. 32.731053310978403,
  1615. 36.123657666448762,
  1616. 39.469206825243883,
  1617. 42.780439265447158,
  1618. 46.06571091157561,
  1619. 49.330780096443524,
  1620. 52.579769064383396,
  1621. 55.815719876305778,
  1622. 59.040934037249271,
  1623. 62.257189393731728,
  1624. 65.465883797232125],
  1625. [18.899997953174024,
  1626. 23.115778347252756,
  1627. 26.907368976182104,
  1628. 30.505950163896036,
  1629. 33.993184984781542,
  1630. 37.408185128639695,
  1631. 40.772827853501868,
  1632. 44.100590565798301,
  1633. 47.400347780543231,
  1634. 50.678236946479898,
  1635. 53.93866620912693,
  1636. 57.184898598119301,
  1637. 60.419409852130297,
  1638. 63.644117508962281,
  1639. 66.860533012260103]]
  1640. jnp_small_zeros = \
  1641. [[0.0,
  1642. 3.8317059702075123,
  1643. 7.0155866698156188,
  1644. 10.173468135062722,
  1645. 13.323691936314223,
  1646. 16.470630050877633,
  1647. 19.615858510468242,
  1648. 22.760084380592772,
  1649. 25.903672087618383,
  1650. 29.046828534916855,
  1651. 32.189679910974404,
  1652. 35.332307550083865,
  1653. 38.474766234771615,
  1654. 41.617094212814451,
  1655. 44.759318997652822],
  1656. [1.8411837813406593,
  1657. 5.3314427735250326,
  1658. 8.5363163663462858,
  1659. 11.706004902592064,
  1660. 14.863588633909033,
  1661. 18.015527862681804,
  1662. 21.16436985918879,
  1663. 24.311326857210776,
  1664. 27.457050571059246,
  1665. 30.601922972669094,
  1666. 33.746182898667383,
  1667. 36.889987409236811,
  1668. 40.033444053350675,
  1669. 43.176628965448822,
  1670. 46.319597561173912],
  1671. [3.0542369282271403,
  1672. 6.7061331941584591,
  1673. 9.9694678230875958,
  1674. 13.170370856016123,
  1675. 16.347522318321783,
  1676. 19.512912782488205,
  1677. 22.671581772477426,
  1678. 25.826037141785263,
  1679. 28.977672772993679,
  1680. 32.127327020443474,
  1681. 35.275535050674691,
  1682. 38.422654817555906,
  1683. 41.568934936074314,
  1684. 44.714553532819734,
  1685. 47.859641607992093],
  1686. [4.2011889412105285,
  1687. 8.0152365983759522,
  1688. 11.345924310743006,
  1689. 14.585848286167028,
  1690. 17.78874786606647,
  1691. 20.9724769365377,
  1692. 24.144897432909265,
  1693. 27.310057930204349,
  1694. 30.470268806290424,
  1695. 33.626949182796679,
  1696. 36.781020675464386,
  1697. 39.933108623659488,
  1698. 43.083652662375079,
  1699. 46.232971081836478,
  1700. 49.381300092370349],
  1701. [5.3175531260839944,
  1702. 9.2823962852416123,
  1703. 12.681908442638891,
  1704. 15.964107037731551,
  1705. 19.196028800048905,
  1706. 22.401032267689004,
  1707. 25.589759681386733,
  1708. 28.767836217666503,
  1709. 31.938539340972783,
  1710. 35.103916677346764,
  1711. 38.265316987088158,
  1712. 41.423666498500732,
  1713. 44.579623137359257,
  1714. 47.733667523865744,
  1715. 50.886159153182682],
  1716. [6.4156163757002403,
  1717. 10.519860873772308,
  1718. 13.9871886301403,
  1719. 17.312842487884625,
  1720. 20.575514521386888,
  1721. 23.803581476593863,
  1722. 27.01030789777772,
  1723. 30.20284907898166,
  1724. 33.385443901010121,
  1725. 36.560777686880356,
  1726. 39.730640230067416,
  1727. 42.896273163494417,
  1728. 46.058566273567043,
  1729. 49.218174614666636,
  1730. 52.375591529563596],
  1731. [7.501266144684147,
  1732. 11.734935953042708,
  1733. 15.268181461097873,
  1734. 18.637443009666202,
  1735. 21.931715017802236,
  1736. 25.183925599499626,
  1737. 28.409776362510085,
  1738. 31.617875716105035,
  1739. 34.81339298429743,
  1740. 37.999640897715301,
  1741. 41.178849474321413,
  1742. 44.352579199070217,
  1743. 47.521956905768113,
  1744. 50.687817781723741,
  1745. 53.85079463676896],
  1746. [8.5778364897140741,
  1747. 12.932386237089576,
  1748. 16.529365884366944,
  1749. 19.941853366527342,
  1750. 23.268052926457571,
  1751. 26.545032061823576,
  1752. 29.790748583196614,
  1753. 33.015178641375142,
  1754. 36.224380548787162,
  1755. 39.422274578939259,
  1756. 42.611522172286684,
  1757. 45.793999658055002,
  1758. 48.971070951900596,
  1759. 52.143752969301988,
  1760. 55.312820330403446],
  1761. [9.6474216519972168,
  1762. 14.115518907894618,
  1763. 17.774012366915256,
  1764. 21.229062622853124,
  1765. 24.587197486317681,
  1766. 27.889269427955092,
  1767. 31.155326556188325,
  1768. 34.39662855427218,
  1769. 37.620078044197086,
  1770. 40.830178681822041,
  1771. 44.030010337966153,
  1772. 47.221758471887113,
  1773. 50.407020967034367,
  1774. 53.586995435398319,
  1775. 56.762598475105272],
  1776. [10.711433970699945,
  1777. 15.28673766733295,
  1778. 19.004593537946053,
  1779. 22.501398726777283,
  1780. 25.891277276839136,
  1781. 29.218563499936081,
  1782. 32.505247352375523,
  1783. 35.763792928808799,
  1784. 39.001902811514218,
  1785. 42.224638430753279,
  1786. 45.435483097475542,
  1787. 48.636922645305525,
  1788. 51.830783925834728,
  1789. 55.01844255063594,
  1790. 58.200955824859509],
  1791. [11.770876674955582,
  1792. 16.447852748486498,
  1793. 20.223031412681701,
  1794. 23.760715860327448,
  1795. 27.182021527190532,
  1796. 30.534504754007074,
  1797. 33.841965775135715,
  1798. 37.118000423665604,
  1799. 40.371068905333891,
  1800. 43.606764901379516,
  1801. 46.828959446564562,
  1802. 50.040428970943456,
  1803. 53.243223214220535,
  1804. 56.438892058982552,
  1805. 59.628631306921512],
  1806. [12.826491228033465,
  1807. 17.600266557468326,
  1808. 21.430854238060294,
  1809. 25.008518704644261,
  1810. 28.460857279654847,
  1811. 31.838424458616998,
  1812. 35.166714427392629,
  1813. 38.460388720328256,
  1814. 41.728625562624312,
  1815. 44.977526250903469,
  1816. 48.211333836373288,
  1817. 51.433105171422278,
  1818. 54.645106240447105,
  1819. 57.849056857839799,
  1820. 61.046288512821078],
  1821. [13.878843069697276,
  1822. 18.745090916814406,
  1823. 22.629300302835503,
  1824. 26.246047773946584,
  1825. 29.72897816891134,
  1826. 33.131449953571661,
  1827. 36.480548302231658,
  1828. 39.791940718940855,
  1829. 43.075486800191012,
  1830. 46.337772104541405,
  1831. 49.583396417633095,
  1832. 52.815686826850452,
  1833. 56.037118687012179,
  1834. 59.249577075517968,
  1835. 62.454525995970462],
  1836. [14.928374492964716,
  1837. 19.88322436109951,
  1838. 23.81938909003628,
  1839. 27.474339750968247,
  1840. 30.987394331665278,
  1841. 34.414545662167183,
  1842. 37.784378506209499,
  1843. 41.113512376883377,
  1844. 44.412454519229281,
  1845. 47.688252845993366,
  1846. 50.945849245830813,
  1847. 54.188831071035124,
  1848. 57.419876154678179,
  1849. 60.641030026538746,
  1850. 63.853885828967512],
  1851. [15.975438807484321,
  1852. 21.015404934568315,
  1853. 25.001971500138194,
  1854. 28.694271223110755,
  1855. 32.236969407878118,
  1856. 35.688544091185301,
  1857. 39.078998185245057,
  1858. 42.425854432866141,
  1859. 45.740236776624833,
  1860. 49.029635055514276,
  1861. 52.299319390331728,
  1862. 55.553127779547459,
  1863. 58.793933759028134,
  1864. 62.02393848337554,
  1865. 65.244860767043859]]
  1866. yn_small_zeros = \
  1867. [[0.89357696627916752,
  1868. 3.9576784193148579,
  1869. 7.0860510603017727,
  1870. 10.222345043496417,
  1871. 13.361097473872763,
  1872. 16.500922441528091,
  1873. 19.64130970088794,
  1874. 22.782028047291559,
  1875. 25.922957653180923,
  1876. 29.064030252728398,
  1877. 32.205204116493281,
  1878. 35.346452305214321,
  1879. 38.487756653081537,
  1880. 41.629104466213808,
  1881. 44.770486607221993],
  1882. [2.197141326031017,
  1883. 5.4296810407941351,
  1884. 8.5960058683311689,
  1885. 11.749154830839881,
  1886. 14.897442128336725,
  1887. 18.043402276727856,
  1888. 21.188068934142213,
  1889. 24.331942571356912,
  1890. 27.475294980449224,
  1891. 30.618286491641115,
  1892. 33.761017796109326,
  1893. 36.90355531614295,
  1894. 40.045944640266876,
  1895. 43.188218097393211,
  1896. 46.330399250701687],
  1897. [3.3842417671495935,
  1898. 6.7938075132682675,
  1899. 10.023477979360038,
  1900. 13.209986710206416,
  1901. 16.378966558947457,
  1902. 19.539039990286384,
  1903. 22.69395593890929,
  1904. 25.845613720902269,
  1905. 28.995080395650151,
  1906. 32.143002257627551,
  1907. 35.289793869635804,
  1908. 38.435733485446343,
  1909. 41.581014867297885,
  1910. 44.725777117640461,
  1911. 47.870122696676504],
  1912. [4.5270246611496439,
  1913. 8.0975537628604907,
  1914. 11.396466739595867,
  1915. 14.623077742393873,
  1916. 17.81845523294552,
  1917. 20.997284754187761,
  1918. 24.166235758581828,
  1919. 27.328799850405162,
  1920. 30.486989604098659,
  1921. 33.642049384702463,
  1922. 36.794791029185579,
  1923. 39.945767226378749,
  1924. 43.095367507846703,
  1925. 46.2438744334407,
  1926. 49.391498015725107],
  1927. [5.6451478942208959,
  1928. 9.3616206152445429,
  1929. 12.730144474090465,
  1930. 15.999627085382479,
  1931. 19.22442895931681,
  1932. 22.424810599698521,
  1933. 25.610267054939328,
  1934. 28.785893657666548,
  1935. 31.954686680031668,
  1936. 35.118529525584828,
  1937. 38.278668089521758,
  1938. 41.435960629910073,
  1939. 44.591018225353424,
  1940. 47.744288086361052,
  1941. 50.896105199722123],
  1942. [6.7471838248710219,
  1943. 10.597176726782031,
  1944. 14.033804104911233,
  1945. 17.347086393228382,
  1946. 20.602899017175335,
  1947. 23.826536030287532,
  1948. 27.030134937138834,
  1949. 30.220335654231385,
  1950. 33.401105611047908,
  1951. 36.574972486670962,
  1952. 39.743627733020277,
  1953. 42.908248189569535,
  1954. 46.069679073215439,
  1955. 49.228543693445843,
  1956. 52.385312123112282],
  1957. [7.8377378223268716,
  1958. 11.811037107609447,
  1959. 15.313615118517857,
  1960. 18.670704965906724,
  1961. 21.958290897126571,
  1962. 25.206207715021249,
  1963. 28.429037095235496,
  1964. 31.634879502950644,
  1965. 34.828638524084437,
  1966. 38.013473399691765,
  1967. 41.19151880917741,
  1968. 44.364272633271975,
  1969. 47.53281875312084,
  1970. 50.697961822183806,
  1971. 53.860312300118388],
  1972. [8.919605734873789,
  1973. 13.007711435388313,
  1974. 16.573915129085334,
  1975. 19.974342312352426,
  1976. 23.293972585596648,
  1977. 26.5667563757203,
  1978. 29.809531451608321,
  1979. 33.031769327150685,
  1980. 36.239265816598239,
  1981. 39.435790312675323,
  1982. 42.623910919472727,
  1983. 45.805442883111651,
  1984. 48.981708325514764,
  1985. 52.153694518185572,
  1986. 55.322154420959698],
  1987. [9.9946283820824834,
  1988. 14.190361295800141,
  1989. 17.817887841179873,
  1990. 21.26093227125945,
  1991. 24.612576377421522,
  1992. 27.910524883974868,
  1993. 31.173701563441602,
  1994. 34.412862242025045,
  1995. 37.634648706110989,
  1996. 40.843415321050884,
  1997. 44.04214994542435,
  1998. 47.232978012841169,
  1999. 50.417456447370186,
  2000. 53.596753874948731,
  2001. 56.771765754432457],
  2002. [11.064090256031013,
  2003. 15.361301343575925,
  2004. 19.047949646361388,
  2005. 22.532765416313869,
  2006. 25.91620496332662,
  2007. 29.2394205079349,
  2008. 32.523270869465881,
  2009. 35.779715464475261,
  2010. 39.016196664616095,
  2011. 42.237627509803703,
  2012. 45.4474001519274,
  2013. 48.647941127433196,
  2014. 51.841036928216499,
  2015. 55.028034667184916,
  2016. 58.209970905250097],
  2017. [12.128927704415439,
  2018. 16.522284394784426,
  2019. 20.265984501212254,
  2020. 23.791669719454272,
  2021. 27.206568881574774,
  2022. 30.555020011020762,
  2023. 33.859683872746356,
  2024. 37.133649760307504,
  2025. 40.385117593813002,
  2026. 43.619533085646856,
  2027. 46.840676630553575,
  2028. 50.051265851897857,
  2029. 53.253310556711732,
  2030. 56.448332488918971,
  2031. 59.637507005589829],
  2032. [13.189846995683845,
  2033. 17.674674253171487,
  2034. 21.473493977824902,
  2035. 25.03913093040942,
  2036. 28.485081336558058,
  2037. 31.858644293774859,
  2038. 35.184165245422787,
  2039. 38.475796636190897,
  2040. 41.742455848758449,
  2041. 44.990096293791186,
  2042. 48.222870660068338,
  2043. 51.443777308699826,
  2044. 54.655042589416311,
  2045. 57.858358441436511,
  2046. 61.055036135780528],
  2047. [14.247395665073945,
  2048. 18.819555894710682,
  2049. 22.671697117872794,
  2050. 26.276375544903892,
  2051. 29.752925495549038,
  2052. 33.151412708998983,
  2053. 36.497763772987645,
  2054. 39.807134090704376,
  2055. 43.089121522203808,
  2056. 46.350163579538652,
  2057. 49.594769786270069,
  2058. 52.82620892320143,
  2059. 56.046916910756961,
  2060. 59.258751140598783,
  2061. 62.463155567737854],
  2062. [15.30200785858925,
  2063. 19.957808654258601,
  2064. 23.861599172945054,
  2065. 27.504429642227545,
  2066. 31.011103429019229,
  2067. 34.434283425782942,
  2068. 37.801385632318459,
  2069. 41.128514139788358,
  2070. 44.425913324440663,
  2071. 47.700482714581842,
  2072. 50.957073905278458,
  2073. 54.199216028087261,
  2074. 57.429547607017405,
  2075. 60.65008661807661,
  2076. 63.862406280068586],
  2077. [16.354034360047551,
  2078. 21.090156519983806,
  2079. 25.044040298785627,
  2080. 28.724161640881914,
  2081. 32.260472459522644,
  2082. 35.708083982611664,
  2083. 39.095820003878235,
  2084. 42.440684315990936,
  2085. 45.75353669045622,
  2086. 49.041718113283529,
  2087. 52.310408280968073,
  2088. 55.56338698149062,
  2089. 58.803488508906895,
  2090. 62.032886550960831,
  2091. 65.253280088312461]]
  2092. ynp_small_zeros = \
  2093. [[2.197141326031017,
  2094. 5.4296810407941351,
  2095. 8.5960058683311689,
  2096. 11.749154830839881,
  2097. 14.897442128336725,
  2098. 18.043402276727856,
  2099. 21.188068934142213,
  2100. 24.331942571356912,
  2101. 27.475294980449224,
  2102. 30.618286491641115,
  2103. 33.761017796109326,
  2104. 36.90355531614295,
  2105. 40.045944640266876,
  2106. 43.188218097393211,
  2107. 46.330399250701687],
  2108. [3.6830228565851777,
  2109. 6.9414999536541757,
  2110. 10.123404655436613,
  2111. 13.285758156782854,
  2112. 16.440058007293282,
  2113. 19.590241756629495,
  2114. 22.738034717396327,
  2115. 25.884314618788867,
  2116. 29.029575819372535,
  2117. 32.174118233366201,
  2118. 35.318134458192094,
  2119. 38.461753870997549,
  2120. 41.605066618873108,
  2121. 44.74813744908079,
  2122. 47.891014070791065],
  2123. [5.0025829314460639,
  2124. 8.3507247014130795,
  2125. 11.574195465217647,
  2126. 14.760909306207676,
  2127. 17.931285939466855,
  2128. 21.092894504412739,
  2129. 24.249231678519058,
  2130. 27.402145837145258,
  2131. 30.552708880564553,
  2132. 33.70158627151572,
  2133. 36.849213419846257,
  2134. 39.995887376143356,
  2135. 43.141817835750686,
  2136. 46.287157097544201,
  2137. 49.432018469138281],
  2138. [6.2536332084598136,
  2139. 9.6987879841487711,
  2140. 12.972409052292216,
  2141. 16.19044719506921,
  2142. 19.38238844973613,
  2143. 22.559791857764261,
  2144. 25.728213194724094,
  2145. 28.890678419054777,
  2146. 32.048984005266337,
  2147. 35.204266606440635,
  2148. 38.357281675961019,
  2149. 41.508551443818436,
  2150. 44.658448731963676,
  2151. 47.807246956681162,
  2152. 50.95515126455207],
  2153. [7.4649217367571329,
  2154. 11.005169149809189,
  2155. 14.3317235192331,
  2156. 17.58443601710272,
  2157. 20.801062338411128,
  2158. 23.997004122902644,
  2159. 27.179886689853435,
  2160. 30.353960608554323,
  2161. 33.521797098666792,
  2162. 36.685048382072301,
  2163. 39.844826969405863,
  2164. 43.001910515625288,
  2165. 46.15685955107263,
  2166. 49.310088614282257,
  2167. 52.461911043685864],
  2168. [8.6495562436971983,
  2169. 12.280868725807848,
  2170. 15.660799304540377,
  2171. 18.949739756016503,
  2172. 22.192841809428241,
  2173. 25.409072788867674,
  2174. 28.608039283077593,
  2175. 31.795195353138159,
  2176. 34.973890634255288,
  2177. 38.14630522169358,
  2178. 41.313923188794905,
  2179. 44.477791768537617,
  2180. 47.638672065035628,
  2181. 50.797131066967842,
  2182. 53.953600129601663],
  2183. [9.8147970120105779,
  2184. 13.532811875789828,
  2185. 16.965526446046053,
  2186. 20.291285512443867,
  2187. 23.56186260680065,
  2188. 26.799499736027237,
  2189. 30.015665481543419,
  2190. 33.216968050039509,
  2191. 36.407516858984748,
  2192. 39.590015243560459,
  2193. 42.766320595957378,
  2194. 45.937754257017323,
  2195. 49.105283450953203,
  2196. 52.269633324547373,
  2197. 55.431358715604255],
  2198. [10.965152105242974,
  2199. 14.765687379508912,
  2200. 18.250123150217555,
  2201. 21.612750053384621,
  2202. 24.911310600813573,
  2203. 28.171051927637585,
  2204. 31.40518108895689,
  2205. 34.621401012564177,
  2206. 37.824552065973114,
  2207. 41.017847386464902,
  2208. 44.203512240871601,
  2209. 47.3831408366063,
  2210. 50.557907466622796,
  2211. 53.728697478957026,
  2212. 56.896191727313342],
  2213. [12.103641941939539,
  2214. 15.982840905145284,
  2215. 19.517731005559611,
  2216. 22.916962141504605,
  2217. 26.243700855690533,
  2218. 29.525960140695407,
  2219. 32.778568197561124,
  2220. 36.010261572392516,
  2221. 39.226578757802172,
  2222. 42.43122493258747,
  2223. 45.626783824134354,
  2224. 48.815117837929515,
  2225. 51.997606404328863,
  2226. 55.175294723956816,
  2227. 58.348990221754937],
  2228. [13.232403808592215,
  2229. 17.186756572616758,
  2230. 20.770762917490496,
  2231. 24.206152448722253,
  2232. 27.561059462697153,
  2233. 30.866053571250639,
  2234. 34.137476603379774,
  2235. 37.385039772270268,
  2236. 40.614946085165892,
  2237. 43.831373184731238,
  2238. 47.037251786726299,
  2239. 50.234705848765229,
  2240. 53.425316228549359,
  2241. 56.610286079882087,
  2242. 59.790548623216652],
  2243. [14.35301374369987,
  2244. 18.379337301642568,
  2245. 22.011118775283494,
  2246. 25.482116178696707,
  2247. 28.865046588695164,
  2248. 32.192853922166294,
  2249. 35.483296655830277,
  2250. 38.747005493021857,
  2251. 41.990815194320955,
  2252. 45.219355876831731,
  2253. 48.435892856078888,
  2254. 51.642803925173029,
  2255. 54.84186659475857,
  2256. 58.034439083840155,
  2257. 61.221578745109862],
  2258. [15.466672066554263,
  2259. 19.562077985759503,
  2260. 23.240325531101082,
  2261. 26.746322986645901,
  2262. 30.157042415639891,
  2263. 33.507642948240263,
  2264. 36.817212798512775,
  2265. 40.097251300178642,
  2266. 43.355193847719752,
  2267. 46.596103410173672,
  2268. 49.823567279972794,
  2269. 53.040208868780832,
  2270. 56.247996968470062,
  2271. 59.448441365714251,
  2272. 62.642721301357187],
  2273. [16.574317035530872,
  2274. 20.73617763753932,
  2275. 24.459631728238804,
  2276. 27.999993668839644,
  2277. 31.438208790267783,
  2278. 34.811512070805535,
  2279. 38.140243708611251,
  2280. 41.436725143893739,
  2281. 44.708963264433333,
  2282. 47.962435051891027,
  2283. 51.201037321915983,
  2284. 54.427630745992975,
  2285. 57.644369734615238,
  2286. 60.852911791989989,
  2287. 64.054555435720397],
  2288. [17.676697936439624,
  2289. 21.9026148697762,
  2290. 25.670073356263225,
  2291. 29.244155124266438,
  2292. 32.709534477396028,
  2293. 36.105399554497548,
  2294. 39.453272918267025,
  2295. 42.766255701958017,
  2296. 46.052899215578358,
  2297. 49.319076602061401,
  2298. 52.568982147952547,
  2299. 55.805705507386287,
  2300. 59.031580956740466,
  2301. 62.248409689597653,
  2302. 65.457606670836759],
  2303. [18.774423978290318,
  2304. 23.06220035979272,
  2305. 26.872520985976736,
  2306. 30.479680663499762,
  2307. 33.971869047372436,
  2308. 37.390118854896324,
  2309. 40.757072537673599,
  2310. 44.086572292170345,
  2311. 47.387688809191869,
  2312. 50.66667461073936,
  2313. 53.928009929563275,
  2314. 57.175005343085052,
  2315. 60.410169281219877,
  2316. 63.635442539153021,
  2317. 66.85235358587768]]
  2318. @pytest.mark.slow
  2319. def test_bessel_zeros_extra():
  2320. mp.dps = 15
  2321. for v in range(V):
  2322. for m in range(1,M+1):
  2323. print(v, m, "of", V, M)
  2324. # Twice to test cache (if used)
  2325. assert besseljzero(v,m).ae(jn_small_zeros[v][m-1])
  2326. assert besseljzero(v,m).ae(jn_small_zeros[v][m-1])
  2327. assert besseljzero(v,m,1).ae(jnp_small_zeros[v][m-1])
  2328. assert besseljzero(v,m,1).ae(jnp_small_zeros[v][m-1])
  2329. assert besselyzero(v,m).ae(yn_small_zeros[v][m-1])
  2330. assert besselyzero(v,m).ae(yn_small_zeros[v][m-1])
  2331. assert besselyzero(v,m,1).ae(ynp_small_zeros[v][m-1])
  2332. assert besselyzero(v,m,1).ae(ynp_small_zeros[v][m-1])