PARAMETER (NMAX=1000,MMAX=40,EPS=0.02)
	DIMENSION X(NMAX*MMAX),R(MMAX*MMAX),RT(MMAX*MMAX)
	DIMENSION XX(NMAX*MMAX),AA(MMAX)
	DIMENSION A(MMAX*MMAX),S(MMAX),Q(MMAX),U(MMAX)
	DIMENSION RR(MMAX),F(MMAX)
	DIMENSION INDEX(MMAX),INDEX1(MMAX),INDEX2(MMAX),IDX(MMAX)
	INTEGER N,M

	INDEX=1
	INDEX1=-1
	INDEX1(1)=1
	INDEX2(1)=1
!	INDEX(1)=0
	MC=1

	OPEN(3,FILE='C:\MASTER\G_TR\Bai Tap TK trong KH\New\KQ.TXT',
     &     STATUS='UNKNOWN')
	OPEN(1,FILE='C:\MASTER\G_TR\Bai Tap TK trong KH\New\SL.DAT')
	READ(1,*)
	READ(1,*)
	READ(1,*) N,M
	DO I=1,N
	   L=(M-1)*N+I
	   READ(1,*)  (X(J),J=I,L,N)
	ENDDO
	
	CALL CORRE(X,N,M,R)
	WRITE(*,'(4F8.4/)') (R(JK),JK=1,M*M)
	WRITE(*,*)
! XAC DINH HE SO TQ MAX (VOI X1)
	KMAX=2
	RMAX=ABS(R(KMAX))
	DO J=3,M
	   IF (ABS(R(J)).GT.RMAX) THEN
	      KMAX=J
	      RMAX=ABS(R(J))
	   ENDIF
	ENDDO
	
	INDEX1(KMAX)=1
	IDX=INDEX1       ! SAO LUU
	MC=MC+1
	INDEX2(MC)=KMAX  ! BIEN DUA VAO DAU TIEN
	INDEX(KMAX)=0    ! LOAI BO BIEN NAY
	
!	PRINT*,INDEX
!	PRINT*,INDEX1
!	PRINT*,INDEX2
!	PRINT*,IDX

	CALL SUBMTX(R,M,RT,IDX)
!	WRITE(*,'(2F8.4/)') (RT(JK),JK=1,MC*MC)

	CALL SUBMTX1(X,N,M,XX,IDX)

	CALL REGRE(XX,N,MC,AA,SS,QQ,UU,RRR,FF)

	DO J=1,MC
	   A(J)=AA(J)
	ENDDO
	S(1)=SS
	Q(1)=QQ
	U(1)=UU
	RR(1)=RRR
	F(1)=FF
	
!	WRITE(*,'(2F12.2/)') (A(JK),JK=1,MC)

!	PAUSE

! BAT DAU VONG LAP
	K=1
  100 K=K+1
	RMAX=0.0
	DO J=2,M
	   IDX=INDEX1
	   IF (INDEX1(J).NE.1) THEN 
	      IDX(J)=1
	      CALL SUBMTX(R,M,RT,IDX)
	      MM=MC+1
		  TQR=HSTQR(RT,MM,1,MM)
	      IF (ABS(TQR).GT.RMAX) THEN
	         RMAX=ABS(TQR)
	         KMAX=J
	      ENDIF
	   ENDIF
	ENDDO
	INDEX1(KMAX)=1
	IDX=INDEX1       ! SAO LUU
	MC=MC+1
	INDEX2(MC)=KMAX  ! BIEN DUA VAO DAU TIEN
	INDEX(KMAX)=0    ! LOAI BO BIEN NAY

	CALL SUBMTX(R,M,RT,IDX)
	CALL SUBMTX1(X,N,M,XX,IDX)
	DO I=1,N
!	   WRITE(*,*)(XX(IJ),IJ=I,N*MC,N)
	ENDDO
	CALL REGRE(XX,N,MC,AA,SS,QQ,UU,RRR,FF)
	
	DO JJ=1,MC
	   A((K-1)*M+JJ)=AA(JJ)
	ENDDO
	S(K)=SS
	Q(K)=QQ
	U(K)=UU
	RR(K)=RRR
	F(K)=FF
	SSO=ABS( (S(K)-S(K-1))/S(K) )
