program main

!*****************************************************************************80
!
!! MAIN is the main program for INT_EXACTNESS.
!
!  Discussion:
!
!    This is an interactive program which allows a user to
!    specify a quadrature rule to be checked for exactness by computing
!    approximate integrals of all polynomials up to a given degree.
!
!    The user specifies:
!    * the "root" name of the R, W and X files that specify the rule;
!    * DEGREE_MAX, the maximum monomial degree to be checked.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    03 July 2008
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer   ( kind = 4 ) arg_num
  integer   ( kind = 4 ) degree
  integer   ( kind = 4 ) degree_max
  integer   ( kind = 4 ) dim
  integer   ( kind = 4 ) dim_num
  integer   ( kind = 4 ) dim_num2
  integer   ( kind = 4 ) exact_unit
  integer   ( kind = 4 ) expon(1)
  integer   ( kind = 4 ) h
  integer   ( kind = 4 ) i
  integer   ( kind = 4 ) iarg
  integer   ( kind = 4 ) iargc
  integer   ( kind = 4 ) ierror
  integer   ( kind = 4 ) ios
  integer   ( kind = 4 ) last
  logical                more
  integer   ( kind = 4 ) order
  integer   ( kind = 4 ) point_num
  integer   ( kind = 4 ) point_num2
  real      ( kind = 8 ) quad_error
  character ( len = 80 ) quad_exact_filename
  character ( len = 80 ) quad_filename
  character ( len = 80 ) quad_r_filename
  character ( len = 80 ) quad_w_filename
  character ( len = 80 ) quad_x_filename
  real      ( kind = 8 ), allocatable, dimension ( :, : ) :: r
  character ( len = 80 ) string
  integer   ( kind = 4 ) t
  real      ( kind = 8 ) volume
  real      ( kind = 8 ), allocatable, dimension ( : ) :: w
  real      ( kind = 8 ), allocatable, dimension ( :, : ) :: x

  call timestamp ( )
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'INT_EXACTNESS'
  write ( *, '(a)' ) '  FORTRAN90 version'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Investigate the polynomial exactness of a quadrature'
  write ( *, '(a)' ) '  rule by integrating all monomials of a given degree'
  write ( *, '(a)' ) '  over the [0,1] interval.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  If necessary, the rule is adjusted to the [0,1] interval.'
!
!  Get the number of command line arguments.
!
  arg_num = iargc ( )
!
!  Get the quadrature file root name:
!
  if ( 1 <= arg_num ) then

    iarg = 1
    call getarg ( iarg, quad_filename )
!
!  A commandline argument can't be blank.  To let a user
!  specify the shortcut names, let "D" indicate them.
!
    if ( quad_filename == 'D' ) then
      quad_filename = ''
    end if

  else

    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'INT_EXACTNESS:'
    write ( *, '(a)' ) '  Enter the "root" name of the quadrature files.'

    read ( *, '(a)' ) quad_filename

  end if
!
!  Create the names of:
!    the quadrature X file;
!    the quadrature W file;
!    the quadrature R file;
!    the output "exactness" file.
!
  if ( 0 < len_trim ( quad_filename ) ) then
    quad_x_filename = trim ( quad_filename ) // '_x.txt'
    quad_w_filename = trim ( quad_filename ) // '_w.txt'
    quad_r_filename = trim ( quad_filename ) // '_r.txt'
    quad_exact_filename = trim ( quad_filename ) // '_exact.txt'
  else
    quad_x_filename = 'x.txt'
    quad_w_filename = 'w.txt'
    quad_r_filename = 'r.txt'
    quad_exact_filename = 'exact.txt'
  end if
!
!  The second command line argument is the maximum degree.
!
  if ( 2 <= arg_num ) then

    iarg = 2
    call getarg ( iarg, string )
    call s_to_i4 ( string, degree_max, ierror, last )

  else

    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'INT_EXACTNESS:'
    write ( *, '(a)' ) '  Please enter DEGREE_MAX, the maximum monomial degree to check.'

    read ( *, * ) degree_max

  end if
!
!  Open the output file.
!
  call get_unit ( exact_unit )

  open ( unit = exact_unit, file = quad_exact_filename, status = 'replace', &
    iostat = ios )

  if ( ios /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'INT_EXACTNESS - Fatal error!'
    write ( *, '(a)' ) '  Could not open the output file "' &
      // trim ( quad_exact_filename ) // '".'
    stop
  end if
!
!  Summarize the input.
!
  call timestring ( string )
  write ( exact_unit, '(a)' ) string
  write ( exact_unit, '(a)' ) ' '
  write ( exact_unit, '(a)' ) 'INT_EXACTNESS'
  write ( exact_unit, '(a)' ) '  FORTRAN90 version'
  write ( exact_unit, '(a)' ) ' '
  write ( exact_unit, '(a)' ) '  Investigate the polynomial exactness of a quadrature'
  write ( exact_unit, '(a)' ) '  rule by integrating all monomials up to a given degree'
  write ( exact_unit, '(a)' ) '  over the [0,1] interval.'
  write ( exact_unit, '(a)' ) ' '
  write ( exact_unit, '(a)' ) '  If necessary, the rule is adjusted to the [0,1] interval.'
  write ( exact_unit, '(a)' ) ' '
  write ( exact_unit, '(a)' ) 'INT_EXACTNESS: User input:'
  write ( exact_unit, '(a)' ) '  Quadrature rule X file = "' &
    // trim ( quad_x_filename ) // '".'
  write ( exact_unit, '(a)' ) '  Quadrature rule W file = "' &
    // trim ( quad_w_filename ) // '".'
  write ( exact_unit, '(a)' ) '  Quadrature rule R file = "' &
    // trim ( quad_r_filename ) // '".'
  write ( exact_unit, '(a,i8)' ) '  Maximum degree to check = ', &
    degree_max
!
!  Read the X file.
!
  call dtable_header_read ( quad_x_filename, dim_num, point_num )

  if ( .false. ) then
    if ( dim_num /= 1 ) then
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'INT_EXACTNESS - Fatal error!'
      write ( *, '(a)' ) '  The spatial dimension should be 1.'
      write ( *, '(a,i8)' ) '  The implicit input dimension was DIM_NUM = ', dim_num
      stop
    end if
  end if

  order = point_num

  write ( exact_unit, '(a)' ) ' '
  write ( exact_unit, '(a,i8)' ) '  Spatial dimension = ', dim_num
  write ( exact_unit, '(a,i8)' ) '  Number of points  = ', point_num

  allocate ( x(dim_num,point_num) )

  call dtable_data_read ( quad_x_filename, dim_num, point_num, x )
!
!  Read the W file.
!
  call dtable_header_read ( quad_w_filename, dim_num2, point_num2 )

  if ( dim_num2 /= 1 ) then
    write ( exact_unit, '(a)' ) ' '
    write ( exact_unit, '(a)' ) 'INT_EXACTNESS - Fatal error!'
    write ( exact_unit, '(a)' ) '  The quadrature weight file should have exactly'
    write ( exact_unit, '(a)' ) '  one value on each line.'
    stop
  end if

  if ( point_num2 /= point_num ) then
    write ( exact_unit, '(a)' ) ' '
    write ( exact_unit, '(a)' ) 'INT_EXACTNESS - Fatal error!'
    write ( exact_unit, '(a)' ) '  The quadrature weight file should have exactly'
    write ( exact_unit, '(a)' ) '  the same number of lines as the abscissa file.'
    stop
  end if

  allocate ( w(point_num) )

  call dtable_data_read ( quad_w_filename, 1, point_num, w )
!
!  Read the R file.
!
  call dtable_header_read ( quad_r_filename, dim_num2, point_num2 )

  if ( dim_num2 /= dim_num ) then
    write ( exact_unit, '(a)' ) ' '
    write ( exact_unit, '(a)' ) 'INT_EXACTNESS - Fatal error!'
    write ( exact_unit, '(a)' ) '  The quadrature region file should have the'
    write ( exact_unit, '(a)' ) '  same number of values on each line as the'
    write ( exact_unit, '(a)' ) '  abscissa file does.'
    stop
  end if

  if ( point_num2 /= 2 ) then
    write ( exact_unit, '(a)' ) ' '
    write ( exact_unit, '(a)' ) 'INT_EXACTNESS - Fatal error!'
    write ( exact_unit, '(a)' ) '  The quadrature region file should have two lines.'
    stop
  end if

  allocate ( r(dim_num,2) )

  call dtable_data_read ( quad_r_filename, dim_num, 2, r )
!
!  Print the input quadrature rule.
!
  write ( exact_unit, '(a)'       ) ' '
  write ( exact_unit, '(a)'       ) '  The quadrature rule to be tested:'
  write ( exact_unit, '(a,i8)'    ) '  ORDER = ', order
  write ( exact_unit, '(a)'       ) ' '
  write ( exact_unit, '(a)'       ) '  Standard rule:'
  write ( exact_unit, '(a)'       ) '    Integral ( R(1) <= x <= R(2) ) f(x) dx'
  write ( exact_unit, '(a)'       ) '    is to be approximated by'
  write ( exact_unit, '(a)'       ) '    sum ( 1 <= I <= ORDER ) w(i) * f(x(i)).'
  write ( exact_unit, '(a)'       ) ' '
  write ( exact_unit, '(a)'       ) '  Weights W:'
  write ( exact_unit, '(a)'       ) ' '
  do i = 1, order
    write ( exact_unit, '(g24.16)' ) w(i)
  end do
  write ( exact_unit, '(a)'       ) ' '
  write ( exact_unit, '(a)'       ) '  Abscissas X:'
  write ( exact_unit, '(a)'       ) ' '
  do i = 1, order
    write ( exact_unit, '(5g24.16)' ) x(1:dim_num,i)
  end do
  write ( exact_unit, '(a)'       ) ' '
  write ( exact_unit, '(a)'       ) '  Region R:'
  write ( exact_unit, '(a)'       ) ' '

  do i = 1, 2
    write ( exact_unit, '(5g24.16)' ) r(1:dim_num,i)
  end do
!
!  Rescale the weights, and translate the abscissas.
!
  volume = abs ( product ( r(1:dim_num,2) - r(1:dim_num,1) ) )
  w(1:order) = w(1:order) / volume

  do dim = 1, dim_num
    x(dim,1:order) = ( x(dim,1:order) - r(dim,1) ) &
                   / ( r(dim,2) - r(dim,1) )
  end do
!
!  Explore the monomials.
!
  write ( exact_unit, '(a)' ) ' '
  write ( exact_unit, '(a)' ) '          Error          Degree  Exponents'
  write ( exact_unit, '(a)' ) ' '

  do degree = 0, degree_max

    more = .false.
    h = 0
    t = 0

    do

      call comp_next ( degree, dim_num, expon, more, h, t )

      if ( .not. more ) then
        exit
      end if

      call monomial_quadrature ( dim_num, expon, order, w, x, &
        quad_error )

      write ( exact_unit, '(2x,f24.16,3x,i2,4x,10i2)' ) &
        quad_error, degree, expon(1:dim_num)

    end do

  end do

  write ( exact_unit, '(a)' ) ' '
  write ( exact_unit, '(a)' ) 'INT_EXACTNESS:'
  write ( exact_unit, '(a)' ) '  Normal end of execution.'
  write ( exact_unit, '(a)' ) ' '
  call timestring ( string )
  write ( exact_unit, '(a)' ) string

  close ( unit = exact_unit )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'INT_EXACTNESS:'
  write ( *, '(a)' ) '  The exactness results were written to "' &
    // trim ( quad_exact_filename ) // '".'

  deallocate ( r )
  deallocate ( w )
  deallocate ( x )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'INT_EXACTNESS:'
  write ( *, '(a)' ) '  Normal end of execution.'
  write ( *, '(a)' ) ' '
  call timestamp ( )

  stop
end
subroutine ch_cap ( c )

!*****************************************************************************80
!
!! CH_CAP capitalizes a single character.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    05 February 2008
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input/output, character C, the character to capitalize.
!
  implicit none

  character c
  integer ( kind = 4 ) itemp

  itemp = ichar ( c )

  if ( 97 <= itemp .and. itemp <= 122 ) then
    c = char ( itemp - 32 )
  end if

  return
end
function ch_eqi ( c1, c2 )

!*****************************************************************************80
!
!! CH_EQI is a case insensitive comparison of two characters for equality.
!
!  Example:
!
!    CH_EQI ( 'A', 'a' ) is .TRUE.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    05 February 2008
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, character C1, C2, the characters to compare.
!
!    Output, logical CH_EQI, the result of the comparison.
!
  implicit none

  logical ch_eqi
  character c1
  character c1_cap
  character c2
  character c2_cap

  c1_cap = c1
  c2_cap = c2

  call ch_cap ( c1_cap )
  call ch_cap ( c2_cap )

  if ( c1_cap == c2_cap ) then
    ch_eqi = .true.
  else
    ch_eqi = .false.
  end if

  return
end
subroutine ch_to_digit ( c, digit )

!*****************************************************************************80
!
!! CH_TO_DIGIT returns the integer value of a base 10 digit.
!
!  Example:
!
!     C   DIGIT
!    ---  -----
!    '0'    0
!    '1'    1
!    ...  ...
!    '9'    9
!    ' '    0
!    'X'   -1
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    05 February 2008
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, character C, the decimal digit, '0' through '9' or blank
!    are legal.
!
!    Output, integer ( kind = 4 ) DIGIT, the corresponding integer value.  
!    If C was 'illegal', then DIGIT is -1.
!
  implicit none

  character c
  integer ( kind = 4 ) digit

  if ( lge ( c, '0' ) .and. lle ( c, '9' ) ) then

    digit = ichar ( c ) - 48

  else if ( c == ' ' ) then

    digit = 0

  else

    digit = -1

  end if

  return
end
subroutine asm_enum ( n, asm_num )

!*****************************************************************************80
!
!! ASM_ENUM returns the number of alternating sign matrices of a given order.
!
!  Discussion:
!
!    N     ASM_NUM
!
!    0       1
!    1       1
!    2       2
!    3       7
!    4      42
!    5     429
!    6    7436
!    7  218348
!
!    A direct formula is
!
!      ASM_NUM ( N ) = product ( 0 <= I <= N-1 ) ( 3 * I + 1 )! / ( N + I )!
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    01 January 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) N, the order of the matrices.
!
!    Output, integer ( kind = 4 ) ASM_NUM, the number of alternating sign 
!    matrices of order N.
!
  implicit none

  integer ( kind = 4 ) n

  integer ( kind = 4 ) a(n+1)
  integer ( kind = 4 ) asm_num
  integer ( kind = 4 ) b(n+1)
  integer ( kind = 4 ) c(n+1)
  integer ( kind = 4 ) i
  integer ( kind = 4 ) nn

  asm_num = 0

  nn = n + 1

  if ( n + 1 <= 0 ) then
    return
  end if
!
!  Row 1
!
  a(1) = 1

  if ( n + 1 == 1 ) then
    asm_num = a(1)
    return
  end if
!
!  Row 2
!
  a(1) = 1

  if ( n + 1 == 2 ) then
    asm_num = a(1)
    return
  end if

  b(1) = 2
  c(1) = 2
  a(2) = 1
!
!  Row 3 and on.
!
  do nn = 3, n

    b(nn-1) = nn
    do i = nn-2, 2, -1
      b(i) = b(i) + b(i-1)
    end do
    b(1) = 2

    c(nn-1) = 2
    do i = nn-2, 2, -1
      c(i) = c(i) + c(i-1)
    end do
    c(1) = nn

    a(1) = sum ( a(1:nn-1) )
    do i = 2, nn
      a(i) = a(i-1) * c(i-1) / b(i-1)
    end do

  end do

  asm_num = sum ( a(1:n) )

  return
end
subroutine asm_triangle ( n, a )

!*****************************************************************************80
!
!! ASM_TRIANGLE returns a row of the alternating sign matrix triangle.
!
!  Discussion:
!
!    The first seven rows of the triangle are as follows:
!
!          1      2      3      4      5      6     7
!
!    0     1
!    1     1      1
!    2     2      3      2
!    3     7     14     14      7
!    4    42    105    135    105     42
!    5   429   1287   2002   2002   1287    429
!    6  7436  26026  47320  56784  47320  26026  7436
!
!    For a given N, the value of A(J) represents entry A(I,J) of
!    the triangular matrix, and gives the number of alternating sign matrices
!    of order N in which the (unique) 1 in row 1 occurs in column J.
!
!    Thus, of alternating sign matrices of order 3, there are
!    2 with a leading 1 in column 1:
!
!      1 0 0  1 0 0
!      0 1 0  0 0 1
!      0 0 1  0 1 0
!
!    3 with a leading 1 in column 2, and
!
!      0 1 0  0 1 0  0 1 0
!      1 0 0  0 0 1  1-1 1
!      0 0 1  1 0 0  0 1 0
!
!    2 with a leading 1 in column 3:
!
!      0 0 1  0 0 1
!      1 0 0  0 1 0
!      0 1 0  1 0 0
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    01 January 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) N, the desired row.
!
!    Output, integer ( kind = 4 ) A(N+1), the entries of the row.
!
  implicit none

  integer ( kind = 4 ) n

  integer ( kind = 4 ) a(n+1)
  integer ( kind = 4 ) b(n+1)
  integer ( kind = 4 ) c(n+1)
  integer ( kind = 4 ) i
  integer ( kind = 4 ) nn

  if ( n+1 <= 0 ) then
    return
  end if
