extratest_gamma.py

from mpmath import *
from mpmath.libmp import ifac

import sys
if "-dps" in sys.argv:
    maxdps = int(sys.argv[sys.argv.index("-dps")+1])
else:
    maxdps = 1000

raise_ = "-raise" in sys.argv

errcount = 0

def check(name, func, z, y):
    global errcount
    try:
        x = func(z)
    except:
        errcount += 1
        if raise_:
            raise
        print()
        print(name)
        print("EXCEPTION")
        import traceback
        traceback.print_tb(sys.exc_info()[2])
        print()
        return
    xre = x.real
    xim = x.imag
    yre = y.real
    yim = y.imag
    tol = eps*8
    err = 0
    if abs(xre-yre) > abs(yre)*tol:
        err = 1
        print()
        print("Error! %s (re = %s, wanted %s, err=%s)" % (name, nstr(xre,10), nstr(yre,10), nstr(abs(xre-yre))))
        errcount += 1
        if raise_:
            raise SystemExit
    if abs(xim-yim) > abs(yim)*tol:
        err = 1
        print()
        print("Error! %s (im = %s, wanted %s, err=%s)" % (name, nstr(xim,10), nstr(yim,10), nstr(abs(xim-yim))))
        errcount += 1
        if raise_:
            raise SystemExit
    if not err:
        sys.stdout.write("%s ok; " % name)

def testcase(case):
    z, result = case
    print("Testing z =", z)
    mp.dps = 1010
    z = eval(z)
    mp.dps = maxdps + 50
    if result is None:
        gamma_val = gamma(z)
        loggamma_val = loggamma(z)
        factorial_val = factorial(z)
        rgamma_val = rgamma(z)
    else:
        loggamma_val = eval(result)
        gamma_val = exp(loggamma_val)
        factorial_val = z * gamma_val
        rgamma_val = 1/gamma_val
    for dps in [5, 10, 15, 25, 40, 60, 90, 120, 250, 600, 1000, 1800, 3600]:
        if dps > maxdps:
            break
        mp.dps = dps
        print("dps = %s" % dps)
        check("gamma", gamma, z, gamma_val)
        check("rgamma", rgamma, z, rgamma_val)
        check("loggamma", loggamma, z, loggamma_val)
        check("factorial", factorial, z, factorial_val)
        print()
        mp.dps = 15

testcases = []

# Basic values
for n in list(range(1,200)) + list(range(201,2000,17)):
    testcases.append(["%s" % n, None])
for n in range(-200,200):
    testcases.append(["%s+0.5" % n, None])
    testcases.append(["%s+0.37" % n, None])

