cossim.py

def expect_cossim(n, approx=True):
    """
    returns gamma(n / 2) / gamma((n + 1) / 2) / sqrt(pi) for positive integers.
    """
    assert n > 0
    
    if approx:
        # abs(error) < 1e-8
        if n >= 128000:
            return np.reciprocal(np.sqrt(np.pi / 2 * n + 1))
        elif n >= 80:
            poly = (2.4674010 * n - 2.4673232) * n + 1.2274046
            return np.reciprocal(np.sqrt(np.sqrt(poly)))
        # fall-through when it's faster just to compute iteratively.

    even = n % 2 == 0
    res = 2.0 if even else 1.0
    for i in range(2 if even else 1, n, 2):
        res *= i / (i + 1)
    return res / np.pi if even else res