// Copyright (c) 2026 Petr Balvín <opensource@petrbalvin.org> (https://petrbalvin.org)
// SPDX-License-Identifier: MIT

package stats

import "sourcedock.dev/petrbalvin/tensor/internal/base"

import (
	"math"
)

// Distribution functions: cumulative distribution functions and their
// quantiles for the distributions the library draws from, so a
// p-value or a confidence interval needs no second library.
//
// The foundations are the regularised incomplete gamma and beta
// functions: every continuous CDF here reduces to one of them, and
// the two discrete CDFs reduce through the exact identities
// P(Poisson ≤ k) = ΓUpper(k+1, λ) and P(Binomial ≤ k) =
// BetaIncomplete(1−p, n−k, k+1), which hold in floating point to the
// accuracy of the incomplete functions themselves. Quantiles invert
// the CDF by a bracketed Newton iteration through the distribution's
// density, falling back to bisection whenever the derivative step
// would leave the bracket, so the convergence stays unconditional on
// every monotone CDF.

// gammaSeriesEps bounds the series and continued-fraction iterations
// of the incomplete functions.
const gammaSeriesEps = 3e-16

// gammaIterations bounds the series and continued-fraction iterations
// for a shape a. Both converge in O(√a) rounds near x ≈ a, so a fixed
// cap would silently truncate the large-shape answers (a shape of
// 50000 at x ≈ a needs about 1800 series rounds); the budget grows
// with √a and every caller errors when it is still not enough.
func gammaIterations(a float64) int {
	return 1000 + int(20*math.Sqrt(a))
}

// GammaLower returns the regularised lower incomplete gamma
// P(a, x) = γ(a, x)/Γ(a), the CDF of a Gamma(shape a, rate 1) draw.
// It evaluates the power series below x < a+1 and the continued
// fraction of the complement above, each to rounding level; for large
// shapes the series budget grows with √a, and either iteration that
// fails to converge inside its budget is an error rather than a
// silently truncated value.
func GammaLower(a, x float64) (float64, error) {
	p, _, err := gammaIncomplete("GammaLower", a, x)
	return p, err
}

// GammaUpper returns the regularised upper incomplete gamma
// Q(a, x) = Γ(a, x)/Γ(a) = 1 − P(a, x), computed on the same
// foundation as GammaLower.
func GammaUpper(a, x float64) (float64, error) {
	_, q, err := gammaIncomplete("GammaUpper", a, x)
	return q, err
}

// gammaIncomplete returns P(a, x) and Q(a, x) together, sharing the
// one exponent factor both need.
func gammaIncomplete(name string, a, x float64) (p, q float64, err error) {
	if !(a > 0) {
		return 0, 0, base.Errf("%s: shape a must be positive, got %g", name, a)
	}
	if !(x >= 0) {
		return 0, 0, base.Errf("%s: x must not be negative, got %g", name, x)
	}
	if x == 0 {
		return 0, 1, nil
	}
	if math.IsInf(x, 1) {
		return 1, 0, nil
	}
	lg, _ := math.Lgamma(a)
	front := math.Exp(-x + a*math.Log(x) - lg)
	// The iteration budget is fixed by the shape alone, so it is
	// computed once: the same bound the inline call evaluated, without
	// re-deriving it on every round.
	rounds := gammaIterations(a)
	if x < a+1 {
		// Power series for P: γ(a, x) = e^{−x}x^a Σ x^n/(a(a+1)…(a+n)).
		// Every term is positive, so the sum carries no cancellation and
		// stays accurate for arbitrarily large shapes given enough
		// rounds.
		sum := 1 / a
		del := sum
		ap := a
		converged := false
		for range rounds {
			ap++
			del *= x / ap
			sum += del
			if math.Abs(del) < math.Abs(sum)*gammaSeriesEps {
				converged = true
				break
			}
		}
		if !converged {
			return 0, 0, base.Errf("%s: the power series did not converge for a = %g, x = %g", name, a, x)
		}
		return sum * front, 1 - sum*front, nil
	}
	// Continued fraction for Q (Numerical Recipes form), Lentz's
	// modified method with the tiny floor.
	const fpmin = 1e-300
	b := x + 1 - a
	c := 1 / fpmin
	d := 1 / b
	h := d
	converged := false
	for i := 1; i <= rounds; i++ {
		an := -float64(i) * (float64(i) - a)
		b += 2
		d = an*d + b
		if math.Abs(d) < fpmin {
			d = fpmin
		}
		c = b + an/c
		if math.Abs(c) < fpmin {
			c = fpmin
		}
		d = 1 / d
		del := d * c
		h *= del
		if math.Abs(del-1) < gammaSeriesEps {
			converged = true
			break
		}
	}
	if !converged {
		return 0, 0, base.Errf("%s: the continued fraction did not converge for a = %g, x = %g", name, a, x)
	}
	return 1 - front*h, front * h, nil
}