testcases += [\
["(0.1+1j)", None],
["(-0.1+1j)", None],
["(0.1-1j)", None],
["(-0.1-1j)", None],
["10j", None],
["-10j", None],
["100j", None],
["10000j", None],
["-10000000j", None],
["(10**100)*j", None],
["125+(10**100)*j", None],
["-125+(10**100)*j", None],
["(10**10)*(1+j)", None],
["(10**10)*(-1+j)", None],
["(10**100)*(1+j)", None],
["(10**100)*(-1+j)", None],
["(1.5-1j)", None],
["(6+4j)", None],
["(4+1j)", None],
["(3.5+2j)", None],
["(1.5-1j)", None],
["(-6-4j)", None],
["(-2-3j)", None],
["(-2.5-2j)", None],
["(4+1j)", None],
["(3+3j)", None],
["(2-2j)", None],
["1", "0"],
["2", "0"],
["3", "log(2)"],
["4", "log(6)"],
["5", "log(24)"],
["0.5", "log(pi)/2"],
["1.5", "log(sqrt(pi)/2)"],
["2.5", "log(3*sqrt(pi)/4)"],
["mpf('0.37')", None],
["0.25", "log(sqrt(2*sqrt(2*pi**3)/agm(1,sqrt(2))))"],
["-0.4", None],
["mpf('-1.9')", None],
["mpf('12.8')", None],
["mpf('33.7')", None],
["mpf('95.2')", None],
["mpf('160.3')", None],
["mpf('2057.8')", None],
["25", "log(ifac(24))"],
["80", "log(ifac(79))"],
["500", "log(ifac(500-1))"],
["8000", "log(ifac(8000-1))"],
["8000.5", None],
["mpf('8000.1')", None],
["mpf('1.37e10')", None],
["mpf('1.37e10')*(1+j)", None],
["mpf('1.37e10')*(-1+j)", None],
["mpf('1.37e10')*(-1-j)", None],
["mpf('1.37e10')*(-1+j)", None],
["mpf('1.37e100')", None],
["mpf('1.37e100')*(1+j)", None],
["mpf('1.37e100')*(-1+j)", None],
["mpf('1.37e100')*(-1-j)", None],
["mpf('1.37e100')*(-1+j)", None],
["3+4j",
"mpc('"
"-1.7566267846037841105306041816232757851567066070613445016197619371316057169"
"4723618263960834804618463052988607348289672535780644470689771115236512106002"
"5970873471563240537307638968509556191696167970488390423963867031934333890838"
"8009531786948197210025029725361069435208930363494971027388382086721660805397"
"9163230643216054580167976201709951509519218635460317367338612500626714783631"
"7498317478048447525674016344322545858832610325861086336204591943822302971823"
"5161814175530618223688296232894588415495615809337292518431903058265147109853"
"1710568942184987827643886816200452860853873815413367529829631430146227470517"
"6579967222200868632179482214312673161276976117132204633283806161971389519137"
"1243359764435612951384238091232760634271570950240717650166551484551654327989"
"9360285030081716934130446150245110557038117075172576825490035434069388648124"
"6678152254554001586736120762641422590778766100376515737713938521275749049949"
"1284143906816424244705094759339932733567910991920631339597278805393743140853"
"391550313363278558195609260225928','"
"4.74266443803465792819488940755002274088830335171164611359052405215840070271"
"5906813009373171139767051863542508136875688550817670379002790304870822775498"
"2809996675877564504192565392367259119610438951593128982646945990372179860613"
"4294436498090428077839141927485901735557543641049637962003652638924845391650"
"9546290137755550107224907606529385248390667634297183361902055842228798984200"
"9591180450211798341715874477629099687609819466457990642030707080894518168924"
"6805549314043258530272479246115112769957368212585759640878745385160943755234"
"9398036774908108204370323896757543121853650025529763655312360354244898913463"
"7115955702828838923393113618205074162812089732064414530813087483533203244056"
"0546577484241423134079056537777170351934430586103623577814746004431994179990"
"5318522939077992613855205801498201930221975721246498720895122345420698451980"
"0051215797310305885845964334761831751370672996984756815410977750799748813563"
"8784405288158432214886648743541773208808731479748217023665577802702269468013"
"673719173759245720489020315779001')"],
]

for z in [4, 14, 34, 64]:
    testcases.append(["(2+j)*%s/3" % z, None])
    testcases.append(["(-2+j)*%s/3" % z, None])
    testcases.append(["(1+2*j)*%s/3" % z, None])
    testcases.append(["(2-j)*%s/3" % z, None])
    testcases.append(["(20+j)*%s/3" % z, None])
    testcases.append(["(-20+j)*%s/3" % z, None])
    testcases.append(["(1+20*j)*%s/3" % z, None])
    testcases.append(["(20-j)*%s/3" % z, None])
    testcases.append(["(200+j)*%s/3" % z, None])
    testcases.append(["(-200+j)*%s/3" % z, None])
    testcases.append(["(1+200*j)*%s/3" % z, None])
    testcases.append(["(200-j)*%s/3" % z, None])

# Poles
for n in [0,1,2,3,4,25,-1,-2,-3,-4,-20,-21,-50,-51,-200,-201,-20000,-20001]:
    for t in ['1e-5', '1e-20', '1e-100', '1e-10000']:
        testcases.append(["fadd(%s,'%s',exact=True)" % (n, t), None])
        testcases.append(["fsub(%s,'%s',exact=True)" % (n, t), None])
        testcases.append(["fadd(%s,'%sj',exact=True)" % (n, t), None])
        testcases.append(["fsub(%s,'%sj',exact=True)" % (n, t), None])

if __name__ == "__main__":
    from timeit import default_timer as clock
    tot_time = 0.0
    for case in testcases:
        t1 = clock()
        testcase(case)
        t2 = clock()
        print("Test time:", t2-t1)
        print()
        tot_time += (t2-t1)
    print("Total time:", tot_time)
    print("Errors:", errcount)