subroutine conp ( matrix, nr, nc, pt, ps, pc )

c*********************************************************************72
c
c  input arguments.
c
c    matrix = specification of the contingency table.
c      this matrix is partioned as follows:
c
c      x(11).....x(ic)  r(1)
c        .  .....  .      .
c        .  .....  .      .
c      x(r1).....x(rc)  r(r)
c       c(1)..... c(c)    n
c
c      where x(ij) are the observed cell frequencies,
c      r(i) are the row totals, c(j) are the column
c      totals, and n is the total sample size.
c      note that the original cell frequencies are
c      destroyed by this subroutine.
c
c    nr = the number of rows in matrix (r=nr-1).
c
c    nc = the number of columns in matrix (c=nc-1).
c
c  output arguments.
c
c    pt = the probability of obtaining the given table.
c
c    ps = the probability of obtaining a table as probable
c      as, or less probable than, the given table.
c
c    pc = the probability of obtaining some of the
c      tables possible within the constraints of the
c      marginal totals.  (this should be 1.0.  deviations
c      from 1.0 reflect the accuracy of the computation.)
c
c  externals
c
c    faclog(n) = function to return the floating point
c      value of log base 10 of n factorial.
c
      dimension matrix(nr,nc)
      integer r,c,temp
c
      r = nr - 1
      c = nc - 1
c
c  compute log of constant numerator.
c
      qxlog = -faclog ( matrix(nr,nc) )
      do i = 1, r
        qxlog = qxlog + faclog ( matrix(i,nc) )
      end do

      do j = 1, c
        qxlog = qxlog + faclog ( matrix(nr,j) )
      end do
c
c  compute probability of given table.
c
      rxlog = 0.0
      do i = 1, r
        do j = 1, c
          rxlog = rxlog + faclog ( matrix(i,j) )
        end do
      end do

      pt = 10.0**( qxlog - rxlog )
c
      ps = 0.0
      pc = 0.0
c
c  fill lower right (r-1) x (c-1) cells with
c  minimimum of row and column totals.
c
      do i = 2, r
        do j = 2, c
          matrix(i,j) = min0 ( matrix(i,nc), matrix(nr,j) )
        end do
      end do

      go to 300
c
c  obtain a new set of frequencies in
c  lower right (r-1) x (c-1) cells.
c
200   continue

      do i = 2, r
        do j = 2, c
          matrix(i,j) = matrix(i,j) - 1
          if ( matrix(i,j) .ge. 0 ) go to 300
          matrix(i,j) = min0 ( matrix(i,nc), matrix(nr,j) )
        end do
      end do

      return
c
c  fill remainder of observed cells
c  ...complete column 1.
c
300   continue

      do i = 2, r
        temp = matrix(i,nc)
        do j = 2, c
          temp = temp - matrix(i,j)
        end do
        if ( temp .lt. 0 ) go to 200
        matrix(i,1) = temp
      end do
c
c  ...complete row 1.
c
      do j = 1, c
        temp = matrix(nr,j)
        do i = 2, r
          temp = temp - matrix(i,j)
        end do
        if ( temp .lt. 0 ) go to 200
        matrix(1,j) = temp
      end do
c
c  compute log of the denominator.
c
      rxlog = 0.0
      do i = 1, r
        do j = 1, c
          rxlog = rxlog + faclog ( matrix(i,j) )
        end do
      end do
c
c  compute px.  add to ps if px .le. pt
c  (allow for round-off error).
c
      px = 10.0** ( qxlog - rxlog )
      pc = pc + px
      if ( ( pt / px ) .gt. 0.99999 ) ps = ps + px
      go to 200
      end
      function faclog ( n )

c*********************************************************************72
c
c  input argument.
c
c    n = an integer greater than or equal to zero.
c
c  function result.
c
c    faclog = the log to the base 10 of n factorial.
c
      dimension table(101)
      save elog,iflag,table,tpilog
      data tpilog / 0.3990899342 /
      data elog / 0.4342944819 /
      data iflag / 0 /
c
c  use stirlings approximation if n gt 100.
c
      if ( n .gt. 100 ) go to 50
c
c  lookup up answer if table was generated.
c
      if ( iflag .eq. 0 ) go to 100
10    faclog = table(n+1)
      return
c
c  here for stirlings approximation.
c
50    x = float ( n )
      faclog = ( x + 0.5 ) * alog10 ( x ) - x * elog + tpilog
     &  + elog / ( 12.0 * x ) - elog / ( 360.0 * x * x * x )
      return
c
c  here to generate log factorial table.
c
100   table(1) = 0.0
      do i = 2, 101
        x = float ( i - 1 )
        table(i) = table(i-1) + alog10 ( x )
      end do
      iflag = 1
      go to 10
      end
      function factorial ( n )

c*********************************************************************72
c
cc FACTORIAL evaluates the factorial function.
c
      real factorial

      factorial = 1.0E+00
      do i = 1, n
        factorial = factorial * i
      end do

      return
      end
      subroutine timestamp ( )

c*********************************************************************72
c
cc TIMESTAMP prints out the current YMDHMS date as a timestamp.
c
c  Modified:
c
c    16 September 2005
c
c  Author:
c
c    John Burkardt
c
      implicit none

      character ( len = 8 ) date
      character ( len = 10 ) time

      call date_and_time ( date, time )

      write ( *, '(a8,2x,a10)' ) date, time

      return
      end