!	PRINT*,K,SSO
!	PAUSE
	IF (SSO.GE.EPS) THEN 
	   IF (K.LT.M-1) GOTO 100
	ENDIF

	DO K=1,MC
	   WRITE(3,'(10F10.3)') (A(J),J=(K-1)*M+1,(K-1)*M+K+1)
	ENDDO
	DO K=1,MC
	   WRITE(3,'(3F10.3)')S(K),RR(K),F(K)
	ENDDO

!	DO I=1,N
!	WRITE(*,'(2F18.4)') (XX(JK),JK=I,N*MC,N)
!	ENDDO
	END

!+++++++++++++++++++++++++++++

	FUNCTION HSTQR(R,M,J,K)
! HAM NAY TINH HE SO TUONG QUAN RIENG GIUA BIEN XJ VA XK
! KHI XET DONG THOI M BIEN X1,X2,...,XM
! INPUT: + R MANG MOT CHIEU CHUA MA TRAN TUONG QUAN CHUAN HOA
!          (KICH THUOC M*M ==> R(M*M)
!        + M KICH THUOC MA TRAN TUONG QUAN, DONG THOI CUNG LA
!          SO BIEN DUOC XET
!        + J,K CHI SO BIEN DUOC TINH TUONG QUAN RIENG
! OUTPUT: HE SO TUONG QUAN RIENG GIUA BIEN THU J VA THU K
!
	DIMENSION R(M*M),RT(M*M)
	RT=R
	RJK=PPDS(RT,M,J,K)
	RT=R
	RJJ=PPDS(RT,M,J,J)
	RT=R
	RKK=PPDS(RT,M,K,K)
	HSTQR=-RJK/(SQRT(RJJ*RKK))
	RETURN
	END

	FUNCTION HSTQB(R,M,J)
! HAM NAY TINH HE SO TUONG QUAN BOI GIUA BIEN XJ VA TAP HOP
!        M-1 BIEN X1,X2,X(J-1),X(J+1),..,XM 
! INPUT: + R MANG MOT CHIEU CHUA MA TRAN TUONG QUAN CHUAN HOA
!          (KICH THUOC M*M ==> R(M*M)
!        + M KICH THUOC MA TRAN TUONG QUAN, DONG THOI CUNG LA
!          SO BIEN DUOC XET
!        + J CHI SO BIEN DUOC TINH TUONG QUAN BOI VOI TAP M-1 
!            BIEN CON LAI
! OUTPUT: HE SO TUONG QUAN BOI GIUA BIEN THU J VA CAC BIEN KHAC
!
	DIMENSION R(M*M),RT(M*M)
	RT=R
	RJJ=PPDS(RT,M,J,J)
	RT=R
	RR=DET(RT,M)
	HSTQB=SQRT(1-RR/RJJ)
	RETURN
	END

	SUBROUTINE SUBMTX1(X,N,M,XC,INDEX)
	DIMENSION X(N*M),XC(N*M),INDEX(M)
	L=0
	DO 10 K=1,M
	KK=0
  1	KK=KK+1
	IF (KK.GT.M) GOTO 10
	IF (INDEX(KK).NE.1) GOTO 1
	IF (KK.NE.K) GOTO 1

	DO J=1,N
	   JK=(KK-1)*N+J
	   L=L+1
	   XC(L)=X(JK)
	ENDDO
  10  CONTINUE
	RETURN
	END	

	SUBROUTINE SUBMTX(R,M,RC,INDEX)
	DIMENSION R(M*M),RC(M*M),INDEX(M)
	L=0
	DO 10 K=1,M
	KK=0
  1	KK=KK+1
	IF (KK.GT.M) GOTO 10
	IF (INDEX(KK).NE.1) GOTO 1
	IF (KK.NE.K) GOTO 1

	DO 11 J=1,M
	JJ=0
  2	JJ=JJ+1
	IF (JJ.GT.M) GOTO 11
	IF (INDEX(JJ).NE.1) GOTO 2
	IF (JJ.NE.J) GOTO 2
	JK=(KK-1)*M+JJ
	L=L+1
	RC(L)=R(JK)
  11	CONTINUE

  10  CONTINUE
	RETURN
	END	


	SUBROUTINE CORRE(X,N,M,R)
! CHUONG TRINH NAY TINH MA TRAN TUONG QUAN CHUAN HOA
! CUA TAP M BIEN TU SO LIEU BAN DAU CO DUNG LUONG N
! INPUT: + X MANG MOT CHIEU KICH THUOC N*M LUU TRU SO
!          LIEU BAN DAU DANG MA TRAN N HANG M COT
!          (X(1,1), X(2,1),...,X(N,1),X(1,2),...)
!        + N DUNG LUONG MAU
!        + M SO BIEN
! OUTPUT: + R MANG MOT CHIEU KICH THUOC M*M LUU TRU MA TRAN
!            TUONG QUAN CHUAN HOA CUA M BIEN
!
	DIMENSION X(N*M),XX(N,M),R(M*M),TB(M),SX(M)
	K=0
	DO J=1,M
	   TB(J)=0.0
	   SX(J)=0.0
	   DO I=1,N
	      K=K+1
	      TB(J)=TB(J)+X(K)
	      SX(J)=SX(J)+X(K)*X(K)
	      XX(I,J)=X(K)
	   ENDDO
	   TB(J)=TB(J)/REAL(N)
	   SX(J)=SQRT(SX(J)/REAL(N)-TB(J)*TB(J))
	ENDDO
	JK=0
	DO J=1,M
	   DO K=1,M
	      JK=JK+1
	      R(JK)=0.0
	      DO I=1,N
	         R(JK)=R(JK)+XX(I,J)*XX(I,K)
	      ENDDO
	      R(JK)=R(JK)/REAL(N)-TB(J)*TB(K)
	      R(JK)=R(JK)/(SX(J)*SX(K))
	   ENDDO
	ENDDO
	RETURN
	END

	SUBROUTINE REGRE(X,N,M,A,S,Q,U,RR,F)
! CHUONG TRINH NAY TINH CAC HE SO HOI QUI CUA PHUONG TRINH
! HOI QUI GIUA X1 VA CAC X2,...,XM VA XAC DINH CAC DAI LUONG
! CHO PHAN TICH PHUONG SAI
! DANG PHUONG TRINH HQ: X1=A1+A2*X2+...+AM*XM
!
! INPUT: + X MANG MOT CHIEU KICH THUOC N*M LUU TRU SO LIEU
!            QUAN TRAC CUA CAC X1..XM (LUU THEO COT)
!        + N DUNG LUONG MAU
!        + M SO BIEN (1 BIEN PHU THUOC VA M-1 BIEN DOC LAP)
!        + A MANG MOT CHIEU KICH THUOC M LUU CAC HE SO
!            A1, A2,..., AM
!        + S SAI SO CHUAN (CHUAN SAI THANG DU)
!        + Q TONG BINH PHUONG CAC BIEN SAI THANG DU
!        + U TONG BINH PHUONG CAC BIEN SAI HOI QUI
!        + RR HE SO TUONG QUAN BOI
!        + F BIEN CO PHAN BO FISHER
!
	DIMENSION X(N*M),A(M)
	DIMENSION R(M*M),RXX(M*M),RY(M),TB(M),RX(M*M)

	CALL COVAR(X,N,M,R,TB)
	L=0
	J=0
	DO I=M+1,M*M
	   IF (MOD(I,M).NE.1) THEN
	      L=L+1
	      RXX(L)=R(I)
	   ELSE
	      J=J+1
	      RY(J)=R(I)
	   ENDIF
	ENDDO
	RX=RXX
	CALL MINVERT(RX,M-1)
	CALL MULTIP(RX,RY,A,M-1,M-1,1)
	TMP=0.0
	DO J=1,M-1
	   TMP=TMP+A(J)*TB(J+1)
	ENDDO
	TMP=TB(1)-TMP
	DO J=M,2,-1
	   A(J)=A(J-1)
	ENDDO
	A(1)=TMP
	
	RX=R
	D=DET(RX,M)
	RX=R
	D11=PPDS(RX,M,1,1)
	TMP=D/D11

	RR=SQRT(1-TMP/R(1))

	Q=0.0
	DO I=1,N
	   YHQ=0.0
	   DO J=2,M
	      IJ=(J-1)*N+I
	      YHQ=YHQ+A(J)*X(IJ)
	   ENDDO
	   YHQ=YHQ+A(1)
	   TMP=X(I)-YHQ
	   Q=Q+TMP*TMP
	ENDDO
	U=REAL(N)*R(1)-Q
	F=U*REAL(N-M)/(Q*REAL(M-1))
	S=SQRT(Q/REAL(N-M))

	RETURN
	END

	SUBROUTINE COVAR(X,N,M,R,TB)
! CHUONG TRINH NAY TINH MA TRAN TUONG QUAN (MA TRAN COVARIANCE)
! CUA TAP M BIEN
! INPUT: + X MANG MOT CHIEU KICH THUOC N*M (N HANG, M COT)
!          LUU TRU SO LIEU CUA X THEO COT
!        + N SO HANG (DUNG LUONG MAU) CUA X
!        + M SO COT (SO BIEN)
! OUTPUT: + R MANG MOT CHIEU KICH THUOC (M*M) CHUA CAC 
!           MOMEN TUONG QUAN (COVARIANCE) GIUA CAC BIEN
!         + TB MANG MOT CHIEU KICH THUOC M CHUA TRUNG BINH CAC BIEN
!
	DIMENSION X(N*M),XX(N,M),R(M*M),RR(M,M),TB(M)

	K=0
	DO J=1,M
	   TB(J)=0.0
	   DO I=1,N
	      K=K+1
	      TB(J)=TB(J)+X(K)
	      XX(I,J)=X(K)
	   ENDDO
	   TB(J)=TB(J)/REAL(N)
	ENDDO
	DO J=1,M
	   DO K=J,M
	      RR(J,K)=0.0
	      DO I=1,N
	         RR(J,K)=RR(J,K)+(XX(I,J)-TB(J))*(XX(I,K)-TB(K))
	      ENDDO
	      RR(J,K)=RR(J,K)/REAL(N)
	      RR(K,J)=RR(J,K)
	   ENDDO
	ENDDO
	L=0
	DO K=1,M
	   DO J=1,M
	      L=L+1
		  R(L)=RR(J,K)
	   ENDDO
	ENDDO	
	RETURN
	END

	SUBROUTINE MINVERT(AA,N)
! CHUONG TRINH NAY TINH MA TRAN NGHICH DAO CUA MA TRAN A
! INPUT: + AA MANG MOT CHIEU KICH THUOC N*N LUU THEO COT CUA A
!        + N SO HANG/COT CUA MA TRAN A
! OUTPUT: AA MA TRAN NGHICH DAO CUA A
!
	DIMENSION AA(N*N), A(N,N)

	K=0
	DO J=1,N
	   DO I=1,N
	      K=K+1
	      A(I,J)=AA(K)
	   ENDDO
	ENDDO
	M=N-1
	DO 10 I=1,N
	   R=A(I,1)
	   DO K=1,M
	      A(I,K)=A(I,K+1)/R
	   ENDDO
	   A(I,N)=1.0/R
	   DO 10 J=1,N
	      IF (I.NE.J) THEN
	         R=A(J,1)
	         DO K=1,M
	            A(J,K)=A(J,K+1)-R*A(I,K)
	         ENDDO
	         A(J,N)=-R*A(I,N)
	      ENDIF
  10  CONTINUE
	K=0
	DO J=1,N
	   DO I=1,N
	      K=K+1
	      AA(K)=A(I,J)
	   ENDDO
	ENDDO
	RETURN
	END

	SUBROUTINE MULTIP(A,B,C,N,M,L)
! CHUONG TRINH NAY TINH TICH CUA HAI MA TRAN
! INPUT: + A MANG MOT CHIEU KICH THUOC N*M (N HANG, M COT)
!            LUU MA TRAN A THEO COT
!        + B MANG MOT CHIEU KICH THUOC M*L (M HANG, L COT)
!            LUU MA TRAN B THEO COT
!        + N SO HANG CUA MA TRAN A
!        + M SO COT CUA A VA LA SO HANG CUA B
!        + L SO COT CUA B
! OUTPUT: + C MANG MOT CHIEU KICH THUOC N*L (N HANG, L COT)
!             LUU MA TRAN C THEO COT

	DIMENSION A(N*M),B(M*L),C(N*L)

	DO I=1,N
	DO K=1,L
	   IK=(K-1)*N+I
	   C(IK)=0.0
	   DO J=1,M
	      IJ=(J-1)*N+I   !(I-1)*M+J
	      JK=(K-1)*M+J
	      C(IK)=C(IK)+A(IJ)*B(JK)
	   ENDDO
	ENDDO
	ENDDO

	RETURN
	END
		
	FUNCTION PPDS(A,N,IR,JC)
!
! HAM NAY TINH PHAN PHU DAI SO CUA PHAN TU B(IR,JC)
! (HANG IR, COT JC) CUA MA TRAN B(NxN) MA CAC PHAN TU
! CUA NO DUOC LUU TRONG MANG A(N*N) THEO QUI CACH COT
! TRUOC DONG SAU
! INPUT: + MANG A(N*N) CHUA MA TRAN DAU VAO CUA B
!        + N KICH THUOC CUA B
!        + IR CHI SO HANG
!        + JC CHI SO COT
! OUTPUT: PHAN PHU DAI SO CUA PHAN TU HANG IR COT JC
!
	DIMENSION A(N*N), B(N,N)
	K=0
	DO J=1,N
	   DO I=1,N
	      K=K+1
	      B(I,J)=A(K)
	   ENDDO
	ENDDO
	
	K=0
	DO J=1,N
	   IF (J.NE.JC) THEN
	   DO I=1,N
	      IF (I.NE.IR) THEN
	         K=K+1
	         A(K)=B(I,J)
	      ENDIF
	   ENDDO
	   ENDIF
	ENDDO
	
	PPDS=DET(A,N-1)*(-1)**(IR+JC)
	RETURN
	END


	FUNCTION DET(X,N)
! HAM NAY TINH DINH THUC CUA MA TRAN A(N,N)
! INPUT: + MANG X DO DAI N*N CHUA CAC PHAN TU CUA MA TRAN A
!        + N KICH THUOC MA TRAN A
!   CHU Y: X(1)=A(1,1), X(2)=A(2,1),..., X(N)=A(N,1),X(N+1)=A(1,2),
!          X(N+2)=A(2,2),...
! OUTPUT: DINH THUC CUA A
!
	PARAMETER (EP=1.0E-6)
	DIMENSION X(N*N),A(N,N)
	K=0
	DO J=1,N
	   DO I=1,N
	      K=K+1
	      A(I,J)=X(K)
	   ENDDO
	ENDDO
	D=1.0
	N1=N-1
	DO 10 K=1,N1
	   AM=0.0
	   DO 11 I=K,N
	      T=A(I,K)
	      IF (ABS(T).GE.ABS(AM)) THEN
	         AM=T
	         J=I
	      ENDIF
  11     CONTINUE
         IF (ABS(AM).LE.EP) THEN
	      DT=0.0
	      DET=DT
	      RETURN
	   ELSE
	      IF (J.NE.K) THEN
	         D=-D
	         DO I=K,N
	            T=A(J,I)
	            A(J,I)=A(K,I)
	            A(K,I)=T
	         ENDDO
	      ENDIF
	   ENDIF
	   M=K+1
	   DO I=M,N
	      T=A(I,K)/AM
	      DO J=M,N
	         A(I,J)=A(I,J)-T*A(K,J)
		  ENDDO
	   ENDDO
	   D=D*A(K,K)
  10  CONTINUE
      DT=D*A(N,N)
	DET=DT
	RETURN
	END

	  FUNCTION FINV(P,N1,N2)
! HAM NAY TINH GIA TRI X0 CUA BIEN NGAU NHIEN X PHAN BO 
! FISHER (F) VOI N1 VA N2 BAC TU DO THOA MAN DIEU KIEN 
! P(X>X0)=P.
! INPUT: + P XAC SUAT DE X>X0 (P(X>X0)=P)
!        + N1 VA N2 LA SO BAC TU DO, PHAI LA NHUNG SO THUC
! OUTPUT: X0
! SUBROUTINE/FUNCTION DUOC GOI TOI: FDIST
	  PARAMETER (EPS=1.0E-6)
	  REAL P,AP,N1,N2, A,B,C,P0,SS,SS1
	  IF (P.LT.0.0.OR.P.GT.1.0) THEN
	     WRITE(*,*)' INVALID NUMERIC INPUT IN FINV FUNCTION'
		 STOP
	  ENDIF
	  IF (P.EQ.0.0) THEN
		 WRITE(*,*)' THE P VALUE TOO SMALL IN FINV FUNCTION'
		 STOP
	  ELSE IF (P.EQ.1) THEN
	     FINV=0.0
		 RETURN
	  ENDIF
!
	  A=0.0
	  B=9999999.0
	  C=A
	  AP=P

 10	  P0=FDIST(B,N1,N2)
	  SS=ABS((P0-AP)/AP)
	  SS1=ABS((B-C)/B)
!	  PRINT*,'P0=',P0,' SS=',SS,' B=',B,' C=',C
	  IF (SS.GE.EPS.AND.SS1.GE.EPS) THEN
	     IF (AP.GT.P0) THEN
		    C=B
			B=(A+B)/2.0
			GOTO 10
		  ELSE IF (AP.LT.P0) THEN
			B=(C+B)/2.0
			GOTO 10
		  ENDIF
	  ENDIF
	  A=(C+B)/2.0
      FINV=A
	  RETURN
	  END


	  FUNCTION FDIST(X0,N1,N2)
! HAM NAY TINH XAC SUAT DE BIEN NGAU NHIEN X CO PHAN BO FISHER
! VOI N1 VA N2 BAC TU DO NHAN GIA TRI >X0:
! P=P(X>X0)=1-P(X<X0)
! INPUT: + X0 MOT GIA TRI NAO DO CUA X
!        + N1,N2 SO BAC TU DO, PHAI LA NHUNG SO THUC
! SUBROUTINE/FUNCTION DUOC GOI TOI: BETAI
!
	  REAL X0,N1,N2,X
	  X=N2/(N2+N1*X0)
	  FDIST=BETAI(0.5*N2,0.5*N1,X)
	  RETURN
	  END 

      FUNCTION betai(a,b,x)
      REAL betai,a,b,x
!CU    USES betacf,gammln
      REAL bt,betacf,gammln
      if(x.lt.0..or.x.gt.1.)pause 'bad argument x in betai'
      if(x.eq.0..or.x.eq.1.)then
        bt=0.
      else
        bt=exp(gammln(a+b)-gammln(a)-gammln(b)+a*log(x)+b*log(1.-x))
      endif
      if(x.lt.(a+1.)/(a+b+2.))then
        betai=bt*betacf(a,b,x)/a
        return
      else
        betai=1.-bt*betacf(b,a,1.-x)/b
        return
      endif
      END

      FUNCTION betacf(a,b,x)
      INTEGER MAXIT
      REAL betacf,a,b,x,EPS,FPMIN
      PARAMETER (MAXIT=100,EPS=3.e-7,FPMIN=1.e-30)
      INTEGER m,m2
      REAL aa,c,d,del,h,qab,qam,qap
      qab=a+b
      qap=a+1.
      qam=a-1.
      c=1.
      d=1.-qab*x/qap
      if(abs(d).lt.FPMIN)d=FPMIN
      d=1./d
      h=d
      do 11 m=1,MAXIT
        m2=2*m
        aa=m*(b-m)*x/((qam+m2)*(a+m2))
        d=1.+aa*d
        if(abs(d).lt.FPMIN)d=FPMIN
        c=1.+aa/c
        if(abs(c).lt.FPMIN)c=FPMIN
        d=1./d
        h=h*d*c
        aa=-(a+m)*(qab+m)*x/((a+m2)*(qap+m2))
        d=1.+aa*d
        if(abs(d).lt.FPMIN)d=FPMIN
        c=1.+aa/c
        if(abs(c).lt.FPMIN)c=FPMIN
        d=1./d
        del=d*c
        h=h*del
        if(abs(del-1.).lt.EPS)goto 1
11    continue
      pause 'a or b too big, or MAXIT too small in betacf'
1     betacf=h
      return
      END

      FUNCTION gammq(a,x)
      REAL a,gammq,x
!CU    USES gcf,gser
      REAL gammcf,gamser,gln
      if(x.lt.0..or.a.le.0.)pause 'bad arguments in gammq'
      if(x.lt.a+1.)then
        call gser(gamser,a,x,gln)
        gammq=1.-gamser
      else
        call gcf(gammcf,a,x,gln)
        gammq=gammcf
      endif
      return
      END

      SUBROUTINE gcf(gammcf,a,x,gln)
      INTEGER ITMAX
      REAL a,gammcf,gln,x,EPS,FPMIN
      PARAMETER (ITMAX=100,EPS=3.e-7,FPMIN=1.e-30)
!CU    USES gammln
      INTEGER i
      REAL an,b,c,d,del,h,gammln
      gln=gammln(a)
      b=x+1.-a
      c=1./FPMIN
      d=1./b
      h=d
      do 11 i=1,ITMAX
        an=-i*(i-a)
        b=b+2.
        d=an*d+b
        if(abs(d).lt.FPMIN)d=FPMIN
        c=b+an/c
        if(abs(c).lt.FPMIN)c=FPMIN
        d=1./d
        del=d*c
        h=h*del
        if(abs(del-1.).lt.EPS)goto 1
11    continue
      pause 'a too large, ITMAX too small in gcf'
1     gammcf=exp(-x+a*log(x)-gln)*h
      return
      END


      SUBROUTINE gser(gamser,a,x,gln)
      INTEGER ITMAX
      REAL a,gamser,gln,x,EPS
      PARAMETER (ITMAX=100,EPS=3.e-7)
!CU    USES gammln
      INTEGER n
      REAL ap,del,sum,gammln
      gln=gammln(a)
      if(x.le.0.)then
        if(x.lt.0.)pause 'x < 0 in gser'
        gamser=0.
        return
      endif
      ap=a
      sum=1./a
      del=sum
      do 11 n=1,ITMAX
        ap=ap+1.
        del=del*x/ap
        sum=sum+del
        if(abs(del).lt.abs(sum)*EPS)goto 1
11    continue
      pause 'a too large, ITMAX too small in gser'
1     gamser=sum*exp(-x+a*log(x)-gln)
      return
      END

      FUNCTION gammln(xx)
      REAL gammln,xx
      INTEGER j
      DOUBLE PRECISION ser,stp,tmp,x,y,cof(6)
      SAVE cof,stp
      DATA cof,stp/76.18009172947146d0,-86.50532032941677d0, 
     & 24.01409824083091d0,-1.231739572450155d0,.1208650973866179d-2, 
     & -.5395239384953d-5,2.5066282746310005d0/
      x=xx
      y=x
      tmp=x+5.5d0
      tmp=(x+0.5d0)*log(tmp)-tmp
      ser=1.000000000190015d0
      do 11 j=1,6
        y=y+1.d0
        ser=ser+cof(j)/y
11    continue
      gammln=tmp+log(stp*ser/x)
      return
      END