      SUBROUTINE PCSHFT(A,B,D,N)
      DIMENSION D(N)
      CONST=2./(B-A)
      FAC=CONST
      DO 11 J=2,N
        D(J)=D(J)*FAC
        FAC=FAC*CONST
11    CONTINUE
      CONST=0.5*(A+B)
      DO 13 J=1,N-1
        DO 12 K=N-1,J,-1
          D(K)=D(K)-CONST*D(K+1)
12      CONTINUE
13    CONTINUE
      RETURN
      END
      SUBROUTINE PEARSN(X,Y,N,R,PROB,Z)
      PARAMETER (TINY=1.E-20)
      DIMENSION X(N),Y(N)
      AX=0.
      AY=0.
      DO 11 J=1,N
        AX=AX+X(J)
        AY=AY+Y(J)
11    CONTINUE
      AX=AX/N
      AY=AY/N
      SXX=0.
      SYY=0.
      SXY=0.
      DO 12 J=1,N
        XT=X(J)-AX
        YT=Y(J)-AY
        SXX=SXX+XT**2
        SYY=SYY+YT**2
        SXY=SXY+XT*YT
12    CONTINUE
      R=SXY/SQRT(SXX*SYY)
      Z=0.5*ALOG(((1.+R)+TINY)/((1.-R)+TINY))
      DF=N-2
      T=R*SQRT(DF/(((1.-R)+TINY)*((1.+R)+TINY)))
      PROB=BETAI(0.5*DF,0.5,DF/(DF+T**2))
C     PROB=ERFCC(ABS(Z*SQRT(N-1.))/1.4142136)
      RETURN
      END
      SUBROUTINE PIKSR2(N,ARR,BRR)
      DIMENSION ARR(N),BRR(N)
      DO 12 J=2,N
        A=ARR(J)
        B=BRR(J)
        DO 11 I=J-1,1,-1
          IF(ARR(I).LE.A)GO TO 10
          ARR(I+1)=ARR(I)
          BRR(I+1)=BRR(I)
11      CONTINUE
        I=0
10      ARR(I+1)=A
        BRR(I+1)=B
12    CONTINUE
      RETURN
      END
      SUBROUTINE PIKSRT(N,ARR)
      DIMENSION ARR(N)
      DO 12 J=2,N
        A=ARR(J)
        DO 11 I=J-1,1,-1
          IF(ARR(I).LE.A)GO TO 10
          ARR(I+1)=ARR(I)
11      CONTINUE
        I=0
10      ARR(I+1)=A
12    CONTINUE
      RETURN
      END
      SUBROUTINE PINVS(IE1,IE2,JE1,JSF,JC1,K,C,NCI,NCJ,NCK,S,NSI,NSJ)
      PARAMETER (ZERO=0.,ONE=1.,NMAX=10)
      DIMENSION C(NCI,NCJ,NCK),S(NSI,NSJ),PSCL(NMAX),INDXR(NMAX)
      JE2=JE1+IE2-IE1
      JS1=JE2+1
      DO 12 I=IE1,IE2
        BIG=ZERO
        DO 11 J=JE1,JE2
          IF(ABS(S(I,J)).GT.BIG) BIG=ABS(S(I,J))
11      CONTINUE
        IF(BIG.EQ.ZERO) PAUSE 'Singular matrix, row all 0'
        PSCL(I)=ONE/BIG
        INDXR(I)=0
12    CONTINUE
      DO 18 ID=IE1,IE2
        PIV=ZERO
        DO 14 I=IE1,IE2
          IF(INDXR(I).EQ.0) THEN
            BIG=ZERO
            DO 13 J=JE1,JE2
              IF(ABS(S(I,J)).GT.BIG) THEN
                JP=J
                BIG=ABS(S(I,J))
              ENDIF
13          CONTINUE
            IF(BIG*PSCL(I).GT.PIV) THEN
              IPIV=I
              JPIV=JP
              PIV=BIG*PSCL(I)
            ENDIF
          ENDIF
