PROGRAM BT5_4
      PARAMETER (NMAX=100,MMAX=6)
	DIMENSION X(NMAX*MMAX)
	DIMENSION R(MMAX*MMAX),RC(MMAX*MMAX)
	CHARACTER VARNAME(6)*2, ANSWER(2)*16,ANS*16
	REAL ALFA
	DATA VARNAME /'X1','X2','X3','Y1','Y2','Y3'/
	DATA ANSWER /'LON RO RET','KHONG LON RO RET'/
	DATA ALFA /0.05/

	OPEN(1,FILE='C:\MASTER\G_TR\BAITAP~1\NEW\BANG2.TXT')
	READ(1,*)
	READ(1,*)
	READ(1,*) N,M
	DO I=1,N
	   L=(M-1)*N+I
	   READ(1,*) TMP, (X(J),J=I,L,N)
	ENDDO

	CALL CORRE(X,N,M,R)

	TALFA=TINV(ALFA,REAL(N-2))
	WRITE(*,*)'KET QUA KIEM NGHIEM DO RO RET CUA HE SO TUONG QUAN'
	WRITE(*,*)
	WRITE(*,'("    MUC Y NGHIA ALFA= ",F6.3)') ALFA
	DO J=1,M-1
	   DO K=(J-1)*M+J+1, J*M
	      T=R(K)*SQRT(REAL(N-2))/SQRT(1.0-R(K)**2)
	      IF (T.GT.TALFA) THEN
	         ANS=ANSWER(1)
	      ELSE
	         ANS=ANSWER(2)
            ENDIF
	      L=MOD(K,M)
	      IF (L.EQ.0) L=6
	      WRITE(*,300) VARNAME(J),VARNAME(L),R(K),ANS
	   ENDDO
	ENDDO
  300 FORMAT(' R(',A2,'-',A2,')=',F7.4,'==>',A17)
	
	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))
!	PRINT*,TB(J),SX(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

	FUNCTION TINV(P,N)
! HAM NAY TINH GIA TRI X0 CUA BIEN NGAU NHIEN X PHAN BO 
! STUDENT VOI N BAC TU DO THOA MAN DIEU KIEN 
! P(|X|>X0)=P.
! INPUT: + P XAC SUAT DE |X|>X0 (P(|X|>X0)=P)
!        + N SO BAC TU DO, PHAI LA MOT SO THUC
! OUTPUT: X0
! SUBROUTINE/FUNCTION DUOC GOI TOI: TDIST
	  PARAMETER (EPS=1.0E-5)
	  REAL P,AP,N, A,B,C,P0
	  IF (P.LT.0.0.OR.P.GT.1.0) THEN
	     WRITE(*,*)' INVALID NUMERIC INPUT IN AKBPINV FUNCTION'
		 STOP
	  ENDIF
	  IF (P.EQ.0.0) THEN
		 WRITE(*,*)' THE P VALUE TOO SMALL IN AKBPINV FUNCTION'
		 STOP
	  ELSE IF (P.EQ.1) THEN
	     WRITE(*,*)' THE P VALUE TOO BIG IN AKBPINV FUNCTION'
		 STOP
	  ENDIF
!
	  A=0.0
	  B=99999.0
!	  C=A
	  AP=P

 10	  P0=TDIST(B,N)
	  SS=ABS((P0-AP)/AP)
	  SS1=ABS((B-C)/B)
	  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
      TINV=A
	  RETURN
	  END


	  FUNCTION TDIST(X0,N)
! HAM NAY TINH XAC SUAT DE BIEN NGAU NHIEN X CO PHAN BO STUDENT
! VOI N BAC TU DO NHAN GIA TRI NGOAI KHOANG (-X0;X0):
! P=P(|X|>X0)=1-P(|X|<X0)
! INPUT: + X0 MOT GIA TRI NAO DO CUA X
!        + N SO BAC TU DO, PHAI LA MOT SO THUC
! SUBROUTINE/FUNCTION DUOC GOI TOI: BETAI
!
	  REAL X0,N,X
	  X=N/(N+X0**2)
	  TDIST=BETAI(N/2.0,0.5,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