!
!  Row 1
!
  a(1) = 1

  if ( n + 1 == 1 ) then
    return
  end if
!
!  Row 2
!
  nn = 2
  b(1) = 2
  c(1) = nn

  a(1) = sum ( a(1:nn-1) )
  do i = 2, nn
    a(i) = a(i-1) * c(i-1) / b(i-1)
  end do

  if ( n + 1 == 2 ) then
    return
  end if
!
!  Row 3 and on.
!
  do nn = 3, n + 1

    b(nn-1) = nn
    do i = nn-2, 2, -1
      b(i) = b(i) + b(i-1)
    end do
    b(1) = 2

    c(nn-1) = 2
    do i = nn-2, 2, -1
      c(i) = c(i) + c(i-1)
    end do
    c(1) = nn

    a(1) = sum ( a(1:nn-1) )
    do i = 2, nn
      a(i) = a(i-1) * c(i-1) / b(i-1)
    end do

  end do

  return
end
subroutine bell ( n, b )

!*****************************************************************************80
!
!! BELL returns the Bell numbers from 0 to N.
!
!  Discussion:
!
!    The Bell number B(N) is defined as the number of partitions (of
!    any size) of a set of N distinguishable objects.
!
!    A partition of a set is a division of the objects of the set into
!    subsets.
!
!    The Bell number B(N) is the number of restricted growth functions
!    on N.
!
!    Note that the Stirling numbers of the second kind, S^m_n, count the
!    number of partitions of N objects into M classes, and so it is
!    true that
!
!      B(N) = S^1_N + S^2_N + ... + S^N_N.
!
!  Examples:
!
!    There are 15 partitions of a set of 4 objects:
!
!      (1234), (123)(4), (124)(3), (12)(34), (12)(3)(4),
!      (134)(2), (13)(24), (13)(2)(4), (14)(23), (1)(234),
!      (1)(23)(4), (14)(2)(3), (1)(24)(3), (1)(2)(34), (1)(2)(3)(4)
!
!    and so B(4) = 15.
!
!  First values:
!
!     N         B(N)
!     0           1
!     1           1
!     2           2
!     3           5
!     4          15
!     5          52
!     6         203
!     7         877
!     8        4140
!     9       21147
!    10      115975
!
!  Recursion:
!
!    B(I) = sum ( 1 <= J <= I ) Binomial ( I-1, J-1 ) * B(I-J)
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    02 June 2007
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) N, the number of Bell numbers desired.
!
!    Output, integer ( kind = 4 ) B(0:N), the Bell numbers from 0 to N.
!
  implicit none

  integer ( kind = 4 ) n

  integer ( kind = 4 ) b(0:n)
  integer ( kind = 4 ) combo
  integer ( kind = 4 ) i
  integer ( kind = 4 ) i4_choose
  integer ( kind = 4 ) j

  b(0) = 1

  do i = 1, n
    b(i) = 0
    do j = 1, i
      combo = i4_choose ( i-1, j-1 )
      b(i) = b(i) + combo * b(i-j)
    end do
  end do

  return
end
subroutine bell_values ( n_data, n, c )

!*****************************************************************************80
!
!! BELL_VALUES returns some values of the Bell numbers for testing.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    06 February 2003
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Milton Abramowitz, Irene Stegun,
!    Handbook of Mathematical Functions,
!    National Bureau of Standards, 1964,
!    ISBN: 0-486-61272-4,
!    LC: QA47.A34.
!
!  Parameters:
!
!    Input/output, integer ( kind = 4 ) N_DATA.
!    On input, if N_DATA is 0, the first test data is returned, and N_DATA
!    is set to 1.  On each subsequent call, the input value of N_DATA is
!    incremented and that test data item is returned, if available.  When
!    there is no more test data, N_DATA is set to 0.
!
!    Output, integer ( kind = 4 ) N, the order of the Bell number.
!
!    Output, integer ( kind = 4 ) C, the value of the Bell number.
!
  implicit none

  integer ( kind = 4 ), parameter :: nmax = 11

  integer ( kind = 4 ) c
  integer ( kind = 4 ), save, dimension ( nmax ) :: c_vec = (/ &
    1, 1, 2, 5, 15, 52, 203, 877, 4140, 21147, 115975 /)
  integer ( kind = 4 ) n
  integer ( kind = 4 ) n_data
  integer ( kind = 4 ), save, dimension ( nmax ) :: n_vec = (/ &
     0,  1,  2,  3,  4, 5,  6,  7,  8,  9,  10 /)

  if ( n_data < 0 ) then
    n_data = 0
  end if

  n_data = n_data + 1

  if ( nmax < n_data ) then
    n_data = 0
    n = 0
    c = 0
  else
    n = n_vec(n_data)
    c = c_vec(n_data)
  end if

  return
end
subroutine bvec_add ( n, bvec1, bvec2, bvec3 )

!*****************************************************************************80
!
!! BVEC_ADD adds two (signed) binary vectors.
!
!  Discussion:
!
!    A BVEC is an integer vector of binary digits, intended to
!    represent an integer.  BVEC(1) is the units digit, BVEC(N-1)
!    is the coefficient of 2**(N-2), and BVEC(N) contains sign
!    information.  It is 0 if the number is positive, and 1 if
!    the number is negative.
!
!  Example:
!
!    N = 5
!
!      BVEC1       +   BVEC2       =   BVEC3
!
!    ( 1 0 0 0 0 ) + ( 1 1 0 0 0 ) = ( 0 0 1 0 0 )
!
!      +1          +   +3         =   + 4
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    30 November 2006
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) N, the length of the vectors.
!
!    Input, integer ( kind = 4 ) BVEC1(N), BVEC2(N), the vectors to be added.
!
!    Output, integer ( kind = 4 ) BVEC3(N), the sum of the two input vectors.
!
  implicit none

  integer ( kind = 4 ) n

  integer ( kind = 4 ), parameter :: base = 2
  integer ( kind = 4 ) bvec1(n)
  integer ( kind = 4 ) bvec2(n)
  integer ( kind = 4 ) bvec3(n)
  integer ( kind = 4 ) i
  logical overflow

  overflow = .false.

  bvec3(1:n) = bvec1(1:n) + bvec2(1:n)

  do i = 1, n

    do while ( base <= bvec3(i) )

      bvec3(i) = bvec3(i) - base

      if ( i < n ) then
        bvec3(i+1) = bvec3(i+1) + 1
      else
        overflow = .true.
      end if

    end do

  end do

  return
end
subroutine bvec_and ( n, bvec1, bvec2, bvec3 )

!*****************************************************************************80
!
!! BVEC_AND computes the AND of two binary vectors.
!
!  Discussion:
!
!    A BVEC is an integer vector of binary digits, intended to
!    represent an integer.  BVEC(1) is the units digit, BVEC(N-1)
!    is the coefficient of 2**(N-2), and BVEC(N) contains sign
!    information.  It is 0 if the number is positive, and 1 if
!    the number is negative.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    30 November 2006
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) N, the length of the vectors.
!
!    Input, integer ( kind = 4 ) BVEC1(N), BVEC2(N), the binary vectors.
!
!    Input, integer ( kind = 4 ) BVEC3(N), the AND of the two vectors.
!
  implicit none

  integer ( kind = 4 ) n

  integer ( kind = 4 ) bvec1(n)
  integer ( kind = 4 ) bvec2(n)
  integer ( kind = 4 ) bvec3(n)

  bvec3(1:n) = min ( bvec1(1:n), bvec2(1:n) )

  return
end
subroutine bvec_check ( n, bvec, ierror )

!*****************************************************************************80
!
!! BVEC_CHECK checks a binary vector.
!
!  Discussion:
!
!    A BVEC is an integer vector of binary digits, intended to
!    represent an integer.  BVEC(1) is the units digit, BVEC(N-1)
!    is the coefficient of 2**(N-2), and BVEC(N) contains sign
!    information.  It is 0 if the number is positive, and 1 if
!    the number is negative.
!
!    The only check made is that the entries are all 0 or 1.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    02 May 2003
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) N, the length of the vectors.
!
!    Input, integer ( kind = 4 ) BVEC(N), the vector to be checked.
!
!    Output, integer ( kind = 4 ) IERROR, is nonzero if an error occurred.
!
  implicit none

  integer ( kind = 4 ) n

  integer ( kind = 4 ), parameter :: base = 2
  integer ( kind = 4 ) bvec(n)
  integer ( kind = 4 ) i
  integer ( kind = 4 ) ierror

  ierror = 0

  do i = 1, n
    if ( bvec(i) < 0 .or. base <= bvec(i) ) then
      ierror = i
      return
    end if
  end do

  return
end
subroutine bvec_complement2 ( n, bvec1, bvec2 )

!*****************************************************************************80
!
!! BVEC_COMPLEMENT2 computes the two's complement of a binary vector.
!
!  Discussion:
!
!    A BVEC is an integer vector of binary digits, intended to
!    represent an integer.  BVEC(1) is the units digit, BVEC(N-1)
!    is the coefficient of 2**(N-2), and BVEC(N) contains sign
!    information.  It is 0 if the number is positive, and 1 if
!    the number is negative.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    30 November 2006
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) N, the length of the vectors.
!
!    Input, integer ( kind = 4 ) BVEC1(N), the vector to be complemented.
!
!    Output, integer ( kind = 4 ) BVEC2(N), the two's complemented vector.
!
  implicit none

  integer ( kind = 4 ) n

  integer ( kind = 4 ), parameter :: base = 2
  integer ( kind = 4 ) bvec1(n)
  integer ( kind = 4 ) bvec2(n)
  integer ( kind = 4 ) bvec3(n)
  integer ( kind = 4 ) bvec4(n)

  bvec3(1:n) = ( base - 1 ) - bvec1(1:n)

  bvec4(1) = 1
  bvec4(2:n) = 0

  call bvec_add ( n, bvec3, bvec4, bvec2 )

  return
end
subroutine bvec_mul ( n, bvec1, bvec2, bvec3 )

!*****************************************************************************80
!
!! BVEC_MUL computes the product of two binary vectors.
!
!  Discussion:
!
!    A BVEC is an integer vector of binary digits, intended to
!    represent an integer.  BVEC(1) is the units digit, BVEC(N-1)
!    is the coefficient of 2**(N-2), and BVEC(N) contains sign
!    information.  It is 0 if the number is positive, and 1 if
!    the number is negative.
!
!    Since the user may want to make calls like
!
!      call bvec_mul ( n, bvec1, bvec1, bvec3 )
!    or even
!      call bvec_mul ( n, bvec1, bvec1, bvec1 )
!
!    we need to copy the arguments, work on them, and then copy out the result.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    02 December 2006
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) N, the length of the vectors.
!
!    Input, integer ( kind = 4 ) BVEC1(N), BVEC2(N), the vectors to 
!    be multiplied.
!
!    Output, integer ( kind = 4 ) BVEC3(N), the product of the two 
!    input vectors.
!
  implicit none

  integer ( kind = 4 ) n

  integer ( kind = 4 ), parameter :: base = 2
  integer ( kind = 4 ) carry
  integer ( kind = 4 ) bvec1(n)
  integer ( kind = 4 ) bvec2(n)
  integer ( kind = 4 ) bvec3(n)
  integer ( kind = 4 ) bveca(n)
  integer ( kind = 4 ) bvecb(n)
  integer ( kind = 4 ) bvecc(n)
  integer ( kind = 4 ) i
  integer ( kind = 4 ) product_sign
!
!  Copy the input.
!
  bveca(1:n) = bvec1(1:n)
  bvecb(1:n) = bvec2(1:n)
!
!  Record the sign of the product.
!  Make the factors positive.
!
  product_sign = 1

  if ( bveca(n) /= 0 ) then
    product_sign = - product_sign
    call bvec_complement2 ( n, bveca, bveca )
  end if

  if ( bvecb(n) /= 0 ) then
    product_sign = - product_sign
    call bvec_complement2 ( n, bvecb, bvecb )
  end if

  bvecc(1:n) = 0
!
!  Multiply.
!
  do i = 1, n-1
    bvecc(i:n-1) = bvecc(i:n-1) + bveca(i) * bvecb(1:n-i)
  end do
!
!  Take care of carries.
!
  do i = 1, n-1

    carry = bvecc(i) / base
    bvecc(i) = bvecc(i) - carry * base
!
!  Unlike the case of BVEC_ADD, we do NOT allow carries into
!  the N-th position when multiplying.
!
    if ( i < n - 1 ) then
      bvecc(i+1) = bvecc(i+1) + carry
    end if

  end do
!
!  Take care of the sign of the product.
!
  if ( product_sign < 0 ) then
    call bvec_complement2 ( n, bvecc, bvecc )
  end if
!
!  Copy the output.
!
  bvec3(1:n) = bvecc(1:n)

  return
end
subroutine bvec_next ( n, bvec )

!*****************************************************************************80
!
!! BVEC_NEXT generates the next binary vector.
!
!  Discussion:
!
!    The vectors have the order
!    (0,0,...,0),
!    (0,0,...,1), 
!    ...
!    (1,1,...,1)
!
!    and the "next" vector after (1,1,...,1) is (0,0,...,0).  That is,
!    we allow wrap around.
!
!  Example:
!
!    N = 3
!
!    Input      Output
!    -----      ------
!    0 0 0  =>  0 0 1
!    0 0 1  =>  0 1 0
!    0 1 0  =>  0 1 1
!    0 1 1  =>  1 0 0
!    1 0 0  =>  1 0 1
!    1 0 1  =>  1 1 0
!    1 1 0  =>  1 1 1
!    1 1 1  =>  0 0 0
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    02 May 2008
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer N, the dimension of the vectors.
!
!    Input/output, integer BVEC(N), on output, the successor to the
!    input vector.  
!
  implicit none

  integer n

  integer bvec(n)
  integer i
  
  do i = n, 1, -1

    if ( bvec(i) == 0 ) then
      bvec(i) = 1
      return
    end if

    bvec(i) = 0

  end do

  return
end
subroutine bvec_not ( n, bvec1, bvec2 )

!*****************************************************************************80
!
!! BVEC_NOT "negates" or takes the 1's complement of a binary vector.
!
!  Discussion:
!
!    A BVEC is an integer vector of binary digits, intended to
!    represent an integer.  BVEC(1) is the units digit, BVEC(N-1)
!    is the coefficient of 2**(N-2), and BVEC(N) contains sign
!    information.  It is 0 if the number is positive, and 1 if
!    the number is negative.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    02 May 2003
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) N, the length of the vectors.
!
!    Input, integer ( kind = 4 ) BVEC1(N), the vector to be negated.
!
!    Output, integer ( kind = 4 ) BVEC2(N), the negated vector.
!
  implicit none

  integer ( kind = 4 ) n

  integer ( kind = 4 ), parameter :: base = 2
  integer ( kind = 4 ) bvec1(n)
  integer ( kind = 4 ) bvec2(n)

  bvec2(1:n) = ( base - 1 ) - bvec1(1:n)

  return
end
subroutine bvec_or ( n, bvec1, bvec2, bvec3 )

!*****************************************************************************80
!
!! BVEC_OR computes the inclusive OR of two binary vectors.
!
!  Discussion:
!
!    A BVEC is an integer vector of binary digits, intended to
!    represent an integer.  BVEC(1) is the units digit, BVEC(N-1)
!    is the coefficient of 2**(N-2), and BVEC(N) contains sign
!    information.  It is 0 if the number is positive, and 1 if
!    the number is negative.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    30 November 2006
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) N, the length of the vectors.
!
!    Input, integer ( kind = 4 ) BVEC1(N), BVEC2(N), the binary vectors.
!
!    Input, integer ( kind = 4 ) BVEC3(N), the inclusive OR of the two vectors.
!
  implicit none

  integer ( kind = 4 ) n

  integer ( kind = 4 ) bvec1(n)
  integer ( kind = 4 ) bvec2(n)
  integer ( kind = 4 ) bvec3(n)

  bvec3(1:n) = max ( bvec1(1:n), bvec2(1:n) )

  return
end
subroutine bvec_print ( n, bvec, title )

!*****************************************************************************80
!
!! BVEC_PRINT prints a binary integer vector, with an optional title.
!
!  Discussion:
!
!    A BVEC is an integer vector of binary digits, intended to
!    represent an integer.  BVEC(1) is the units digit, BVEC(N-1)
!    is the coefficient of 2**(N-2), and BVEC(N) contains sign
!    information.  It is 0 if the number is positive, and 1 if
!    the number is negative.
!
!    The vector is printed "backwards", that is, the first entry
!    printed is BVEC(N).
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    30 November 2006
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) N, the number of components of the vector.
!
!    Input, integer ( kind = 4 ) BVEC(N), the vector to be printed.
!
!    Input, character ( len = * ) TITLE, a title to be printed first.
!    TITLE may be blank.
!
  implicit none

  integer ( kind = 4 ) n

  integer ( kind = 4 ) bvec(n)
  integer ( kind = 4 ) ihi
  integer ( kind = 4 ) ilo
  character ( len = * ) title

  if ( title /= ' ' ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) trim ( title )
    write ( *, '(a)' ) ' '
  end if

  do ihi = n, 1, -80
    ilo = max ( ihi - 80 + 1, 1 )
    write ( *, '(2x,80i1)' ) bvec(ihi:ilo:-1)
  end do

  return