// BetaIncomplete returns the regularised incomplete beta function
// I_x(a, b), the CDF of a Beta(a, b) draw, by the continued fraction
// with the standard symmetry switch at x > (a+1)/(a+b+2).
func BetaIncomplete(x, a, b float64) (float64, error) {
	if !(a > 0 && b > 0) {
		return 0, base.Errf("BetaIncomplete: shape parameters must be positive, got %g, %g", a, b)
	}
	if !(x >= 0 && x <= 1) {
		return 0, base.Errf("BetaIncomplete: x must lie in [0, 1], got %g", x)
	}
	if x == 0 || x == 1 {
		return x, nil
	}
	lab, _ := math.Lgamma(a + b)
	la, _ := math.Lgamma(a)
	lb, _ := math.Lgamma(b)
	front := math.Exp(lab - la - lb + a*math.Log(x) + b*math.Log(1-x))
	// The continued fraction converges fastest away from the skew side.
	if x < (a+1)/(a+b+2) {
		h, err := betaContinue(a, b, x)
		if err != nil {
			return 0, base.Errf("BetaIncomplete: %w", err)
		}
		return front * h / a, nil
	}
	h, err := betaContinue(b, a, 1-x)
	if err != nil {
		return 0, base.Errf("BetaIncomplete: %w", err)
	}
	return 1 - front*h/b, nil
}

// betaContinue evaluates the incomplete-beta continued fraction
// (Numerical Recipes form) by Lentz's modified method. The round
// budget grows with the shapes exactly as the incomplete gamma's does,
// and an iteration that fails to converge inside it is an error rather
// than a silently truncated value.
func betaContinue(a, b, x float64) (float64, error) {
	const fpmin = 1e-300
	rounds := max(1000, gammaIterations(a)+gammaIterations(b))
	qab := a + b
	qap := a + 1
	qam := a - 1
	c := 1.0
	d := 1 - qab*x/qap
	if math.Abs(d) < fpmin {
		d = fpmin
	}
	d = 1 / d
	h := d
	converged := false
	for i := 1; i <= rounds; i++ {
		m := float64(i)
		m2 := 2 * m
		aa := m * (b - m) * x / ((qam + m2) * (a + m2))
		d = 1 + aa*d
		if math.Abs(d) < fpmin {
			d = fpmin
		}
		c = 1 + aa/c
		if math.Abs(c) < fpmin {
			c = fpmin
		}
		d = 1 / d
		h *= d * c
		aa = -(a + m) * (qab + m) * x / ((a + m2) * (qap + m2))
		d = 1 + aa*d
		if math.Abs(d) < fpmin {
			d = fpmin
		}
		c = 1 + aa/c
		if math.Abs(c) < fpmin {
			c = fpmin
		}
		d = 1 / d
		del := d * c
		h *= del
		if math.Abs(del-1) < gammaSeriesEps {
			converged = true
			break
		}
	}
	if !converged {
		return 0, base.Errf("the continued fraction did not converge for a = %g, b = %g, x = %g", a, b, x)
	}
	return h, nil
}

