cephes_incbeta.jule

// Copyright 2025 mertcandav.
// Use of this source code is governed by a BSD 3-Clause
// license that can be found in the LICENSE file.

// Implementation derived from the gonum.
//
// Copyright ©2016 The Gonum Authors. All rights reserved.
// Use of this source code is governed by a BSD-style
// license that can be found in the LICENSE file.

/*
 * Cephes Math Library, Release 2.3:  March, 1995
 * Copyright 1984, 1995 by Stephen L. Moshier
 */

use "std/math"

const (
	maxGam = 171.624376956302725
	big    = 4.503599627370496e15
	biginv = 2.22044604925031308085e-16
)

// Computes the regularized incomplete beta function.
fn incbet(aa: f64, bb: f64, xx: f64): f64 {
	if aa <= 0 || bb <= 0 {
		panic(paramOutOfBounds)
	}
	if xx <= 0 || xx >= 1 {
		if xx == 0 {
			return 0
		}
		if xx == 1 {
			return 1
		}
		panic(paramOutOfBounds)
	}

	let mut flag: int
	if bb*xx <= 1 && xx <= 0.95 {
		t := pseries(aa, bb, xx)
		return transformT(t, flag)
	}

	mut w := 1 - xx

	// Reverse a and b if x is greater than the mean.
	let mut a: f64
	let mut b: f64
	let mut xc: f64
	let mut x: f64
	if xx > aa/(aa+bb) {
		flag = 1
		a = bb
		b = aa
		xc = xx
		x = w
	} else {
		a = aa
		b = bb
		xc = w
		x = xx
	}

	if flag == 1 && (b*x) <= 1.0 && x <= 0.95 {
		t := pseries(a, b, x)
		return transformT(t, flag)
	}

	// Choose expansion for better convergence.
	mut y := x*(a+b-2.0) - (a - 1.0)
	if y < 0.0 {
		w = incbcf(a, b, x)
	} else {
		w = incbd(a, b, x) / xc
	}

	// Multiply w by the factor
	// x^a * (1-x)^b * Γ(a+b) / (a*Γ(a)*Γ(b))
	let mut t: f64
	y = a * math::Log(x)
	t = b * math::Log(xc)
	if (a+b) < maxGam && math::Abs(y) < maxLog && math::Abs(t) < maxLog {
		t = math::Pow(xc, b)
		t *= math::Pow(x, a)
		t /= a
		t *= w
		t *= 1.0 / Beta(a, b)
		return transformT(t, flag)
	}

	// Resort to logarithms.
	y += t - Lbeta(a, b)
	y += math::Log(w / a)
	if y < minLog {
		t = 0.0
	} else {
		t = math::Exp(y)
	}

	return transformT(t, flag)
}

fn transformT(mut t: f64, flag: int): f64 {
	if flag == 1 {
		if t <= machEp {
			t = 1.0 - machEp
		} else {
			t = 1.0 - t
		}
	}
	return t
}

// Returns the incomplete beta integral evaluated by a continued fraction expansion.
fn incbcf(a: f64, b: f64, x: f64): f64 {
	let mut xk: f64
	let mut pk: f64
	let mut pkm1: f64
	let mut pkm2: f64
	let mut qk: f64
	let mut qkm1: f64
	let mut qkm2: f64
	let mut k1: f64
	let mut k2: f64
	let mut k3: f64
	let mut k4: f64
	let mut k5: f64
	let mut k6: f64
	let mut k7: f64
	let mut k8: f64
	let mut r: f64
	let mut t: f64
	let mut ans: f64
	let mut thresh: f64
	let mut n: int

	k1 = a
	k2 = a + b
	k3 = a
	k4 = a + 1.0
	k5 = 1.0
	k6 = b - 1.0
	k7 = k4
	k8 = a + 2.0

	pkm2 = 0.0
	qkm2 = 1.0
	pkm1 = 1.0
	qkm1 = 1.0
	ans = 1.0
	r = 1.0
	thresh = 3.0 * machEp

	n = 0
	for n <= 300; n++ {
		xk = -(x * k1 * k2) / (k3 * k4)
		pk = pkm1 + pkm2*xk
		qk = qkm1 + qkm2*xk
		pkm2 = pkm1
		pkm1 = pk
		qkm2 = qkm1
		qkm1 = qk

		xk = (x * k5 * k6) / (k7 * k8)
		pk = pkm1 + pkm2*xk
		qk = qkm1 + qkm2*xk
		pkm2 = pkm1
		pkm1 = pk
		qkm2 = qkm1
		qkm1 = qk

		if qk != 0 {
			r = pk / qk
		}
		if r != 0 {
			t = math::Abs((ans - r) / r)
			ans = r
		} else {
			t = 1.0
		}

		if t < thresh {
			return ans
		}

		k1 += 1.0
		k2 += 1.0
		k3 += 2.0
		k4 += 2.0
		k5 += 1.0
		k6 -= 1.0
		k7 += 2.0
		k8 += 2.0

		if (math::Abs(qk)+math::Abs(pk)) > big {
			pkm2 *= biginv
			pkm1 *= biginv
			qkm2 *= biginv
			qkm1 *= biginv
		}
		if (math::Abs(qk) < biginv) || (math::Abs(pk) < biginv) {
			pkm2 *= big
			pkm1 *= big
			qkm2 *= big
			qkm1 *= big
		}
	}

	return ans
}

