# Various 'special' functions of mathematics and statistics
# Most are based on details and algorithms given in
# "Numerical Recipes - Art Of Scientific Computing" by Press, et al.
#
#	gamma, beta, factorial functions
#	incomplete gamma functions
#	error function and it's complement
#	standard and generalised univariate normal distribution
#
# nb: routines assume array base is 1
# mws, February, 1992.

echo "Defining gamma, factorial and related functions"
silent

ROOT2PI = sqrt(2*PI)
array gammcoef[6] = {
	76.18009173, -86.50532033,
	24.01409822, -1.231739516,
	0.120858003e-2, -0.536382e-5
}

# ln(gamma(z))
func gammaln(z) = {
	local tmp, ser, j
	z -= 1
	tmp = z+5.5; tmp = (z+0.5)*ln(tmp)-tmp
	ser = 1
	for (j = 1; j < 7; j += 1) {
		z += 1
		ser += gammcoef[j]/z
	}
	tmp+ln(ROOT2PI*ser)
}

func gamma(z) = exp(gammaln(z))

array facttab[33] = {1}		#table for speed
facttop = 0			#filled as required
func fact(n) = {
	local j
	if (n < 0) error(1)
	if (n <= facttop) {
		facttab[n+1]	# implicit return (last expr evaluated)
	} else if (n <= 32) {
		for (j = facttop+1; j <= n; j += 1)
			facttab[j+1] = j*facttab[j]
		facttop = n
		facttab[n+1]	# implicit return
	} else gamma(n+1)
}

func beta(z,w) = exp(gammaln(z) + gammaln(w) - gammaln(z+w))


GITMAX = 100	# max iterations to calculate gamma functions
GEPS = 3e-7	#

func gamser(a,x) = {		# evaluate P(a,x) by series
	local gln,ap,sum,del,n

	gln = gammaln(a)
	if (x == 0) return 0 else {
		ap = a
		sum = 1/a		
		del = sum
		for (n = 1; n <= GITMAX; n += 1) {
			ap += 1
			del *= x/ap
			sum += del
			if (abs(del) < abs(sum)*GEPS)
				return sum*exp(-x+a*ln(x)-gln)
		}
		error(2)	# can't reach desired accuracy
	}
}

func gamcf(a,x) = {		# evaluate Q(a,x) by continued fraction
	local gln,gold,a0,a1,b0,b1,factor,n,na,nf,g

 	gln = gammaln(a)
	gold = 0
	a0 = 1; a1 = x; b0 = 0; b1 = 1
	factor = 1
	for (n = 1; n <= GITMAX; n += 1) {
		na = n-a
		a0 = (a1+a0*na)*factor
		b0 = (b1+b0*na)*factor
		nf = n*factor
		a1 = x*a0+nf*a1
		b1 = x*b0+nf*b1
		if (a1 != 0) {
			factor = 1/a1
			g = b1*factor
			if ((abs(g-gold)/g) < GEPS)
				return exp(-x+a*ln(x)-gln)*g
			gold = g
		}
	}
	error(3)	# can't reach desired accuracy
}

# incomplete gamma function P(a,x) (by standard symbology)
func gammap(a,x) = {
	if (x < 0 || a <= 0) error(4)	# invalid parameters
	if (x < a+1) gamser(a,x) else 1-gamcf(a,x)
}

# and the other, Q(a,x) = 1 - P(a,x)
func gammaq(a,x) = 1-gammap(a,x)

# and the error function erf(x) and its complement erfc(x)
func erf(x) = x < 0 ? -gammap(0.5,sqr(x)) : gammap(0.5,sqr(x))
func erfc(x) = x < 0 ? 1+gammap(0.5,sqr(x)) : gammaq(0.5,sqr(x))

# BUT a (much) quicker calculation of erfc (and hence erf), with fractional
# error less than 1.2e-7 is:
func Erfc(x) = {
	local z,t,rv
	z = abs(x)
	t = 1/(1+0.5*z)
	rv = t*exp(-z*z-1.26551223+t*(1.00002368+t*(0.37409196+ \
		t*(0.09678418+t*(-0.18628806+t*(0.27886807+t*(-1.13520398+ \
		t*(1.48851587+t*(-0.82215223+t*0.17087277)))))))))

	x >= 0 ? rv : 2-rv
}

func Erf(x) = 1-Erfc(x)

# cumulative normal distribution functions
ROOT2 = sqrt(2)

# standard normal - mean 0, variance 1
func snorm(x) = 1 - 0.5*Erfc(x/ROOT2)

# generalised normal - mean mu, variance sigma
func gnorm(x,mu,sigma) = 1 - 0.5*Erfc((x-mu)/(sigma*ROOT2))

verbose		# turn messages back on
echo "Done."