cscript
mata:
/*
        ALGORITHM AS 5  APPL. STATIST. (1968) VOL.17, P.193

        COMPUTES LOWER TAIL AREA OF NON-CENTRAL T-DISTRIBUTION
*/

real scalar prncst(real scalar st, real scalar idf, 
		   real scalar d, real scalar ifault)
{
	real scalar     g1, g2, g3
	real scalar     zero, one, two, half
	real scalar     eps, emin
	real scalar     f, a, b, rb, da, drb, fmkm1, fmkm2
	real scalar     sum, ak, fk, fkm1
/*
alnorm() is 66
tfn() is ??  may be 27?
alogam is ??  may be 32?
*/

/*
 
         CONSTANTS - G1 IS 1.0 / SQRT(2.0 * PI)
                     G2 IS 1.0 / (2.0 * PI)
                     G3 IS SQRT(2.0 * PI)
*/
	g1 = 0.3989422804
	g2 = 0.1591549431
	g3 = 2.5066282746

	zero = 0
	one = 1
	two = 2
	half = .5
	eps = 1.0e-6
	emin = 12.5
 
/*
      ZSQRT(A) = SQRT(A)
      ZEXP(A) = EXP(A)
*/
	f = idf
	if (idf>100) goto L50
	ifault = 0
	ioe = mod(idf, 2)
	a = st/sqrt(f)
	b = f / (f+st^2)
	rb = sqrt(b)
	da = d*a
	drb = d*rb
	sum = zero
	if (idf==1) goto L30
	fmkm2 = zero
	if (abs(drb)<emin) fmkm2 = a*rb*exp(-half*drb^2) * 
		alnorm(a*drb, 0)*g1
	fmkm1 = b * da * fmkm2
	if (abs(d)<emin) {
		fmkm1 = fmkm1 + b * a * g2 * exp(-half * d^2)
	}
	if (ioe==0) sum = fmkm2
	else if (ioe==1) sum = fmkm1
	goto L20
	ifm2 = idf - 2
	ak = one
	fk = two
	for (k=2; k<=ifm2; k=k+2) {
		fkm1 = fk - one
		fmkm2 = b * (da * ak * fmkm1 + fmkm2) * fkm1 / fk
		ak = one / (ak * fkm1)
		fmkm1 = b * (da * ak * fmkm2 + fmkm1) * fk / (fk + one)
		if (ioe == 0) sum = sum + fmkm2
		else if (ioe == 1) sum = sum + fmkm1
		ak = one / (ak * fk)
		fk = fk + two
	}
	L20: if (ioe==0) goto L40
	L30: return(alnorm(drb, 1) + two * (sum + tfn(drb, a)))
	L40: return(alnorm(d, 1) + sum*g3)
/*
         NORMAL APPROXIMATION - K IS NOT TESTED AFTER THE TWO CALLS
         OF ALOGAM, BECAUSE A FAULT IS IMPOSSIBLE WHEN F EXCEEDS 100
*/
	L50:
	ifault = 1
	a = sqrt(half * f) * exp(alogam(half * (f - one), k) - 
		- alogam(half * f, k)) * d
	return( alnorm((st - a) / sqrt(f * (one + d^2) /
		(f - two) - a^2), 0) )
}
end