from mpmath import * from mpmath.libmp import round_up, from_float, mpf_zeta_int def test_zeta_int_bug(): assert mpf_zeta_int(0, 10) == from_float(-0.5) def test_bernoulli(): assert bernfrac(0) == (1,1) assert bernfrac(1) == (-1,2) assert bernfrac(2) == (1,6) assert bernfrac(3) == (0,1) assert bernfrac(4) == (-1,30) assert bernfrac(5) == (0,1) assert bernfrac(6) == (1,42) assert bernfrac(8) == (-1,30) assert bernfrac(10) == (5,66) assert bernfrac(12) == (-691,2730) assert bernfrac(18) == (43867,798) p, q = bernfrac(228) assert p % 10**10 == 164918161 assert q == 625170 p, q = bernfrac(1000) assert p % 10**10 == 7950421099 assert q == 342999030 mp.dps = 15 assert bernoulli(0) == 1 assert bernoulli(1) == -0.5 assert bernoulli(2).ae(1./6) assert bernoulli(3) == 0 assert bernoulli(4).ae(-1./30) assert bernoulli(5) == 0 assert bernoulli(6).ae(1./42) assert str(bernoulli(10)) == '0.0757575757575758' assert str(bernoulli(234)) == '7.62772793964344e+267' assert str(bernoulli(10**5)) == '-5.82229431461335e+376755' assert str(bernoulli(10**8+2)) == '1.19570355039953e+676752584' mp.dps = 50 assert str(bernoulli(10)) == '0.075757575757575757575757575757575757575757575757576' assert str(bernoulli(234)) == '7.6277279396434392486994969020496121553385863373331e+267' assert str(bernoulli(10**5)) == '-5.8222943146133508236497045360612887555320691004308e+376755' assert str(bernoulli(10**8+2)) == '1.1957035503995297272263047884604346914602088317782e+676752584' mp.dps = 1000 assert bernoulli(10).ae(mpf(5)/66) mp.dps = 50000 assert bernoulli(10).ae(mpf(5)/66) mp.dps = 15 def test_bernpoly_eulerpoly(): mp.dps = 15 assert bernpoly(0,-1).ae(1) assert bernpoly(0,0).ae(1) assert bernpoly(0,'1/2').ae(1) assert bernpoly(0,'3/4').ae(1) assert bernpoly(0,1).ae(1) assert bernpoly(0,2).ae(1) assert bernpoly(1,-1).ae('-3/2') assert bernpoly(1,0).ae('-1/2') assert bernpoly(1,'1/2').ae(0) assert bernpoly(1,'3/4').ae('1/4') assert bernpoly(1,1).ae('1/2') assert bernpoly(1,2).ae('3/2') assert bernpoly(2,-1).ae('13/6') assert bernpoly(2,0).ae('1/6') assert bernpoly(2,'1/2').ae('-1/12') assert bernpoly(2,'3/4').ae('-1/48') assert bernpoly(2,1).ae('1/6') assert bernpoly(2,2).ae('13/6') assert bernpoly(3,-1).ae(-3) assert bernpoly(3,0).ae(0) assert bernpoly(3,'1/2').ae(0) assert bernpoly(3,'3/4').ae('-3/64') assert bernpoly(3,1).ae(0) assert bernpoly(3,2).ae(3) assert bernpoly(4,-1).ae('119/30') assert bernpoly(4,0).ae('-1/30') assert bernpoly(4,'1/2').ae('7/240') assert bernpoly(4,'3/4').ae('7/3840') assert bernpoly(4,1).ae('-1/30') assert bernpoly(4,2).ae('119/30') assert bernpoly(5,-1).ae(-5) assert bernpoly(5,0).ae(0) assert bernpoly(5,'1/2').ae(0) assert bernpoly(5,'3/4').ae('25/1024') assert bernpoly(5,1).ae(0) assert bernpoly(5,2).ae(5) assert bernpoly(10,-1).ae('665/66') assert bernpoly(10,0).ae('5/66') assert bernpoly(10,'1/2').ae('-2555/33792') assert bernpoly(10,'3/4').ae('-2555/34603008') assert bernpoly(10,1).ae('5/66') assert bernpoly(10,2).ae('665/66') assert bernpoly(11,-1).ae(-11) assert bernpoly(11,0).ae(0) assert bernpoly(11,'1/2').ae(0) assert bernpoly(11,'3/4').ae('-555731/4194304') assert bernpoly(11,1).ae(0) assert bernpoly(11,2).ae(11) assert eulerpoly(0,-1).ae(1) assert eulerpoly(0,0).ae(1) assert eulerpoly(0,'1/2').ae(1) assert eulerpoly(0,'3/4').ae(1) assert eulerpoly(0,1).ae(1) assert eulerpoly(0,2).ae(1) assert eulerpoly(1,-1).ae('-3/2') assert eulerpoly(1,0).ae('-1/2') assert eulerpoly(1,'1/2').ae(0) assert eulerpoly(1,'3/4').ae('1/4') assert eulerpoly(1,1).ae('1/2') assert eulerpoly(1,2).ae('3/2') assert eulerpoly(2,-1).ae(2) assert eulerpoly(2,0).ae(0) assert eulerpoly(2,'1/2').ae('-1/4') assert eulerpoly(2,'3/4').ae('-3/16') assert eulerpoly(2,1).ae(0) assert eulerpoly(2,2).ae(2) assert eulerpoly(3,-1).ae('-9/4') assert eulerpoly(3,0).ae('1/4') assert eulerpoly(3,'1/2').ae(0) assert eulerpoly(3,'3/4').ae('-11/64') assert eulerpoly(3,1).ae('-1/4') assert eulerpoly(3,2).ae('9/4') assert eulerpoly(4,-1).ae(2) assert eulerpoly(4,0).ae(0) assert eulerpoly(4,'1/2').ae('5/16') assert eulerpoly(4,'3/4').ae('57/256') assert eulerpoly(4,1).ae(0) assert eulerpoly(4,2).ae(2) assert eulerpoly(5,-1).ae('-3/2') assert eulerpoly(5,0).ae('-1/2') assert eulerpoly(5,'1/2').ae(0) assert eulerpoly(5,'3/4').ae('361/1024') assert eulerpoly(5,1).ae('1/2') assert eulerpoly(5,2).ae('3/2') assert eulerpoly(10,-1).ae(2) assert eulerpoly(10,0).ae(0) assert eulerpoly(10,'1/2').ae('-50521/1024') assert eulerpoly(10,'3/4').ae('-36581523/1048576') assert eulerpoly(10,1).ae(0) assert eulerpoly(10,2).ae(2) assert eulerpoly(11,-1).ae('-699/4') assert eulerpoly(11,0).ae('691/4') assert eulerpoly(11,'1/2').ae(0) assert eulerpoly(11,'3/4').ae('-512343611/4194304') assert eulerpoly(11,1).ae('-691/4') assert eulerpoly(11,2).ae('699/4') # Potential accuracy issues assert bernpoly(10000,10000).ae('5.8196915936323387117e+39999') assert bernpoly(200,17.5).ae(3.8048418524583064909e244) assert eulerpoly(200,17.5).ae(-3.7309911582655785929e275) def test_gamma(): mp.dps = 15 assert gamma(0.25).ae(3.6256099082219083119) assert gamma(0.0001).ae(9999.4228832316241908) assert gamma(300).ae('1.0201917073881354535e612') assert gamma(-0.5).ae(-3.5449077018110320546) assert gamma(-7.43).ae(0.00026524416464197007186) #assert gamma(Rational(1,2)) == gamma(0.5) #assert gamma(Rational(-7,3)).ae(gamma(mpf(-7)/3)) assert gamma(1+1j).ae(0.49801566811835604271 - 0.15494982830181068512j) assert gamma(-1+0.01j).ae(-0.422733904013474115 + 99.985883082635367436j) assert gamma(20+30j).ae(-1453876687.5534810 + 1163777777.8031573j) # Should always give exact factorials when they can # be represented as mpfs under the current working precision fact = 1 for i in range(1, 18): assert gamma(i) == fact fact *= i for dps in [170, 600]: fact = 1 mp.dps = dps for i in range(1, 105): assert gamma(i) == fact fact *= i mp.dps = 100 assert gamma(0.5).ae(sqrt(pi)) mp.dps = 15 assert factorial(0) == fac(0) == 1 assert factorial(3) == 6 assert isnan(gamma(nan)) assert gamma(1100).ae('4.8579168073569433667e2866') assert rgamma(0) == 0 assert rgamma(-1) == 0 assert rgamma(2) == 1.0 assert rgamma(3) == 0.5 assert loggamma(2+8j).ae(-8.5205176753667636926 + 10.8569497125597429366j) assert loggamma('1e10000').ae('2.302485092994045684017991e10004') assert loggamma('1e10000j').ae(mpc('-1.570796326794896619231322e10000','2.302485092994045684017991e10004')) def test_fac2(): mp.dps = 15 assert [fac2(n) for n in range(10)] == [1,1,2,3,8,15,48,105,384,945] assert fac2(-5).ae(1./3) assert fac2(-11).ae(-1./945) assert fac2(50).ae(5.20469842636666623e32) assert fac2(0.5+0.75j).ae(0.81546769394688069176-0.34901016085573266889j) assert fac2(inf) == inf assert isnan(fac2(-inf)) def test_gamma_quotients(): mp.dps = 15 h = 1e-8 ep = 1e-4 G = gamma assert gammaprod([-1],[-3,-4]) == 0 assert gammaprod([-1,0],[-5]) == inf assert abs(gammaprod([-1],[-2]) - G(-1+h)/G(-2+h)) < 1e-4 assert abs(gammaprod([-4,-3],[-2,0]) - G(-4+h)*G(-3+h)/G(-2+h)/G(0+h)) < 1e-4 assert rf(3,0) == 1 assert rf(2.5,1) == 2.5 assert rf(-5,2) == 20 assert rf(j,j).ae(gamma(2*j)/gamma(j)) assert rf('-255.5815971722918','-0.5119253100282322').ae('-0.1952720278805729485') # issue 421 assert ff(-2,0) == 1 assert ff(-2,1) == -2 assert ff(4,3) == 24 assert ff(3,4) == 0 assert binomial(0,0) == 1 assert binomial(1,0) == 1 assert binomial(0,-1) == 0 assert binomial(3,2) == 3 assert binomial(5,2) == 10 assert binomial(5,3) == 10 assert binomial(5,5) == 1 assert binomial(-1,0) == 1 assert binomial(-2,-4) == 3 assert binomial(4.5, 1.5) == 6.5625 assert binomial(1100,1) == 1100 assert binomial(1100,2) == 604450 assert beta(1,1) == 1 assert beta(0,0) == inf assert beta(3,0) == inf assert beta(-1,-1) == inf assert beta(1.5,1).ae(2/3.) assert beta(1.5,2.5).ae(pi/16) assert (10**15*beta(10,100)).ae(2.3455339739604649879) assert beta(inf,inf) == 0 assert isnan(beta(-inf,inf)) assert isnan(beta(-3,inf)) assert isnan(beta(0,inf)) assert beta(inf,0.5) == beta(0.5,inf) == 0 assert beta(inf,-1.5) == inf assert beta(inf,-0.5) == -inf assert beta(1+2j,-1-j/2).ae(1.16396542451069943086+0.08511695947832914640j) assert beta(-0.5,0.5) == 0 assert beta(-3,3).ae(-1/3.) assert beta('-255.5815971722918','-0.5119253100282322').ae('18.157330562703710339') # issue 421 def test_zeta(): mp.dps = 15 assert zeta(2).ae(pi**2 / 6) assert zeta(2.0).ae(pi**2 / 6) assert zeta(mpc(2)).ae(pi**2 / 6) assert zeta(100).ae(1) assert zeta(0).ae(-0.5) assert zeta(0.5).ae(-1.46035450880958681) assert zeta(-1).ae(-mpf(1)/12) assert zeta(-2) == 0 assert zeta(-3).ae(mpf(1)/120) assert zeta(-4) == 0 assert zeta(-100) == 0 assert isnan(zeta(nan)) assert zeta(1e-30).ae(-0.5) assert zeta(-1e-30).ae(-0.5) # Zeros in the critical strip assert zeta(mpc(0.5, 14.1347251417346937904)).ae(0) assert zeta(mpc(0.5, 21.0220396387715549926)).ae(0) assert zeta(mpc(0.5, 25.0108575801456887632)).ae(0) assert zeta(mpc(1e-30,1e-40)).ae(-0.5) assert zeta(mpc(-1e-30,1e-40)).ae(-0.5) mp.dps = 50 im = '236.5242296658162058024755079556629786895294952121891237' assert zeta(mpc(0.5, im)).ae(0, 1e-46) mp.dps = 15 # Complex reflection formula assert (zeta(-60+3j) / 10**34).ae(8.6270183987866146+15.337398548226238j) # issue #358 assert zeta(0,0.5) == 0 assert zeta(0,0) == 0.5 assert zeta(0,0.5,1).ae(-0.34657359027997265) # see issue #390 assert zeta(-1.5,0.5j).ae(-0.13671400162512768475 + 0.11411333638426559139j) def test_altzeta(): mp.dps = 15 assert altzeta(-2) == 0 assert altzeta(-4) == 0 assert altzeta(-100) == 0 assert altzeta(0) == 0.5 assert altzeta(-1) == 0.25 assert altzeta(-3) == -0.125 assert altzeta(-5) == 0.25 assert altzeta(-21) == 1180529130.25 assert altzeta(1).ae(log(2)) assert altzeta(2).ae(pi**2/12) assert altzeta(10).ae(73*pi**10/6842880) assert altzeta(50) < 1 assert altzeta(60, rounding='d') < 1 assert altzeta(60, rounding='u') == 1 assert altzeta(10000, rounding='d') < 1 assert altzeta(10000, rounding='u') == 1 assert altzeta(3+0j) == altzeta(3) s = 3+4j assert altzeta(s).ae((1-2**(1-s))*zeta(s)) s = -3+4j assert altzeta(s).ae((1-2**(1-s))*zeta(s)) assert altzeta(-100.5).ae(4.58595480083585913e+108) assert altzeta(1.3).ae(0.73821404216623045) assert altzeta(1e-30).ae(0.5) assert altzeta(-1e-30).ae(0.5) assert altzeta(mpc(1e-30,1e-40)).ae(0.5) assert altzeta(mpc(-1e-30,1e-40)).ae(0.5) def test_zeta_huge(): mp.dps = 15 assert zeta(inf) == 1 mp.dps = 50 assert zeta(100).ae('1.0000000000000000000000000000007888609052210118073522') assert zeta(40*pi).ae('1.0000000000000000000000000000000000000148407238666182') mp.dps = 10000 v = zeta(33000) mp.dps = 15 assert str(v-1) == '1.02363019598118e-9934' assert zeta(pi*1000, rounding=round_up) > 1 assert zeta(3000, rounding=round_up) > 1 assert zeta(pi*1000) == 1 assert zeta(3000) == 1 def test_zeta_negative(): mp.dps = 150 a = -pi*10**40 mp.dps = 15 assert str(zeta(a)) == '2.55880492708712e+1233536161668617575553892558646631323374078' mp.dps = 50 assert str(zeta(a)) == '2.5588049270871154960875033337384432038436330847333e+1233536161668617575553892558646631323374078' mp.dps = 15 def test_polygamma(): mp.dps = 15 psi0 = lambda z: psi(0,z) psi1 = lambda z: psi(1,z) assert psi0(3) == psi(0,3) == digamma(3) #assert psi2(3) == psi(2,3) == tetragamma(3) #assert psi3(3) == psi(3,3) == pentagamma(3) assert psi0(pi).ae(0.97721330794200673) assert psi0(-pi).ae(7.8859523853854902) assert psi0(-pi+1).ae(7.5676424992016996) assert psi0(pi+j).ae(1.04224048313859376 + 0.35853686544063749j) assert psi0(-pi-j).ae(1.3404026194821986 - 2.8824392476809402j) assert findroot(psi0, 1).ae(1.4616321449683622) assert psi0(1e-10).ae(-10000000000.57722) assert psi0(1e-40).ae(-1.000000000000000e+40) assert psi0(1e-10+1e-10j).ae(-5000000000.577215 + 5000000000.000000j) assert psi0(1e-40+1e-40j).ae(-5.000000000000000e+39 + 5.000000000000000e+39j) assert psi0(inf) == inf assert psi1(inf) == 0 assert psi(2,inf) == 0 assert psi1(pi).ae(0.37424376965420049) assert psi1(-pi).ae(53.030438740085385) assert psi1(pi+j).ae(0.32935710377142464 - 0.12222163911221135j) assert psi1(-pi-j).ae(-0.30065008356019703 + 0.01149892486928227j) assert (10**6*psi(4,1+10*pi*j)).ae(-6.1491803479004446 - 0.3921316371664063j) assert psi0(1+10*pi*j).ae(3.4473994217222650 + 1.5548808324857071j) assert isnan(psi0(nan)) assert isnan(psi0(-inf)) assert psi0(-100.5).ae(4.615124601338064) assert psi0(3+0j).ae(psi0(3)) assert psi0(-100+3j).ae(4.6106071768714086321+3.1117510556817394626j) assert isnan(psi(2,mpc(0,inf))) assert isnan(psi(2,mpc(0,nan))) assert isnan(psi(2,mpc(0,-inf))) assert isnan(psi(2,mpc(1,inf))) assert isnan(psi(2,mpc(1,nan))) assert isnan(psi(2,mpc(1,-inf))) assert isnan(psi(2,mpc(inf,inf))) assert isnan(psi(2,mpc(nan,nan))) assert isnan(psi(2,mpc(-inf,-inf))) def test_polygamma_high_prec(): mp.dps = 100 assert str(psi(0,pi)) == "0.9772133079420067332920694864061823436408346099943256380095232865318105924777141317302075654362928734" assert str(psi(10,pi)) == "-12.98876181434889529310283769414222588307175962213707170773803550518307617769657562747174101900659238" def test_polygamma_identities(): mp.dps = 15 psi0 = lambda z: psi(0,z) psi1 = lambda z: psi(1,z) psi2 = lambda z: psi(2,z) assert psi0(0.5).ae(-euler-2*log(2)) assert psi0(1).ae(-euler) assert psi1(0.5).ae(0.5*pi**2) assert psi1(1).ae(pi**2/6) assert psi1(0.25).ae(pi**2 + 8*catalan) assert psi2(1).ae(-2*apery) mp.dps = 20 u = -182*apery+4*sqrt(3)*pi**3 mp.dps = 15 assert psi(2,5/6.).ae(u) assert psi(3,0.5).ae(pi**4) def test_foxtrot_identity(): # A test of the complex digamma function. # See http://mathworld.wolfram.com/FoxTrotSeries.html and # http://mathworld.wolfram.com/DigammaFunction.html psi0 = lambda z: psi(0,z) mp.dps = 50 a = (-1)**fraction(1,3) b = (-1)**fraction(2,3) x = -psi0(0.5*a) - psi0(-0.5*b) + psi0(0.5*(1+a)) + psi0(0.5*(1-b)) y = 2*pi*sech(0.5*sqrt(3)*pi) assert x.ae(y) mp.dps = 15 def test_polygamma_high_order(): mp.dps = 100 assert str(psi(50, pi)) == "-1344100348958402765749252447726432491812.641985273160531055707095989227897753035823152397679626136483" assert str(psi(50, pi + 14*e)) == "-0.00000000000000000189793739550804321623512073101895801993019919886375952881053090844591920308111549337295143780341396" assert str(psi(50, pi + 14*e*j)) == ("(-0.0000000000000000522516941152169248975225472155683565752375889510631513244785" "9377385233700094871256507814151956624433 - 0.00000000000000001813157041407010184" "702414110218205348527862196327980417757665282244728963891298080199341480881811613j)") mp.dps = 15 assert str(psi(50, pi)) == "-1.34410034895841e+39" assert str(psi(50, pi + 14*e)) == "-1.89793739550804e-18" assert str(psi(50, pi + 14*e*j)) == "(-5.2251694115217e-17 - 1.81315704140701e-17j)" def test_harmonic(): mp.dps = 15 assert harmonic(0) == 0 assert harmonic(1) == 1 assert harmonic(2) == 1.5 assert harmonic(3).ae(1. + 1./2 + 1./3) assert harmonic(10**10).ae(23.603066594891989701) assert harmonic(10**1000).ae(2303.162308658947) assert harmonic(0.5).ae(2-2*log(2)) assert harmonic(inf) == inf assert harmonic(2+0j) == 1.5+0j assert harmonic(1+2j).ae(1.4918071802755104+0.92080728264223022j) def test_gamma_huge_1(): mp.dps = 500 x = mpf(10**10) / 7 mp.dps = 15 assert str(gamma(x)) == "6.26075321389519e+12458010678" mp.dps = 50 assert str(gamma(x)) == "6.2607532138951929201303779291707455874010420783933e+12458010678" mp.dps = 15 def test_gamma_huge_2(): mp.dps = 500 x = mpf(10**100) / 19 mp.dps = 15 assert str(gamma(x)) == (\ "1.82341134776679e+5172997469323364168990133558175077136829182824042201886051511" "9656908623426021308685461258226190190661") mp.dps = 50 assert str(gamma(x)) == (\ "1.82341134776678875374414910350027596939980412984e+5172997469323364168990133558" "1750771368291828240422018860515119656908623426021308685461258226190190661") def test_gamma_huge_3(): mp.dps = 500 x = 10**80 // 3 + 10**70*j / 7 mp.dps = 15 y = gamma(x) assert str(y.real) == (\ "-6.82925203918106e+2636286142112569524501781477865238132302397236429627932441916" "056964386399485392600") assert str(y.imag) == (\ "8.54647143678418e+26362861421125695245017814778652381323023972364296279324419160" "56964386399485392600") mp.dps = 50 y = gamma(x) assert str(y.real) == (\ "-6.8292520391810548460682736226799637356016538421817e+26362861421125695245017814" "77865238132302397236429627932441916056964386399485392600") assert str(y.imag) == (\ "8.5464714367841748507479306948130687511711420234015e+263628614211256952450178147" "7865238132302397236429627932441916056964386399485392600") def test_gamma_huge_4(): x = 3200+11500j mp.dps = 15 assert str(gamma(x)) == \ "(8.95783268539713e+5164 - 1.94678798329735e+5164j)" mp.dps = 50 assert str(gamma(x)) == (\ "(8.9578326853971339570292952697675570822206567327092e+5164" " - 1.9467879832973509568895402139429643650329524144794e+51" "64j)") mp.dps = 15 def test_gamma_huge_5(): mp.dps = 500 x = 10**60 * j / 3 mp.dps = 15 y = gamma(x) assert str(y.real) == "-3.27753899634941e-227396058973640224580963937571892628368354580620654233316839" assert str(y.imag) == "-7.1519888950416e-227396058973640224580963937571892628368354580620654233316841" mp.dps = 50 y = gamma(x) assert str(y.real) == (\ "-3.2775389963494132168950056995974690946983219123935e-22739605897364022458096393" "7571892628368354580620654233316839") assert str(y.imag) == (\ "-7.1519888950415979749736749222530209713136588885897e-22739605897364022458096393" "7571892628368354580620654233316841") mp.dps = 15 def test_gamma_huge_7(): mp.dps = 100 a = 3 + j/mpf(10)**1000 mp.dps = 15 y = gamma(a) assert str(y.real) == "2.0" # wrong #assert str(y.imag) == "2.16735365342606e-1000" assert str(y.imag) == "1.84556867019693e-1000" mp.dps = 50 y = gamma(a) assert str(y.real) == "2.0" #assert str(y.imag) == "2.1673536534260596065418805612488708028522563689298e-1000" assert str(y.imag) == "1.8455686701969342787869758198351951379156813281202e-1000" def test_stieltjes(): mp.dps = 15 assert stieltjes(0).ae(+euler) mp.dps = 25 assert stieltjes(1).ae('-0.07281584548367672486058637587') assert stieltjes(2).ae('-0.009690363192872318484530386035') assert stieltjes(3).ae('0.002053834420303345866160046543') assert stieltjes(4).ae('0.002325370065467300057468170178') mp.dps = 15 assert stieltjes(1).ae(-0.07281584548367672486058637587) assert stieltjes(2).ae(-0.009690363192872318484530386035) assert stieltjes(3).ae(0.002053834420303345866160046543) assert stieltjes(4).ae(0.0023253700654673000574681701775) def test_barnesg(): mp.dps = 15 assert barnesg(0) == barnesg(-1) == 0 assert [superfac(i) for i in range(8)] == [1, 1, 2, 12, 288, 34560, 24883200, 125411328000] assert str(superfac(1000)) == '3.24570818422368e+1177245' assert isnan(barnesg(nan)) assert isnan(superfac(nan)) assert isnan(hyperfac(nan)) assert barnesg(inf) == inf assert superfac(inf) == inf assert hyperfac(inf) == inf assert isnan(superfac(-inf)) assert barnesg(0.7).ae(0.8068722730141471) assert barnesg(2+3j).ae(-0.17810213864082169+0.04504542715447838j) assert [hyperfac(n) for n in range(7)] == [1, 1, 4, 108, 27648, 86400000, 4031078400000] assert [hyperfac(n) for n in range(0,-7,-1)] == [1,1,-1,-4,108,27648,-86400000] a = barnesg(-3+0j) assert a == 0 and isinstance(a, mpc) a = hyperfac(-3+0j) assert a == -4 and isinstance(a, mpc) def test_polylog(): mp.dps = 15 zs = [mpmathify(z) for z in [0, 0.5, 0.99, 4, -0.5, -4, 1j, 3+4j]] for z in zs: assert polylog(1, z).ae(-log(1-z)) for z in zs: assert polylog(0, z).ae(z/(1-z)) for z in zs: assert polylog(-1, z).ae(z/(1-z)**2) for z in zs: assert polylog(-2, z).ae(z*(1+z)/(1-z)**3) for z in zs: assert polylog(-3, z).ae(z*(1+4*z+z**2)/(1-z)**4) assert polylog(3, 7).ae(5.3192579921456754382-5.9479244480803301023j) assert polylog(3, -7).ae(-4.5693548977219423182) assert polylog(2, 0.9).ae(1.2997147230049587252) assert polylog(2, -0.9).ae(-0.75216317921726162037) assert polylog(2, 0.9j).ae(-0.17177943786580149299+0.83598828572550503226j) assert polylog(2, 1.1).ae(1.9619991013055685931-0.2994257606855892575j) assert polylog(2, -1.1).ae(-0.89083809026228260587) assert polylog(2, 1.1*sqrt(j)).ae(0.58841571107611387722+1.09962542118827026011j) assert polylog(-2, 0.9).ae(1710) assert polylog(-2, -0.9).ae(-90/6859.) assert polylog(3, 0.9).ae(1.0496589501864398696) assert polylog(-3, 0.9).ae(48690) assert polylog(-3, -4).ae(-0.0064) assert polylog(0.5+j/3, 0.5+j/2).ae(0.31739144796565650535 + 0.99255390416556261437j) assert polylog(3+4j,1).ae(zeta(3+4j)) assert polylog(3+4j,-1).ae(-altzeta(3+4j)) # issue 390 assert polylog(1.5, -48.910886523731889).ae(-6.272992229311817) assert polylog(1.5, 200).ae(-8.349608319033686529 - 8.159694826434266042j) def test_bell_polyexp(): mp.dps = 15 # TODO: more tests for polyexp assert (polyexp(0,1e-10)*10**10).ae(1.00000000005) assert (polyexp(1,1e-10)*10**10).ae(1.0000000001) assert polyexp(5,3j).ae(-607.7044517476176454+519.962786482001476087j) assert polyexp(-1,3.5).ae(12.09537536175543444) # bell(0,x) = 1 assert bell(0,0) == 1 assert bell(0,1) == 1 assert bell(0,2) == 1 assert bell(0,inf) == 1 assert bell(0,-inf) == 1 assert isnan(bell(0,nan)) # bell(1,x) = x assert bell(1,4) == 4 assert bell(1,0) == 0 assert bell(1,inf) == inf assert bell(1,-inf) == -inf assert isnan(bell(1,nan)) # bell(2,x) = x*(1+x) assert bell(2,-1) == 0 assert bell(2,0) == 0 # large orders / arguments assert bell(10) == 115975 assert bell(10,1) == 115975 assert bell(10, -8) == 11054008 assert bell(5,-50) == -253087550 assert bell(50,-50).ae('3.4746902914629720259e74') mp.dps = 80 assert bell(50,-50) == 347469029146297202586097646631767227177164818163463279814268368579055777450 assert bell(40,50) == 5575520134721105844739265207408344706846955281965031698187656176321717550 assert bell(74) == 5006908024247925379707076470957722220463116781409659160159536981161298714301202 mp.dps = 15 assert bell(10,20j) == 7504528595600+15649605360020j # continuity of the generalization assert bell(0.5,0).ae(sinc(pi*0.5)) def test_primezeta(): mp.dps = 15 assert primezeta(0.9).ae(1.8388316154446882243 + 3.1415926535897932385j) assert primezeta(4).ae(0.076993139764246844943) assert primezeta(1) == inf assert primezeta(inf) == 0 assert isnan(primezeta(nan)) def test_rs_zeta(): mp.dps = 15 assert zeta(0.5+100000j).ae(1.0730320148577531321 + 5.7808485443635039843j) assert zeta(0.75+100000j).ae(1.837852337251873704 + 1.9988492668661145358j) assert zeta(0.5+1000000j, derivative=3).ae(1647.7744105852674733 - 1423.1270943036622097j) assert zeta(1+1000000j, derivative=3).ae(3.4085866124523582894 - 18.179184721525947301j) assert zeta(1+1000000j, derivative=1).ae(-0.10423479366985452134 - 0.74728992803359056244j) assert zeta(0.5-1000000j, derivative=1).ae(11.636804066002521459 + 17.127254072212996004j) # Additional sanity tests using fp arithmetic. # Some more high-precision tests are found in the docstrings def ae(x, y, tol=1e-6): return abs(x-y) < tol*abs(y) assert ae(fp.zeta(0.5-100000j), 1.0730320148577531321 - 5.7808485443635039843j) assert ae(fp.zeta(0.75-100000j), 1.837852337251873704 - 1.9988492668661145358j) assert ae(fp.zeta(0.5+1e6j), 0.076089069738227100006 + 2.8051021010192989554j) assert ae(fp.zeta(0.5+1e6j, derivative=1), 11.636804066002521459 - 17.127254072212996004j) assert ae(fp.zeta(1+1e6j), 0.94738726251047891048 + 0.59421999312091832833j) assert ae(fp.zeta(1+1e6j, derivative=1), -0.10423479366985452134 - 0.74728992803359056244j) assert ae(fp.zeta(0.5+100000j, derivative=1), 10.766962036817482375 - 30.92705282105996714j) assert ae(fp.zeta(0.5+100000j, derivative=2), -119.40515625740538429 + 217.14780631141830251j) assert ae(fp.zeta(0.5+100000j, derivative=3), 1129.7550282628460881 - 1685.4736895169690346j) assert ae(fp.zeta(0.5+100000j, derivative=4), -10407.160819314958615 + 13777.786698628045085j) assert ae(fp.zeta(0.75+100000j, derivative=1), -0.41742276699594321475 - 6.4453816275049955949j) assert ae(fp.zeta(0.75+100000j, derivative=2), -9.214314279161977266 + 35.07290795337967899j) assert ae(fp.zeta(0.75+100000j, derivative=3), 110.61331857820103469 - 236.87847130518129926j) assert ae(fp.zeta(0.75+100000j, derivative=4), -1054.334275898559401 + 1769.9177890161596383j) def test_siegelz(): mp.dps = 15 assert siegelz(100000).ae(5.87959246868176504171) assert siegelz(100000, derivative=2).ae(-54.1172711010126452832) assert siegelz(100000, derivative=3).ae(-278.930831343966552538) assert siegelz(100000+j,derivative=1).ae(678.214511857070283307-379.742160779916375413j) def test_zeta_near_1(): # Test for a former bug in mpf_zeta and mpc_zeta mp.dps = 15 s1 = fadd(1, '1e-10', exact=True) s2 = fadd(1, '-1e-10', exact=True) s3 = fadd(1, '1e-10j', exact=True) assert zeta(s1).ae(1.000000000057721566490881444e10) assert zeta(s2).ae(-9.99999999942278433510574872e9) z = zeta(s3) assert z.real.ae(0.57721566490153286060) assert z.imag.ae(-9.9999999999999999999927184e9) mp.dps = 30 s1 = fadd(1, '1e-50', exact=True) s2 = fadd(1, '-1e-50', exact=True) s3 = fadd(1, '1e-50j', exact=True) assert zeta(s1).ae('1e50') assert zeta(s2).ae('-1e50') z = zeta(s3) assert z.real.ae('0.57721566490153286060651209008240243104215933593992') assert z.imag.ae('-1e50')