FUNCTION rngpi(t, n) RESULT(fn_val) ! N.B. Argument IFAULT has been removed. ! Code converted using TO_F90 by Alan Miller ! Date: 2003-02-09 Time: 14:24:28 ! Algorithm AS 126 Appl. Statist. (1978) Vol. 27, No. 2 ! Computes the probability of the normal range given t, the ! upper limit of integration, and n, the sample size. ! Auxiliary function required: ALNORM = algorithm AS66 IMPLICIT NONE INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(15, 100) REAL (dp), INTENT(IN) :: t INTEGER, INTENT(IN) :: n REAL (dp) :: fn_val ! Local variables REAL (dp) :: a, b, c, xl, y INTEGER :: i REAL (dp), PARAMETER :: g(8) = (/ & 0.4947004675D0, 0.4722875115D0, 0.4328156012D0, 0.3777022042D0, & 0.3089381222D0, 0.2290083888D0, 0.1408017754D0, 0.04750625492D0 /) REAL (dp), PARAMETER :: h(8) = (/ & 0.01357622971D0, 0.03112676197D0, 0.04757925584D0, 0.06231448563D0, & 0.07479799441D0, 0.08457825969D0, 0.09130170752D0, 0.09472530523D0 /) REAL (dp), PARAMETER :: zero = 0.0_dp, half = 0.5_dp, two = 2.0_dp, & eight = 8.0_dp INTERFACE FUNCTION alnorm( x, upper ) RESULT( fn_val ) IMPLICIT NONE INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(15, 100) REAL(DP), INTENT(IN) :: x LOGICAL, INTENT(IN) :: upper REAL(DP) :: fn_val END FUNCTION alnorm END INTERFACE fn_val = zero IF (t <= zero .OR. n <= 1) RETURN xl = half * t a = half * (eight+xl) b = eight - xl y = zero DO i = 1, 8 c = b * g(i) y = y + h(i) * (risf(a+c) + risf(a-c)) END DO fn_val = (two*(alnorm(xl,.false.)-half)) ** n + two * b * y * n RETURN CONTAINS FUNCTION risf(x) RESULT(fn_val) REAL (dp), INTENT(IN) :: x REAL (dp) :: fn_val fn_val = 0.3989422804D0 * EXP(-half*x*x) * (alnorm(x, .false.) - & alnorm(x-t, .false.)) ** (n-1) RETURN END FUNCTION risf END FUNCTION rngpi