"""High-precision reference values for reproduce.py.

P(T <= x) for a noncentral t variable T with df degrees of freedom and noncentrality nc equals
  integral over v > 0 of Phi(x * sqrt(v / df) - nc) * f_chi2(v; df) dv.
The integral is evaluated with mpmath at 50 and 60 digits using different interval splits; the two runs
agree to all digits printed.
Run: python reference.py   (mpmath 1.4.1)
"""
import mpmath as mp

CASES = [(3, 10, -4.0), (5, 10, -4.0), (5, 10, -2.0), (30, 6, -4.0), (100, 6, -4.0), (4, 12.25, -2.776)]


def nct_cdf(x, df, nc, dps, splits):
    mp.mp.dps = dps
    x, df, nc = mp.mpf(x), mp.mpf(df), mp.mpf(nc)
    f = lambda v: mp.ncdf(x * mp.sqrt(v / df) - nc) * mp.exp((df / 2 - 1) * mp.log(v) - v / 2 - (df / 2) * mp.log(2) - mp.loggamma(df / 2))
    return mp.quad(f, [0] + [s * df for s in splits] + [mp.inf])


for df, nc, x in CASES:
    a = nct_cdf(x, df, nc, 50, [1, 10])
    b = nct_cdf(x, df, nc, 60, [0.25, 1, 4, 20])
    print(f"df={df:<4} nc={nc:<6} x={x:<7} reference={mp.nstr(b, 12)}  (50-digit run: {mp.nstr(a, 12)})")