from math import lgamma, exp, log
 
# Computes P(a, x)
# Regularized lower incomplete gamma function.
def lowergamma(a, x, itmax=100, eps=3.0E-7, fpmin=1.0E-30):
    return 1 - uppergamma(a, x, itmax, eps, fpmin)
 
# Computes Q(a, x)
# Regularized upper incomplete gamma function.
# @see Numerical Recipes in Fortran77 (Ch. 6.2 Special Functions).
def uppergamma(a, x, itmax=100, eps=3.0E-7, fpmin=1.0E-30):
    assert x >= 0 and a > 0, 'arguments zero or negative'
    b=x+1-a
    c=1/fpmin
    d=1/b
    h=d
    for i in range(1, 1+itmax):
        an=-i*(i-a)
        b=b+2
        d=an*d+b
        if abs(d) < fpmin: d=fpmin
        c=b+an/c
        if abs(c) < fpmin: c=fpmin
        d=1/d
        dl=d*c
        h=h*dl
        if abs(dl-1) < eps:
            return h*exp(-x+a*log(x)-lgamma(a))
    raise ValueError('a too large, itmax too small')
 
def chisqcdf(x, k):
    return lowergamma(k/2, x/2)
 
def pvalue(x, k):
    return 1 - chisqcdf(x, k)
 
print(round(pvalue(0.71, 4), 2)) # 0.95
print(round(pvalue(3.36, 4), 2)) # 0.5
print(round(pvalue(9.49, 4), 2)) # 0.05
 
print(round(pvalue(3.66, 3), 2)) # 0.3
print(round(pvalue(6.06, 5), 2)) # 0.3
print(round(pvalue(3.83, 6), 2)) # 0.7
 
print(round(pvalue( 124.352,  100), 2)) # 0.05
print(round(pvalue( 969.484, 1000), 2)) # 0.75
print(round(pvalue( 999.333, 1000), 2)) # 0.5
print(round(pvalue(1029.790, 1000), 2)) # 0.25
print(round(pvalue(1057.724, 1000), 2)) # 0.1

ZnJvbSBtYXRoIGltcG9ydCBsZ2FtbWEsIGV4cCwgbG9nCgojIENvbXB1dGVzIFAoYSwgeCkKIyBSZWd1bGFyaXplZCBsb3dlciBpbmNvbXBsZXRlIGdhbW1hIGZ1bmN0aW9uLgpkZWYgbG93ZXJnYW1tYShhLCB4LCBpdG1heD0xMDAsIGVwcz0zLjBFLTcsIGZwbWluPTEuMEUtMzApOgogICAgcmV0dXJuIDEgLSB1cHBlcmdhbW1hKGEsIHgsIGl0bWF4LCBlcHMsIGZwbWluKQoKIyBDb21wdXRlcyBRKGEsIHgpCiMgUmVndWxhcml6ZWQgdXBwZXIgaW5jb21wbGV0ZSBnYW1tYSBmdW5jdGlvbi4KIyBAc2VlIE51bWVyaWNhbCBSZWNpcGVzIGluIEZvcnRyYW43NyAoQ2guIDYuMiBTcGVjaWFsIEZ1bmN0aW9ucykuCmRlZiB1cHBlcmdhbW1hKGEsIHgsIGl0bWF4PTEwMCwgZXBzPTMuMEUtNywgZnBtaW49MS4wRS0zMCk6CiAgICBhc3NlcnQgeCA+PSAwIGFuZCBhID4gMCwgJ2FyZ3VtZW50cyB6ZXJvIG9yIG5lZ2F0aXZlJwogICAgYj14KzEtYQogICAgYz0xL2ZwbWluCiAgICBkPTEvYgogICAgaD1kCiAgICBmb3IgaSBpbiByYW5nZSgxLCAxK2l0bWF4KToKICAgICAgICBhbj0taSooaS1hKQogICAgICAgIGI9YisyCiAgICAgICAgZD1hbipkK2IKICAgICAgICBpZiBhYnMoZCkgPCBmcG1pbjogZD1mcG1pbgogICAgICAgIGM9Yithbi9jCiAgICAgICAgaWYgYWJzKGMpIDwgZnBtaW46IGM9ZnBtaW4KICAgICAgICBkPTEvZAogICAgICAgIGRsPWQqYwogICAgICAgIGg9aCpkbAogICAgICAgIGlmIGFicyhkbC0xKSA8IGVwczoKICAgICAgICAgICAgcmV0dXJuIGgqZXhwKC14K2EqbG9nKHgpLWxnYW1tYShhKSkKICAgIHJhaXNlIFZhbHVlRXJyb3IoJ2EgdG9vIGxhcmdlLCBpdG1heCB0b28gc21hbGwnKQoKZGVmIGNoaXNxY2RmKHgsIGspOgogICAgcmV0dXJuIGxvd2VyZ2FtbWEoay8yLCB4LzIpCgpkZWYgcHZhbHVlKHgsIGspOgogICAgcmV0dXJuIDEgLSBjaGlzcWNkZih4LCBrKQoKcHJpbnQocm91bmQocHZhbHVlKDAuNzEsIDQpLCAyKSkgIyAwLjk1CnByaW50KHJvdW5kKHB2YWx1ZSgzLjM2LCA0KSwgMikpICMgMC41CnByaW50KHJvdW5kKHB2YWx1ZSg5LjQ5LCA0KSwgMikpICMgMC4wNQoKcHJpbnQocm91bmQocHZhbHVlKDMuNjYsIDMpLCAyKSkgIyAwLjMKcHJpbnQocm91bmQocHZhbHVlKDYuMDYsIDUpLCAyKSkgIyAwLjMKcHJpbnQocm91bmQocHZhbHVlKDMuODMsIDYpLCAyKSkgIyAwLjcKCnByaW50KHJvdW5kKHB2YWx1ZSggMTI0LjM1MiwgIDEwMCksIDIpKSAjIDAuMDUKcHJpbnQocm91bmQocHZhbHVlKCA5NjkuNDg0LCAxMDAwKSwgMikpICMgMC43NQpwcmludChyb3VuZChwdmFsdWUoIDk5OS4zMzMsIDEwMDApLCAyKSkgIyAwLjUKcHJpbnQocm91bmQocHZhbHVlKDEwMjkuNzkwLCAxMDAwKSwgMikpICMgMC4yNQpwcmludChyb3VuZChwdmFsdWUoMTA1Ny43MjQsIDEwMDApLCAyKSkgIyAwLjE=