14      CONTINUE
        IF(S(IPIV,JPIV).EQ.ZERO) PAUSE 'Singular matrix'
        INDXR(IPIV)=JPIV
        PIVINV=ONE/S(IPIV,JPIV)
        DO 15 J=JE1,JSF
          S(IPIV,J)=S(IPIV,J)*PIVINV
15      CONTINUE
        S(IPIV,JPIV)=ONE
        DO 17 I=IE1,IE2
          IF(INDXR(I).NE.JPIV) THEN
            IF(S(I,JPIV).NE.ZERO) THEN
              DUM=S(I,JPIV)
              DO 16 J=JE1,JSF
                S(I,J)=S(I,J)-DUM*S(IPIV,J)
16            CONTINUE
              S(I,JPIV)=ZERO
            ENDIF
          ENDIF
17      CONTINUE
18    CONTINUE
      JCOFF=JC1-JS1
      ICOFF=IE1-JE1
      DO 21 I=IE1,IE2
        IROW=INDXR(I)+ICOFF
        DO 19 J=JS1,JSF
          C(IROW,J+JCOFF,K)=S(I,J)
19      CONTINUE
21    CONTINUE
      RETURN
      END
      FUNCTION PLGNDR(L,M,X)
      IF(M.LT.0.OR.M.GT.L.OR.ABS(X).GT.1.)PAUSE 'bad arguments'
      PMM=1.
      IF(M.GT.0) THEN
        SOMX2=SQRT((1.-X)*(1.+X))
        FACT=1.
        DO 11 I=1,M
          PMM=-PMM*FACT*SOMX2
          FACT=FACT+2.
11      CONTINUE
      ENDIF
      IF(L.EQ.M) THEN
        PLGNDR=PMM
      ELSE
        PMMP1=X*(2*M+1)*PMM
        IF(L.EQ.M+1) THEN
          PLGNDR=PMMP1
        ELSE
          DO 12 LL=M+2,L
            PLL=(X*(2*LL-1)*PMMP1-(LL+M-1)*PMM)/(LL-M)
            PMM=PMMP1
            PMMP1=PLL
12        CONTINUE
          PLGNDR=PLL
        ENDIF
      ENDIF
      RETURN
      END
      FUNCTION POIDEV(XM,IDUM)
      PARAMETER (PI=3.141592654)
      DATA OLDM /-1./
      IF (XM.LT.12.)THEN
        IF (XM.NE.OLDM) THEN
          OLDM=XM
          G=EXP(-XM)
        ENDIF
        EM=-1
        T=1.
2       EM=EM+1.
        T=T*RAN1(IDUM)
        IF (T.GT.G) GO TO 2
      ELSE
        IF (XM.NE.OLDM) THEN
          OLDM=XM
          SQ=SQRT(2.*XM)
          ALXM=ALOG(XM)
          G=XM*ALXM-GAMMLN(XM+1.)
        ENDIF
1       Y=TAN(PI*RAN1(IDUM))
        EM=SQ*Y+XM
        IF (EM.LT.0.) GO TO 1
        EM=INT(EM)
        T=0.9*(1.+Y**2)*EXP(EM*ALXM-GAMMLN(EM+1.)-G)
        IF (RAN1(IDUM).GT.T) GO TO 1
      ENDIF
      POIDEV=EM
      RETURN
      END
      SUBROUTINE POLCOE(X,Y,N,COF)
      PARAMETER (NMAX=15)
      DIMENSION X(N),Y(N),COF(N),S(NMAX)
      DO 11 I=1,N
        S(I)=0.
        COF(I)=0.
11    CONTINUE
      S(N)=-X(1)
      DO 13 I=2,N
        DO 12 J=N+1-I,N-1
          S(J)=S(J)-X(I)*S(J+1)
12      CONTINUE
        S(N)=S(N)-X(I)
13    CONTINUE
      DO 16 J=1,N
        PHI=N
        DO 14 K=N-1,1,-1
          PHI=K*S(K+1)+X(J)*PHI
14      CONTINUE
        FF=Y(J)/PHI
        B=1.
        DO 15 K=N,1,-1
          COF(K)=COF(K)+B*FF
          B=S(K)+X(J)*B