// NormalCDF returns Φ(x), the standard normal cumulative distribution.
func NormalCDF(x float64) float64 {
	return 0.5 * math.Erfc(-x/math.Sqrt2)
}

// ExponentialCDF returns P(X ≤ x) for X ~ Exponential(rate).
func ExponentialCDF(x, rate float64) (float64, error) {
	if !(rate > 0) {
		return 0, base.Errf("ExponentialCDF: rate must be positive, got %g", rate)
	}
	if math.IsNaN(x) {
		// The x ≤ 0 branch below is the left tail, not a rejection, so
		// NaN needs its own guard: without it the answer is a silent NaN.
		return 0, base.Errf("ExponentialCDF: x must be a number, got %g", x)
	}
	if x <= 0 {
		return 0, nil
	}
	// −Expm1(−rate·x) is 1 − e^{−rate·x} without the cancellation: the
	// literal form loses the left tail, an 11 % relative error already
	// at rate·x = 1e-16 and a total loss not far below, while Expm1
	// keeps every bit of the small answer.
	return -math.Expm1(-rate * x), nil
}

// GammaCDF returns P(X ≤ x) for X ~ Gamma(shape, rate), the same
// parametrisation GammaDraws samples.
func GammaCDF(x, shape, rate float64) (float64, error) {
	if !(shape > 0 && rate > 0) {
		return 0, base.Errf("GammaCDF: shape and rate must be positive, got %g, %g", shape, rate)
	}
	return GammaLower(shape, x*rate)
}

// ChiSquareCDF returns P(X ≤ x) for X ~ χ²(df).
func ChiSquareCDF(x float64, df int) (float64, error) {
	if df < 1 {
		return 0, base.Errf("ChiSquareCDF: df must be ≥ 1, got %d", df)
	}
	return GammaLower(float64(df)/2, x/2)
}

// StudentTCDF returns P(T ≤ t) for T ~ Student t(df).
func StudentTCDF(t float64, df int) (float64, error) {
	if df < 1 {
		return 0, base.Errf("StudentTCDF: df must be ≥ 1, got %d", df)
	}
	// One house tail: the upper-tail helper carries the asymptotic forms
	// the heavy df ≤ 2 laws keep past t²'s overflow, where the closed
	// form's z = df/(df+t²) collapses to 0 and the lower tail answered a
	// silent 0 for a tail the format still holds.
	upper, err := studentTUpperTail(math.Abs(t), df)
	if err != nil {
		return 0, base.Errf("StudentTCDF: %w", err)
	}
	if t >= 0 {
		return 1 - upper, nil
	}
	return upper, nil
}

// PoissonCDF returns P(N ≤ k) for N ~ Poisson(lambda), through the
// identity P(N ≤ k) = ΓUpper(k+1, λ).
func PoissonCDF(k int, lambda float64) (float64, error) {
	if !(lambda > 0) {
		return 0, base.Errf("PoissonCDF: lambda must be positive, got %g", lambda)
	}
	if k < 0 {
		return 0, nil
	}
	return GammaUpper(float64(k+1), lambda)
}

// BinomialCDF returns P(X ≤ k) for X ~ Binomial(trials, p), through
// the identity P(X ≤ k) = I_{1−p}(n−k, k+1).
func BinomialCDF(k, trials int, p float64) (float64, error) {
	if trials < 1 {
		return 0, base.Errf("BinomialCDF: trials must be ≥ 1, got %d", trials)
	}
	if !(p > 0 && p < 1) {
		return 0, base.Errf("BinomialCDF: p must lie in (0, 1), got %g", p)
	}
	if k < 0 {
		return 0, nil
	}
	if k >= trials {
		return 1, nil
	}
	return BetaIncomplete(1-p, float64(trials-k), float64(k+1))
}

