PROGRAM test_as241 IMPLICIT NONE INTEGER, PARAMETER :: sp = SELECTED_REAL_KIND(6, 30), & dp = SELECTED_REAL_KIND(12, 60) REAL (sp) :: p_single, z_single REAL (dp) :: p_dble, z_dble INTEGER :: ifault INTERFACE SUBROUTINE ppnd7 (p, normal_dev, ifault) IMPLICIT NONE INTEGER, PARAMETER :: sp = SELECTED_REAL_KIND(6, 30) REAL (sp), INTENT(IN) :: p INTEGER, INTENT(OUT) :: ifault REAL (sp), INTENT(OUT) :: normal_dev END SUBROUTINE ppnd7 SUBROUTINE ppnd16 (p, normal_dev, ifault) IMPLICIT NONE INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(12, 60) REAL (dp), INTENT(IN) :: p INTEGER, INTENT(OUT) :: ifault REAL (dp), INTENT(OUT) :: normal_dev END SUBROUTINE ppnd16 END INTERFACE DO WRITE(*, *)'Enter area under normal curve: ' READ(*, *) p_dble p_single = p_dble CALL ppnd7(p_single, z_single, ifault) IF (ifault /= 0) THEN WRITE(*, '(1x, "IFAULT =", i3)') ifault CYCLE END IF CALL ppnd16(p_dble, z_dble, ifault) IF (ifault /= 0) THEN WRITE(*, '(1x, "IFAULT =", i3)') ifault CYCLE END IF WRITE(*, '(1x, a, f11.6, 2x, a, f16.11)') & 'Low prec. result =', z_single, ' High prec. result =', z_dble END DO STOP END PROGRAM test_as241 SUBROUTINE ppnd7 (p, normal_dev, ifault) ! ALGORITHM AS241 APPL. STATIST. (1988) VOL. 37, NO. 3, 477- 484. ! Produces the normal deviate Z corresponding to a given lower tail area of P; ! Z is accurate to about 1 part in 10**7. ! 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 :: sp = SELECTED_REAL_KIND(6, 30) REAL (sp), INTENT(IN) :: p INTEGER, INTENT(OUT) :: ifault REAL (sp), INTENT(OUT) :: normal_dev ! Local variables REAL (sp) :: zero = 0.0, one = 1.0, half = 0.5, split1 = 0.425, & split2 = 5.0, const1 = 0.180625, const2 = 1.6, q, r ! Coefficients for P close to 0.5 REAL (sp) :: a0 = 3.3871327179E+00, a1 = 5.0434271938E+01, & a2 = 1.5929113202E+02, a3 = 5.9109374720E+01, & b1 = 1.7895169469E+01, b2 = 7.8757757664E+01, & b3 = 6.7187563600E+01 ! HASH SUM AB 32.3184577772 ! Coefficients for P not close to 0, 0.5 or 1. REAL (sp) :: c0 = 1.4234372777E+00, c1 = 2.7568153900E+00, & c2 = 1.3067284816E+00, c3 = 1.7023821103E-01, & d1 = 7.3700164250E-01, d2 = 1.2021132975E-01 ! HASH SUM CD 15.7614929821 ! Coefficients for P near 0 or 1. REAL (sp) :: e0 = 6.6579051150E+00, e1 = 3.0812263860E+00, & e2 = 4.2868294337E-01, e3 = 1.7337203997E-02, & f1 = 2.4197894225E-01, f2 = 1.2258202635E-02 ! HASH SUM EF 19.4052910204 ifault = 0 q = p - half IF (ABS(q) <= split1) THEN r = const1 - q * q normal_dev = q * (((a3 * r + a2) * r + a1) * r + a0) / & (((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 = (((c3 * r + c2) * r + c1) * r + c0) / ((d2 * r + d1) * r + one) ELSE r = r - split2 normal_dev = (((e3 * r + e2) * r + e1) * r + e0) / ((f2 * r + f1) * r + one) END IF IF (q < zero) normal_dev = - normal_dev RETURN END IF END SUBROUTINE ppnd7 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 :: dp = SELECTED_REAL_KIND(12, 60) REAL (dp), INTENT(IN) :: p INTEGER, INTENT(OUT) :: ifault REAL (dp), INTENT(OUT) :: normal_dev ! Local variables REAL (dp) :: 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 (dp) :: 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 (dp) :: 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 (dp) :: 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 END SUBROUTINE ppnd16