end
subroutine bvec_reverse ( n, bvec1, bvec2 )

!*****************************************************************************80
!
!! BVEC_REVERSE reverses a binary vector.
!
!  Discussion:
!
!    A BVEC is an integer vector of binary digits, intended to
!    represent an integer.  BVEC(1) is the units digit, BVEC(N-1)
!    is the coefficient of 2**(N-2), and BVEC(N) contains sign
!    information.  It is 0 if the number is positive, and 1 if
!    the number is negative.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    03 May 2003
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) N, the length of the vectors.
!
!    Input, integer ( kind = 4 ) BVEC1(N), the vector to be reversed.
!
!    Output, integer ( kind = 4 ) BVEC2(N), the reversed vector.
!
  implicit none

  integer ( kind = 4 ) n

  integer ( kind = 4 ) bvec1(n)
  integer ( kind = 4 ) bvec2(n)

  bvec2(1:n) = bvec1(n:1:-1)

  return
end
subroutine bvec_sub ( n, bvec1, bvec2, bvec3 )

!*****************************************************************************80
!
!! BVEC_SUB subtracts two binary vectors.
!
!  Discussion:
!
!    A BVEC is an integer vector of binary digits, intended to
!    represent an integer.  BVEC(1) is the units digit, BVEC(N-1)
!    is the coefficient of 2**(N-2), and BVEC(N) contains sign
!    information.  It is 0 if the number is positive, and 1 if
!    the number is negative.
!
!  Example:
!
!    N = 4
!
!    BVEC1    dec  BVEC2    dec  BVEC3    dec
!    -------  ---  -------  ---  -------  ---
!    0 0 1 0   4   1 0 0 0   1   1 1 0 0   3
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    08 May 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) N, the length of the vectors.
!
!    Input, integer ( kind = 4 ) BVEC1(N), BVEC2(N), the vectors to 
!    be subtracted.
!
!    Output, integer ( kind = 4 ) BVEC3(N), the value of BVEC1 - BVEC2.
!
  implicit none

  integer ( kind = 4 ) n

  integer ( kind = 4 ) bvec1(n)
  integer ( kind = 4 ) bvec2(n)
  integer ( kind = 4 ) bvec3(n)
  integer ( kind = 4 ) bvec4(n)

  call bvec_complement2 ( n, bvec2, bvec4 )

  call bvec_add ( n, bvec1, bvec4, bvec3 )

  return
end
subroutine bvec_to_i4 ( n, bvec, i4 )

!*****************************************************************************80
!
!! BVEC_TO_I4 makes an integer from a (signed) binary vector.
!
!  Discussion:
!
!    A BVEC is an integer vector of binary digits, intended to
!    represent an integer.  BVEC(1) is the units digit, BVEC(N-1)
!    is the coefficient of 2**(N-2), and BVEC(N) contains sign
!    information.  It is 0 if the number is positive, and 1 if
!    the number is negative.
!
!  Example:
!
!         BVEC   binary  I
!    ----------  -----  --
!    1  2  3  4
!    ----------
!    1, 0, 0, 0       1  1
!    0, 1, 0, 0      10  2
!    0, 0, 1, 1    -100 -4
!    0, 0, 1, 0     100  4
!    1, 0, 0, 1    -111 -9
!    1, 1, 1, 1      -0  0
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    30 November 2006
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) N, the dimension of the vector.
!
!    Input, integer ( kind = 4 ) BVEC(N), the binary representation.
!
!    Output, integer ( kind = 4 ) I4, the integer.
!
  implicit none

  integer ( kind = 4 ) n

  integer ( kind = 4 ), parameter :: base = 2
  integer ( kind = 4 ) bvec(n)
  integer ( kind = 4 ) bvec2(n)
  integer ( kind = 4 ) i
  integer ( kind = 4 ) i_sign
  integer ( kind = 4 ) i4
  integer ( kind = 4 ), parameter :: i4_1 = 1

  bvec2(1:n) = bvec(1:n)

  if ( bvec2(n) == base - 1 ) then
    i_sign = -1
    call bvec_complement2 ( n-i4_1, bvec2, bvec2 )
  else 
    i_sign = 1
  end if

  i4 = 0
  do i = n-1, 1, -1
    i4 = base * i4 + bvec2(i)
  end do

  i4 = i_sign * i4

  return
end
subroutine bvec_xor ( n, bvec1, bvec2, bvec3 )

!*****************************************************************************80
!
!! BVEC_XOR computes the exclusive OR of two binary vectors.
!
!  Discussion:
!
!    A BVEC is an integer vector of binary digits, intended to
!    represent an integer.  BVEC(1) is the units digit, BVEC(N-1)
!    is the coefficient of 2**(N-2), and BVEC(N) contains sign
!    information.  It is 0 if the number is positive, and 1 if
!    the number is negative.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    08 May 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) N, the length of the vectors.
!
!    Input, integer ( kind = 4 ) BVEC1(N), BVEC2(N), the binary vectors 
!    to be XOR'ed.
!
!    Input, integer ( kind = 4 ) BVEC3(N), the exclusive OR of the two vectors.
!
  implicit none

  integer ( kind = 4 ) n

  integer ( kind = 4 ) bvec1(n)
  integer ( kind = 4 ) bvec2(n)
  integer ( kind = 4 ) bvec3(n)
  integer ( kind = 4 ), parameter :: i4_2 = 2

  bvec3(1:n) = mod ( bvec1(1:n) + bvec2(1:n), i4_2 )

  return
end
subroutine catalan ( n, c )

!*****************************************************************************80
!
!! CATALAN computes the Catalan numbers, from C(0) to C(N).
!
!  First values:
!
!     C(0)     1
!     C(1)     1
!     C(2)     2
!     C(3)     5
!     C(4)    14
!     C(5)    42
!     C(6)   132
!     C(7)   429
!     C(8)  1430
!     C(9)  4862
!    C(10) 16796
!
!  Formula:
!
!    C(N) = (2*N)! / ( (N+1) * (N!) * (N!) )
!         = 1 / (N+1) * COMB ( 2N, N )
!         = 1 / (2N+1) * COMB ( 2N+1, N+1).
!
!  Recursion:
!
!    C(N) = 2 * (2*N-1) * C(N-1) / (N+1)
!    C(N) = sum ( 1 <= I <= N-1 ) C(I) * C(N-I)
!
!  Discussion:
!
!    The Catalan number C(N) counts:
!
!    1) the number of binary trees on N vertices;
!    2) the number of ordered trees on N+1 vertices;
!    3) the number of full binary trees on 2N+1 vertices;
!    4) the number of well formed sequences of 2N parentheses;
!    5) number of ways 2N ballots can be counted, in order,
!       with N positive and N negative, so that the running sum
!       is never negative;
!    6) the number of standard tables in a 2 by N rectangular Ferrers diagram;
!    7) the number of monotone functions from [1..N} to [1..N} which
!       satisfy f(i) <= i for all i;
!    8) the number of ways to triangulate a polygon with N+2 vertices.
!
!  Example:
!
!    N = 3
!
!    ()()()
!    ()(())
!    (()())
!    (())()
!    ((()))
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    14 August 1998
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Dennis Stanton, Dennis White,
!    Constructive Combinatorics,
!    Springer, 1986,
!    ISBN: 0387963472,
!    LC: QA164.S79.
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) N, the number of Catalan numbers desired.
!
!    Output, integer ( kind = 4 ) C(0:N), the Catalan numbers from C(0) to C(N).
!
  implicit none

  integer ( kind = 4 ) n

  integer ( kind = 4 ) c(0:n)
  integer ( kind = 4 ) i

  if ( n < 0 ) then
    return
  end if

  c(0) = 1
!
!  The extra parentheses ensure that the integer division is
!  done AFTER the integer multiplication.
!
  do i = 1, n
    c(i) = ( c(i-1) * 2 * ( 2 * i - 1 ) ) / ( i + 1 )
  end do

  return
end
subroutine catalan_row_next ( ido, n, irow )

!*****************************************************************************80
!
!! CATALAN_ROW_NEXT computes row N of Catalan's triangle.
!
!  Example:
!
!    I\J 0   1   2   3   4   5   6
!
!    0   1
!    1   1   1
!    2   1   2   2
!    3   1   3   5   5
!    4   1   4   9  14  14
!    5   1   5  14  28  42  42
!    6   1   6  20  48  90 132 132
!
!  Recursion:
!
!    C(0,0) = 1
!    C(I,0) = 1
!    C(I,J) = 0 for I < J
!    C(I,J) = C(I,J-1) + C(I-1,J)
!    C(I,I) is the I-th Catalan number.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    12 December 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) IDO, indicates whether this is a call for
!    the 'next' row of the triangle.
!    IDO = 0, this is a startup call.  Row N is desired, but
!    presumably this is a first call, or row N-1 was not computed
!    on the previous call.
!    IDO = 1, this is not the first call, and row N-1 was computed
!    on the previous call.  In this case, much work can be saved
!    by using the information from the previous values of IROW
!    to build the next values.
!
!    Input, integer ( kind = 4 ) N, the index of the row of the triangle 
!    desired.  
!
!    Input/output, integer ( kind = 4 ) IROW(0:N), the row of coefficients.
!    If IDO = 0, then IROW is not required to be set on input.
!    If IDO = 1, then IROW must be set on input to the value of
!    row N-1.
!
  implicit none

  integer ( kind = 4 ) n

  integer ( kind = 4 ) i
  integer ( kind = 4 ) ido
  integer ( kind = 4 ) irow(0:n)
  integer ( kind = 4 ) j

  if ( n < 0 ) then
    return
  end if

  if ( ido == 0 ) then
 
    irow(0) = 1
    irow(1:n) = 0
 
    do i = 1, n

      irow(0) = 1

      do j = 1, i-1
        irow(j) = irow(j) + irow(j-1)
      end do

      irow(i) = irow(i-1)

    end do
 
  else
 
    irow(0) = 1

    do j = 1, n-1
      irow(j) = irow(j) + irow(j-1)
    end do

    if ( 1 <= n ) then
      irow(n) = irow(n-1)
    end if
 
  end if
 
  return
end
subroutine catalan_values ( n_data, n, c )

!*****************************************************************************80
!
!! CATALAN_VALUES returns some values of the Catalan numbers for testing.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    03 February 2003
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Milton Abramowitz, Irene Stegun,
!    Handbook of Mathematical Functions,
!    National Bureau of Standards, 1964,
!    ISBN: 0-486-61272-4,
!    LC: QA47.A34.
!
!  Parameters:
!
!    Input/output, integer ( kind = 4 ) N_DATA.
!    On input, if N_DATA is 0, the first test data is returned, and N_DATA
!    is set to 1.  On each subsequent call, the input value of N_DATA is
!    incremented and that test data item is returned, if available.  When 
!    there is no more test data, N_DATA is set to 0.
!
!    Output, integer ( kind = 4 ) N, the order of the Catalan number.
!
!    Output, integer ( kind = 4 ) C, the value of the Catalan number.
!
  implicit none

  integer ( kind = 4 ), parameter :: nmax = 11

  integer ( kind = 4 ) c
  integer ( kind = 4 ), save, dimension ( nmax ) :: c_vec = (/ &
    1, 1, 2, 5, 14, 42, 132, 429, 1430, 4862, 16796 /)
  integer ( kind = 4 ) n
  integer ( kind = 4 ) n_data
  integer ( kind = 4 ), save, dimension ( nmax ) :: n_vec = (/ &
     0,  1,  2,  3,  4, 5,  6,  7,  8,  9,  10 /)

  if ( n_data < 0 ) then
    n_data = 0
  end if

  n_data = n_data + 1

  if ( nmax < n_data ) then
    n_data = 0
    n = 0
    c = 0
  else
    n = n_vec(n_data)
    c = c_vec(n_data)
  end if

  return
end
subroutine cfrac_to_rat ( n, a, p, q )

!*****************************************************************************80
!
!! CFRAC_TO_RAT converts a monic continued fraction to an ordinary fraction.
!
!  Discussion:
!
!    The routine is given the monic or "simple" continued fraction with
!    integer coefficients:
!
!      A(1) + 1 / ( A(2) + 1 / ( A(3) ... + 1 / A(N) ) )
!
!    and returns the N successive approximants P(I)/Q(I)
!    to the value of the rational number represented by the continued
!    fraction, with the value exactly equal to the final ratio P(N)/Q(N).
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    12 June 2004
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    John Hart, Ward Cheney, Charles Lawson, Hans Maehly, Charles Mesztenyi, 
!    John Rice, Henry Thatcher, Christoph Witzgall,
!    Computer Approximations,
!    Wiley, 1968.
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) N, the number of continued fraction 
!    coefficients.
!
!    Input, integer ( kind = 4 ) A(N), the continued fraction coefficients.
!
!    Output, integer ( kind = 4 ) P(N), Q(N), the N successive approximations
!    to the value of the continued fraction.
!
  implicit none

  integer ( kind = 4 ) n

  integer ( kind = 4 ) a(n)
  integer ( kind = 4 ) i
  integer ( kind = 4 ) p(n)
  integer ( kind = 4 ) q(n)

  do i = 1, n

    if ( i == 1 ) then
      p(i) = a(i) * 1 + 0
      q(i) = a(i) * 0 + 1
    else if ( i == 2 ) then
      p(i) = a(i) * p(i-1) + 1
      q(i) = a(i) * q(i-1) + 0
    else
      p(i) = a(i) * p(i-1) + p(i-2)
      q(i) = a(i) * q(i-1) + q(i-2)
    end if

  end do

  return
end
subroutine cfrac_to_rfrac ( m, g, h, p, q )

!*****************************************************************************80
!
!! CFRAC_TO_RFRAC converts polynomial fractions from continued to rational form.
!
!  Discussion:
!
!    The routine accepts a continued polynomial fraction:
!
!      G(1)     / ( H(1) +
!      G(2) * X / ( H(2) +
!      G(3) * X / ( H(3) + ...
!      G(M) * X / ( H(M) )...) ) )
!
!    and returns the equivalent rational polynomial fraction:
!
!      P(1) + P(2) * X + ... + P(L1) * X**(L1)
!      -------------------------------------------------------
!      Q(1) + Q(2) * X + ... + Q(L2) * X**(L2-1)
!
!    where
!
!      L1 = (M+1)/2
!      L2 = (M+2)/2.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    12 June 2004
!
!  Author:
!
!    FORTRAN77 original version by John Hart, Ward Cheney, Charles Lawson, 
!    Hans Maehly, Charles Mesztenyi, John Rice, Henry Thatcher, 
!    Christoph Witzgall.
!    FORTRAN90 version by John Burkardt
!
!  Reference:
!
!    John Hart, Ward Cheney, Charles Lawson, Hans Maehly, Charles Mesztenyi,
!    John Rice, Henry Thatcher, Christoph Witzgall,
!    Computer Approximations,
!    Wiley, 1968.
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) M, the number of continued fraction 
!    polynomial coefficients.
!
!    Input, real ( kind = 8 ) G(M), H(M), the continued polynomial 
!    fraction coefficients.
!
!    Output, real ( kind = 8 ) P((M+1)/2), Q((M+2)/2), the rational 
!    polynomial fraction coefficients.
!
  implicit none

  integer ( kind = 4 ) m

  real    ( kind = 8 ) a(m,(m+2)/2)
  real    ( kind = 8 ) g(m)
  real    ( kind = 8 ) h(m)
  integer ( kind = 4 ) i
  integer ( kind = 4 ) j
  real    ( kind = 8 ) p((m+1)/2)
  real    ( kind = 8 ) q((m+2)/2)

  if ( m == 1 ) then
    p(1) = g(1)
    q(1) = h(1)
    return
  end if

  do i = 1, m
    do j = 1, (m+2)/2
      a(i,j) = 0.0D+00
    end do
  end do
!
!  Solve for P's.
!
  a(1,1) = g(1)
  a(2,1) = g(1) * h(2)

  do i = 3, m
    a(i,1) = h(i) * a(i-1,1)
    do j = 2, (i+1)/2
      a(i,j) = h(i) * a(i-1,j) + g(i) * a(i-2,j-1)
    end do
  end do

  do j = 1, (m+1)/2
    p(j) = a(m,j)
  end do
!
!  Solve for Q's.
!
  a(1,1) = h(1)
  a(2,1) = h(1) * h(2)
  a(2,2) = g(2)

  do i = 3, m
    a(i,1) = h(i) * a(i-1,1)
    do j = 2, (i+2) / 2
      a(i,j) = h(i) * a(i-1,j) + g(i) * a(i-2,j-1)
    end do
  end do

  do j = 1, (m+2)/2
    q(j) = a(m,j)
  end do

  return
end
subroutine ch_cap ( c )