// continuousQuantile inverts a continuous CDF on the positive axis:
// the bracket starts at seed and doubles outward until the CDF
// straddles q, then the crossing is refined by a bracketed Newton
// iteration through the CDF's derivative pdf when one is given, and
// by plain bisection when it is not. The Newton walk carries the
// bracket with it and bisects whenever the derivative step would
// leave the bracket or would not cut it fast enough, so the
// convergence stays unconditional on every monotone CDF either way.
func continuousQuantile(name string, q float64, seed float64,
	cdf func(float64) (float64, error), pdf func(float64) float64) (float64, error) {
	// The guard is NaN-rejecting on purpose: NaN compares false against
	// both bounds, so the written-out form would let it past and every
	// bracket comparison below would then also be false.
	if !(q >= 0 && q <= 1) {
		return 0, base.Errf("%s: q must lie in [0, 1], got %g", name, q)
	}
	if q == 0 || q == 1 {
		return 0, base.Errf("%s: q = %g has no finite quantile", name, q)
	}
	lo, hi := seed, seed
	fLo, err := cdf(lo)
	if err != nil {
		return 0, base.Errf("%s: %w", name, err)
	}
	fHi, err := cdf(hi)
	if err != nil {
		return 0, base.Errf("%s: %w", name, err)
	}
	// Establish a bracket where the CDF crosses q.
	for fLo > q {
		lo /= 2
		// Halving from a finite seed underflows through the denormals
		// to exactly zero, which is the whole representable range: no
		// CDF value can sit below a probability of that size.
		if lo == 0 {
			return 0, base.Errf("%s: failed to bracket q = %g from below", name, q)
		}
		if fLo, err = cdf(lo); err != nil {
			return 0, base.Errf("%s: %w", name, err)
		}
	}
	for fHi < q {
		hi *= 2
		if math.IsInf(hi, 0) {
			return 0, base.Errf("%s: failed to bracket q = %g from above", name, q)
		}
		if fHi, err = cdf(hi); err != nil {
			return 0, base.Errf("%s: %w", name, err)
		}
	}
	if pdf == nil {
		// Halve to rounding level. The interval at least halves every pass
		// and the loop leaves on mid == lo || mid == hi, which a finite
		// bracket always reaches (about 1100 passes from the widest one),
		// so the iteration count is a guard against a non-monotone CDF,
		// not a working part of the convergence: hitting it is an error,
		// never a quietly unconverged midpoint.
		converged := false
		for range 4096 {
			mid := (lo + hi) / 2
			if mid == lo || mid == hi {
				converged = true
				break
			}
			f, err := cdf(mid)
			if err != nil {
				return 0, base.Errf("%s: %w", name, err)
			}
			if f < q {
				lo = mid
			} else {
				hi = mid
			}
		}
		if !converged {
			return 0, base.Errf("%s: the bisection for q = %g did not converge", name, q)
		}
		return (lo + hi) / 2, nil
	}
	return continuousQuantileNewton(name, q, lo, hi, cdf, pdf)
}

