C* FORTRAN CODE WRITTEN FOR INCLUSION IN IBM RESEARCH REPORT RC20525, *
C* 'FORTRAN ROUTINES FOR USE WITH THE METHOD OF L-MOMENTS, VERSION 3' *
C* *
C* J. R. M. HOSKING *
C* IBM RESEARCH DIVISION *
C* T. J. WATSON RESEARCH CENTER *
C* YORKTOWN HEIGHTS *
C* NEW YORK 10598, U.S.A. *
C* *
C* VERSION 3 AUGUST 1996 *
C* *
C* VERSION 3.03 JUNE 2000 *
C* * Fixed: WRITE(6,7000) and FORMAT statement 7000 incompatible *
C* *
C***********************************************************************
C
C QUANTILE FUNCTION OF THE STANDARD NORMAL DISTRIBUTION
C
C BASED ON ALGORITHM AS241, APPL. STATIST. (1988) VOL.37 NO.3
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DATA ZERO/0D0/,HALF/0.5D0/,ONE/1D0/
DATA SPLIT1/0.425D0/,SPLIT2/5D0/,CONST1/0.180625D0/,CONST2/1.6D0/
C
C COEFFICIENTS OF RATIONAL-FUNCTION APPROXIMATIONS
C
DATA A0,A1,A2,A3,A4,A5,A6,A7,B1,B2,B3,B4,B5,B6,B7/
* 0.33871 32872 79636 661D 1,
* 0.13314 16678 91784 377D 3, 0.19715 90950 30655 144D 4,
* 0.13731 69376 55094 611D 5, 0.45921 95393 15498 715D 5,
* 0.67265 77092 70087 009D 5, 0.33430 57558 35881 281D 5,
* 0.25090 80928 73012 267D 4, 0.42313 33070 16009 113D 2,
* 0.68718 70074 92057 908D 3, 0.53941 96021 42475 111D 4,
* 0.21213 79430 15865 959D 5, 0.39307 89580 00927 106D 5,
* 0.28729 08573 57219 427D 5, 0.52264 95278 85285 456D 4/
DATA C0,C1,C2,C3,C4,C5,C6,C7,D1,D2,D3,D4,D5,D6,D7/
* 0.14234 37110 74968 358D 1,
* 0.46303 37846 15654 530D 1, 0.57694 97221 46069 141D 1,
* 0.36478 48324 76320 461D 1, 0.12704 58252 45236 838D 1,
* 0.24178 07251 77450 612D 0, 0.22723 84498 92691 846D -1,
* 0.77454 50142 78341 408D -3, 0.20531 91626 63775 882D 1,
* 0.16763 84830 18380 385D 1, 0.68976 73349 85100 005D 0,
* 0.14810 39764 27480 075D 0, 0.15198 66656 36164 572D -1,
* 0.54759 38084 99534 495D -3, 0.10507 50071 64441 684D -8/
DATA E0,E1,E2,E3,E4,E5,E6,E7,F1,F2,F3,F4,F5,F6,F7/
* 0.66579 04643 50110 378D 1,
* 0.54637 84911 16411 437D 1, 0.17848 26539 91729 133D 1,
* 0.29656 05718 28504 891D 0, 0.26532 18952 65761 230D -1,
* 0.12426 60947 38807 844D -2, 0.27115 55568 74348 758D -4,
* 0.20103 34399 29228 813D -6, 0.59983 22065 55887 938D 0,
* 0.13692 98809 22735 805D 0, 0.14875 36129 08506 149D -1,
* 0.78686 91311 45613 259D -3, 0.18463 18317 51005 468D -4,
* 0.14215 11758 31644 589D -6, 0.20442 63103 38993 979D-14/
C
Q=F-HALF
IF(DABS(Q).GT.SPLIT1)GOTO 10
R=CONST1-Q*Q
QUASTN=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
10 R=F
IF(Q.GE.ZERO)R=ONE-F
IF(R.LE.ZERO)GOTO 1000
R=DSQRT(-DLOG(R))
IF(R.GT.SPLIT2)GOTO 20
R=R-CONST2
QUASTN=(((((((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)
GOTO 30
20 R=R-SPLIT2
QUASTN=(((((((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)
30 IF(Q.LT.ZERO)QUASTN=-QUASTN
177