15      CONTINUE
16    CONTINUE
      RETURN
      END
      SUBROUTINE POLCOF(XA,YA,N,COF)
      PARAMETER (NMAX=15)
      DIMENSION XA(N),YA(N),COF(N),X(NMAX),Y(NMAX)
      DO 11 J=1,N
        X(J)=XA(J)
        Y(J)=YA(J)
11    CONTINUE
      DO 14 J=1,N
        CALL POLINT(X,Y,N+1-J,0.,COF(J),DY)
        XMIN=1.E38
        K=0
        DO 12 I=1,N+1-J
          IF (ABS(X(I)).LT.XMIN)THEN
            XMIN=ABS(X(I))
            K=I
          ENDIF
          IF(X(I).NE.0.)Y(I)=(Y(I)-COF(J))/X(I)
12      CONTINUE
        IF (K.LT.N+1-J) THEN
          DO 13 I=K+1,N+1-J
            Y(I-1)=Y(I)
            X(I-1)=X(I)
13        CONTINUE
        ENDIF
14    CONTINUE
      RETURN
      END
      SUBROUTINE POLDIV(U,N,V,NV,Q,R)
      DIMENSION U(N),V(NV),Q(N),R(N)
      DO 11 J=1,N
        R(J)=U(J)
        Q(J)=0.
11    CONTINUE
      DO 13 K=N-NV,0,-1
        Q(K+1)=R(NV+K)/V(NV)
        DO 12 J=NV+K-1,K+1,-1
          R(J)=R(J)-Q(K+1)*V(J-K)
12      CONTINUE
13    CONTINUE
      R(NV)=0.
      RETURN
      END
      SUBROUTINE POLIN2(X1A,X2A,YA,M,N,X1,X2,Y,DY)
      PARAMETER (NMAX=20,MMAX=20)
      DIMENSION X1A(M),X2A(N),YA(M,N),YNTMP(NMAX),YMTMP(MMAX)
      DO 12 J=1,M
        DO 11 K=1,N
          YNTMP(K)=YA(J,K)
11      CONTINUE
        CALL POLINT(X2A,YNTMP,N,X2,YMTMP(J),DY)
12    CONTINUE
      CALL POLINT(X1A,YMTMP,M,X1,Y,DY)
      RETURN
      END
      SUBROUTINE POLINT(XA,YA,N,X,Y,DY)
      PARAMETER (NMAX=10) 
      DIMENSION XA(N),YA(N),C(NMAX),D(NMAX)
      NS=1
      DIF=ABS(X-XA(1))
      DO 11 I=1,N 
        DIFT=ABS(X-XA(I))
        IF (DIFT.LT.DIF) THEN
          NS=I
          DIF=DIFT
        ENDIF
        C(I)=YA(I)
        D(I)=YA(I)
11    CONTINUE
      Y=YA(NS)
      NS=NS-1
      DO 13 M=1,N-1
        DO 12 I=1,N-M
          HO=XA(I)-X
          HP=XA(I+M)-X
          W=C(I+1)-D(I)
          DEN=HO-HP
          IF(DEN.EQ.0.)PAUSE
          DEN=W/DEN
          D(I)=HP*DEN
          C(I)=HO*DEN
12      CONTINUE
        IF (2*NS.LT.N-M)THEN
          DY=C(NS+1)
        ELSE
          DY=D(NS)
          NS=NS-1
        ENDIF
        Y=Y+DY
13    CONTINUE
      RETURN
      END
      SUBROUTINE POWELL(P,XI,N,NP,FTOL,ITER,FRET)
      PARAMETER (NMAX=20,ITMAX=200)
      DIMENSION P(NP),XI(NP,NP),PT(NMAX),PTT(NMAX),XIT(NMAX)
      FRET=FUNC(P)
      DO 11 J=1,N
        PT(J)=P(J)
11    CONTINUE
      ITER=0
1     ITER=ITER+1
      FP=FRET
      IBIG=0
      DEL=0.
      DO 13 I=1,N
        DO 12 J=1,N
          XIT(J)=XI(J,I)