// continuousQuantileNewton refines a bracket where the CDF crosses q
// by a safeguarded Newton walk: a Newton step on cdf(x) = q through
// the pdf derivative whenever that step lands inside the bracket and
// cuts it at least as fast as the bisection it replaces, a bisection
// step otherwise. The bisection guarantee survives untouched: the
// bracket at least halves every second round, the walk leaves when the
// bracket has no float strictly inside it, the same rounding-level
// exit the plain bisection takes, and a vanishing or non-finite
// derivative fails both Newton guards and bisects.
func continuousQuantileNewton(name string, q float64, lo, hi float64,
	cdf func(float64) (float64, error), pdf func(float64) float64) (float64, error) {
	x := 0.5 * (lo + hi)
	f, err := cdf(x)
	if err != nil {
		return 0, base.Errf("%s: %w", name, err)
	}
	if f == q {
		return x, nil
	}
	if f < q {
		lo = x
	} else {
		hi = x
	}
	// stepSize carries the previous round's step, the yardstick the
	// second Newton guard measures against: a derivative step that does
	// not halve it is stalling, and the round bisects instead.
	stepSize := hi - lo
	converged := false
	for range 4096 {
		d := pdf(x)
		g := f - q
		// The two rtsafe guards: the derivative step is taken only when
		// it lands strictly inside the bracket, x - g/d in (lo, hi),
		// and cuts the last step at least in half. A vanishing or
		// non-finite derivative fails the first guard and bisects.
		newton := d > 0 && !math.IsInf(d, 1) &&
			(x-lo)*d > g && g > (x-hi)*d &&
			2*math.Abs(g) <= math.Abs(stepSize*d)
		if newton {
			stepSize = g / d
			x -= stepSize
		} else {
			stepSize = 0.5 * (hi - lo)
			x = 0.5 * (lo + hi)
		}
		if x == lo || x == hi {
			converged = true
			break
		}
		if f, err = cdf(x); err != nil {
			return 0, base.Errf("%s: %w", name, err)
		}
		if f == q {
			converged = true
			break
		}
		if f < q {
			lo = x
		} else {
			hi = x
		}
	}
	if !converged {
		return 0, base.Errf("%s: the quantile iteration for q = %g did not converge", name, q)
	}
	return x, nil
}

// continuousQuantileUpper inverts a monotone decreasing upper-tail
// probability on the positive axis: it returns the z > 0 with
// upper(z) = q, for q below what the reflected CDF can represent
// (below 2⁻⁵³, where 1−q rounds to 1). The bracket starts at [0, 1]
// and doubles outward; upper(0) = 0.5 covers every such q.
func continuousQuantileUpper(name string, q float64,
	upper func(float64) (float64, error)) (float64, error) {
	lo, hi := 0.0, 1.0
	fHi, err := upper(hi)
	if err != nil {
		return 0, base.Errf("%s: %w", name, err)
	}
	for fHi > q {
		lo = hi
		hi *= 2
		if math.IsInf(hi, 0) {
			return 0, base.Errf("%s: failed to bracket q = %g from above", name, q)
		}
		if fHi, err = upper(hi); err != nil {
			return 0, base.Errf("%s: %w", name, err)
		}
	}
	converged := false
	for range 4096 {
		mid := (lo + hi) / 2
		if mid == lo || mid == hi {
			converged = true
			break
		}
		f, err := upper(mid)
		if err != nil {
			return 0, base.Errf("%s: %w", name, err)
		}
		if f > q {
			lo = mid
		} else {
			hi = mid
		}
	}
	if !converged {
		return 0, base.Errf("%s: the bisection for q = %g did not converge", name, q)
	}
	return (lo + hi) / 2, nil
}

// discreteQuantile returns the smallest integer k whose CDF reaches
// q, found by doubling the upper bound then bisecting the integer
// grid. The result is a whole number in float form.
func discreteQuantile(name string, q float64, cdf func(int) (float64, error)) (float64, error) {
	// NaN-rejecting for the same reason continuousQuantile's guard is:
	// a NaN q makes every comparison below false, and the doubling
	// bracket would then run without end.
	if !(q >= 0 && q <= 1) {
		return 0, base.Errf("%s: q must lie in [0, 1], got %g", name, q)
	}
	if q == 0 {
		return 0, nil
	}
	hi := 1
	for {
		f, err := cdf(hi)
		if err != nil {
			return 0, base.Errf("%s: %w", name, err)
		}
		if f >= q {
			break
		}
		// A monotone CDF reaches any q > 0 at some finite k, so a bound
		// this high means the CDF never will: refuse rather than let hi
		// wrap to zero and the doubling loop spin for ever.
		if hi > math.MaxInt/2 {
			return 0, base.Errf("%s: failed to bracket q = %g", name, q)
		}
		hi *= 2
	}
	lo := 0
	for lo < hi {
		mid := (lo + hi) / 2
		f, err := cdf(mid)
		if err != nil {
			return 0, base.Errf("%s: %w", name, err)
		}
		if f >= q {
			hi = mid
		} else {
			lo = mid + 1
		}
	}
	return float64(lo), nil
}

