from Numeric import *
def cdf(y, n):
    if n % 2:
        print "cdf: n must be even"
        return 0
    m = n/2
    x = y/2.0
    sum = 0.0
    for i in xrange(0, m):
        tmp = -x - logFactor(i) + i * log(x)
        if tmp > -200.0:
            sum += exp(tmp)
    # factorial version overflowed
    return 1.0 - sum 
    
def pdf(y, n):
    """
    chi_squre's pdf of a degree of freedome of n
    """
    if n % 2:
        print "cdf: n must be even"
        return 0
        
    m = n/2
    tmp = (m-1)*log(y) - (y/2) - m * log(2.0) - logFactor(m - 1)
    return exp(tmp)
    
def logFactor(n):
    if n == 0:
        return 0.0
    else:
        sum = 0.0
        for i in xrange(n):
            sum += log(i + 1)
        return sum