!*****************************************************************************80
!
!! CH_CAP capitalizes a single character.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    19 July 1998
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input/output, character C, the character to capitalize.
!
  implicit none

  character c
  integer ( kind = 4 ) itemp

  itemp = ichar ( c )

  if ( 97 <= itemp .and. itemp <= 122 ) then
    c = char ( itemp - 32 )
  end if

  return
end
subroutine change_greedy ( total, coin_num, coin_value, change_num, change )

!*****************************************************************************80
!
!! CHANGE_GREEDY makes change for a given total using the biggest coins first.
!
!  Discussion:
!
!    The algorithm is simply to use as many of the largest coin first,
!    then the next largest, and so on.
!
!    It is assumed that there is always a coin of value 1.  The
!    algorithm will otherwise fail!
!
!  Examples:
!
!    Total = 17
!    COIN_NUM = 3
!    COIN_VALUE = (/ 1, 5, 10 /)
!
!
!    #  CHANGE              COIN_VALUE(CHANGE)
!
!    4  3 2 1 1             10 5 1 1
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    20 November 2002
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) TOTAL, the total for which change is to be 
!    made.
!
!    Input, integer ( kind = 4 ) COIN_NUM, the number of types of coins.
!
!    Input, integer ( kind = 4 ) COIN_VALUE(COIN_NUM), the value of each coin.
!    The values should be in ascending order, and if they are not,
!    they will be sorted.
!
!    Output, integer ( kind = 4 ) CHANGE_NUM, the number of coins given in 
!    change.
!
!    Output, integer ( kind = 4 ) CHANGE(TOTAL), the indices of the coins 
!    will be in entries 1 through CHANGE_NUM.
!
  implicit none

  integer ( kind = 4 ) coin_num
  integer ( kind = 4 ) total

  integer ( kind = 4 ) change(total)
  integer ( kind = 4 ) change_num
  integer ( kind = 4 ) coin_value(coin_num)
  integer ( kind = 4 ) j
  integer ( kind = 4 ) total_copy

  change_num = 0
!
!  Find the largest coin smaller than the total.
!
  j = coin_num

  do while ( 0 < j )
    if ( coin_value(j) <= total ) then
      exit
    end if
    j = j - 1
  end do

  if ( j <= 0 ) then
    return
  end if
!
!  Subtract the current coin from the total.
!  Once that coin is too big, use the next coin.
!
  total_copy = total

  do while ( 0 < total_copy )

    if ( coin_value(j) <= total_copy ) then

      total_copy = total_copy - coin_value(j)
      change_num = change_num + 1
      change(change_num) = j

    else

      j = j - 1
      if ( j <= 0 ) then
        exit
      end if

    end if

  end do

  return
end
subroutine change_next ( total, coin_num, coin_value, change_num, change, &
  done  )

!*****************************************************************************80
!
!! CHANGE_NEXT computes the next set of change for a given sum.
!
!  Examples:
!
!    Total = 17
!    COIN_NUM = 3
!    COIN_VALUE = (/ 1, 5, 10 /)
!
!
!        #  CHANGE              COIN_VALUE(CHANGE)
!
!    1   4  3 2 1 1             10 5 1 1
!    2   8  3 1 1 1 1 1 1 1     10 1 1 1 1 1 1 1
!    3   5  2 2 2 1 1            5 5 5 1 1
!    4   9  2 2 1 1 1 1 1 1 1    5 5 1 1 1 1 1 1 1
!    5  13  2 1 1 1 1 1 1 1 1 1  5 1 1 1 1 1 1 1 1 1
!           1 1 1                1 1 1
!    6  17  1 1 1 1 1 1 1 1 1 1  1 1 1 1 1 1 1 1 1 1 1
!           1 1 1 1 1 1 1        1 1 1 1 1 1
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    21 June 2004
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) TOTAL, the total for which change is to be 
!    made.
!
!    Input, integer ( kind = 4 ) COIN_NUM, the number of types of coins.
!
!    Input, integer ( kind = 4 ) COIN_VALUE(COIN_NUM), the value of each coin.
!    The values must be in ascending order.
!
!    Input/output, integer ( kind = 4 ) CHANGE_NUM, the number of coins given
!    in change for this form of the change.
!
!    Input/output, integer ( kind = 4 ) CHANGE(CHANGE_NUM), the indices of the 
!    coins.  The user must dimension this array to have dimension TOTAL!
!
!    Input/output, logical DONE.  The user sets DONE = TRUE on
!    first call to tell the routine this is the beginning of a computation.
!    The program resets DONE to FALSE and it stays that way until
!    the last possible change combination is made, at which point the
!    program sets DONE to TRUE again.
!
  implicit none

  integer ( kind = 4 ) coin_num
  integer ( kind = 4 ) total

  integer ( kind = 4 ) change(total)
  integer ( kind = 4 ) change_num
  integer ( kind = 4 ) change_num2
  integer ( kind = 4 ) coin_num2
  integer ( kind = 4 ) coin_value(coin_num)
  logical done
  integer ( kind = 4 ) i
  logical i4vec_ascends
  integer ( kind = 4 ) last
  integer ( kind = 4 ) total2

  if ( done ) then
!
!  Make sure the coin values are sorted.
!
    if ( .not. i4vec_ascends ( coin_num, coin_value ) ) then
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'CHANGE_NEXT - Fatal error!'
      write ( *, '(a)' ) '  The COIN_VALUE array is not in ascending order.'
      stop
    end if
!
!  Start with the greedy change.
!
    call change_greedy ( total, coin_num, coin_value, change_num, change )
!
!  In a few cases, like change for 4 cents, we're done after the first call.
!
    if ( change_num == total ) then
      done = .true.
    else
      done = .false.
    end if

  else
!
!  Find the last location in the input change which is NOT a penny.
!
    last = 0

    do i = change_num, 1, -1
      if ( change(i) /= 1 ) then
        last = i
        exit
      end if
    end do
!
!  If that location is still 0, an error was made.
!
    if ( last == 0 ) then
      done = .true.
      return
    end if
!
!  Sum the entries from that point to the end.
!
    total2 = sum ( coin_value ( change(last:change_num) ) )
!
!  Make greedy change for the partial sum using coins smaller than that one.
!
    coin_num2 = change(last) - 1

    call change_greedy ( total2, coin_num2, coin_value, change_num2, &
      change(last:total) )

    change_num = ( last - 1 ) + change_num2

  end if

  return
end
subroutine chinese_check ( n, m, ierror )

!*****************************************************************************80
!
!! CHINESE_CHECK checks the Chinese remainder moduluses.
!
!  Discussion:
!
!    For a Chinese remainder representation, the moduluses M(I) must
!    be positive and pairwise prime.  Also, in case this is not obvious,
!    no more than one of the moduluses may be 1.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    05 July 2000
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) N, the number of moduluses.
!
!    Input, integer ( kind = 4 ) M(N), the moduluses.  These should be positive
!    and pairwise prime.
!
!    Output, integer ( kind = 4 ) IERROR, an error flag.
!    0, no error was detected.
!    nonzero, an error was detected.
!
  implicit none

  integer ( kind = 4 ) n

  integer ( kind = 4 ) i
  integer ( kind = 4 ) ierror
  logical i4vec_pairwise_prime
  integer ( kind = 4 ) j
  integer ( kind = 4 ) m(n)

  ierror = 0
!
!  Do not allow nonpositive entries.
!
  if ( any ( m(1:n) <= 0 ) ) then
    ierror = 1
    return
  end if
!
!  Allow one entry to be 1, but not two entries.
!
  do i = 1, n
    do j = i+1, n
      if ( m(i) == 1 .and. m(j) == 1 ) then
        ierror = 2
        return
      end if
    end do
  end do
!
!  Now check pairwise primeness.
!
  if ( .not. i4vec_pairwise_prime ( n, m ) ) then
    ierror = 3
    return
  end if

  return
end
subroutine chinese_to_i4 ( n, m, r, j )

!*****************************************************************************80
!
!! CHINESE_TO_I4 converts a set of Chinese remainders to an equivalent integer.
!
!  Discussion:
!
!    Given a set of N pairwise prime, positive moduluses M(I), and
!    a corresponding set of remainders R(I), this routine finds an
!    integer J such that, for all I,
!
!      J = R(I) mod M(I)
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    10 July 2000
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) N, the number of moduluses.
!
!    Input, integer ( kind = 4 ) M(N), the moduluses.  These should be 
!    positive and pairwise prime.
!
!    Input, integer ( kind = 4 ) R(N), the Chinese remainder representation 
!    of the integer.
!
!    Output, integer ( kind = 4 ) J, the corresponding integer.
!
  implicit none

  integer ( kind = 4 ) n

  integer ( kind = 4 ) a
  integer ( kind = 4 ) b(n)
  integer ( kind = 4 ) big_m
  integer ( kind = 4 ) c
  integer ( kind = 4 ) i
  integer ( kind = 4 ) ierror
  integer ( kind = 4 ) j
  integer ( kind = 4 ) m(n)
  integer ( kind = 4 ) r(n)

  call chinese_check ( n, m, ierror )

  if ( ierror /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'CHINESE_TO_I4 - Fatal error!'
    write ( *, '(a)' ) '  The moduluses are not legal.'
    stop
  end if
!
!  Set BIG_M.
!
  big_m = product ( m(1:n) )
!
!  Solve BIG_M / M(I) * B(I) = 1, mod M(I)
!
  do i = 1, n
    a = big_m / m(i)
    c = 1
    call congruence ( a, m(i), c, ierror, b(i) )
  end do
!
!  Set J = sum ( 1 <= I <= N ) ( R(I) * B(I) * BIG_M / M(I) ) mod M
!
  j = 0
  do i = 1, n
    j = mod ( j + r(i) * b(i) * ( big_m / m(i) ), big_m )
  end do

  return
end
subroutine comb_next ( n, k, a, done )

!*****************************************************************************80
!
!! COMB_NEXT computes combinations of K things out of N.
!
!  Discussion:
!
!    The combinations are computed one at a time, in lexicographical order.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    21 June 2004
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Charles Mifsud,
!    Algorithm 154:
!    Combination in Lexicographic Order,
!    Communications of the ACM, 
!    March 1963.
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) N, the total number of things.
!
!    Input, integer ( kind = 4 ) K, the number of things in each combination.
!
!    Input/output, integer ( kind = 4 ) A(K), contains the list of elements in
!    the current combination.
!
!    Input/output, logical DONE.  On first call, set DONE to FALSE,
!    and therafter, its input value should be the output value from
!    the previous call.  The output value will normally be TRUE,
!    indicating that there are further combinations that can be
!    returned.  When DONE is returned FALSE, the sequence is exhausted.
!
  implicit none

  integer ( kind = 4 ) k

  integer ( kind = 4 ) a(k)
  logical done
  integer ( kind = 4 ) i
  integer ( kind = 4 ) j
  integer ( kind = 4 ) n

  if ( done ) then

    call i4vec_indicator ( k, a )

    if ( 1 < k ) then
      done = .false.
    else
      done = .true.
    end if

  else

    if ( a(k) < n ) then
      a(k) = a(k) + 1
      return
    end if

    do i = k, 2, -1

      if ( a(i-1) < n-k+i-1 ) then

        a(i-1) = a(i-1) + 1

        do j = i, k
          a(j) = a(i-1) + j - ( i-1 )
        end do

        return

      end if

    end do

    done = .true.

  end if

  return
end
subroutine comb_row ( ido, n, irow )

!*****************************************************************************80
!
!! COMB_ROW computes row N of Pascal's triangle.
!
!  Discussion:
!
!    Row N contains the combinatorial coefficients
!
!      C(N,0), C(N,1), C(N,2), ... C(N,N)
!
!    The sum of the elements of row N is equal to 2**N.
!
!  Formula:
!
!    C(N,K) = N! / ( K! * (N-K)! )
!
!  First terms:
!
!     N K:0  1   2   3   4   5   6   7  8  9 10
!
!     0   1
!     1   1  1
!     2   1  2   1
!     3   1  3   3   1
!     4   1  4   6   4   1
!     5   1  5  10  10   5   1
!     6   1  6  15  20  15   6   1
!     7   1  7  21  35  35  21   7   1
!     8   1  8  28  56  70  56  28   8  1
!     9   1  9  36  84 126 126  84  36  9  1
!    10   1 10  45 120 210 252 210 120 45 10  1
!
!  Recursion:
!
!    C(N,K) = C(N-1,K-1)+C(N-1,K)
!
!  Special values:
!
!    C(N,0) = C(N,N) = 1
!    C(N,1) = C(N,N-1) = N
!    C(N,N-2) = sum ( 1 <= I <= N ) N
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    09 May 2003
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) IDO, indicates whether this is a call for
!    the 'next' row of the triangle.
!
!    0 means this is a startup call.  Row N is desired, but
!    presumably this is a first call, or row N-1 was not computed
!    on the previous call.
!
!    1 means this is not the first call, and row N-1 was computed
!    on the previous call.  In this case, much work can be saved
!    by using the information from the previous values of IROW
!    to build the next values.
!
!    Input, integer ( kind = 4 ) N, the row of the triangle desired.  The 
!    triangle begins with row 0.
!
!    Output, integer ( kind = 4 ) IROW(N+1), the row of coefficients.
!    IROW(I) = C(N,I-1).
!
  implicit none

  integer ( kind = 4 ) n

  integer ( kind = 4 ) i
  integer ( kind = 4 ) ido
  integer ( kind = 4 ) irow(n+1)
  integer ( kind = 4 ) j

  if ( n < 0 ) then
    return
  end if

  if ( ido == 1 ) then

    do i = n, 2, -1
      irow(i) = irow(i) + irow(i-1)
    end do

    irow(n+1) = 1

  else

    irow(1) = 1
    irow(2:n+1) = 0

    do j = 1, n
      do i = j+1, 2, -1
        irow(i) = irow(i) + irow(i-1)
      end do
    end do

  end if

  return
end
subroutine comb_unrank ( m, n, rank, a )

!*****************************************************************************80
!
!! COMB_UNRANK returns the RANK-th combination of N things out of M.
!
!  Discussion:
!
!    Going from a rank to a thing is called "unranking".
!
!    The combinations are ordered lexically.
!
!    Lexical order can be illustrated for the general case of N and M as
!    follows:
!
!    1:       1,     2,     3,     ..., N-2, N-1, N
!    2:       1,     2,     3,     ..., N-2, N-1, N+1
!    3:       1,     2,     3,     ..., N-2, N-1, N+2
!    ...
!    M-N+1:   1,     2,     3,     ..., N-2, N-1, M
!    M-N+2:   1,     2,     3,     ..., N-2, N,   N+1
!    M-N+3:   1,     2,     3,     ..., N-2, N,   N+2
!    ...
!    LAST-2:  M-N,   M-N+1, M-N+3, ..., M-2, M-1, M
!    LAST-1:  M-N,   M-N+2, M-N+3, ..., M-2, M-1, M
!    LAST:    M-N+1, M-N+2, M-N+3, ..., M-2, M-1, M
!
!    There are a total of M!/(N!*(M-N)!) combinations of M
!    things taken N at a time.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    02 June 2007
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Bill Buckles, Matthew Lybanon,
!    Algorithm 515:
!    Generation of a Vector from the Lexicographical Index,
!    ACM Transactions on Mathematical Software,
!    Volume 3, Number 2, pages 180-182, June 1977.
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) M, the size of the set.
!
!    Input, integer ( kind = 4 ) N, the number of things in the combination.
!    N must be greater than 0, and no greater than M.
!
!    Input, integer ( kind = 4 ) RANK, the lexicographical rank of the 
!    combination sought.  RANK must be at least 1, and no greater 
!    than M!/(N!*(M-N)!).
!
!    Output, integer ( kind = 4 ) A(N), the combination.
!
  implicit none

  integer ( kind = 4 ) n

  integer ( kind = 4 ) a(n)
  integer ( kind = 4 ) i
  integer ( kind = 4 ) i4_choose
  integer ( kind = 4 ) j
  integer ( kind = 4 ) k
  integer ( kind = 4 ) m
  integer ( kind = 4 ) rank
!
!  Initialize the lower bound index.
!
  k = 0
!
!  Select elements in ascending order.
!
  do i = 1, n - 1
!
!  Set the lower bound element number for next element value.
!
    a(i) = 0

    if ( 1 < i ) then
      a(i) = a(i-1)
    end if
!
!  Check each element value.
!
    do

      a(i) = a(i) + 1
      j = i4_choose ( m-a(i), n-i )
      k = k + j

      if ( rank <= k ) then
        exit
      end if

    end do

    k = k - j

  end do

  a(n) = a(n-1) + rank - k

  return
end
subroutine comp_enum ( n, k, number )

