MODULE rsquared
IMPLICIT NONE

INTEGER, SAVE :: ifault

PUBLIC :: ifault, qr2

CONTAINS

!------------------------------------------------------------

! The following code for calculating quantiles of R^2 (qr2.for)
! is by the author of AS 261 and uses an improved algorithm.

FUNCTION qr2(m, size, rho2, p) RESULT(quantile)

!  Computes the quantile of the distribution of the square of the
!  sample multiple correlation coefficient for given number of
!  random variables M, sample size SIZE, square of the population
!  multiple correlation coefficient RHO2, and lower tail area P

!  Reference:
!  Ding, C.G. (1996) `On the computation of the distribution of
!  the square of the sample multiple correlation coefficient',
!  Comput. Statist. & Data Analysis, vol. 22, 345-350.

!  IFAULT is a fault indicator:
!  = 1 if there is no convergence after 10 Newton's iterations
!  = 2 if any of the input values is illegal
!  = 0 otherwise

!  No auxiliary algorithm is required

!  ELF90 compatible version by Alan Miller
!  N.B. The error indicator (ifault) has been put into the module contents
!       as functions can only return one result in ELF90.
!  Latest revision - 20 May 1997

INTEGER, INTENT(IN) :: m, size
REAL, INTENT(IN)    :: rho2, p
REAL                :: quantile

! Local variables
REAL               :: n, a, b, ab, ga, q0, yp, y, coeff, q, s, t, cdf, pdf,  &
                      gb, gab, v, bndcdf, bndpdf, ynew, diff
INTEGER            :: na, i, nab, nb, iter
REAL, PARAMETER    :: zero = 0.0, half = 0.5, one = 1.0, two = 2.0,   &
                      eps = 1.E-06, delta = 1.E-04, rp = 1.772453850905516028
INTEGER, PARAMETER :: itrmax = 10

quantile = p

!        Test for admissibility of arguments

ifault = 2
IF (m <= 1 .OR. size <= m .OR. rho2 < zero .OR.  &
    rho2 > one .OR. p < zero .OR. p > one) RETURN
ifault = 0
IF (p == zero .OR. p == one) RETURN

!        Calculate the constants needed for each Newton's iteration

a = (m - 1) / two
b = (size - m) / two
ab = (size - 1) / two
IF (MOD(m + 1, 2) == 0) THEN
  na = a + half
  ga = one
  DO i = 1, na
    ga = ga * i
  END DO
ELSE
  na = a + one
  ga = rp
  DO i = 1, na
    ga = ga * (i - half)
  END DO
END IF
IF (MOD(size - m, 2) == 0) THEN
  nb = b - half
  gb = one
  DO i = 1, nb
    gb = gb * i
  END DO
ELSE
  nb = b
  gb = rp
  DO i = 1, nb
    gb = gb * (i - half)
  END DO
END IF
IF (MOD(size - 1, 2) == 0) THEN
  nab = ab - half
  gab = one
  DO i = 1, nab
    gab = gab * i
  END DO
ELSE
  nab = ab
  gab = rp
  DO i = 1, nab
    gab = gab * (i - half)
  END DO
END IF
q0 = (one - rho2) ** ab
coeff = gab / ga / gb

!        Use 0.5 as a starting value for Newton's iterations

y = half

!        Perform Newton's iterations

DO iter = 1, itrmax
  
!        Evaluate the first terms of the series for CDF (distribution
!        function) and PDF (density)
  
  n = one
  yp = one - y
  t = coeff * y ** a * yp ** b
  s = a * t / y / yp
  q = q0
  v = q
  cdf = v * t
  pdf = q * s
  
!        Check if a + n > (a + b + n)y
  
  80 IF (a + n > (a + b + n) * y) GO TO 90
  
!        Evaluate the next terms of two series and then the
!        partial sums
  
  q = q * (a + b + n - one) * rho2 / n
  v = v + q
  s = t * (a + b + n - one) / yp
  t = t * y * (a + b + n - one) / (a + n)
  cdf = cdf + v * t
  pdf = pdf + q * s
  n = n + one
  GO TO 80
  
!        Find the error bounds and check for convergence for both
!        series
  
  90 bndcdf = t * y * (a + b + n - one) / (a + n - (a + b + n) * y)
  bndpdf = t * (a + b + n - one) * (one - v) / yp
  100 IF (bndcdf <= eps .AND. bndpdf <= eps) GO TO 110
  
!        Continue to update the terms and then accumulate
  
  q = q * (a + b + n - one) * rho2 / n
  v = v + q
  IF (bndcdf <= eps) THEN
    s = s * y * (a + b + n - one) / (a + n - one)
    pdf = pdf + q * s
    n = n + one
    bndpdf = s * y * (a + b + n - one) * (one - v) / (a + n - one)
    GO TO 100
  ELSE IF (bndpdf <= eps) THEN
    t = t * y * (a + b + n - one) / (a + n)
    cdf = cdf + v * t
    n = n + one
    bndcdf = t * y * (a+b+n-one) / (a+n - (a+b+n)*y)
    GO TO 100
  ELSE
    s = t * (a + b + n - one) / yp
    t = t * y * (a + b + n - one) / (a + n)
    cdf = cdf + v * t
    pdf = pdf + q * s
    n = n + one
    GO TO 90
  END IF
  
!        Obtain a new Y and make changes if it is illegal
  
  110 diff = (cdf - p) / pdf
  ynew = y - diff
  IF (ynew <= zero) THEN
    y = y / two
  ELSE IF (ynew >= one) THEN
    y = (y + one) / two
  ELSE
    y = ynew
  END IF
  
!        Check for convergence of Newton's iterations
  
  IF (ABS(diff) <= delta * y) THEN
    quantile = y
    RETURN
  END IF
END DO
ifault = 1
RETURN
END FUNCTION qr2

END MODULE rsquared



PROGRAM main

!        This is a driver program that calls QR2 and produces output

USE rsquared
IMPLICIT NONE

INTEGER :: m, size, icode, k
REAL    :: rho2, p, y

10 WRITE (*, 11)
11 FORMAT (/' ENTER M (>1), N (>M), RHO2 (BETWEEN 0 AND 1), and',  &
           /' P (BETWEEN 0 AND 1) FOR QR2 ==> ')
READ (*,*) m, size, rho2, p
y = qr2(m, size, rho2, p)
icode = ifault + 1
SELECT CASE (icode)
  CASE (1)
    WRITE (*, 21) m, size, rho2, p, y
    21 FORMAT (/' qr2(', i3, ',', i3, ',', g13.4, ',', g13.4, ') = ', g13.5)
  CASE (2)
    WRITE (*, 31)
    31 FORMAT (/' NO convergence after 10 Newton ITERATIONS')
  CASE (3)
    WRITE (*, 41)
    41 FORMAT (/' THE INPUT VALUE IS ILLEGAL!')
END SELECT
WRITE (*, 61)
61 FORMAT (/ ' ENTER 1 TO CONTINUE OR 0 TO QUIT ==> ')
READ (*,*) k
IF (k == 1) GO TO 10

STOP
END PROGRAM main