PROGRAM BT4_12
      PARAMETER (NMAX=100,MMAX=12)
	DIMENSION X(NMAX,MMAX),Y(NMAX,MMAX)
	REAL MUY1,MUY2,XICMA, ALFA,T,TALFA
	REAL TBX(NMAX),TBY(NMAX),SX,SY
	INTEGER NX,NY,M,I,J
	CHARACTER*15 H0,H1,ANSWER1,ANSWER2
! XICMA: DO LECH CHUAN CUA X VA Y
! ALFA: XAC SUAT PHAM SAI LAM LOAI I
	DATA XICMA,ALFA /0.4,0.05/
	DATA H0,H1 /'CHAP NHAN H0','BAC BO H0'/

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

	OPEN(1,FILE='C:\MASTER\G_TR\BAITAP~1\NEW\BANG11-C.TXT')
	READ(1,*)
	READ(1,*)
	READ(1,*) NY,M
	DO I=1,NY
         READ(1,*) TMP, (Y(I,J),J=1,M)
	ENDDO

! TINH TRUNG BINH SO HOC
	DO I=1,NX
	TBX(I)=0.0
	DO J=1,M
	   TBX(I)=TBX(I)+X(I,J)
	ENDDO
	TBX(I)=TBX(I)/REAL(M)
	ENDDO
	MUY1=AMMGOC(TBX,NX,1)
	SX=AMMTAM(TBX,NX,2)
	SX=SQRT(SX)	

	DO I=1,NY
	TBY(I)=0.0
	DO J=1,M
	   TBY(I)=TBY(I)+Y(I,J)
	ENDDO
	TBY(I)=TBY(I)/REAL(M)
	ENDDO
	MUY2=AMMGOC(TBY,NY,1)
	SY=AMMTAM(TBY,NY,2)
	SY=SQRT(SY)	

! TINH T VA TALFA
	
	A=(REAL(NX-1)*SX+REAL(NY-1)*SY)/REAL(NX+NY-2)
	A=SQRT(A)
	B=(MUY1-MUY2)/SQRT(1.0/REAL(NX)+1.0/REAL(NY))
	T=B/A

	TALFA=TINV(ALFA,REAL(NX+NY-2))
	IF (ABS(T).LT.TALFA) THEN 
	   ANSWER1=H0
	ELSE
	   ANSWER1=H1
	ENDIF
	WRITE(*,*)' KET QUA KIEM NGHIEM T - H0: MUY1=MUY2'
	WRITE(*,'("   MUY1  = ",F10.3)') MUY1
	WRITE(*,'("   MUY2  = ",F10.3)') MUY2
	WRITE(*,'("     N1  = ",I10)') NX
	WRITE(*,'("     N2  = ",I10)') NY
	WRITE(*,'("      T  = ",F10.3)') T
	WRITE(*,'("   ALFA  = ",F10.3)') ALFA
	WRITE(*,'("  TALFA  = ",F10.3)') TALFA
	WRITE(*,'(" KET LUAN: ",A15)') ANSWER1

	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-6)
	  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
!	  PRINT*,'1- B=',B,'  C=',C
			GOTO 10
		 ELSE IF (AP.LT.P0) THEN
			B=(C+B)/2.0
!	  PRINT*,'2- B=',B,'  C=',C
			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 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

	  FUNCTION AMMGOC(X,N,K)
!
! HAM NAY TINH MOMEN GOC BAC K CUA CHUOI SO LIEU
! INPUT: + X MANG SO LIEU BAN DAU DO DAI N
!        + N DUNG LUONG MAU
!        + K BAC CUA MOMEN
! OUTPUT: MOMEN GOC BAC K
!
      DIMENSION X(N)
	  INTEGER N,K,I
	  REAL TMP
	  TMP=0.0
	  DO I=1,N
	     TMP=TMP+X(I)**K
	  ENDDO
	  AMMGOC=TMP/REAL(N)
	  RETURN
	  END

	  FUNCTION FAC(N)
! TINH N! (N GIAI THUA)
	  REAL TMP
	  IF (N.LT.0) THEN
	     WRITE(*,*)' INVALID NUMERIC INPUT IN FAC FUNCTION'
		 STOP
	  ELSE IF (N.EQ.0.OR.N.EQ.1) THEN
	     FAC=1.0
		 RETURN
	  ELSE
	     TMP=1.0
		 DO I=2,N
	        TMP=TMP*REAL(I)
		 ENDDO
		 FAC=TMP
		 RETURN
	  ENDIF
	  END

      FUNCTION COMBIN(N,K)
! TINH TO HOP CHAP K CUA N
! FUNCTION/SUBROUTINE DUOC GOI TOI: FAC(N)
	  IF (N.LT.0.OR.K.LT.0.OR.N.LT.K) THEN
	     WRITE(*,*)' INVALID NUMERIC INPUT IN COMBIN FUNCTION'
		 STOP
      ELSE 
         COMBIN=FAC(N)/FAC(K)/FAC(N-K)
	     RETURN
	  ENDIF
	  END

	  FUNCTION AMMTAM(X,N,K)
! HAM NAY TINH MOMEN TRUNG TAM BAC K TU CHUOI SO LIEU X
! INPUT: + X MANG SO LIEU BAN DAU DO DAI N
!        + N DUNG LUONG MAU
!        + K BAC CUA MOMEN
! OUTPUT: MOMEN TRUNG TAM BAC K
! FUNCTION/SUBROUTINE DUOC GOI TOI: AMMGOC(X,N,K), COMBIN(N,K)
	  DIMENSION X(N)
	  INTEGER N,K,I
	  REAL TMP,TB
	  TB=AMMGOC(X,N,1)
	  TMP=0
	  DO I=0,K
	     TMP=TMP+(-1)**I*COMBIN(K,I)*TB**I*AMMGOC(X,N,K-I)
	  ENDDO
	  AMMTAM=TMP
	  RETURN
	  END