!*****************************************************************************80
!
!! COMP_ENUM returns the number of compositions of the integer N into K parts.
!
!  Discussion:
!
!    A composition of the integer N into K parts is an ordered sequence
!    of K nonnegative integers which sum to N.  The compositions (1,2,1)
!    and (1,1,2) are considered to be distinct.
!
!    The 28 compositions of 6 into three parts are:
!
!      6 0 0,  5 1 0,  5 0 1,  4 2 0,  4 1 1,  4 0 2,
!      3 3 0,  3 2 1,  3 1 2,  3 0 3,  2 4 0,  2 3 1,
!      2 2 2,  2 1 3,  2 0 4,  1 5 0,  1 4 1,  1 3 2,
!      1 2 3,  1 1 4,  1 0 5,  0 6 0,  0 5 1,  0 4 2,
!      0 3 3,  0 2 4,  0 1 5,  0 0 6.
!
!    The formula for the number of compositions of N into K parts is
!
!      Number = ( N + K - 1 )! / ( N! * ( K - 1 )! )
!
!    Describe the composition using N '1's and K-1 dividing lines '|'.
!    The number of distinct permutations of these symbols is the number
!    of compositions.  This is equal to the number of permutations of 
!    N+K-1 things, with N identical of one kind and K-1 identical of another.
!
!    Thus, for the above example, we have:
!
!      Number = ( 6 + 3 - 1 )! / ( 6! * (3-1)! ) = 28
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    02 June 2007
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Albert Nijenhuis, Herbert Wilf,
!    Combinatorial Algorithms for Computers and Calculators,
!    Second Edition,
!    Academic Press, 1978,
!    ISBN: 0-12-519260-6,
!    LC: QA164.N54.
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) N, the integer whose compositions are desired.
!
!    Input, integer ( kind = 4 ) K, the number of parts in the composition.
!
!    Output, integer ( kind = 4 ) NUMBER, the number of compositions of N 
!    into K parts.
!
  implicit none

  integer ( kind = 4 ) i4_choose
  integer ( kind = 4 ) k
  integer ( kind = 4 ) n
  integer ( kind = 4 ) number

  number = i4_choose ( n + k - 1, n )

  return
end
subroutine comp_next ( n, k, a, more, h, t )

!*****************************************************************************80
!
!! COMP_NEXT computes the compositions of the integer N into K parts.
!
!  Discussion:
!
!    A composition of the integer N into K parts is an ordered sequence
!    of K nonnegative integers which sum to N.  The compositions (1,2,1)
!    and (1,1,2) are considered to be distinct.
!
!    The routine computes one composition on each call until there are no more.
!    For instance, one composition of 6 into 3 parts is
!    3+2+1, another would be 6+0+0.
!
!    On the first call to this routine, set MORE = FALSE.  The routine
!    will compute the first element in the sequence of compositions, and
!    return it, as well as setting MORE = TRUE.  If more compositions
!    are desired, call again, and again.  Each time, the routine will
!    return with a new composition.
!
!    However, when the LAST composition in the sequence is computed 
!    and returned, the routine will reset MORE to FALSE, signaling that
!    the end of the sequence has been reached.
!
!    This routine originally used a SAVE statement to maintain the
!    variables H and T.  I have decided that it is safer
!    to pass these variables as arguments, even though the user should
!    never alter them.  This allows this routine to safely shuffle
!    between several ongoing calculations.
!
!
!    There are 28 compositions of 6 into three parts.  This routine will
!    produce those compositions in the following order:
!
!     I         A
!     -     ---------
!     1     6   0   0
!     2     5   1   0
!     3     4   2   0
!     4     3   3   0
!     5     2   4   0
!     6     1   5   0
!     7     0   6   0
!     8     5   0   1
!     9     4   1   1
!    10     3   2   1
!    11     2   3   1
!    12     1   4   1
!    13     0   5   1
!    14     4   0   2
!    15     3   1   2
!    16     2   2   2
!    17     1   3   2
!    18     0   4   2
!    19     3   0   3
!    20     2   1   3
!    21     1   2   3
!    22     0   3   3
!    23     2   0   4
!    24     1   1   4
!    25     0   2   4
!    26     1   0   5
!    27     0   1   5
!    28     0   0   6
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    02 July 2008
!
!  Author:
!
!    FORTRAN77 original version by Albert Nijenhuis, Herbert Wilf.
!    FORTRAN90 version by John Burkardt
!
!  Reference:
!
!    Albert Nijenhuis, Herbert Wilf,
!    Combinatorial Algorithms for Computers and Calculators,
!    Second Edition,
!    Academic Press, 1978,
!    ISBN: 0-12-519260-6,
!    LC: QA164.N54.
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) N, the integer whose compositions are desired.
!
!    Input, integer ( kind = 4 ) K, the number of parts in the composition.
!
!    Input/output, integer ( kind = 4 ) A(K), the parts of the composition.
!
!    Input/output, logical MORE, set by the user to start the computation,
!    and by the routine to terminate it.
!
!    Input/output, integer ( kind = 4 )  H, T, two internal parameters needed for the
!    computation.  The user should allocate space for these in the calling
!    program, include them in the calling sequence, but never alter them!
!
  implicit none

  integer ( kind = 4 ) k

  integer ( kind = 4 ) a(k)
  integer ( kind = 4 ) h
  logical more
  integer ( kind = 4 ) n
  integer ( kind = 4 ) t
!
!  The first computation.
!
  if ( .not. more ) then

    t = n
    h = 0
    a(1) = n
    a(2:k) = 0
!
!  The next computation.
!
  else

    if ( 1 < t ) then
      h = 0
    end if

    h = h + 1
    t = a(h)
    a(h) = 0
    a(1) = t - 1
    a(h+1) = a(h+1) + 1

  end if
!
!  This is the last element of the sequence if all the
!  items are in the last slot.
!
  more = ( a(k) /= n )

  return
end
subroutine comp_random ( n, k, seed, a )

!*****************************************************************************80
!
!! COMP_RANDOM selects a random composition of the integer N into K parts.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    30 April 2003
!
!  Author:
!
!    FORTRAN77 original version by Albert Nijenhuis, Herbert Wilf.
!    FORTRAN90 version by John Burkardt
!
!  Reference:
!
!    Albert Nijenhuis, Herbert Wilf,
!    Combinatorial Algorithms for Computers and Calculators,
!    Second Edition,
!    Academic Press, 1978,
!    ISBN: 0-12-519260-6,
!    LC: QA164.N54.
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) N, the integer to be decomposed.
!
!    Input, integer ( kind = 4 ) K, the number of parts in the composition.
!
!    Input/output, integer ( kind = 4 ) SEED, a seed for the random number 
!    generator.
!
!    Output, integer ( kind = 4 ) A(K), the parts of the composition.
!
  implicit none

  integer ( kind = 4 ) k

  integer ( kind = 4 ) a(k)
  integer ( kind = 4 ) i
  integer ( kind = 4 ), parameter :: i4_1 = 1
  integer ( kind = 4 ) l
  integer ( kind = 4 ) m
  integer ( kind = 4 ) n
  integer ( kind = 4 ) seed

  call ksub_random ( n+k-i4_1, k-i4_1, seed, a )

  a(k) = n + k
  l = 0

  do i = 1, k
    m = a(i)
    a(i) = a(i) - l - 1
    l = m
  end do

  return
end
subroutine compnz_enum ( n, k, number )

!*****************************************************************************80
!
!! COMPNZ_ENUM returns the number of nonzero compositions of the N into K parts.
!
!  Discussion:
!
!    A composition of the integer N into K nonzero parts is an ordered sequence
!    of K positive integers which sum to N.  The compositions (1,2,1)
!    and (1,1,2) are considered to be distinct.
!
!    The routine computes one composition on each call until there are no more.
!    For instance, one composition of 6 into 3 parts is 3+2+1, another would
!    be 4+1+1 but 5+1+0 is not allowed since it includes a zero part.
!
!    On the first call to this routine, set MORE = FALSE.  The routine
!    will compute the first element in the sequence of compositions, and
!    return it, as well as setting MORE = TRUE.  If more compositions
!    are desired, call again, and again.  Each time, the routine will
!    return with a new composition.
!
!    However, when the LAST composition in the sequence is computed 
!    and returned, the routine will reset MORE to FALSE, signaling that
!    the end of the sequence has been reached.
!
!    The 10 compositions of 6 into three nonzero parts are:
!
!      4 1 1,  3 2 1,  3 1 2,  2 3 1,  2 2 2,  2 1 3,  
!      1 4 1,  1 3 2,  1 2 3,  1 1 4.
!
!    The formula for the number of compositions of N into K nonzero
!    parts is
!
!      Number = ( N - 1 )! / ( ( N - K )! * ( K - 1 )! )
!
!    (Describe the composition using N-K '1's and K-1 dividing lines '|'.
!    The number of distinct permutations of these symbols is the number
!    of compositions into nonzero parts.  This is equal to the number of 
!    permutations of  N-1 things, with N-K identical of one kind 
!    and K-1 identical of another.)
!
!    Thus, for the above example, we have:
!
!      Number = ( 6 - 1 )! / ( ( 6 - 3 )! * ( 3 - 1 )! ) = 10
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    02 June 2007
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Albert Nijenhuis, Herbert Wilf,
!    Combinatorial Algorithms for Computers and Calculators,
!    Second Edition,
!    Academic Press, 1978,
!    ISBN: 0-12-519260-6,
!    LC: QA164.N54.
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) N, the integer whose compositions are desired.
!
!    Input, integer ( kind = 4 ) K, the number of parts in the composition.
!
!    Output, integer ( kind = 4 ) NUMBER, the number of compositions of N into 
!    K nonzero parts.
!
  implicit none

  integer ( kind = 4 ) i4_choose
  integer ( kind = 4 ) k
  integer ( kind = 4 ) n
  integer ( kind = 4 ) number

  number = i4_choose ( n - 1, n - k )

  return
end
subroutine compnz_next ( n, k, a, more )

!*****************************************************************************80
!
!! COMPNZ_NEXT computes the compositions of the integer N into K nonzero parts.
!
!  Discussion:
!
!    A composition of the integer N into K nonzero parts is an ordered sequence
!    of K positive integers which sum to N.  The compositions (1,2,1)
!    and (1,1,2) are considered to be distinct.
!
!    The routine computes one composition on each call until there are no more.
!    For instance, one composition of 6 into 3 parts is 3+2+1, another would
!    be 4+1+1 but 5+1+0 is not allowed since it includes a zero part.
!
!    On the first call to this routine, set MORE = FALSE.  The routine
!    will compute the first element in the sequence of compositions, and
!    return it, as well as setting MORE = TRUE.  If more compositions
!    are desired, call again, and again.  Each time, the routine will
!    return with a new composition.
!
!    However, when the LAST composition in the sequence is computed 
!    and returned, the routine will reset MORE to FALSE, signaling that
!    the end of the sequence has been reached.
!
!    The 10 compositions of 6 into three nonzero parts are:
!
!      4 1 1,  3 2 1,  3 1 2,  2 3 1,  2 2 2,  2 1 3,  
!      1 4 1,  1 3 2,  1 2 3,  1 1 4.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    01 December 2005
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Albert Nijenhuis, Herbert Wilf,
!    Combinatorial Algorithms for Computers and Calculators,
!    Second Edition,
!    Academic Press, 1978,
!    ISBN: 0-12-519260-6,
!    LC: QA164.N54.
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) N, the integer whose compositions are desired.
!
!    Input, integer ( kind = 4 ) K, the number of parts in the composition.  
!    K must be less than or equal to N.
!
!    Input/output, integer ( kind = 4 ) A(K), the parts of the composition.
!
!    Input/output, logical MORE, set by the user to start the computation,
!    and by the routine to terminate it.
!
  implicit none

  integer ( kind = 4 ) k

  integer ( kind = 4 ) a(k)
  integer ( kind = 4 ), save :: h = 0
  logical more
  integer ( kind = 4 ) n
  integer ( kind = 4 ), save :: t = 0
!
!  We use the trick of computing ordinary compositions of (N-K)
!  into K parts, and adding 1 to each part.
!
  if ( n < k ) then
    more = .false.
    a(1:k) = -1
    return
  end if
!
!  The first computation.
!
  if ( .not. more ) then

    t = n - k
    h = 0
    a(1) = n - k
    a(2:k) = 0
!
!  The next computation.
!
  else

    a(1:k) = a(1:k) - 1

    if ( 1 < t ) then
      h = 0
    end if

    h = h + 1
    t = a(h)
    a(h) = 0
    a(1) = t - 1
    a(h+1) = a(h+1) + 1

  end if
!
!  This is the last element of the sequence if all the
!  items are in the last slot.
!
  more = ( a(k) /= ( n - k ) )

  a(1:k) = a(1:k) + 1

  return
end
subroutine compnz_random ( n, k, seed, a )

!*****************************************************************************80
!
!! COMPNZ_RANDOM selects a random composition of N into K nonzero parts.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    01 December 2005
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Albert Nijenhuis, Herbert Wilf,
!    Combinatorial Algorithms for Computers and Calculators,
!    Second Edition,
!    Academic Press, 1978,
!    ISBN: 0-12-519260-6,
!    LC: QA164.N54.
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) N, the integer to be decomposed.
!
!    Input, integer ( kind = 4 ) K, the number of parts in the composition.  
!    K must be no greater than N.
!
!    Input/output, integer ( kind = 4 ) SEED, a seed for the random number 
!    generator.
!
!    Output, integer ( kind = 4 ) A(K), the parts of the composition.
!
  implicit none

  integer ( kind = 4 ) k

  integer ( kind = 4 ) a(k)
  integer ( kind = 4 ) i
  integer ( kind = 4 ), parameter :: i4_1 = 1
  integer ( kind = 4 ) l
  integer ( kind = 4 ) m
  integer ( kind = 4 ) n
  integer ( kind = 4 ) seed

  if ( n < k ) then
    a(1:k) = -1
    return
  end if

  call ksub_random ( n-i4_1, k-i4_1, seed, a )

  a(k) = n
  l = 0

  do i = 1, k
    m = a(i)
    a(i) = a(i) - l - 1
    l = m
  end do

  a(1:k) = a(1:k) + 1

  return
end
subroutine congruence ( a, b, c, ierror, x )

!*****************************************************************************80
!
!! CONGRUENCE solves a congruence of the form A * X = C ( mod B ).
!
!  Discussion:
!
!    A, B and C are given integers.  The equation is solvable if and only
!    if the greatest common divisor of A and B also divides C.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    15 November 2004
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Eric Weisstein, editor,
!    CRC Concise Encylopedia of Mathematics,
!    CRC Press, 1998, page 446.
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) A, B, C, the coefficients of the Diophantine
!    equation.
!
!    Output, integer ( kind = 4 ) IERROR, error flag.
!    0, no error, X was computed.
!    1, A = B = 0, C is nonzero.
!    2, A = 0, B and C nonzero, but C is not a multiple of B.
!    3, A nonzero, B zero, C nonzero, but C is not a multiple of A.
!    4, A, B, C nonzero, but GCD of A and B does not divide C.
!    5, algorithm ran out of internal space.
!
!    Output, integer ( kind = 4 ) X, the solution of the Diophantine equation.
!    X will be between 0 and B-1.
!
  implicit none

  integer ( kind = 4 ), parameter :: nmax = 100

  integer ( kind = 4 ) a
  integer ( kind = 4 ) a_copy
  integer ( kind = 4 ) a_mag
  integer ( kind = 4 ) a_sign
  integer ( kind = 4 ) b
  integer ( kind = 4 ) b_copy
  integer ( kind = 4 ) b_mag
  integer ( kind = 4 ) b_sign
  integer ( kind = 4 ) c
  integer ( kind = 4 ) c_copy
  integer ( kind = 4 ) g
  integer ( kind = 4 ), parameter :: i4_1 = 1
  integer ( kind = 4 ) i4_gcd
  integer ( kind = 4 ) ierror
  integer ( kind = 4 ) k
  integer ( kind = 4 ) n
  integer ( kind = 4 ) q(nmax)
  logical swap
  integer ( kind = 4 ) x
  integer ( kind = 4 ) y
  integer ( kind = 4 ) z
!
!  Defaults for output parameters.
!
  ierror = 0
  x = 0
  y = 0
!
!  Special cases.
!
  if ( a == 0 .and. b == 0 .and. c == 0 ) then
    x = 0
    return
  else if ( a == 0 .and. b == 0 .and. c /= 0 ) then
    ierror = 1
    x = 0
    return
  else if ( a == 0 .and. b /= 0 .and. c == 0 ) then
    x = 0
    return
  else if ( a == 0 .and. b /= 0 .and. c /= 0 ) then
    x = 0
    if ( mod ( c, b ) /= 0 ) then
      ierror = 2
    end if
    return
  else if ( a /= 0 .and. b == 0 .and. c == 0 ) then
    x = 0
    return
  else if ( a /= 0 .and. b == 0 .and. c /= 0 ) then
    x = c / a
    if ( mod ( c, a ) /= 0 ) then
      ierror = 3
    end if
    return
  else if ( a /= 0 .and. b /= 0 .and. c == 0 ) then
    g = i4_gcd ( a, b )
    x = b / g
    return
  end if
!
!  Now handle the "general" case: A, B and C are nonzero.
!
!  Step 1: Compute the GCD of A and B, which must also divide C.
!
  g = i4_gcd ( a, b )

  if ( mod ( c, g ) /= 0 ) then
    ierror = 4
    return
  end if

  a_copy = a / g
  b_copy = b / g
  c_copy = c / g
!
!  Step 2: Split A and B into sign and magnitude.
!
  a_mag = abs ( a_copy )
  a_sign = sign ( i4_1, a_copy )
  b_mag = abs ( b_copy )
  b_sign = sign ( i4_1, b_copy )
