!	Program NormalQuantile
!	QUAD precision FORTRAN implementation
!	Uses non-linear iteration of central power series
!	to compute Normal quantile function.
!	Relative error < 1.0 E-29 on 0.0007 < u < 0.9993
!	Relative error < 1.0 E-24 on 0.0005 < u < 0.9995
!	with ABSOFT MacOS on Intel FORTRAN 95
!	Version 0.6 (simple termination, no convergence accln and no tail management)
!	William Shaw, April 2007
!	william.shaw@kcl.ac.uk

SUBROUTINE RECURSIVEQUANTILE(QUANTILE,U,RAT)
IMPLICIT NONE
INTEGER, PARAMETER :: QP=SELECTED_REAL_KIND(31, 300)
REAL(QP) QUANTILE
INTEGER ARRSIZE, I
PARAMETER (ARRSIZE=25000)
REAL(QP) U
REAL(QP) RAT(ARRSIZE-1)
REAL(QP) B,S,T,Q
REAL(QP) ROOTPIOTWO
ROOTPIOTWO = 1.2533141373155002512078826424055_QP
B=2*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*ROOTPIOTWO
END SUBROUTINE RECURSIVEQUANTILE

PROGRAM NormalQuantile
IMPLICIT NONE
INTEGER, PARAMETER :: QP=SELECTED_REAL_KIND(31, 300)
INTEGER ARRSIZE
PARAMETER (ARRSIZE=25000)
INTEGER K, M
REAL(QP) U, QUANTILE
REAL(QP) PIOVERFOUR
REAL(QP) QARR(ARRSIZE)
REAL(QP) RATIOS(ARRSIZE-1)
PIOVERFOUR = 0.78539816339744830961566084581988_QP
U = 0.5_QP
QUANTILE = 0.0_QP

WRITE(6,*) "Steinbrecher's recursive normal quantile"
WRITE(6,*) "Initializing series....."
QARR(1) = 1.0_QP
QARR(2) = PIOVERFOUR
DO K=2, ARRSIZE-1
	QARR(K+1)=0.0_QP
	DO M=0, K/2-1, 1
	QARR(K+1)=QARR(K+1)+QARR(M+1)*QARR(K-M)*(1.0_QP/(M+1)/(2*M+1)+1.0_QP/(K-M)/(2*K-2*M-1)) 
	END DO
	IF (MOD(K+1,2)==0) THEN
	QARR(K+1) = QARR(K+1)+2*QARR((K+1)/2)*QARR((K+1)/2)/K/(K+1)
	END IF 
	QARR(K+1) = QARR(K+1)*PIOVERFOUR
END DO 

DO K=1, ARRSIZE-1
QARR(K+1) = QARR(K+1)/(2*K+1)
RATIOS(K) = QARR(K+1)/QARR(K)
END DO

WRITE(6,*) "Series ratios computed in quad precision to order 25000"
OPEN(UNIT=5, FILE="fquantiles.txt", STATUS="UNKNOWN")
DO K=5, 9995
U = K/10000.0_QP
CALL RECURSIVEQUANTILE(QUANTILE,U, RATIOS)
WRITE(5,FMT= "(F36.32)") QUANTILE
END DO
CLOSE(5, STATUS="KEEP")

END PROGRAM NormalRecursion