// normalPdf is the standard normal density, the derivative
// NormalQuantile's Newton walk steps through.
func normalPdf(x float64) float64 {
	return math.Exp(-0.5*x*x) / math.Sqrt(2*math.Pi)
}

// gammaShapeRatePdf is the Gamma(shape, rate) density at x, evaluated
// from its logarithm so a large shape's normalising factor never
// overflows on its way to a density the format cannot hold. x at or
// below zero answers 0: the quantile walk treats a vanishing
// derivative as a signal to bisect.
func gammaShapeRatePdf(x, shape, rate float64) float64 {
	if x <= 0 {
		return 0
	}
	return math.Exp((shape-1)*math.Log(x) - rate*x + shape*math.Log(rate) - logGamma(shape))
}

// studentTPdf is the Student t(df) density at x, from its logarithm
// like the gamma density, so the wide tails' underflow lands on the 0
// the quantile walk reads as a bisection signal.
func studentTPdf(x float64, df int) float64 {
	d := float64(df)
	return math.Exp(logGamma(0.5*(d+1)) - logGamma(0.5*d) -
		0.5*math.Log(d*math.Pi) - 0.5*(d+1)*math.Log1p(x*x/d))
}

// NormalQuantile returns the q-quantile of the standard normal
// distribution, the inverse of NormalCDF.
func NormalQuantile(q float64) (float64, error) {
	if !(q >= 0 && q <= 1) {
		return 0, base.Errf("NormalQuantile: q must lie in [0, 1], got %g", q)
	}
	if q == 0 || q == 1 {
		return 0, base.Errf("NormalQuantile: q = %g has no finite quantile", q)
	}
	// The bracket expansion below only walks up from the seed, and
	// Phi(0) = 0.5 is the infimum it can reach, so the whole lower half
	// of the domain is bracketed on the reflected probability instead:
	// Phi is symmetric, so the lower-tail quantile is the mirror of the
	// upper one. Below 2⁻⁵³ the reflection 1−q rounds to exactly 1 and
	// would be read as the q = 1 refusal, so the extreme tail bisects
	// the accurate upper tail 0.5·Erfc directly, on the positive axis.
	if q < 0.5 {
		if 1-q == 1 {
			z, err := continuousQuantileUpper("NormalQuantile", q, func(x float64) (float64, error) {
				return 0.5 * math.Erfc(x/math.Sqrt2), nil
			})
			if err != nil {
				return 0, err
			}
			return -z, nil
		}
		z, err := continuousQuantile("NormalQuantile", 1-q, 1, func(x float64) (float64, error) {
			return NormalCDF(x), nil
		}, normalPdf)
		if err != nil {
			return 0, err
		}
		return -z, nil
	}
	return continuousQuantile("NormalQuantile", q, 1, func(x float64) (float64, error) {
		return NormalCDF(x), nil
	}, normalPdf)
}

// ExponentialQuantile returns the q-quantile of Exponential(rate).
func ExponentialQuantile(q, rate float64) (float64, error) {
	return continuousQuantile("ExponentialQuantile", q, 1, func(x float64) (float64, error) {
		return ExponentialCDF(x, rate)
	}, func(x float64) float64 {
		return rate * math.Exp(-rate*x)
	})
}

// GammaQuantile returns the q-quantile of Gamma(shape, rate).
func GammaQuantile(q, shape, rate float64) (float64, error) {
	return continuousQuantile("GammaQuantile", q, shape/rate, func(x float64) (float64, error) {
		return GammaCDF(x, shape, rate)
	}, func(x float64) float64 {
		return gammaShapeRatePdf(x, shape, rate)
	})
}

