!	Program StudentRecursion
!	QUAD precision FORTRAN implementation
!	Uses non-linear iteration of pair of power series
!	to compute Student quantile function.
!	with ABSOFT MacOS on Intel FORTRAN 95
!	Version 0.2 (early prototype!)
!	William Shaw, April 2007
!	william.shaw@kcl.ac.uk

SUBROUTINE ISTUDENTSCALE(N, SCALE)
! This subroutine computes the normalization coefficient for the series
! for integer cases. Tabulated up to dof = 10 then recursion used. Avoids calls to 
! gamma functions when not needed. 
IMPLICIT NONE
INTEGER, PARAMETER :: QP=SELECTED_REAL_KIND(31, 300)
REAL(QP) SCALE
INTEGER*4 N, K
REAL(QP) MYPI
MYPI = 3.141592653589793238462643383279503_QP
SELECT CASE(N)
	CASE(1) 
	SCALE = MYPI/2.0_QP
	CASE(2) 
	SCALE = SQRT(2.0_QP)
	CASE(3) 
	SCALE = MYPI*SQRT(3.0_QP)/4.0_QP
	CASE(4) 
	SCALE = 4.0_QP/3.0_QP
	CASE(5) 
	SCALE = 3.0_QP*SQRT(5.0_QP)*MYPI/16.0_QP
	CASE(6) 
	SCALE = 8.0_QP*SQRT(2.0_QP/3.0_QP)/5.0_QP
	CASE(7) 
	SCALE = 5.0_QP*SQRT(7.0_QP)*MYPI/32.0_QP
	CASE(8) 
	SCALE =32.0_QP*SQRT(2.0_QP)/35.0_QP
	CASE(9) 
	SCALE = 105.0_QP*MYPI/256.0_QP
	CASE(10) 
	SCALE = 128.0_QP*SQRT(2.0_QP/5.0_QP)/63.0_QP
	CASE DEFAULT 
	IF (MOD(N,2)==0) THEN
		SCALE = 128.0_QP*SQRT(2.0_QP/5.0_QP)/63.0_QP
		DO K=12, N, 2
			SCALE = SCALE*SQRT(K*(K-2.0_QP))/(K-1.0_QP)
		END DO
	ELSE
		SCALE = 105.0_QP*MYPI/256.0_QP
		DO K=11, N, 2
			SCALE = SCALE*SQRT(K*(K-2.0_QP))/(K-1.0_QP)
		END DO
	END IF
END SELECT
END SUBROUTINE ISTUDENTSCALE

SUBROUTINE RSTUDENTSCALE(DOF, SCALE)
! This variation for real dof uses upward recursion and then
! an asymptotic series
IMPLICIT NONE
INTEGER, PARAMETER :: QP=SELECTED_REAL_KIND(31, 300)
REAL(QP) SCALE, DOF, ASYMPTOTE, INVDOF, CUPOWER
INTEGER N
REAL(QP) ROOTPIOTWO
REAL(QP) AARR(12)
ROOTPIOTWO = 1.2533141373155002512078826424055_QP
SCALE = ROOTPIOTWO
AARR(1)=1.0_QP/4.0_QP
AARR(2)=1.0_QP/3.2Q1_QP
AARR(3)=-5.0_QP/1.28Q2_QP 
AARR(4)=-2.1Q1_QP/2.048Q3_QP 
AARR(5)=3.99Q2_QP/8.192Q3_QP
AARR(6)=8.69Q2_QP/6.5536Q4_QP
AARR(7)=-3.9325Q4_QP/2.62144Q5_QP
AARR(8)=-3.34477Q5_QP/8.388608Q6_QP
AARR(9)=2.8717403Q7_QP/3.3554432Q7_QP
AARR(10)=5.9697183Q7_QP/2.68435456Q8_QP 
!AARR(11)=-8.400372435Q9_QP/1.073741824Q9_QP 
!AARR(12)=-3.4429291905Q10_QP/1.7179869184Q10_QP 
DO WHILE(DOF<1000)
	SCALE = (DOF+1.0_QP)/SQRT(DOF*(DOF+2.0_QP))*SCALE
	DOF = DOF+2.0_QP