// Returns the incomplete beta integral evaluated by a continued fraction expansion.
fn incbd(a: f64, b: f64, x: f64): f64 {
	let mut xk: f64
	let mut pk: f64
	let mut pkm1: f64
	let mut pkm2: f64
	let mut qk: f64
	let mut qkm1: f64
	let mut qkm2: f64
	let mut k1: f64
	let mut k2: f64
	let mut k3: f64
	let mut k4: f64
	let mut k5: f64
	let mut k6: f64
	let mut k7: f64
	let mut k8: f64
	let mut r: f64
	let mut t: f64
	let mut ans: f64
	let mut z: f64
	let mut thresh: f64
	let mut n: int

	k1 = a
	k2 = b - 1.0
	k3 = a
	k4 = a + 1.0
	k5 = 1.0
	k6 = a + b
	k7 = a + 1.0
	k8 = a + 2.0

	pkm2 = 0.0
	qkm2 = 1.0
	pkm1 = 1.0
	qkm1 = 1.0
	z = x / (1.0 - x)
	ans = 1.0
	r = 1.0
	thresh = 3.0 * machEp
	n = 0
	for n <= 300; n++ {
		xk = -(z * k1 * k2) / (k3 * k4)
		pk = pkm1 + pkm2*xk
		qk = qkm1 + qkm2*xk
		pkm2 = pkm1
		pkm1 = pk
		qkm2 = qkm1
		qkm1 = qk

		xk = (z * k5 * k6) / (k7 * k8)
		pk = pkm1 + pkm2*xk
		qk = qkm1 + qkm2*xk
		pkm2 = pkm1
		pkm1 = pk
		qkm2 = qkm1
		qkm1 = qk

		if qk != 0 {
			r = pk / qk
		}
		if r != 0 {
			t = math::Abs((ans - r) / r)
			ans = r
		} else {
			t = 1.0
		}

		if t < thresh {
			return ans
		}

		k1 += 1.0
		k2 -= 1.0
		k3 += 2.0
		k4 += 2.0
		k5 += 1.0
		k6 += 1.0
		k7 += 2.0
		k8 += 2.0

		if (math::Abs(qk)+math::Abs(pk)) > big {
			pkm2 *= biginv
			pkm1 *= biginv
			qkm2 *= biginv
			qkm1 *= biginv
		}
		if (math::Abs(qk) < biginv) || (math::Abs(pk) < biginv) {
			pkm2 *= big
			pkm1 *= big
			qkm2 *= big
			qkm1 *= big
		}
	}
	return ans
}

// Returns the incomplete beta integral evaluated by a power series. Use
// when b*x is small and x not too close to 1.
fn pseries(a: f64, b: f64, x: f64): f64 {
	let mut s: f64
	let mut t: f64
	let mut u: f64
	let mut v: f64
	let mut n: f64
	let mut t1: f64
	let mut z: f64
	let mut ai: f64
	ai = 1.0 / a
	u = (1.0 - b) * x
	v = u / (a + 1.0)
	t1 = v
	t = u
	n = 2.0
	s = 0.0
	z = machEp * ai
	for math::Abs(v) > z {
		u = (n - b) * x / n
		t *= u
		v = t / (a + n)
		s += v
		n += 1.0
	}
	s += t1
	s += ai

	u = a * math::Log(x)
	if (a+b) < maxGam && math::Abs(u) < maxLog {
		t = 1.0 / Beta(a, b)
		s = s * t * math::Pow(x, a)
	} else {
		t = -Lbeta(a, b) + u + math::Log(s)
		if t < minLog {
			s = 0.0
		} else {
			s = math::Exp(t)
		}
	}
	return s
}