!
!  Another special case, A_MAG = 1 or B_MAG = 1.
!
  if ( a_mag == 1 ) then
    x = a_sign * c_copy
    return
  else if ( b_mag == 1 ) then
    x = 0
    return
  end if
!
!  Step 3: Produce the Euclidean remainder sequence.
!
  if ( b_mag <= a_mag ) then

    swap = .false.
    q(1) = a_mag
    q(2) = b_mag

  else

    swap = .true.
    q(1) = b_mag
    q(2) = a_mag

  end if

  n = 3

  do

    q(n) = mod ( q(n-2), q(n-1) )

    if ( q(n) == 1 ) then
      exit
    end if

    n = n + 1

    if ( nmax < n ) then
      ierror = 5
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'CONGRUENCE - Fatal error!'
      write ( *, '(a)' ) '  Exceeded number of iterations.'
      stop
    end if

  end do
!
!  Step 4: Now go backwards to solve X * A_MAG + Y * B_MAG = 1.
!
  y = 0
  do k = n, 2, -1
    x = y
    y = ( 1 - x * q(k-1) ) / q(k)
  end do
!
!  Step 5: Undo the swapping.
!
  if ( swap ) then
    z = x
    x = y
    y = z
  end if
!
!  Step 6: Now apply signs to X and Y so that X * A + Y * B = 1.
!
  x = x * a_sign
!
!  Step 7: Multiply by C, so that X * A + Y * B = C.
!
  x = x * c_copy
!
!  Step 8: Now force 0 <= X < B.
!
  x = mod ( x, b )
!
!  Force positivity.
!
  if ( x < 0 ) then
    x = x + b
  end if

  return
end
subroutine count_pose_random ( seed, blocks, goal )

!*****************************************************************************80
!
!! COUNT_POSE_RANDOM poses a problem for the game "The Count is Good"
!
!  Discussion:
!
!    The French television show "The Count is Good" has a game that goes
!    as follows:
!
!      A number is chosen at random between 100 and 999.  This is the GOAL.
!
!      Six numbers are randomly chosen from the set 1, 2, 3, 4, 5, 6, 7, 8,
!      9, 10, 25, 50, 75, 100.  These numbers are the BLOCKS.
!
!      The player must construct a formula, using some or all of the blocks,
!      (but not more than once), and the operations of addition, subtraction,
!      multiplication and division.  Parentheses should be used to remove
!      all ambiguity.  However, it is forbidden to use subtraction in a
!      way that produces a negative result, and all division must come out
!      exactly, with no remainder.
!
!    This routine poses a sample problem from the show.  The point is,
!    to determine how to write a program that can solve such a problem.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    30 April 2003
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Raymond Seroul,
!    Programming for Mathematicians,
!    Springer Verlag, 2000, pages 355-357.
!
!  Parameters:
!
!    Input/output, integer ( kind = 4 ) SEED, a seed for the random number 
!    generator.
!
!    Output, integer ( kind = 4 ) BLOCKS(6), the six numbers available for 
!    the formula.
!
!    Output, integer ( kind = 4 ) GOAL, the goal number.
!
  implicit none

  integer ( kind = 4 ) blocks(6)
  integer ( kind = 4 ) goal
  integer ( kind = 4 ), parameter :: i4_6 = 6
  integer ( kind = 4 ), parameter :: i4_14 = 14
  integer ( kind = 4 ), parameter :: i4_100 = 100
  integer ( kind = 4 ), parameter :: i4_999 = 999
  integer ( kind = 4 ) i4_uniform
  integer ( kind = 4 ) ind(6)
  integer ( kind = 4 ) seed
  integer ( kind = 4 ), dimension ( 14 ) :: stuff = &
    (/ 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 25, 50, 75, 100 /)

  goal = i4_uniform ( i4_100, i4_999, seed )

  call ksub_random ( i4_14, i4_6, seed, ind )

  blocks(1:6) = stuff(ind(1:6))

  return
end
subroutine debruijn ( m, n, string )

!*****************************************************************************80
!
!! DEBRUIJN constructs a de Bruijn sequence.
!
!  Discussion:
!
!    Suppose we have an alphabet of M letters, and we are interested in
!    all possible strings of length N.  If M = 2 and N = 3, then we are
!    interested in the M**N strings:
!
!      000
!      001
!      010
!      011
!      100
!      101
!      110
!      111
!
!    Now, instead of making a list like this, we prefer, if possible, to
!    write a string of letters, such that every consecutive sequence of
!    N letters is one of the strings, and every string occurs once, if
!    we allow wraparound.
!
!    For the above example, a suitable sequence would be the 8 characters:
!
!      00011101(00...
!
!    where we have suggested the wraparound feature by repeating the first
!    two characters at the end.
!
!    Such a sequence is called a de Bruijn sequence.  It can easily be
!    constructed by considering a directed graph, whose nodes are all
!    M**(N-1) strings of length N-1.  A node I has a directed edge to
!    node J (labeled with character K) if the string at node J can
!    be constructed by beheading the string at node I and adding character K.
!
!    In this setting, a de Bruijn sequence is simply an Eulerian circuit
!    of the directed graph, with the edge labels being the entries of the
!    sequence.  In general, there are many distinct de Bruijn sequences
!    for the same parameter M and N.  This program will only find one
!    of them.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    06 July 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) M, the number of letters in the alphabet.
!
!    Input, integer ( kind = 4 ) N, the number of letters in a codeword.
!
!    Output, integer ( kind = 4 ) STRING(M**N), a deBruijn string.
!
  implicit none

  integer ( kind = 4 ) m
  integer ( kind = 4 ) n

  integer ( kind = 4 ) i
  integer ( kind = 4 ), parameter :: i4_1 = 1
  integer ( kind = 4 ) iedge
  integer ( kind = 4 ) inode(m**n)
  integer ( kind = 4 ) ivec(n-1)
  integer ( kind = 4 ) j
  integer ( kind = 4 ) jnode(m**n)
  integer ( kind = 4 ) jvec(n-1)
  integer ( kind = 4 ) k
  integer ( kind = 4 ) knode(m**n)
  integer ( kind = 4 ) nedge
  integer ( kind = 4 ) nnode
  logical success
  integer ( kind = 4 ) string(m**n)
  integer ( kind = 4 ) trail(m**n)
!
!  Construct the adjacency information.
!
  nnode = m**(n-1)
  nedge = m**n

  iedge = 0

  do i = 1, nnode

    call index_unrank0 ( n-i4_1, m, i, ivec )

    do k = 1, m
      jvec(1:n-2) = ivec(2:n-1)
      jvec(n-1) = k
      call index_rank0 ( n-i4_1, m, jvec, j )
      iedge = iedge + 1
      inode(iedge) = i
      jnode(iedge) = j
      knode(iedge) = k
    end do

  end do
!
!  Determine a circuit.
!
  call digraph_arc_euler ( nnode, nedge, inode, jnode, success, trail )
!
!  The string is constructed from the labels of the edges in the circuit.
!
  string(1:nedge) = knode(trail(1:nedge))

  return
end
subroutine dec_add ( mantissa1, exponent1, mantissa2, exponent2, &
  dec_digit, mantissa, exponent )

!*****************************************************************************80
!
!! DEC_ADD adds two decimal quantities.
!
!  Discussion:
!
!    A decimal value is represented by MANTISSA * 10**EXPONENT.
!
!    The routine computes
!
!      MANTISSA * 10**EXPONENT = 
!        MANTISSA1 * 10**EXPONENT1 
!      + MANTISSA2 * 10**EXPONENT2
!
!    while trying to avoid integer overflow.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    22 July 2004
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) MANTISSA1, EXPONENT1, the first number to 
!    be added.
!
!    Input, integer ( kind = 4 ) MANTISSA2, EXPONENT2, the second number to 
!    be added.
!
!    Input, integer ( kind = 4 ) DEC_DIGIT, the number of decimal digits.
!
!    Output, integer ( kind = 4 ) MANTISSA, EXPONENT, the sum.
!
  implicit none

  integer ( kind = 4 ) dec_digit
  integer ( kind = 4 ) exponent
  integer ( kind = 4 ) exponent1
  integer ( kind = 4 ) exponent2
  integer ( kind = 4 ) mantissa
  integer ( kind = 4 ) mantissa1
  integer ( kind = 4 ) mantissa2
  integer ( kind = 4 ) mantissa3
  integer ( kind = 4 ) mantissa4

  if ( mantissa1 == 0 ) then
    mantissa = mantissa2
    exponent = exponent2
    return
  else if ( mantissa2 == 0 ) then
    mantissa = mantissa1
    exponent = exponent1
    return
  else if ( exponent1 == exponent2 ) then
    mantissa = mantissa1 + mantissa2
    exponent = exponent1
    call dec_round ( mantissa, exponent, dec_digit, mantissa, exponent )
    return
  end if
!
!  Line up the exponents.
!
  mantissa3 = mantissa1
  mantissa4 = mantissa2

  if ( exponent1 < exponent2 ) then
    mantissa4 = mantissa4 * 10**( exponent2 - exponent1 )
  else
    mantissa3 = mantissa3 * 10**( exponent1 - exponent2 )
  end if
!
!  Add the coefficients.
!
  mantissa = mantissa3 + mantissa4
  exponent = min ( exponent1, exponent2 )
!
!  Clean up the result.
!
  call dec_round ( mantissa, exponent, dec_digit, mantissa, exponent )

  return
end
subroutine dec_div ( mantissa1, exponent1, mantissa2, exponent2, &
  dec_digit, mantissa, exponent, ierror )

!*****************************************************************************80
!
!! DEC_DIV divides two decimal values.
!
!  Discussion:
!
!    A decimal value is represented by MANTISSA * 10**EXPONENT.
!
!    The routine computes
!
!      MANTISSA * 10**EXPONENT
!      = ( MANTISSA1 * 10**EXPONENT1 ) / ( MANTISSA2 * 10**EXPONENT2 )
!      = ( MANTISSA1 / MANTISSA2 ) * 10**( EXPONENT1 - EXPONENT2 )
!
!    while avoiding integer overflow.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    22 July 2004
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) MANTISSA1, EXPONENT1, the numerator.
!
!    Input, integer ( kind = 4 ) MANTISSA2, EXPONENT2, the denominator.
!
!    Input, integer ( kind = 4 ) DEC_DIGIT, the number of decimal digits.
!
!    Output, integer ( kind = 4 ) MANTISSA, EXPONENT, the result.
!
!    Output, integer ( kind = 4 ) IERROR.
!    0, no error occurred.
!    1, an error occurred.
!
  implicit none

  integer ( kind = 4 ) dec_digit
  real    ( kind = 8 ) dval
  integer ( kind = 4 ) exponent
  integer ( kind = 4 ) exponent1
  integer ( kind = 4 ) exponent2
  integer ( kind = 4 ) exponent3
  integer ( kind = 4 ) ierror
  integer ( kind = 4 ) mantissa
  integer ( kind = 4 ) mantissa1
  integer ( kind = 4 ) mantissa2
  integer ( kind = 4 ) mantissa3
!
!  First special case, top fraction is 0.
!
  if ( mantissa1 == 0 ) then
    mantissa = 0
    exponent = 0
    return
  end if
!
!  First error, bottom of fraction is 0.
!
  if ( mantissa2 == 0 ) then
    ierror = 1
    mantissa = 0
    exponent = 0
    return
  end if
!
!  Second special case, result is 1.
!
  if ( mantissa1 == mantissa2 .and. exponent1 == exponent2 ) then
    mantissa = 1
    exponent = 0
    return
  end if
!
!  Third special case, result is power of 10.
!
  if ( mantissa1 == mantissa2 ) then
    mantissa = 1
    exponent = exponent1 - exponent2
    return
  end if
!
!  Fourth special case: MANTISSA1/MANTISSA2 is exact.
!
  if ( ( mantissa1 / mantissa2 ) * mantissa2 == mantissa1 ) then
    mantissa = mantissa1 / mantissa2
    exponent = exponent1 - exponent2
    return
  end if
!
!  General case.
!
  dval = real ( mantissa1, kind = 8 ) / real ( mantissa2, kind = 8 )

  call r8_to_dec ( dval, dec_digit, mantissa3, exponent3 )

  mantissa = mantissa3
  exponent = exponent3 + exponent1 - exponent2

  return
end
subroutine dec_mul ( mantissa1, exponent1, mantissa2, exponent2, &
  dec_digit, mantissa, exponent )

!*****************************************************************************80
!
!! DEC_MUL multiplies two decimals.
!
!  Discussion:
!
!    A decimal value is represented by MANTISSA * 10**EXPONENT.
!
!    The routine computes
!
!      MANTISSA * 10**EXPONENT 
!      = ( MANTISSA1 * 10**EXPONENT1) * (MANTISSA2 * 10**EXPONENT2)
!      = ( MANTISSA1 * MANTISSA2 ) * 10**( EXPONENT1 + EXPONENT2 )
!
!    while avoiding integer overflow.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    22 July 2004
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) MANTISSA1, EXPONENT1, the first multiplier.
!
!    Input, integer ( kind = 4 ) MANTISSA2, EXPONENT2, the second multiplier.
!
!    Input, integer ( kind = 4 ) DEC_DIGIT, the number of decimal digits.
!
!    Output, integer ( kind = 4 ) MANTISSA, EXPONENT, the product.
!
  implicit none

  integer ( kind = 4 ) dec_digit
  real    ( kind = 8 ) dval
  integer ( kind = 4 ) exponent
  integer ( kind = 4 ) exponent1
  integer ( kind = 4 ) exponent2
  integer ( kind = 4 ) exponent3
  integer ( kind = 4 ) i_max
  integer ( kind = 4 ) i4_huge
  integer ( kind = 4 ) mantissa
  integer ( kind = 4 ) mantissa1
  integer ( kind = 4 ) mantissa2
  integer ( kind = 4 ) mantissa3
  real    ( kind = 8 ) temp

  i_max = i4_huge ( )
!
!  The result is zero if either MANTISSA1 or MANTISSA2 is zero.
!
  if ( mantissa1 == 0 .or. mantissa2 == 0 ) then
    mantissa = 0
    exponent = 0
    return
  end if
!
!  The result is simple if either MANTISSA1 or MANTISSA2 is one.
!
  if ( abs ( mantissa1 ) == 1 .or. abs ( mantissa2 ) == 1 ) then
    mantissa = mantissa1 * mantissa2
    exponent = exponent1 + exponent2
    return
  end if

  temp = log ( real ( abs ( mantissa1 ), kind = 8 ) ) &
       + log ( real ( abs ( mantissa2 ), kind = 8 ) )

  if ( temp < log ( real ( i_max, kind = 8 ) ) ) then

    mantissa = mantissa1 * mantissa2
    exponent = exponent1 + exponent2

  else

    dval = real ( mantissa1, kind = 8 ) * real ( mantissa2, kind = 8 )

    call r8_to_dec ( dval, dec_digit, mantissa3, exponent3 )

    mantissa = mantissa3
    exponent = exponent3 + ( exponent1 + exponent2 )

  end if

  call dec_round ( mantissa, exponent, dec_digit, mantissa, exponent )

  return
end
subroutine dec_round ( mantissa1, exponent1, dec_digit, mantissa2, exponent2 )

!*****************************************************************************80
!
!! DEC_ROUND rounds a decimal fraction to a given number of digits.
!
!  Discussion:
!
!    A decimal value is represented by MANTISSA * 10**EXPONENT.
!
!    The routine takes an arbitrary decimal fraction makes sure that 
!    MANTISSA has no more than DEC_DIGIT digits.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    22 July 2004
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) MANTISSA1, EXPONENT1, the coefficient and
!    exponent of a decimal fraction to be rounded.
!
!    Input, integer ( kind = 4 ) DEC_DIGIT, the number of decimal digits.
!
!    Output, integer ( kind = 4 ) MANTISSA2, EXPONENT2, the rounded coefficient
!    and  exponent of a decimal fraction.  MANTISSA2 has no more than
!    DEC_DIGIT decimal digits.
!
  implicit none

  integer ( kind = 4 ) dec_digit
  integer ( kind = 4 ) exponent1
  integer ( kind = 4 ) exponent2
  integer ( kind = 4 ) mantissa1
  integer ( kind = 4 ) mantissa2

  mantissa2 = mantissa1
  exponent2 = exponent1

  if ( mantissa2 == 0 ) then
    exponent2 = 0
    return
  end if

  do while ( 10**dec_digit <= abs ( mantissa2 ) )
    mantissa2 = nint ( real ( mantissa2, kind = 8 ) / 10.0D+00 )
    exponent2 = exponent2 + 1
  end do
!
!  Absorb trailing 0's into the exponent.
!
  do while ( ( mantissa2 / 10 ) * 10 == mantissa2 )
    mantissa2 = mantissa2 / 10
    exponent2 = exponent2 + 1
  end do

  return
end
subroutine dec_to_r8 ( mantissa, exponent, r )

