*
* Kiso enshu Sample program for Student t-test
**************************************
      function betai( df, t0 )
*
*     df: degree of freedom, N-2
*     t0: t-statistics, t=r*sqrt(df/(1-r**2))
*     significance level s=1-betai
*
      real  betai, df, t0
      real  a, b, x
      real  bt, betacf, gammln

      a = df / 2.
      b = 0.5
      x = df / ( df + t0**2 )
      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( a, b, 1.-x ) / b
         return
      endif
      end
      function gammln( xx )

      real     gammln, xx 
      integer  j
      real*8   ser, stp, tmp, x, y, cof( 6 )
      save     cof, stp
      data     cof, stp 
     &     / 76.18009172947146d0, -86.50532032941677d0,
     &       24.01409824083091d0, -1.231739572450155d0,
     &       0.1208650973866179d-2, -0.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
      function betacf( a, b, x )

      integer maxit
      real    betacf, a, b, x, eps, fpmin
      parameter ( maxit = 500, 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 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 / ( ( 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
         del = d * c
         h = h * del
         if( abs( del-1. ) .lt. eps ) goto 1
      enddo

      pause 'a or b too big, or MAXIT too small in betacf'
 1    betacf = h
      return
      end