*! version 1.0.0  07jul2005
program prho
	version 9

	syntax varlist(min=2 max=2) [if] [in]
	marksample touse

	spearman `varlist' if `touse'
	mata: my_prho(`r(N)', `=abs(r(rho))')
	di as txt "      New test = " as res %12.4f 2*(1-r(p2))
end

mata:
version 9
mata set matastrict on

void my_prho(real scalar n, real scalar r)
{
	st_numscalar("r(p2)", scalar_prho(n, r))
}

real scalar scalar_prho(real scalar n, real scalar r)
{
	real scalar	ifault
	real scalar	res

	res = sl89_prho(n, (n^3-n)*(1-r)/6, ifault)
	if (ifault) res = . 
	return(res) 
}



real scalar sl89_prho(real scalar n, real scalar is, real scalar ifault)
{
/*
      double precision function prho(n, is, ifault)
c
c        Algorithm AS 89   Appl. Statist. (1975) Vol.24, No. 3, P377.
c       
c        To evaluate the probability of obtaining a value greater than or
c        equal to is, where is=(n**3-n)*(1-r)/6, r=Spearman's rho and n
c        must be greater than 1
c
c     Auxiliary function required: ALNORM = algorithm AS66
c
*/
	real vector l /*(16)*/
	real scalar b, x, y, u
	real scalar c1, c2, c3, c4, c5, c6, c7, c8, c9, c10, c11, c12
	real scalar prho
	real scalar js, nfac, i, ifr, m, ise, n1, mt, nn

	c1  = 0.2274d0 
	c2  = 0.2531d0 
	c3  = 0.1745d0 
	c4  = 0.0758d0 
	c5  = 0.1033d0 
	c6  = 0.3932d0
	c7  = 0.0879d0 
	c8  = 0.0151d0 
	c9  = 0.0072d0 
	c10 = 0.0831d0 
	c11 = 0.0131d0 
	c12 = 0.00046d0
/*
c
c        Test admissibility of arguments and initialize
c
*/
	l = J(1,16,.)

	prho = 1
	ifault = 1
	if (n < 1) return(prho)
	ifault = 0
	if (is < 0) return(prho)
	prho = 0
	if (is > n * (n * n -1) / 3) return(prho)
	js = is
	if (js != 2 * (js / 2)) js++
	if (n > 6) goto L6

/*
c
c        Exact evaluation of probability
c
*/
	nfac = 1
	for (i=1; i<=n; i++) {
		nfac = nfac * i
		l[i] = i
	}
	prho = 1 / nfac
	if (js == n * (n * n -1) / 3) return(prho)
	ifr = 0

	for (m=1; m<=nfac; m++) {
		ise = 0
		for (i=1; i<=n; i++) {
			ise = ise + (i - l[i]) ** 2
			L2:
		}
		if (js < ise) ifr = ifr++
		n1 = n
		L3:  mt = l[1]
		nn = n1 - 1
		for (i=1; i<=nn; i++) {
			l[i] = l[i + 1]
		}
		l[n1] = mt
		if (l[n1]!=n1 | n1== 2) goto L5
		n1--
		if (m != nfac) goto L3
		L5:
	}
	prho = ifr / nfac
	return(prho)
/*
c
c        Evaluation by Edgeworth series expansion
c
*/
	L6:
	b = 1 / n
	x = (6 * (js - 1) * b / (1 / (b * b) -1) -
						1) * sqrt(1 / b - 1)
	y = x * x
	u = x * b * (c1 + b * (c2 + c3 * b) + y * (-c4 +
     		b * (c5 + c6 * b) - y * b * (c7 + c8 * b -
     		y * (c9 - c10 * b + y * b * (c11 - c12 * y)))))
/*
c
c      Call to algorithm AS 66
	Original line was 
		prho = u / exp(y / 2) + alnorm(x, /* .true.*/ 1)

	alnorm() says area from x to infinity if .true.
	Thus
		alnorm(x, .true.)  <=>   normal(-x)

c
*/
	prho = u / exp(y / 2) + normal(-x)
	if (prho < 0) prho = 0
	else if (prho > 1) prho = 1
	return(prho)
}

scalar_prho(20, .2)
end