END DO 
INVDOF = 1.0_QP/DOF
CUPOWER = INVDOF
ASYMPTOTE = 1.0_QP
DO N = 1, 10
	ASYMPTOTE = ASYMPTOTE+CUPOWER*AARR(N)
	CUPOWER=INVDOF*CUPOWER
END DO
SCALE=SCALE*ASYMPTOTE
END SUBROUTINE RSTUDENTSCALE


SUBROUTINE CENTRALRECURSIVEQUANTILE(QUANTILE,U,RAT,SCALE)
IMPLICIT NONE
INTEGER, PARAMETER :: QP=SELECTED_REAL_KIND(31, 300)
REAL(QP) QUANTILE
INTEGER ARRSIZE, I
PARAMETER (ARRSIZE=200)
REAL(QP) U, SCALE
REAL(QP) RAT(ARRSIZE-1)
REAL(QP) B,S,T,Q
REAL(QP) ROOTPIOTWO
B=SCALE*(2.0_QP*U-1.0_QP)
I=-1
S=B
Q=B*B
DO WHILE((ABS(B)>1.0Q-36.AND.I<= ARRSIZE-3).OR.I<10)
	T=S
	I=I+1
	B=B*Q*RAT(I+1)
	S=T+B
END DO
IF (I.EQ.ARRSIZE-2) THEN 
WRITE(6,*) "End of series reached - switch to tail model"
END IF
QUANTILE = S
END SUBROUTINE CENTRALRECURSIVEQUANTILE

SUBROUTINE TAILRECURSIVEQUANTILE(QUANTILE,X,RAT)
IMPLICIT NONE
INTEGER, PARAMETER :: QP=SELECTED_REAL_KIND(31, 300)
REAL(QP) QUANTILE,X
INTEGER TARRSIZE, I
PARAMETER (TARRSIZE=400)
REAL(QP) RAT(TARRSIZE-1)
REAL(QP) B,S,T
B=X
I=-1
S=B
DO WHILE((ABS(B)>1.0Q-36.AND.I<= TARRSIZE-3).OR.I<10)
	T=S
	I=I+1
	B=B*X*RAT(I+1)
	S=T+B
END DO
IF (I.EQ.TARRSIZE-2) THEN 
WRITE(6,*) "End of series reached - switch to tail model"
END IF
QUANTILE = S
END SUBROUTINE TAILRECURSIVEQUANTILE



PROGRAM StudentRecursion
IMPLICIT NONE
INTEGER, PARAMETER :: QP=SELECTED_REAL_KIND(31, 300)
INTEGER ARRSIZE, TARRSIZE
PARAMETER (ARRSIZE=200)
PARAMETER (TARRSIZE=400)
INTEGER K, L, M, N, P, IDOFQ
REAL(QP) DOF
REAL(QP) SCALE, TEMPA, TEMPB, X
REAL(QP) U, QUANTILE
REAL(QP) PIOVERFOUR
REAL(QP) QARR(ARRSIZE)
REAL(QP) TQARR(TARRSIZE)
REAL(QP) RATIOS(ARRSIZE-1)
REAL(QP) TRATIOS(TARRSIZE-1)
PIOVERFOUR = 0.78539816339744830961566084581988_QP
U = 0.5_QP
QUANTILE = 0.0_QP
SCALE = 0.0_QP
IDOFQ = 0
WRITE(6,*) "Student quantile function"
WRITE(6,*) "Enter degrees of freedom as a.b, where a>= 1 and b may be empty"
WRITE(6,*) "e.g. 3., 4.1, 20.0"
WRITE(6,*) "Program accepts REAL degrees of freedom >= 1.0"
10 READ(*,*) DOF
IF (DOF<1) THEN 
WRITE(6,*) "Must be >= 1., try again!"
GOTO 10
END IF

WRITE(6,*) DOF
N = NINT(DOF)
IF (ABS(N-DOF)<1.0E-6_QP) THEN 
IDOFQ=1
WRITE(6,*) "Treated as exact integer degree of freedom"
ELSE
WRITE(6,*) "Treated as non-integer degree of freedom"
END IF
IF (IDOFQ==1) THEN
CALL ISTUDENTSCALE(N, SCALE)
ELSE
CALL RSTUDENTSCALE(DOF, SCALE)
END IF

