import mpmath as mp

mp.dps = 256  # increase if needed

def foxh(z, a, A, b, B, m, n):
    """
    Compute the Fox H-function H_{p,q}^{m,n}(z) with integer scale parameters Ai and Bj
    based on formula (28) using mpmath's meijer G-function.

    Parameters:
    - z : argument at which to evaluate H
    - a : list of length p containing ai
    - A : list of length p containing integer Ai
    - b : list of length q containing bj
    - B : list of length q containing integer Bj
    - m : integer, number of bj terms in numerator
    - n : integer, number of ai terms in denominator

    Returns:
    - value of H_{p,q}^{m,n}(z)
    """
    # Validate integer scale parameters
    if not all(isinstance(Ai, int) and Ai > 0 for Ai in A):
        raise ValueError("All Ai must be positive integers")
    if not all(isinstance(Bj, int) and Bj > 0 for Bj in B):
        raise ValueError("All Bj must be positive integers")

    p = len(a)
    q = len(b)

    # Build the expanded parameter lists for the Meijer G-function
    a_tilde = []
    for ai, Ai in zip(a, A):
        for k in range(Ai):
            a_tilde.append(ai + k / Ai)

    b_tilde = []
    for bj, Bj in zip(b, B):
        for k in range(Bj):
            b_tilde.append(bj + k / Bj)

    # Compute beta = prod(B_j^B_j) / prod(A_i^A_i)
    beta = mp.mpf(1)
    for Ai in A:
        beta /= Ai**Ai
    for Bj in B:
        beta *= Bj**Bj

    # Compute m_tilde and n_tilde (lengths for numerator/denominator lists)
    m_tilde = sum(B[:m])
    n_tilde = sum(A[:n])

    # Compute a_star and c_star
    a_star = sum(A[:n]) - sum(A[n:]) + sum(B[:m]) - sum(B[m:])
    c_star = m + n - (p + q) / 2

    # Compute M factor = prod(B_j^(b_j - 1/2)) / prod(A_i^(a_i - 1/2))
    M = mp.mpf(1)
    for bj, Bj in zip(b, B):
        M *= Bj**(bj - mp.mpf(1)/2)
    for ai, Ai in zip(a, A):
        M /= Ai**(ai - mp.mpf(1)/2)

    # Prefactor from formula (28)
    prefactor = (2 * mp.pi)**(c_star - a_star/2) * M

    # Evaluate the Meijer G-function
    return prefactor * mp.meijerg(
        [a_tilde[:n_tilde], a_tilde[n_tilde:]],
        [b_tilde[:m_tilde], b_tilde[m_tilde:]],
        z / beta
    )

def chisq1(z):
    prefactor = mp.sqrt(1 / (2 * mp.pi)) * z**(-1/2)
    return prefactor * mp.exp(-z / 2)

def chisq1_foxH(z):
    b = [-1/2]
    B = [2]
    w = z * z / (4)
    prefactor = mp.mpf(1) / (mp.sqrt(mp.pi))
    return prefactor * foxh(w, [], [], b, B, m=1, n=0)

print(chisq1_foxH(4.0))  # Example for chi-squared distribution with 1 degree of freedom
print(chisq1(4.0))