C     ALGORITHM 487 COLLECTED ALGORITHMS FROM ACM.
C     ALGORITHM APPEARED IN COMM. ACM, VOL. 17, NO. 12,
C     P. 703.
      REAL FUNCTION PKS2(N, D)                                          PKS   10
      INTEGER N
C N IS THE SAMPLE SIZE USED.
      REAL D
C D IS THE MAXIMUM MAGNITUDE (OF THE DISCREPANCY
C BETWEEN THE EMPIRICAL AND PROPOSED DISTRIBUTIONS)
C IN EITHER THE POSITIVE OR NEGATIVE DIRECTION.
C PKS2 IS THE EXACT PROBABILITY OF OBTAINING A
C DEVIATION NO LARGER THAN D.
C THESE FORMULAS APPEAR AS (23) AND (24) IN
C J. DURBIN.  THE PROBABILITY THAT THE SAMPLE
C DISTRIBUTION FUNCTION LIES BETWEEN TWO PARALLEL
C STRAIGHT LINES. ANNALS OF MATHEMATICAL STATISTICS
C 39, 2(APRIL 1968),398-411.
      DOUBLE PRECISION Q(141), FACT(141), SUM, CI,
     * FT, FU, FV
      IF (N.EQ.1) GO TO 90
      FN = FLOAT(N)
      FND = FN*D
      NDT = IFIX(2.*FND)
      IF (NDT.LT.1) GO TO 100
      ND = IFIX(FND)
      NDD = MIN0(2*ND,N)
      NDP = ND + 1
      NDDP = NDD + 1
      FACT(1) = 1.
      CI = 1.
      DO 10 I=1,N
        FACT(I+1) = FACT(I)*CI
        CI = CI + 1.
   10 CONTINUE
      Q(1) = 1.
      IF (NDD.EQ.0) GO TO 50
      CI = 1.
      DO 20 I=1,NDD
        Q(I+1) = CI**I/FACT(I+1)
        CI = CI + 1.
   20 CONTINUE
      IF (NDP.GT.N) GO TO 80
      FV = FLOAT(NDP) - FND
      JMAX = IDINT(FV) + 1
      DO 40 I=NDP,NDD
        SUM = 0.
        FT = FND
        K = I
        FU = FV
        DO 30 J=1,JMAX
          SUM = SUM + FT**(J-2)/FACT(J)*FU**K/
     *     FACT(K+1)
          FT = FT + 1.
          FU = FU - 1.
          K = K - 1
   30   CONTINUE
        Q(I+1) = Q(I+1) - 2.*FND*SUM
        JMAX = JMAX + 1
        FV = FV + 1.
   40 CONTINUE
      IF (NDD.EQ.N) GO TO 80
   50 DO 70 I=NDDP,N
        SUM = 0.
        SIGN = 1.
        FT = 2.*FND
        DO 60 J=1,NDT
          FT = FT - 1.
          K = I - J + 1
          SUM = SUM + SIGN*FT**J/FACT(J+1)*Q(K)
          SIGN = -SIGN
   60   CONTINUE
        Q(I+1) = SUM
   70 CONTINUE
   80 PKS2 = Q(N+1)*FACT(N+1)/FN**N
      RETURN
   90 PKS2 = 2.*D - 1.
      RETURN
  100 PKS2 = 0.
      RETURN
      END
      SUBROUTINE PRFAC                                                  PRF   10
      DOUBLE PRECISION PF(4,40)
      DIMENSION DXA(4)
      COMMON DX, DXA, PF, J
      DATA I /1/
      DO 10 J=1,4
        IF (DXA(J).EQ.DX) RETURN
   10 CONTINUE
      J = I
      I = I + 1
      IF (I.EQ.5) I = 1
      DXA(J) = DX
      PF(J,1) = 1.
      DO 20 K=2,38
        PF(J,K) = (PF(J,K-1)*DX)/FLOAT(K-1)
   20 CONTINUE
      RETURN
      END
      FUNCTION CEIL(X)                                                  CEI   10
      IF (X.GE.0.) GO TO 10
      I = -X
      CEIL = -I
      RETURN
   10 I = X + .99999999
      CEIL = I
      RETURN
      END
      FUNCTION PKS(N, EPS)                                              PKS   10
C CALCULATE THE CUMULATIVE DISTRIBUTION OF THE
C KOLMOGOROV-SMIRNOV STATISTIC USING THE FORMULAS OF
C JOHN POMERANZ. EXACT VALUES OF THE TWO-SIDED
C KOLMOGOROV- SMIRNOV CUMULATIVE DISTRIBUTION FOR
C FINITE SAMPLE SIZE. TECHNICAL REPORT NUMBER 88,
C COMPUTER SCIENCES DEPARTMENT, PURDUE UNIVERSITY,
C FEBRUARY 1973.
      DOUBLE PRECISION PF(4,40), U(40), V(40)
      DOUBLE PRECISION SUM
      DIMENSION DXA(4)
      COMMON DX, DXA, PF, L
      DATA MNP /40/
      FN = N
      RN = 1./FN
      K = EPS*FN + .00000001
      FK = K
      IF (ABS(FK-EPS*FN).GT..00000001) GO TO 10
      K = K - 1
      FK = K
   10 CONTINUE
      DEL = EPS - FK*RN
      XUP = RN - DEL
      XLO = DEL
      IF (ABS(XUP-XLO).LT..00000001) XUP = XLO
      XPREV = 0.
      DO 20 I=1,MNP
        U(I) = 0.
   20 CONTINUE
      U(K+1) = 1.
      IMIN = -K
   30 X = AMIN1(XUP,XLO)
      IF (X.GT..999999) X = 1.
      DX = X - XPREV
      JMIN = CEIL((X-EPS)*FN-.00000001)
      IF (ABS(FLOAT(JMIN)-(X-EPS)*FN).LT..00000001)
     * JMIN = JMIN + 1
      JMAX = (X+EPS)*FN + .00000001
      IF (ABS(FLOAT(JMAX)-(X+EPS)*FN).LT..00000001)
     * JMAX = JMAX - 1
      JMAX = JMAX - JMIN + 1
      CALL PRFAC
      DO 60 J=1,MNP
        SUM = 0.
        IF (J.GT.JMAX) GO TO 50
        I = 1
   40   IP = J - I + 1 + JMIN - IMIN
        SUM = SUM + U(I)*PF(L,IP)
        I = I + 1
        IF ((IMIN+I).LE.(JMIN+J)) GO TO 40
   50   V(J) = SUM
   60 CONTINUE
      DO 70 I=1,MNP
        U(I) = V(I)
   70 CONTINUE
      IMIN = JMIN
      XPREV = X
      IF (X.EQ.XUP) XUP = XUP + RN
      IF (X.EQ.XLO) XLO = XLO + RN
      IF (X.LT.1.) GO TO 30
      DO 80 I=1,N
        U(K+1) = U(K+1)*FLOAT(I)
   80 CONTINUE
      PKS = U(K+1)
      RETURN
      END