WRITE(6,*) SCALE
WRITE(6,*) "Initializing series....."
WRITE(6,*) "Coefficients are:"
! Do the cubic recursion for the central Student
QARR(1) = 1.0_QP
DO P=1, ARRSIZE-1
	QARR(P+1)=0.0_QP
	DO K=0, P-1
		DO L=0, P-K-1
		TEMPA = QARR(K+1)*QARR(L+1)*QARR(P-K-L)
		TEMPB = (2.0_QP*L+1.0_QP)*(2.0_QP*(P-K-L)-1.0_QP)*(1.0_QP + 1.0_QP/DOF)
		TEMPB = TEMPB-2.0_QP*K*(2.0_QP*K+1)/DOF
		QARR(P+1)=QARR(P+1)+TEMPA*TEMPB
		END DO
	END DO
	QARR(P+1) = 0.5_QP*QARR(P+1)/P/(2.0_QP*P+1.0_QP)
END DO 

DO P=1, ARRSIZE-1
RATIOS(P) = QARR(P+1)/QARR(P)
END DO

WRITE(6,*) "Tail coefficients are"
! Do the mixed quadratic-cubic recursion for the tail regions
TQARR(1) = 1.0_QP
DO P=2, TARRSIZE
	TQARR(P)=0.0_QP
	DO K=1, P-1
		DO L=1, P-K
		TEMPA = TQARR(K)*TQARR(L)*TQARR(P+1-K-L)
		TEMPB = K*(K-DOF/2.0_QP)+(DOF/2.0_QP-1.5_QP)*L*(P+1.0_QP-K-L)
		TQARR(P)=TQARR(P)+TEMPA*TEMPB
		END DO
	END DO
	IF (P>2) THEN
		DO K=2,P-1
		TEMPA=TQARR(K)*TQARR(P+1-K)
		TEMPB=(1.0_QP-DOF/2.0_QP)*K*(P-K)-K*(K-1.0_QP)
		TQARR(P)=TQARR(P)+TEMPA*TEMPB
		END DO 
	END IF
	TQARR(P) = TQARR(P)/(P*P+P*(DOF/2.0_QP-2.0_QP)+1.0_QP-DOF/2.0_QP)
END DO 

DO P=1, TARRSIZE-1
WRITE(6,*) TQARR(P)
IF (ABS(TQARR(P))>0.0_QP) THEN
TRATIOS(P) = TQARR(P+1)/TQARR(P)
ELSE
TRATIOS(P)=0.0_QP
END IF
END DO




WRITE(6,*) "Series ratios computed in quad precision to order 200"
WRITE(6,*) "Writing quantiles to fstuquantiles.txt"
OPEN(UNIT=5, FILE="fstuquantiles.txt", STATUS="UNKNOWN")
DO K=1, 99
U = K/100.0_QP
IF ((U>0.1_QP).AND.(U<0.9_QP)) CALL CENTRALRECURSIVEQUANTILE(QUANTILE,U, RATIOS,SCALE)
IF (U>=0.9_QP) THEN
X = (1.0_QP-U)*SCALE*2.0_QP*SQRT(DOF)
X = X**(2.0_QP/DOF)
CALL TAILRECURSIVEQUANTILE(QUANTILE,X,TRATIOS)
QUANTILE=SQRT(DOF*(1.0_QP/QUANTILE-1.0_QP))
END IF
IF (U<=0.1_QP) THEN
X = U*SCALE*2.0_QP*SQRT(DOF)
X = X**(2.0_QP/DOF)
CALL TAILRECURSIVEQUANTILE(QUANTILE,X,TRATIOS)
QUANTILE=-SQRT(DOF*(1.0_QP/QUANTILE-1.0_QP))
END IF


WRITE(5,*) U
WRITE(5,FMT= "(F36.32)") QUANTILE
END DO
CLOSE(5, STATUS="KEEP")
READ(*,*) DOF

END PROGRAM StudentRecursion