// ChiSquareQuantile returns the q-quantile of χ²(df), the critical
// value tables quote.
func ChiSquareQuantile(q float64, df int) (float64, error) {
	shape := float64(df) / 2
	return continuousQuantile("ChiSquareQuantile", q, float64(df), func(x float64) (float64, error) {
		return ChiSquareCDF(x, df)
	}, func(x float64) float64 {
		return gammaShapeRatePdf(x, shape, 0.5)
	})
}

// StudentTQuantile returns the q-quantile of Student t(df). The signed
// axis is bracketed through a two-sided bracket expansion, the crossing
// refined by the bracketed Newton walk through the t density.
func StudentTQuantile(q float64, df int) (float64, error) {
	if !(q >= 0 && q <= 1) {
		return 0, base.Errf("StudentTQuantile: q must lie in [0, 1], got %g", q)
	}
	if q == 0 || q == 1 {
		return 0, base.Errf("StudentTQuantile: q = %g has no finite quantile", q)
	}
	// Bracket on the positive side, then mirror for q < 0.5. Below 2⁻⁵³
	// the reflection 1−q rounds to exactly 1 and would be read as the
	// q = 1 refusal, so the extreme tail bisects the one-sided upper
	// tail directly, on the positive axis.
	if q < 0.5 {
		if 1-q == 1 {
			t, err := continuousQuantileUpper("StudentTQuantile", q, func(x float64) (float64, error) {
				return studentTUpperTail(x, df)
			})
			if err != nil {
				return 0, err
			}
			return -t, nil
		}
		t, err := continuousQuantile("StudentTQuantile", 1-q, 1, func(x float64) (float64, error) {
			return StudentTCDF(x, df)
		}, func(x float64) float64 {
			return studentTPdf(x, df)
		})
		if err != nil {
			return 0, err
		}
		return -t, nil
	}
	return continuousQuantile("StudentTQuantile", q, 1, func(x float64) (float64, error) {
		return StudentTCDF(x, df)
	}, func(x float64) float64 {
		return studentTPdf(x, df)
	})
}

// studentTUpperTail returns P(T > t) for Student t(df), one tail of
// the symmetric two-sided form twoSidedT uses. Beyond t²'s overflow
// the incomplete-beta argument saturates to 0 and would report a
// phantom zero tail exactly where the heavy df ≤ 2 tails stay
// representable, so those degrees answer by their asymptotic forms.
func studentTUpperTail(t float64, df int) (float64, error) {
	if t > 1.3e154 {
		switch {
		case df == 1:
			// Cauchy: the exact tail 1/2 − atan(t)/π, in a form the
			// huge t keeps accurate.
			return 1 / (math.Pi * t), nil
		case df == 2:
			// The tail 1/(2t²) divided one t at a time: the literal
			// denominator overflows past √MaxFloat64, and the quotient
			// would flush the still-representable tail to zero.
			return 0.5 / t / t, nil
		}
	}
	z := float64(df) / (float64(df) + t*t)
	p, err := BetaIncomplete(z, float64(df)/2, 0.5)
	if err != nil {
		return 0, err
	}
	return p / 2, nil
}

// PoissonQuantile returns the smallest k with P(N ≤ k) ≥ q for
// N ~ Poisson(lambda).
func PoissonQuantile(q, lambda float64) (float64, error) {
	if !(lambda > 0) {
		return 0, base.Errf("PoissonQuantile: lambda must be positive, got %g", lambda)
	}
	return discreteQuantile("PoissonQuantile", q, func(k int) (float64, error) {
		return PoissonCDF(k, lambda)
	})
}

// BinomialQuantile returns the smallest k with P(X ≤ k) ≥ q for
// X ~ Binomial(trials, p).
func BinomialQuantile(q, p float64, trials int) (float64, error) {
	if !(p > 0 && p < 1) {
		return 0, base.Errf("BinomialQuantile: p must lie in (0, 1), got %g", p)
	}
	return discreteQuantile("BinomialQuantile", q, func(k int) (float64, error) {
		return BinomialCDF(k, trials, p)
	})
}