!*****************************************************************************80
!
!! DEC_TO_R8 converts a decimal to an R8.
!
!  Discussion:
!
!    A decimal value is represented by MANTISSA * 10**EXPONENT.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    12 May 2003
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) MANTISSA, EXPONENT, the coefficient and 
!    exponent of the decimal value.
!
!    Output, real ( kind = 8 ) R, the equivalent real value.
!
  implicit none

  integer ( kind = 4 ) exponent
  integer ( kind = 4 ) mantissa
  real    ( kind = 8 ) r

  r = mantissa * 10.0D+00**exponent

  return
end
subroutine dec_to_rat ( mantissa, exponent, rat_top, rat_bot )

!*****************************************************************************80
!
!! DEC_TO_RAT converts a decimal to a rational representation.
!
!  Discussion:
!
!    A decimal value is represented by MANTISSA * 10**EXPONENT.
!
!    A rational value is represented by RAT_TOP / RAT_BOT.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    07 May 2003
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) MANTISSA, EXPONENT, the decimal number.
!
!    Output, integer ( kind = 4 ) RAT_TOP, RAT_BOT, the rational value.
!
  implicit none

  integer ( kind = 4 ) gcd
  integer ( kind = 4 ) exponent
  integer ( kind = 4 ) i4_gcd
  integer ( kind = 4 ) mantissa
  integer ( kind = 4 ) rat_bot
  integer ( kind = 4 ) rat_top

  if ( exponent == 0 ) then
    rat_top = mantissa
    rat_bot = 1
  else if ( 0 < exponent ) then
    rat_top = mantissa * 10**exponent
    rat_bot = 1
  else
    rat_top = mantissa
    rat_bot = 10**( - exponent )
    gcd = i4_gcd ( rat_top, rat_bot )
    rat_top = rat_top / gcd
    rat_bot = rat_bot / gcd
  end if

  return
end
subroutine dec_to_s ( mantissa, exponent, s )

!*****************************************************************************80
!
!! DEC_TO_S returns a string representation of a decimal.
!
!  Discussion:
!
!    A decimal value is represented by MANTISSA * 10**EXPONENT.
!
!  Examples:
!
!    MANTISSA EXPONENT   S
!    ----     ----       ------
!       0        0       0
!      21        3       21000
!      -3        0       -3
!     147       -2       14.7
!      16       -5       0.00016
!      34       30       Inf
!     123      -21       0.0000000000000000012
!      34      -30       0.0
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    13 September 1998
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) MANTISSA, EXPONENT, integers which represent
!    the decimal.
!
!    Output, character ( len = * ) S, the representation of the value.
!    The string is 'Inf' or '0.0' if the value was too large
!    or small to represent with a fixed point format.
!
  implicit none

  character ( len = 22 ) chrrep
  integer ( kind = 4 ) exponent
  integer ( kind = 4 ) i
  integer ( kind = 4 ) iget1
  integer ( kind = 4 ) iget2
  integer ( kind = 4 ) iput1
  integer ( kind = 4 ) iput2
  integer ( kind = 4 ) mantissa
  integer ( kind = 4 ) maxdigit
  integer ( kind = 4 ) ndigit
  integer ( kind = 4 ) nleft
  character ( len = * ) s

  s = ' '

  if ( mantissa == 0 ) then
    s = '0'
    return
  end if

  maxdigit = len ( s )
!
!  Store a representation of MANTISSA in CHRREP.
!
  write ( chrrep, '(i22)' ) mantissa
  call s_blank_delete ( chrrep )
  ndigit = len_trim ( chrrep )
!
!  Overflow if EXPONENT is positive, and MAXDIGIT < NDIGIT + EXPONENT.
!
  if ( 0 < exponent ) then
    if ( maxdigit < ndigit + exponent ) then
      s = 'Inf'
      return
    end if
  end if
!
!  Underflow if EXPONENT is negative, and MAXDIGIT < 3 + NDIGIT - EXPONENT.
!
  if ( exponent < 0 ) then
    if ( 0 < mantissa ) then
      if ( maxdigit < 3 - ndigit - exponent ) then
        s = '0.0'
        return
      end if
    else
      if ( maxdigit < 5 - ndigit - exponent ) then
        s = '0.0'
        return
      end if
    end if
  end if
!
!  If EXPONENT is nonnegative, insert trailing zeros.
!
  if ( 0 <= exponent ) then

    s(1:ndigit) = chrrep(1:ndigit)

    do i = ndigit + 1, ndigit + exponent
      s(i:i) = '0'
    end do

  else if ( exponent < 0 ) then

    iput2 = 0
    iget2 = 0
!
!  Sign.
!
    if ( mantissa < 0 ) then
      iput1 = 1
      iput2 = 1
      iget2 = 1
      s(iput1:iput2) = '-'
      ndigit = ndigit - 1
    end if
!
!  Digits of the integral part.
!
    if ( 0 < ndigit + exponent ) then
      iput1 = iput2 + 1
      iput2 = iput1 + ndigit + exponent -1
      iget1 = iget2 + 1
      iget2 = iget1 + ndigit + exponent - 1
      s(iput1:iput2) = chrrep(iget1:iget2)
    else
      iput1 = iput2 + 1
      iput2 = iput1
      s(iput1:iput2) = '0'
    end if
!
!  Decimal point.
!
    iput1 = iput2 + 1
    iput2 = iput1
    s(iput1:iput2) = '.'
!
!  Leading zeroes.
!
    do i = 1, - exponent - ndigit
      iput1 = iput2 + 1
      iput2 = iput1
      s(iput1:iput2) = '0'
    end do

    nleft = min ( -exponent, ndigit )
    nleft = min ( nleft, maxdigit - iput2 )
    iput1 = iput2 + 1
    iput2 = iput1 + nleft - 1
    iget1 = iget2 + 1
    iget2 = iget1 + nleft - 1
    s(iput1:iput2) = chrrep(iget1:iget2)

  end if

  return
end
function dec_width ( mantissa, exponent )

!*****************************************************************************80
!
!! DEC_WIDTH returns the "width" of a decimal number.
!
!  Discussion:
!
!    A decimal value is represented by MANTISSA * 10**EXPONENT.
!
!    The "width" of a decimal number is the number of characters
!    required to print it.
!
!  Example:
!
!    Mantissa  Exponent Width  Representation:
!
!         523      -1       4           5.23
!         134       2       5       13400
!           0      10       1           0
!      123456     -10      12           0.0000123456
!      123456      -3       7         123.456
!      123456       0       6      123456
!      123456       3       9   123456000
! 
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    01 June 2003
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) MANTISSA, EXPONENT, the decimal number.
!
!    Output, integer ( kind = 4 ) DEC_WIDTH, the "width" of the decimal number.
!
  implicit none

  integer ( kind = 4 ) dec_width
  integer ( kind = 4 ) exponent
  integer ( kind = 4 ) mantissa
  integer ( kind = 4 ) mantissa_abs
  integer ( kind = 4 ) ten_pow
  integer ( kind = 4 ) value

  value = 1
  ten_pow = 10

  if ( mantissa == 0 ) then
    dec_width = value
    return
  end if
  
  mantissa_abs = abs ( mantissa )

  do while ( ten_pow <= mantissa_abs )
    value = value + 1
    ten_pow = ten_pow * 10
  end do

  if ( 0 < exponent ) then
    value = value + exponent
  else if ( exponent < 0 ) then
    value = max ( value, 1-exponent )
!
!  An internal decimal point adds one position.
!
    if ( 0 < value ) then
      value = value + 1
!
!  A leading "0." adds two positions.
!
    else
      value = 2 - value
    end if
  end if

  if ( mantissa < 0 ) then
    value = value + 1
  end if

  dec_width = value

  return
end
subroutine decmat_det ( n, atop, abot, dec_digit, dtop, dbot, ierror )

!*****************************************************************************80
!
!! DECMAT_DET finds the determinant of an N by N matrix of decimal entries.
!
!  Discussion:
!
!    The brute force method is used.  The routine should only be used for
!    small matrices, since this calculation requires the summation of N!
!    products of N numbers.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    22 July 2004
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) N, the number of rows and columns of A.
!
!    Input, integer ( kind = 4 ) ATOP(N,N), ABOT(N,N), the decimal
!    representation of the matrix.
!
!    Output, integer ( kind = 4 ) DTOP, DBOT, the decimal determinant of 
!    the matrix.
!
!    Output, integer ( kind = 4 ) IERROR.
!    0, no error occurred.
!    1, an error occurred (probably overflow).
!
  implicit none

  integer ( kind = 4 ) n

  integer ( kind = 4 ) dec_digit
  logical even
  integer ( kind = 4 ) i
  integer ( kind = 4 ) abot(n,n)
  integer ( kind = 4 ) atop(n,n)
  integer ( kind = 4 ) iarray(n)
  integer ( kind = 4 ) ibot
  integer ( kind = 4 ) ibot1
  integer ( kind = 4 ) ibot2
  integer ( kind = 4 ) dbot
  integer ( kind = 4 ) dtop
  integer ( kind = 4 ) ierror
  integer ( kind = 4 ) itop
  integer ( kind = 4 ) itop1
  integer ( kind = 4 ) itop2
  logical more

  ierror = 0
  more = .false.
  dtop = 0
  dbot = 1
!
!  Compute the next permutation.
!
  do

    call perm_next ( n, iarray, more, even )
!
!  The sign of this term depends on the sign of the permutation.
!
    if ( even ) then
      itop = 1
    else
      itop = -1
    end if
!
!  Choose one item from each row, as specified by the permutation,
!  and multiply them
!
    ibot = 0

    do i = 1, n

      itop1 = itop
      ibot1 = ibot
      itop2 = atop(i,iarray(i))
      ibot2 = abot(i,iarray(i))

      call dec_mul ( itop1, ibot1, itop2, ibot2, dec_digit, itop, ibot )

    end do
!
!  Add this term to the total.
!
    itop1 = itop
    ibot1 = ibot

    call dec_add ( itop1, ibot1, dtop, dbot, dec_digit, itop, ibot )

    dtop = itop
    dbot = ibot

    if ( .not. more ) then
      exit
    end if

  end do

  return
end
subroutine decmat_print ( m, n, a, b, title )

!*****************************************************************************80
!
!! DECMAT_PRINT prints out decimal vectors and matrices.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    12 July 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) M, N, the number of rows and columns in 
!    the matrix.
!
!    Input, integer ( kind = 4 ) A(M,N), B(M,N), the decimal matrix.
!
!    Input, character ( len = * ) TITLE, a label for the object being printed.
!
  implicit none

  integer ( kind = 4 ) m
  integer ( kind = 4 ) n

  integer ( kind = 4 ) a(m,n)
  integer ( kind = 4 ) b(m,n)
  character ( len = 22 ) chrtmp
  character ( len = 10 ) chrtmp2
  character ( len = 10 ) chrtmp3
  character ( len = 40 ) format2
  integer ( kind = 4 ) i
  integer ( kind = 4 ) imax
  integer ( kind = 4 ) j
  integer ( kind = 4 ) jmax
  integer ( kind = 4 ) jmin
  integer ( kind = 4 ) khi
  integer ( kind = 4 ) klo
  integer ( kind = 4 ) kmax
  integer ( kind = 4 ) lenc
  integer ( kind = 4 ), parameter :: ncolum = 80
  integer ( kind = 4 ) npline
  character ( len = 100 ) output
  character ( len = * ) title
!
!  Figure out how wide we must make each column.
!
  imax = 0
  jmax = 0

  do i = 1, m
    do j = 1, n

      call dec_to_s ( a(i,j), b(i,j), chrtmp )
      lenc = len_trim ( chrtmp )
      jmax = max ( jmax, lenc )

    end do
  end do

  kmax = 2 + imax + 1 + jmax
  npline = ncolum / kmax
!
!  Set up the format for the heading.
!
    call i4_to_s_left ( npline, chrtmp2 )
    call i4_to_s_left ( kmax, chrtmp3 )
    format2 = '(' // chrtmp2 // 'i' // chrtmp3 // ')'

  call s_blank_delete ( format2 )

  do jmin = 1, n, npline

    jmax = min ( jmin + npline - 1, n )

    if ( 0 < len_trim ( title ) ) then
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) trim ( title )
    end if

    write ( *, '(a)' ) ' '

    if ( 1 < jmin .or. jmax < n ) then
      write ( output, * ) 'Columns ', jmin, ' to ', jmax
      call s_blanks_delete ( output )
      write ( *, '(a)' ) trim ( output )
      write ( *, '(a)' ) ' '
    end if

    do i = 1, m

      output = ' '

      do j = jmin, jmax
        klo = 4 + ( j - jmin ) * kmax + 1
        khi = 4 + ( j - jmin ) * kmax + kmax
        call dec_to_s ( a(i,j), b(i,j), chrtmp )
        output(klo:khi) = adjustr ( chrtmp(1:kmax) )
      end do

      write ( *, '(a)' ) trim ( output )

    end do

  end do

  return
end
subroutine derange_back_candidate ( n, maxstack, a, k, nstack, stack, ncan )

!*****************************************************************************80
!
!! DERANGE_BACK_CANDIDATE finds values for the K-th entry of a derangement.
!
!  Discussion:
!
!    A derangement of N objects is a permutation which leaves no object
!    unchanged.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    10 August 2000
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) N, the order of the derangement.
!
!    Input, integer ( kind = 4 ) MAXSTACK, the maximum stack length.
!
!    Input, integer ( kind = 4 ) A(N).  The first K-1 entries of A
!    record the currently set values of the derangement.
!
!    Input, integer ( kind = 4 ) K, the entry of the derangement for which 
!    candidates are to be found.
!
!    Input/output, integer ( kind = 4 ) NSTACK, the length of the stack.
!
!    Input/output, integer ( kind = 4 ) STACK(MAXSTACK).  On output, we have 
!    added the candidates for entry K to the end of the stack.
!
!    Input/output, integer ( kind = 4 ) NCAN(N), the number of candidates 
!    for each level.
!
  implicit none

  integer ( kind = 4 ) maxstack
  integer ( kind = 4 ) n

  integer ( kind = 4 ) a(n)
  integer ( kind = 4 ), parameter :: i4_1 = 1
  integer ( kind = 4 ) ican
  integer ( kind = 4 ) ifree(n)
  integer ( kind = 4 ) k
  integer ( kind = 4 ) ncan(n)
  integer ( kind = 4 ) nfree
  integer ( kind = 4 ) nstack
  integer ( kind = 4 ) stack(maxstack)
!
!  Consider all the integers from 1 through N that have not been used yet.
!
  nfree = n - k + 1

  call perm_free ( a, k-i4_1, nfree, ifree )
!
!  Everything but K is a legitimate candidate for the K-th entry.
!
  ncan(k) = 0

  do ican = 1, nfree

    if ( ifree(ican) /= k ) then

      ncan(k) = ncan(k) + 1
      nstack = nstack + 1

      if ( maxstack < nstack ) then
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'DERANGE_BACK_CANDIDATE - Fatal error!'
        write ( *, '(a,i8)' ) '  Exceeding stacksize limit of ', maxstack
        stop
      end if

      stack(nstack) = ifree(ican)

    end if

  end do

  return
end
subroutine derange_back_next ( n, a, more )

!*****************************************************************************80
!
!! DERANGE_BACK_NEXT returns the next derangement of N items.
!
!  Discussion:
!
!    A derangement of N objects is a permutation which leaves no object
!    unchanged.
!
!    This routine uses backtracking.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    16 July 2004
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) N, the number of items to be deranged.  
!    N should be at least 2.
!
!    Input/output, integer ( kind = 4 ) A(N).
!    On the first call, the input value of A is not important.
!    On return with MORE = TRUE, A contains the next derangement.
!    On subsequent input, A should not be changed.
!
!    Input/output, logical MORE.
!    On first call, set MORE to FALSE, and do not alter it after.
!    On return, MORE is TRUE if another derangement is being treturned in A, 
!    and FALSE if no more can be found.
!
  implicit none

  integer ( kind = 4 ) n

  integer ( kind = 4 ) a(n)
  integer ( kind = 4 ), save :: indx = 0
  integer ( kind = 4 ), save :: k = 0
  integer ( kind = 4 ), save :: maxstack = 0
  logical more
  integer ( kind = 4 ), save, allocatable, dimension (:) :: ncan
  integer ( kind = 4 ), save :: nstack = 0
  integer ( kind = 4 ), save, allocatable, dimension (:) :: stack

  if ( .not. more ) then

    if ( n < 2 ) then
      more = .false.
      return
    end if

    indx = 0
    k = 0
    maxstack = ( n * ( n + 1 ) ) / 2
    nstack = 0

    if ( allocated ( stack ) ) then
      deallocate ( stack )
    end if

    allocate ( stack(1:maxstack) )
    stack(1:maxstack) = 0

    if ( allocated ( ncan ) ) then
      deallocate ( ncan )
    end if

    allocate ( ncan(1:n) )
    ncan(1:n) = 0

    more = .true.

  end if

  do

    call i4vec_backtrack ( n, maxstack, stack, a, indx, k, nstack, ncan )

    if ( indx == 1 ) then

      exit

    else if ( indx == 2 ) then

      call derange_back_candidate ( n, maxstack, a, k, nstack, stack, ncan )
 
    else

      more = .false.
      exit

    end if

  end do

  return
end
subroutine derange_check ( n, a, deranged )

