      program fprobadr
      real data1(3),data2(3)
      n1=3
      n2=3
      data data1/290,311,284/
      data data2/271,304,260/
      call ftest(data1,n1,data2,n2,f,prob)

      write(*,*)f,prob

      stop
      end

      SUBROUTINE ftest(data1,n1,data2,n2,f,prob)
      INTEGER n1,n2
      REAL f,prob,data1(n1),data2(n2)
C     USES avevar,betai
C     Given the arrays data1(1:n1) and data2(1:n2),
C     this routine returns the value of f, and its signicance as prob.
C     Small values of prob indicate that the two arrays have signicantly
C     different variances.

      REAL ave1,ave2,df1,df2,var1,var2,betai
      call avevar(data1,n1,ave1,var1)
      call avevar(data2,n2,ave2,var2)
      if(var1.gt.var2)then
C      Make F the ratio of the larger variance to the smaller one.
      f=var1/var2
      df1=n1-1
      df2=n2-1
      else
      f=var2/var1
      df1=n2-1
      df2=n1-1
      endif
      prob=2.*betai(0.5*df2,0.5*df1,df2/(df2+df1*f))
      if(prob.gt.1.)prob=2.-prob
      return
      END

      SUBROUTINE avevar(data,n,ave,var)
      INTEGER n
      REAL ave,var,data(n)
C      Given array data(1:n), returns its mean as ave and its variance
C     as var.
      INTEGER j
      REAL s,ep
      ave=0.0
      do j=1,n
      ave=ave+data(j)
      enddo
      ave=ave/n
      var=0.0
      ep=0.0
      do j=1,n
      s=data(j)-ave
      ep=ep+s
      var=var+s*s
      enddo
      var=(var-ep**2/n)/(n-1)

      return
      END

      FUNCTION betai(a,b,x)
      REAL betai,a,b,x
C     USES betacf,gammln
C     Returns the incomplete beta function Ix(a,b).
      REAL bt,betacf,gammln
      if(x.lt.0..or.x.gt.1.)pause
C      'bad argument x in betai'
      if(x.lt.0..or.x.eq.1.) then
	    bt=0.
      else
C     Factors in front of the continued fraction.
	    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
C     Use continued fraction directly.
	    betai=bt*betacf(a,b,x)/a
	    return
      else
      betai=1.-bt*betacf(b,a,1.-x)/b

C     Use continued fraction after making the symmetry transformation.
      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)
C     Used by betai: Evaluates continued fraction for
C     incomplete beta function by modied
C     Lentz's method ( x 5.2).
      INTEGER m,m2
      REAL aa,c,d,del,h,qab,qam,qap
      qab=a+b
C     These q's will be used in factors that occur in the coecients
C     (6.4.6).
      qap=a+1.
      qam=a-1.
      c=1.
C     First step of Lentz's method.
      d=1.-qab*x/qap
      if(abs(d).lt.FPMIN)d=FPMIN
      d=1./d
      h=d
      do m=1,MAXIT
      m2=2*m
      aa=m*(b-m)*x/((qam+m2)*(a+m2))
      d=1.+aa*d
C     One step (the even one) of the recurrence.
      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
C     Next step of the recurrence (the odd one).
      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
C     Are we done?
      enddo
C      pause
C     'a or b too big, or MAXIT too small in betacf'
    1 betacf=h
      return
      END

      FUNCTION gammln(xx)
      REAL gammln,xx
C     Returns the value ln((xx)) for xx > 0.
      INTEGER j
      DOUBLE PRECISION ser,stp,tmp,x,y,cof(6)
C     Internal arithmetic will be done in double precision, a nicety
C     that you can omit if five-figure accuracy is good enough.
      SAVE cof,stp
      DATA cof,stp/76.18009172947146d0,-86.50532032941677d0,
     X24.01409824083091d0,-1.231739572450155d0,.1208650973866179d-2,
     X-.5395239384953d-5,2.5066282746310005d0/
      x=xx
      y=x
      tmp=x+5.5d0
      tmp=(x+0.5d0)*log(tmp)-tmp
      ser=1.000000000190015d0
      do j=1,6
	    y=y+1.d0
	    ser=ser+cof(j)/y
      enddo
      gammln=tmp+log(stp*ser/x)
      return
      END