12      CONTINUE
        CALL LINMIN(P,XIT,N,FRET)
        IF(ABS(FP-FRET).GT.DEL)THEN
          DEL=ABS(FP-FRET)
          IBIG=I
        ENDIF
13    CONTINUE
      IF(2.*ABS(FP-FRET).LE.FTOL*(ABS(FP)+ABS(FRET)))RETURN
      IF(ITER.EQ.ITMAX) PAUSE 'Powell exceeding maximum iterations.'
      DO 14 J=1,N
        PTT(J)=2.*P(J)-PT(J)
        XIT(J)=P(J)-PT(J)
        PT(J)=P(J)
14    CONTINUE
      FPTT=FUNC(PTT)
      IF(FPTT.GE.FP)GO TO 1
      T=2.*(FP-2.*FRET+FPTT)*(FP-FRET-DEL)**2-DEL*(FP-FPTT)**2
      IF(T.GE.0.)GO TO 1
      CALL LINMIN(P,XIT,N,FRET)
      DO 15 J=1,N
        XI(J,IBIG)=XIT(J)
15    CONTINUE
      GO TO 1
      END
      SUBROUTINE PREDIC(DATA,NDATA,D,NPOLES,FUTURE,NFUT)
      PARAMETER (NPMAX=100)
      DIMENSION DATA(NDATA),D(NPOLES),FUTURE(NFUT),REG(NPMAX)
      DO 11 J=1,NPOLES
        REG(J)=DATA(NDATA+1-J)
11    CONTINUE
      DO 14 J=1,NFUT
        DISCRP=0.
        SUM=DISCRP
        DO 12 K=1,NPOLES
          SUM=SUM+D(K)*REG(K)
12      CONTINUE
        DO 13 K=NPOLES,2,-1
          REG(K)=REG(K-1)
13      CONTINUE
        REG(1)=SUM
        FUTURE(J)=SUM
14    CONTINUE
      RETURN
      END
      FUNCTION PROBKS(ALAM)
      PARAMETER (EPS1=0.001, EPS2=1.E-8)
      A2=-2.*ALAM**2
      FAC=2.
      PROBKS=0.
      TERMBF=0.
      DO 11 J=1,100
        TERM=FAC*EXP(A2*J**2)
        PROBKS=PROBKS+TERM
        IF(ABS(TERM).LT.EPS1*TERMBF.OR.ABS(TERM).LT.EPS2*PROBKS)RETURN
        FAC=-FAC
        TERMBF=ABS(TERM)
11    CONTINUE
      PROBKS=1.
      RETURN
      END
      SUBROUTINE PZEXTR(IEST,XEST,YEST,YZ,DY,NV,NUSE)
      PARAMETER (IMAX=11,NCOL=7,NMAX=10)
      DIMENSION X(IMAX),YEST(NV),YZ(NV),DY(NV),QCOL(NMAX,NCOL),D(NMAX)
      X(IEST)=XEST
      DO 11 J=1,NV
        DY(J)=YEST(J)
        YZ(J)=YEST(J)
11    CONTINUE
      IF(IEST.EQ.1) THEN
        DO 12 J=1,NV
          QCOL(J,1)=YEST(J)
12      CONTINUE
      ELSE
        M1=MIN(IEST,NUSE)
        DO 13 J=1,NV
          D(J)=YEST(J)
13      CONTINUE
        DO 15 K1=1,M1-1
          DELTA=1./(X(IEST-K1)-XEST)
          F1=XEST*DELTA
          F2=X(IEST-K1)*DELTA
          DO 14 J=1,NV
            Q=QCOL(J,K1)
            QCOL(J,K1)=DY(J)
            DELTA=D(J)-Q
            DY(J)=F1*DELTA
            D(J)=F2*DELTA
            YZ(J)=YZ(J)+DY(J)
14        CONTINUE
15      CONTINUE
        DO 16 J=1,NV
          QCOL(J,M1)=DY(J)
16      CONTINUE
      ENDIF
      RETURN
      END