!*****************************************************************************80
!
!! DERANGE_CHECK determines whether a permutation is a derangement.
!
!  Discussion:
!
!    A derangement of N objects is a permutation which leaves no object
!    unchanged.
!
!    A derangement of the integers 1 through N is a permutation of the
!    integers such that the first value is not 1, the second is not 2,
!    and so on.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    22 September 1998
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) N, the number of objects permuted.
!
!    Input, integer ( kind = 4 ) A(N), a permutation of the integers 1 
!    through N.
!
!    Output, logical DERANGED, is TRUE if A is a derangement, and
!    FALSE otherwise.
!
  implicit none

  integer ( kind = 4 ) n

  integer ( kind = 4 ) a(n)
  logical deranged
  integer ( kind = 4 ) i

  do i = 1, n
    if ( a(i) == i ) then
      deranged = .false.
      return
    end if
  end do

  deranged = .true.

  return
end
function derange_enum ( n )

!*****************************************************************************80
!
!! DERANGE_ENUM returns the number of derangements of N objects.
!
!  Discussion:
!
!    A derangement of N objects is a permutation which leaves no object
!    unchanged.
!
!    A derangement of N objects is a permutation with no fixed
!    points.  If we symbolize the permutation operation by "P",
!    then for a derangment, P(I) is never equal to I.
!
!    The number of derangements of N objects is sometimes called
!    the subfactorial function, or the derangement number D(N).
!
!    D(N) is the number of ways of placing N non-attacking rooks on
!    an N by N chessboard with one diagonal deleted.
!
!    Limit ( N -> Infinity ) D(N)/N! = 1 / e.
!
!    The number of permutations with exactly K items in the right
!    place is COMB(N,K) * D(N-K).
!
!  Recursion:
!
!      D(0) = 1
!      D(1) = 0
!      D(2) = 1
!      D(N) = (N-1) * ( D(N-1) + D(N-2) )
!
!    or
!
!      D(0) = 1
!      D(1) = 0
!      D(N) = N * D(N-1) + (-1)**N
!
!  Formula:
!
!    D(N) = N! * ( 1 - 1/1! + 1/2! - 1/3! ... 1/N! )
!
!    Based on the inclusion/exclusion law.
!
!    D(N) = nint ( N! / E )
!
!  First values:
!
!     N         D(N)
!     0           1
!     1           0
!     2           1
!     3           2
!     4           9
!     5          44
!     6         265
!     7        1854
!     8       14833
!     9      133496
!    10     1334961
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    01 February 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) N, the number of objects to be permuted.
!
!    Output, integer ( kind = 4 ) DERANGE_ENUM, the number of derangements 
!    of N objects.
!
  implicit none

  integer ( kind = 4 ) derange_enum
  integer ( kind = 4 ) dn
  integer ( kind = 4 ) dnm1
  integer ( kind = 4 ) dnm2
  integer ( kind = 4 ) i
  integer ( kind = 4 ) n

  if ( n < 0 ) then

    dn = 0

  else if ( n == 0 ) then

    dn = 1

  else if ( n == 1 ) then

    dn = 0

  else if ( n == 2 ) then

    dn = 1

  else

    dnm1 = 0
    dn = 1

    do i = 3, n
      dnm2 = dnm1
      dnm1 = dn
      dn = ( i - 1 ) * ( dnm1 + dnm2 )
    end do

  end if

  derange_enum = dn

  return
end
subroutine derange_enum2 ( n, d )

!*****************************************************************************80
!
!! DERANGE_ENUM2 returns the number of derangements of 0 through N objects.
!
!  Discussion:
!
!    A derangement of N objects is a permutation which leaves no object
!    unchanged.
!
!    A derangement of N objects is a permutation with no fixed
!    points.  If we symbolize the permutation operation by "P",
!    then for a derangment, P(I) is never equal to I.
!
!    The number of derangements of N objects is sometimes called
!    the subfactorial function, or the derangement number D(N).
!
!    D(N) is the number of ways of placing N non-attacking rooks on
!    an N by N chessboard with one diagonal deleted.
!
!    Limit ( N -> Infinity ) D(N)/N! = 1 / e.
!
!    The number of permutations with exactly K items in the right
!    place is COMB(N,K) * D(N-K).
!
!  Recursion:
!
!      D(0) = 1
!      D(1) = 0
!      D(2) = 1
!      D(N) = (N-1) * ( D(N-1) + D(N-2) )
!
!    or
!
!      D(0) = 1
!      D(1) = 0
!      D(N) = N * D(N-1) + (-1)**N
!
!  Formula:
!
!    D(N) = N! * ( 1 - 1/1! + 1/2! - 1/3! ... 1/N! )
!
!    Based on the inclusion/exclusion law.
!
!    D(N) = nint ( N! / E )
!
!  First values:
!
!     N         D(N)
!     0           1
!     1           0
!     2           1
!     3           2
!     4           9
!     5          44
!     6         265
!     7        1854
!     8       14833
!     9      133496
!    10     1334961
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    28 April 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) N, the maximum number of objects to be 
!    permuted.
!
!    Output, integer ( kind = 4 ) D(0:N); D(I) is the number of derangements of
!    I objects.
!
  implicit none

  integer ( kind = 4 ) n

  integer ( kind = 4 ) d(0:n)
  integer ( kind = 4 ) i

  d(0) = 1
  d(1) = 0

  do i = 2, n
    d(i) = ( i - 1 ) * ( d(i-1) + d(i-2) )
  end do

  return
end
function derange_enum3 ( n )

!*****************************************************************************80
!
!! DERANGE_ENUM3 returns the number of derangements of 0 through N objects.
!
!  Discussion:
!
!    A derangement of N objects is a permutation which leaves no object
!    unchanged.
!
!    A derangement of N objects is a permutation with no fixed
!    points.  If we symbolize the permutation operation by "P",
!    then for a derangment, P(I) is never equal to I.
!
!    The number of derangements of N objects is sometimes called
!    the subfactorial function, or the derangement number D(N).
!
!    D(N) is the number of ways of placing N non-attacking rooks on
!    an N by N chessboard with one diagonal deleted.
!
!    Limit ( N -> Infinity ) D(N)/N! = 1 / e.
!
!    The number of permutations with exactly K items in the right
!    place is COMB(N,K) * D(N-K).
!
!  Recursion:
!
!      D(0) = 1
!      D(1) = 0
!      D(2) = 1
!      D(N) = (N-1) * ( D(N-1) + D(N-2) )
!
!    or
!
!      D(0) = 1
!      D(1) = 0
!      D(N) = N * D(N-1) + (-1)**N
!
!  Formula:
!
!    D(N) = N! * ( 1 - 1/1! + 1/2! - 1/3! ... 1/N! )
!
!    Based on the inclusion/exclusion law.
!
!    D(N) = nint ( N! / E )
!
!  First values:
!
!     N         D(N)
!     0           1
!     1           0
!     2           1
!     3           2
!     4           9
!     5          44
!     6         265
!     7        1854
!     8       14833
!     9      133496
!    10     1334961
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    12 August 2004
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) N, the maximum number of objects to be 
!    permuted.
!
!    Output, integer ( kind = 4 ) DERANGE_ENUM3, the number of derangements 
!    of N objects.
!
  implicit none

  integer ( kind = 4 ) derange_enum3
  real    ( kind = 8 ), parameter :: e = 2.718281828459045D+00
  integer ( kind = 4 ) n
  real    ( kind = 8 ) r8_factorial

  if ( n < 0 ) then
    derange_enum3 = -1
  else if ( n == 0 ) then
    derange_enum3 = 1
  else if ( n == 1 ) then
    derange_enum3 = 0
  else
    derange_enum3 = nint ( r8_factorial ( n ) / e )
  end if

  return
end
subroutine derange_weed_next ( n, a, more )

!*****************************************************************************80
!
!! DERANGE_WEED_NEXT computes all of the derangements of N objects, one at a time.
!
!  Discussion:
!
!    A derangement of N objects is a permutation which leaves no object
!    unchanged.
!
!    A derangement of N objects is a permutation with no fixed
!    points.  If we symbolize the permutation operation by "P",
!    then for a derangment, P(I) is never equal to I.
!
!    The number of derangements of N objects is sometimes called
!    the subfactorial function, or the derangement number D(N).
!
!    This routine simply generates all permutations, one at a time,
!    and weeds out those that are not derangements.
!
!  Examples:
!
!    Here are the derangements when N = 4:
!
!    2143  3142  4123
!    2341  3412  4312
!    2413  3421  4321
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    18 July 2000
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) N, the number of objects being permuted.
!
!    Input/output, integer ( kind = 4 ) A(N).
!    On first call, the input contents of A are unimportant.  But
!    on the second and later calls, the input value of A should be
!    the output value returned on the previous call.
!    On output, A contains the next derangement.
!
!    Input/output, logical MORE.
!    Set MORE = FALSE before the first call.
!    MORE will be reset to TRUE and a derangement will be returned.
!    Each new call produces a new derangement until MORE is returned FALSE.
!
  implicit none

  integer ( kind = 4 ) n

  integer ( kind = 4 ) a(n)
  logical deranged
  integer ( kind = 4 ) derange_enum
  integer ( kind = 4 ), save :: maxder = 0
  logical more
  integer ( kind = 4 ), save :: numder = 0
!
!  Initialization on call with MORE = FALSE.
!
  if ( .not. more ) then

    maxder = derange_enum ( n )
    numder = 0

  end if
!
!  Watch out for cases where there are no derangements.
!
  if ( maxder == 0 ) then
    more = .false.
    return
  end if
!
!  Get the next permutation.
!
  do

    call perm_lex_next ( n, a, more )
!
!  See if it is a derangment.
!
    call derange_check ( n, a, deranged )

    if ( deranged ) then
      exit
    end if

  end do

  numder = numder + 1

  if ( maxder <= numder ) then
    more = .false.
  end if

  return
end
subroutine digit_to_ch ( digit, c )

!*****************************************************************************80
!
!! DIGIT_TO_CH returns the character representation of a decimal digit.
!
!  Example:
!
!    DIGIT   C
!    -----  ---
!      0    '0'
!      1    '1'
!    ...    ...
!      9    '9'
!     17    '*'
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    04 August 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) DIGIT, the digit value between 0 and 9.
!
!    Output, character C, the corresponding character, or '*' if DIGIT
!    was illegal.
!
  implicit none

  character c
  integer ( kind = 4 ) digit

  if ( 0 <= digit .and. digit <= 9 ) then

    c = char ( digit + 48 )

  else

    c = '*'

  end if

  return
end
subroutine digraph_arc_euler ( nnode, nedge, inode, jnode, success, trail )

!*****************************************************************************80
!
!! DIGRAPH_ARC_EULER returns an Euler circuit in a digraph.
!
!  Discussion:
!
!    An Euler circuit of a digraph is a path which starts and ends at
!    the same node and uses each directed edge exactly once.  A digraph is
!    eulerian if it has an Euler circuit.  The problem is to decide whether
!    a given digraph is eulerian and to find an Euler circuit if the
!    answer is affirmative.
!
!
!    A digraph has an Euler circuit if and only if the number of incoming
!    edges is equal to the number of outgoing edges at each node.
!
!    This characterization gives a straightforward procedure to decide whether
!    a digraph is eulerian.  Furthermore, an Euler circuit in an eulerian
!    digraph G of NEDGE edges can be determined by the following method:
!
!      STEP 1: Choose any node U as the starting node, and traverse any edge
!      ( U, V ) incident to node U, and than traverse any unused edge incident
!      to node U.  Repeat this process of traversing unused edges until the
!      starting node U is reached.  Let P be the resulting walk consisting of
!      all used edges.  If all edges of G are in P, than stop.
!
!      STEP 2: Choose any unused edge ( X,  Y) in G such that X is
!      in P and Y is not in P.  Use node X as the starting node and
!      find another walk Q using all unused edges as in step 1.
!
!      STEP 3: Walk P and walk Q share a common node X, they can be merged
!      to form a walk R by starting at any node S of P and to traverse P
!      until node X is reached; than, detour and traverse all edges of Q
!      until node X is reached and continue to traverse the edges of P until
!      the starting node S is reached.  Set P = R.
!
!      STEP 4: Repeat steps 2 and 3 until all edges are used.
!
!    The running time of the algorithm is O ( NEDGE ).
!
!    The digraph is assumed to be connected.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    25 July 2000
!
!  Author:
!
!    FORTRAN77 original version by Hang Tong Lau.
!    FORTRAN90 version by John Burkardt
!
!  Reference:
!
!    Hang Tong Lau,
!    Algorithms on Graphs,
!    Tab Books, 1989.
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) NNODE, the number of nodes.
!
!    Input, integer ( kind = 4 ) NEDGE, the number of edges.
!
!    Input, integer ( kind = 4 ) INODE(NEDGE), JNODE(NEDGE); the I-th edge 
!    starts at node INODE(I) and ends at node JNODE(I).
!
!    Output, logical SUCCESS, is TRUE if an Euler circuit was found,
!    and FALSE otherwise.
!
!    Output, integer ( kind = 4 ) TRAIL(NEDGE).  TRAIL(I) is the edge number 
!    of the I-th edge in the Euler circuit.
!
  implicit none

  integer ( kind = 4 ) nedge

  logical candid(nedge)
  integer ( kind = 4 ) endnod(nedge)
  integer ( kind = 4 ) i
  integer ( kind = 4 ) inode(nedge)
  integer ( kind = 4 ) istak
  integer ( kind = 4 ) j
  integer ( kind = 4 ) jnode(nedge)
  integer ( kind = 4 ) k
  integer ( kind = 4 ) l
  integer ( kind = 4 ) len
  integer ( kind = 4 ) lensol
  integer ( kind = 4 ) lenstk
  integer ( kind = 4 ) nnode
  integer ( kind = 4 ) stack(2*nedge)
  logical success
  integer ( kind = 4 ) trail(nedge)
!
!  Check if the digraph is eulerian.
!
  trail(1:nedge) = 0
  endnod(1:nedge) = 0

  do i = 1, nedge
    j = inode(i)
    trail(j) = trail(j) + 1
    j = jnode(i)
    endnod(j) = endnod(j) + 1
  end do

  do i = 1, nnode
    if ( trail(i) /= endnod(i) ) then
      success = .false.
      return
    end if
  end do
!
!  The digraph is eulerian; find an Euler circuit.
!
  success = .true.
  lensol = 1
  lenstk = 0
!
!  Find the next edge.
!
  do

    if ( lensol == 1 ) then

      endnod(1) = inode(1)
      stack(1) = 1
      stack(2) = 1
      lenstk = 2

    else

      l = lensol - 1

      if ( lensol /= 2 ) then
        endnod(l) = inode(trail(l)) + jnode(trail(l)) - endnod(l-1)
      end if

      k = endnod(l)

      do i = 1, nedge
        candid(i) = ( k == jnode(i) )
      end do

      do i = 1, l
        candid(trail(i)) = .false.
      end do

      len = lenstk

      do i = 1, nedge

        if ( candid(i) ) then
          len = len + 1
          stack(len) = i
        end if

      end do

      stack(len+1) = len - lenstk
      lenstk = len + 1

    end if

    do

      istak = stack(lenstk)
      lenstk = lenstk - 1

      if ( istak /= 0 ) then
        exit
      end if

      lensol = lensol - 1

      if ( lensol == 0 ) then
        call i4vec_reverse ( nedge, trail )
        return
      end if

    end do

    trail(lensol) = stack(lenstk)
    stack(lenstk) = istak - 1

    if ( lensol == nedge ) then
      exit
    end if

    lensol = lensol + 1

  end do

  call i4vec_reverse ( nedge, trail )

  return
end
subroutine digraph_arc_print ( nedge, inode, jnode, title )

!*****************************************************************************80
!
!! DIGRAPH_ARC_PRINT prints out a digraph from an edge list.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    04 July 2000
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) NEDGE, the number of edges.
!
!    Input, integer ( kind = 4 ) INODE(NEDGE), JNODE(NEDGE), the beginning 
!    and end nodes of the edges.
!
!    Input, character ( len = * ) TITLE, a title.
!
  implicit none

  integer ( kind = 4 ) nedge

  integer ( kind = 4 ) i
  integer ( kind = 4 ) inode(nedge)
  integer ( kind = 4 ) jnode(nedge)
  character ( len = * ) title

  if ( len_trim ( title ) /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) trim ( title )
  end if

  write ( *, '(a)' ) ' '

  do i = 1, nedge
    write ( *, '(i8,4x,2i8)' ) i, inode(i), jnode(i)
  end do

  return
end
subroutine diophantine ( a, b, c, ierror, x, y )

!*****************************************************************************80
!
!! DIOPHANTINE solves a Diophantine equation A * X + B * Y = C.
!
!  Discussion:
!
!    Given integers A, B and C, produce X and Y so that
!
!      A * X + B * Y = C.
!
!    In general, the equation is solvable if and only if the
!    greatest common divisor of A and B also divides C.
!
!    A solution (X,Y) of the Diophantine equation also gives the solution
!    X to the congruence equation:
!
!      A * X = C mod ( B ).
!
!    Generally, if there is one nontrivial solution, there are an infinite
!    number of solutions to a Diophantine prob