PROGRAM test395 ! Test translated versions of CACM algorithms 395 & 396 IMPLICIT NONE INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(12, 60) REAL (dp) :: ndf, prob, t INTERFACE FUNCTION student(t, ndf, normal) RESULT(fn_val) IMPLICIT NONE INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(12, 60) REAL (dp), INTENT(IN) :: t, ndf REAL (dp) :: fn_val INTERFACE FUNCTION normal(x) RESULT(fx) IMPLICIT NONE INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(12, 60) REAL (dp), INTENT(IN) :: x REAL (dp) :: fx END FUNCTION normal END INTERFACE END FUNCTION student FUNCTION t_quantile(p, ndf, normdev) RESULT(fn_val) IMPLICIT NONE INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(12, 60) REAL (dp), INTENT(IN) :: p, ndf REAL (dp) :: fn_val INTERFACE FUNCTION normdev(p) RESULT(x) IMPLICIT NONE INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(12, 60) REAL (dp), INTENT(IN) :: p REAL (dp) :: x END FUNCTION normdev END INTERFACE END FUNCTION t_quantile FUNCTION normal(x) RESULT(fx) IMPLICIT NONE INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(12, 60) REAL (dp), INTENT(IN) :: x REAL (dp) :: fx END FUNCTION normal FUNCTION normdev(p) RESULT(x) IMPLICIT NONE INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(12, 60) REAL (dp), INTENT(IN) :: p REAL (dp) :: x END FUNCTION normdev END INTERFACE 10 WRITE(*, *) 'Enter no. of deg. of freedom & t-value: ' READ(*, *) ndf, t prob = student(t, ndf, normal) WRITE(*, 900) prob 900 FORMAT(' Prob = ', g15.7) t = t_quantile(prob, ndf, normdev) WRITE(*, 910) t 910 FORMAT(' t-value = ', f10.6/) GO TO 10 END PROGRAM test395 FUNCTION normal(z) RESULT(prob) IMPLICIT NONE INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(12, 60) REAL (dp), INTENT(IN) :: z REAL (dp) :: prob CALL nprob(z, p=prob) RETURN CONTAINS SUBROUTINE nprob(z, p, q, pdf) ! For a given number (z) of standard deviations from the mean, the ! probabilities to the left (p) and right (q) of z are calculated, ! and the probability density (pdf). ! Programmer - Alan J. Miller ! Latest revision - 7 August 1997 ! This version is compatible with the ELF90 subset of Fortran 90. ! REFERENCE: ADAMS,A.G. AREAS UNDER THE NORMAL CURVE, ! ALGORITHM 39, COMPUTER J., VOL. 12, 197-8, 1969. IMPLICIT NONE INTEGER, PARAMETER :: doubleprec = SELECTED_REAL_KIND(15, 60) REAL(doubleprec), INTENT(IN) :: z REAL(doubleprec), INTENT(OUT), OPTIONAL :: p, q, pdf ! Local variables REAL(doubleprec) :: a0 = 0.5D0, a1 = 0.398942280444D0, a2 = 0.399903438504D0, & a3 = 5.75885480458D0, a4 = 29.8213557808D0, & a5 = 2.62433121679D0, a6 = 48.6959930692D0, & a7 = 5.92885724438D0, b0 = 0.398942280385D0, & b1 = 3.8052D-8, b2 = 1.00000615302D0, b3 = 3.98064794D-4, & b4 = 1.98615381364D0, b5 = 0.151679116635D0, & b6 = 5.29330324926D0, b7 = 4.8385912808D0, & b8 = 15.1508972451D0, b9 = 0.742380924027D0, & b10 = 30.789933034D0, b11 = 3.99019417011D0, & zero = 0.D0, one = 1.D0, zabs, y, pp, qq, ppdf zabs = ABS(z) IF (zabs < 12.7D0) THEN y = a0*z*z ppdf = EXP(-y)*b0 IF (PRESENT(pdf)) pdf = ppdf ! Z BETWEEN -1.28 AND +1.28 IF (zabs < 1.28) THEN qq = a0 - zabs*(a1-a2*y/(y+a3-a4/(y+a5+a6/(y+a7)))) IF(z < zero) THEN pp = qq qq = one - pp ELSE pp = one - qq END IF IF (PRESENT(p)) p = pp IF (PRESENT(q)) q = qq RETURN END IF ! ZABS BETWEEN 1.28 AND 12.7 qq = ppdf/(zabs-b1+b2/(zabs+b3+b4/(zabs-b5+b6/(zabs+b7-b8/ & (zabs+b9+b10/(zabs+b11)))))) IF(z < zero) THEN pp = qq qq = one - pp ELSE pp = one - qq END IF IF (PRESENT(p)) p = pp IF (PRESENT(q)) q = qq RETURN ! Z FAR OUT IN TAIL ELSE ppdf = zero IF(z < zero) THEN pp = zero qq = one ELSE pp = one qq = zero END IF IF (PRESENT(p)) p = pp IF (PRESENT(q)) q = qq IF (PRESENT(pdf)) pdf = ppdf RETURN END IF RETURN END SUBROUTINE nprob END FUNCTION normal FUNCTION normdev(p) RESULT(fn_val) IMPLICIT NONE INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(12, 60) REAL (dp), INTENT(IN) :: p REAL (dp) :: fn_val ! Local variable INTEGER :: ifault CALL ppnd16(p, fn_val, ifault) IF (ifault /= 0) WRITE(*, *) 'Error in ppnd16: ifault =', ifault RETURN CONTAINS SUBROUTINE ppnd16 (p, normal_dev, ifault) ! ALGORITHM AS241 APPL. STATIST. (1988) VOL. 37, NO. 3 ! Produces the normal deviate Z corresponding to a given lower ! tail area of P; Z is accurate to about 1 part in 10**16. ! The hash sums below are the sums of the mantissas of the ! coefficients. They are included for use in checking ! transcription. ! This ELF90-compatible version by Alan Miller - 20 August 1996 ! N.B. The original algorithm is as a function; this is a subroutine IMPLICIT NONE INTEGER, PARAMETER :: doub_prec = SELECTED_REAL_KIND(15, 60) REAL (doub_prec), INTENT(IN) :: p INTEGER, INTENT(OUT) :: ifault REAL (doub_prec), INTENT(OUT) :: normal_dev ! Local variables REAL (doub_prec) :: zero = 0.d0, one = 1.d0, half = 0.5d0, & split1 = 0.425d0, split2 = 5.d0, const1 = 0.180625d0, & const2 = 1.6d0, q, r ! Coefficients for P close to 0.5 REAL (doub_prec) :: a0 = 3.3871328727963666080D0, & a1 = 1.3314166789178437745D+2, & a2 = 1.9715909503065514427D+3, & a3 = 1.3731693765509461125D+4, & a4 = 4.5921953931549871457D+4, & a5 = 6.7265770927008700853D+4, & a6 = 3.3430575583588128105D+4, & a7 = 2.5090809287301226727D+3, & b1 = 4.2313330701600911252D+1, & b2 = 6.8718700749205790830D+2, & b3 = 5.3941960214247511077D+3, & b4 = 2.1213794301586595867D+4, & b5 = 3.9307895800092710610D+4, & b6 = 2.8729085735721942674D+4, & b7 = 5.2264952788528545610D+3 ! HASH SUM AB 55.8831928806149014439 ! Coefficients for P not close to 0, 0.5 or 1. REAL (doub_prec) :: c0 = 1.42343711074968357734D0, & c1 = 4.63033784615654529590D0, & c2 = 5.76949722146069140550D0, & c3 = 3.64784832476320460504D0, & c4 = 1.27045825245236838258D0, & c5 = 2.41780725177450611770D-1, & c6 = 2.27238449892691845833D-2, & c7 = 7.74545014278341407640D-4, & d1 = 2.05319162663775882187D0, & d2 = 1.67638483018380384940D0, & d3 = 6.89767334985100004550D-1, & d4 = 1.48103976427480074590D-1, & d5 = 1.51986665636164571966D-2, & d6 = 5.47593808499534494600D-4, & d7 = 1.05075007164441684324D-9 ! HASH SUM CD 49.33206503301610289036 ! Coefficients for P near 0 or 1. REAL (doub_prec) :: e0 = 6.65790464350110377720D0, & e1 = 5.46378491116411436990D0, & e2 = 1.78482653991729133580D0, & e3 = 2.96560571828504891230D-1, & e4 = 2.65321895265761230930D-2, & e5 = 1.24266094738807843860D-3, & e6 = 2.71155556874348757815D-5, & e7 = 2.01033439929228813265D-7, & f1 = 5.99832206555887937690D-1, & f2 = 1.36929880922735805310D-1, & f3 = 1.48753612908506148525D-2, & f4 = 7.86869131145613259100D-4, & f5 = 1.84631831751005468180D-5, & f6 = 1.42151175831644588870D-7, & f7 = 2.04426310338993978564D-15 ! HASH SUM EF 47.52583317549289671629 ifault = 0 q = p - half IF (ABS(q) <= split1) THEN r = const1 - q * q normal_dev = q * (((((((a7*r + a6)*r + a5)*r + a4)*r + a3)*r + a2)*r + a1)*r + a0) / & (((((((b7*r + b6)*r + b5)*r + b4)*r + b3)*r + b2)*r + b1)*r + one) RETURN ELSE IF (q < zero) THEN r = p ELSE r = one - p END IF IF (r <= zero) THEN ifault = 1 normal_dev = zero RETURN END IF r = SQRT(-LOG(r)) IF (r <= split2) THEN r = r - const2 normal_dev = (((((((c7*r + c6)*r + c5)*r + c4)*r + c3)*r + c2)*r + c1)*r + c0) / & (((((((d7*r + d6)*r + d5)*r + d4)*r + d3)*r + d2)*r + d1)*r + one) ELSE r = r - split2 normal_dev = (((((((e7*r + e6)*r + e5)*r + e4)*r + e3)*r + e2)*r + e1)*r + e0) / & (((((((f7*r + f6)*r + f5)*r + f4)*r + f3)*r + f2)*r + f1)*r + one) END IF IF (q < zero) normal_dev = - normal_dev RETURN END IF RETURN END SUBROUTINE ppnd16 END FUNCTION normdev