# include <stdlib.h>
# include <stdio.h>
# include <math.h>
# include <time.h>
# include <string.h>

# include "treepack.h"

/******************************************************************************/

int *catalan ( int n )

/******************************************************************************/
/*
  Purpose:

    CATALAN computes the Catalan numbers, from C(0) to C(N).

  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) the 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 tableaus 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.

  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)

  Example:

    N = 3

    ()()()
    ()(())
    (()())
    (())()
    ((()))

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    08 May 2003

  Author:

    John Burkardt

  Reference:

    Dennis Stanton, Dennis White,
    Constructive Combinatorics,
    Springer, 1986,
    ISBN: 0387963472,
    LC: QA164.S79.

  Parameters:

    Input, int N, the number of Catalan numbers desired.

    Output, int CATALAN[N+1], the Catalan numbers from C(0) to C(N).
*/
{
  int *c;
  int i;

  if ( n < 0 )
  {
    return NULL;
  }

  c = ( int * ) malloc ( ( n + 1 ) * sizeof ( int ) );

  c[0] = 1;
/*
  The extra parentheses ensure that the integer division is
  done AFTER the integer multiplication.
*/
  for ( i = 1; i <= n; i++ )
  {
    c[i] = ( c[i-1] * 2 * ( 2 * i - 1 ) ) / ( i + 1 );
  }

  return c;
}
/******************************************************************************/

void catalan_values ( int *n_data, int *n, int *c )

/******************************************************************************/
/*
  Purpose:

    CATALAN_VALUES returns some values of the Catalan numbers.

  Discussion:

    In Mathematica, the function can be evaluated by:

      Binomial[2*n,n] / ( n + 1 )

  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) the 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 tableaus 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

    ()()()
    ()(())
    (()())
    (())()
    ((()))

  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.

    Stephen Wolfram,
    The Mathematica Book,
    Fourth Edition,
    Cambridge University Press, 1999,
    ISBN: 0-521-64314-7,
    LC: QA76.95.W65.

  Parameters:

    Input/output, int *N_DATA.  The user sets N_DATA to 0 before the
    first call.  On each call, the routine increments N_DATA by 1, and
    returns the corresponding data; when there is no more data, the
    output value of N_DATA will be 0 again.

    Output, int *N, the order of the Catalan number.

    Output, int *C, the value of the Catalan number.
*/
{
# define N_MAX 11

  int c_vec[N_MAX] = { 
    1, 1, 2, 5, 14, 42, 132, 429, 1430, 4862, 16796 };

  int n_vec[N_MAX] = { 
     0,  1,  2,  3,  4, 5,  6,  7,  8,  9,  10 };

  if ( *n_data < 0 )
  {
    *n_data = 0;
  }

  *n_data = *n_data + 1;

  if ( N_MAX < *n_data )
  {
    *n_data = 0;
    *n = 0;
    *c = 0;
  }
  else
  {
    *n = n_vec[*n_data-1];
    *c = c_vec[*n_data-1];
  }

  return;
# undef N_MAX
}
/******************************************************************************/

void cbt_traverse ( int depth )

/******************************************************************************/
/*
  Purpose:

    CBT_TRAVERSE traverses a complete binary tree of given depth.

  Discussion:

    There will be 2^DEPTH terminal nodes of the complete binary tree.

    This function traverses the tree, and prints out a binary code of 0's
    and 1's each time it encounters a terminal node.  This results in a
    printout of the binary digits from 0 to 2^DEPTH - 1.

    The function is intended as a framework to be used to traverse a binary
    tree.  Thus, in practice, a user would insert some action when a terminal
    node is encountered.

    Another use would occur when a combinatorial search is being made, for
    example in a knapsack problem.  Each binary string then represents which
    objects are to be included in the knapsack.  In that case, the traversal
    could be speeded up by noticing cases where a nonterminal node has been
    reached, but the knapsack is already full, in which case the only solution
    uses none of the succeeding items, or overfull, in which case no solutions
    exist that include this initial path segment.

  Licensing:

    This code is distributed under the GNU LGPL license.

  Modified:

    14 July 2013

  Author:

    John Burkardt

  Parameters:

    Input, int DEPTH, the depth of the tree.
*/
{
  int *b;
  int direction;
  int DOWNLEFT = 1;
  int i;
  int k;
  int p;
  int UP = 3;
  int UPDOWNRIGHT = 2;

  if ( depth < 1 )
  {
    return;
  }

  b = ( int * ) malloc ( ( depth + 1 ) * sizeof ( int ) );

  for ( i = 0; i <= depth; i++ )
  {
    b[i] = 0;
  }
  p = 0;
  direction = DOWNLEFT;
  k = 0;

  for ( ; ; )
  {
/*
  Try going in direction DOWNLEFT.
*/
    if ( direction == DOWNLEFT )
    {
      p = p + 1;
      b[p-1] = 0;
      if ( p < depth )
      {
        printf ( "           " );
        for ( i = 0; i < p; i++ )
        {
          printf ( "%d", b[i] );
        }
        printf ( "\n" );
      }
      else
      {
        printf ( "  (  %4d  ", k );
        for ( i = 0; i < p; i++ )
        {
          printf ( "%d", b[i] );
        }
        printf ( "\n" );
        k = k + 1;
        direction = UPDOWNRIGHT;
      }
    }
/*
  Try going in direction UPDOWNRIGHT.
*/
    if ( direction == UPDOWNRIGHT )
    {
      b[p-1] = + 1;
      if ( p < depth )
      {
        printf ( "           " );
        for ( i = 0; i < p; i++ )
        {
          printf ( "%d", b[i] );
        }
        printf ( "\n" );
        direction = DOWNLEFT;
      }
      else
      {
        printf ( "  )  %4d  ", k );
        for ( i = 0; i < p; i++ )
        {
          printf ( "%d", b[i] );
        }
        printf ( "\n" );
        k = k + 1;
        direction = UP;
      }
    }
/*
  Try going in direction UP.
*/
    if ( direction == UP )
    {
      p = p - 1;
      if ( 1 <= p )
      {
        printf ( "           " );
        for ( i = 0; i < p; i++ )
        {
          printf ( "%d", b[i] );
        }
        printf ( "\n" );
        if ( b[p-1] == 0 )
        {
          direction = UPDOWNRIGHT;
        }
      }
      else
      {
        break;
      }
    }
  }

  free ( b );

  return;
}
/******************************************************************************/

int graph_adj_edge_count ( int adj[], int nnode )

/******************************************************************************/
/*
  Purpose:

    GRAPH_ADJ_EDGE_COUNT counts the number of edges in a graph.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    14 July 2013

  Author:

    John Burkardt

  Parameters:

    Input, int ADJ[NNODE*NNODE], the adjacency information.
    ADJ(I,J) is 1 if there is an edge from node I to node J.

    Input, int NNODE, the number of nodes.

    Output, int GRAPH_ADJ_EDGE_COUNT, the number of edges in the graph.
*/
{
  int i;
  int j;
  int nedge;

  nedge = 0;

  for ( i = 0; i < nnode; i++ )
  {
    for ( j = 0; j < nnode; j++ )
    {
      if ( i == j )
      {
        nedge = nedge + 2 * adj[i+j*nnode];
      }
      else
      {
        nedge = nedge + adj[i+j*nnode];
      }
    }
  }

  nedge = nedge / 2;

  return nedge;
}
/******************************************************************************/

int graph_adj_is_node_connected ( int adj[], int nnode )

/******************************************************************************/
/*
  Purpose:

    GRAPH_ADJ_IS_NODE_CONNECTED determines if a graph is nodewise connected.

  Discussion:

    A graph is nodewise connected if, from every node, there is a path
    to any other node.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    14 July 2013

  Author:

    John Burkardt

  Parameters:

    Input, int ADJ[NNODE*NNODE], the adjacency matrix for the 
    graph.  ADJ(I,J) is nonzero if there is an edge from node I to node J.

    Input, int NNODE, the number of nodes.

    Output, int GRAPH_ADJ_IS_NODE_CONNECTED.
    0, the graph is not nodewise connected.
    1, the graph is nodewise connected.
*/
{
  int *found;
  int i;
  int ihi;
  int ii;
  int ilo;
  int j;
  int jhi;
  int jlo;
  int *list;
  int result;
/*
  FOUND(I) is 1 if node I has been reached.
  LIST(I) contains a list of the nodes as they are reached.
*/
  list = ( int * ) malloc ( nnode * sizeof ( int ) );
  found = ( int * ) malloc ( nnode * sizeof ( int ) );

  for ( i = 0; i < nnode; i++ )
  {
    list[i] = 0;
    found[i] = 0;
  }
/*
  Start at node 1.
*/
  found[1-1] = 1;
  list[1-1] = 1;
  ilo = 1;
  ihi = 1;
/*
  From the batch of nodes found last time, LIST(ILO:IHI),
  look for unfound neighbors, and store their indices in LIST(JLO:JHI).
*/
  for ( ; ; )
  {
    jlo = ihi + 1;
    jhi = ihi;

    for ( ii = ilo; ii <= ihi; ii++ )
    {
      i = list[ii-1];

      for ( j = 1; j <= nnode; j++ )
      {
        if ( adj[i-1+(j-1)*nnode] != 0 || adj[j-1+(i-1)*nnode] != 0 )
        {
          if ( found[j-1] == 0 )
          {
            jhi = jhi + 1;
            list[jhi-1] = j;
            found[j-1] = 1;
          }
        }
      }
    }
/*
  If no neighbors were found, exit.
*/
    if ( jhi < jlo )
    {
      break;
    }
/*
  If neighbors were found, then go back and find THEIR neighbors.
*/
    ilo = jlo;
    ihi = jhi; 
  }
/*
  No more neighbors were found.  Have we reached all nodes?
*/
  if ( ihi == nnode )
  {
    result = 1;
  }
  else
  {
    result = 0;
  }

  free ( found );
  free ( list );

  return result;
}
/******************************************************************************/

int graph_adj_is_tree ( int adj[], int nnode )

/******************************************************************************/
/*
  Purpose:

    GRAPH_ADJ_IS_TREE determines whether a graph is a tree.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    14 July 2013

  Author:

    John Burkardt

  Parameters:

    Input, int ADJ[NNODE*NNODE], the adjacency matrix for the 
    graph.  ADJ(I,J) is nonzero if there is an edge from node I to node J.

    Input, int NNODE, the number of nodes.

    Output, int GRAPH_ADJ_IS_TREE.
    0, the graph is not a tree.
    1, the graph is a tree.
*/
{
  int nedge;
  int result;

  if ( nnode <= 1 )
  {
    result = 1;
    return result;
  }
/*
  Every node must be connected to every other node.
*/
  result = graph_adj_is_node_connected ( adj, nnode );

  if ( result == 0 )
  {
    return result;
  }
/*
  There must be exactly NNODE-1 edges.
*/
  nedge = graph_adj_edge_count ( adj, nnode );

  if ( nedge == nnode - 1 )
  {
    result = 1;
  }
  else
  {
    result = 0;
  }

  return result;
}
/******************************************************************************/

int *graph_arc_degree ( int nnode, int nedge, int inode[], int jnode[] )

/******************************************************************************/
/*
  Purpose:

    GRAPH_ARC_DEGREE determines the degree of the nodes of a graph.

  Discussion:

    The degree of a node is the number of edges that include the node.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    14 July 2013

  Author:

    John Burkardt

  Parameters:

    Input, int NNODE, the number of nodes.

    Input, int NEDGE, the number of edges.

    Input, int INODE[NEDGE], JNODE[NEDGE], the pairs of nodes
    that form the edges.

    Output, int GRAPH_ARC_DEGREE[NNODE], the degree of each node,
    that is, the number of edges that include the node.
*/
{
  int *degree;
  int i;
  int n;

  degree = ( int * ) malloc ( nnode * sizeof ( int ) );

  for ( i = 0; i < nnode; i++ )
  {
    degree[i] = 0;
  }

  for ( i = 0; i < nedge; i++ )
  {
    n = inode[i];

    if ( 1 <= n && n <= nnode )
    {
      degree[n-1] = degree[n-1] + 1;
    }
    else
    {
      fprintf ( stderr, "\n" );
      fprintf ( stderr, "GRAPH_ARC_DEGREE - Fatal error!\n" );
      fprintf ( stderr, "  Out-of-range node value = %d\n", n );
      exit ( 1 );
    }

    n = jnode[i];
    if ( 1 <= n && n <= nnode )
    {
      degree[n-1] = degree[n-1] + 1;
    }
    else
    {
      fprintf ( stderr, "\n" );
      fprintf ( stderr, "GRAPH_ARC_DEGREE - Fatal error!\n" );
      fprintf ( stderr, "  Out-of-range node value = %d\n", n );
      exit ( 1 );
    }
  }
  return degree;
}
/******************************************************************************/

int graph_arc_is_tree ( int nedge, int inode[], int jnode[], int nnode )

/******************************************************************************/
/*
  Purpose:

    GRAPH_ARC_IS_TREE determines whether a graph is a tree.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    14 July 2013

  Author:

    John Burkardt

  Parameters:

    Input, int INODE[NEDGE], JNODE[NEDGE].  INODE(I) and 
    JNODE(I) are the start and end nodes of the I-th edge of the graph G.

    Input, int NEDGE, the number of edges in the graph G.

    Input, int NNODE, the number of nodes.

    Output, int GRAPH_ARC_IS_TREE.
    0, the graph is not a tree.
    1, the graph is a tree.
*/
{
  int *adj;
  int result;

  adj = graph_arc_to_graph_adj ( nedge, inode, jnode );

  result = graph_adj_is_tree ( adj, nnode );

  free ( adj );

  return result;
}
/******************************************************************************/

int graph_arc_node_count ( int nedge, int inode[], int jnode[] )

/******************************************************************************/
/*
  Purpose:

    GRAPH_ARC_NODE_COUNT counts the number of nodes in a graph.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    14 July 2013

  Author:

    John Burkardt

  Parameters:

    Input, int NEDGE, the number of edges.

    Input, int INODE[NEDGE], JNODE[NEDGE].  INODE(I) and 
    JNODE(I) are the start and end nodes of the I-th edge.

    Output, int GRAPH_ARC_NODE_COUNT, the number of distinct nodes.
*/
{
  int i;
  int *knode;
  int nnode;
/*
  Copy all the node labels into KNODE,
  sort KNODE,
  count the unique entries.  

  That's NNODE.
*/
  knode = ( int * ) malloc ( 2 * nedge * sizeof ( int ) );

  for ( i = 0; i < nedge; i++ )
  {
    knode[i] = inode[i];
    knode[i+nedge] = jnode[i];
  }

  i4vec_sort_heap_a ( 2*nedge, knode );

  nnode = i4vec_sorted_unique_count ( 2*nedge, knode );

  free ( knode );

  return nnode;
}
/******************************************************************************/

int graph_arc_node_max ( int nedge, int inode[], int jnode[] )

/******************************************************************************/
/*
  Purpose:

    GRAPH_ARC_NODE_MAX determines the maximum node label.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    14 July 2013

  Author:

    John Burkardt

  Parameters:

    Input, int NEDGE, the number of edges.

    Input, int INODE[NEDGE], JNODE[NEDGE].  INODE(I) and 
    JNODE(I) are the start and end nodes of the I-th edge.

    Output, int GRAPH_ARC_NODE_MAX, the maximum node index.
*/
{
  int i;
  int node_max;

  node_max = 0;
  for ( i = 0; i < nedge; i++ )
  {
    if ( node_max < inode[i] )
    {
      node_max = inode[i];
    }
    if ( node_max < jnode[i] )
    {
      node_max = jnode[i];
    }
  }
  return node_max;
}
/******************************************************************************/

void graph_arc_print ( int nedge, int inode[], int jnode[], char *title )

/******************************************************************************/
/*
  Purpose:

    GRAPH_ARC_PRINT prints out a graph from an edge list.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    14 July 2013

  Author:

    John Burkardt

  Parameters:

    Input, int NEDGE, the number of edges.

    Input, int INODE[NEDGE], JNODE[NEDGE], the beginning 
    and end nodes of the edges.

    Input, char *TITLE, a title.
*/
{
  int i;

  printf ( "\n" );
  printf ( "%s\n", title );
  printf ( "\n" );

  for ( i = 0; i < nedge; i++ )
  {
    printf ( "  %6d    %6d  %6d\n", i, inode[i], jnode[i] );
  }

  return;
}
/******************************************************************************/

int *graph_arc_to_graph_adj ( int nedge, int inode[], int jnode[] )

/******************************************************************************/
/*
  Purpose:

    GRAPH_ARC_TO_GRAPH_ADJ converts an arc list graph to an adjacency graph.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    14 July 2013

  Author:

    John Burkardt

  Parameters:

    Input, int NEDGE, the number of edges.

    Input, int INODE[NEDGE], JNODE[NEDGE], the edge array for 
    an undirected graph.  The I-th edge connects nodes INODE(I) and JNODE(I).

    Output, int GRAPH_ARC_TO_GRAPH_ADJ[NNODE*NNODE], the adjacency information.

*/
{
  int *adj;
  int i;
  int j;
  int k;
  int nnode;
/*
  Determine the number of nodes.
*/
  nnode = graph_arc_node_count ( nedge, inode, jnode );

  adj = ( int * ) malloc ( nnode * nnode * sizeof ( int ) );

  for ( j = 0; j < nnode; j++ )
  {
    for ( i = 0; i < nnode; i++ )
    {
      adj[i+j*nnode] = 0;
    }
  }

  for ( k = 0; k < nedge; k++ )
  {
    i = inode[k] - 1;
    j = jnode[k] - 1;
    adj[i+j*nnode] = 1;
    adj[j+i*nnode] = 1;
  }

  return adj;
}
/******************************************************************************/

int i4_max ( int i1, int i2 )

/******************************************************************************/
/*
  Purpose:

    I4_MAX returns the maximum of two I4's.

  Licensing:

    This code is distributed under the GNU LGPL license.

  Modified:

    29 August 2006

  Author:

    John Burkardt

  Parameters:

    Input, int I1, I2, are two integers to be compared.

    Output, int I4_MAX, the larger of I1 and I2.
*/
{
  int value;

  if ( i2 < i1 )
  {
    value = i1;
  }
  else
  {
    value = i2;
  }
  return value;
}
/******************************************************************************/

int i4_min ( int i1, int i2 )

/******************************************************************************/
/*
  Purpose:

    I4_MIN returns the smaller of two I4's.

  Licensing:

    This code is distributed under the GNU LGPL license.

  Modified:

    29 August 2006

  Author:

    John Burkardt

  Parameters:

    Input, int I1, I2, two integers to be compared.

    Output, int I4_MIN, the smaller of I1 and I2.
*/
{
  int value;

  if ( i1 < i2 )
  {
    value = i1;
  }
  else
  {
    value = i2;
  }
  return value;
}
/******************************************************************************/

int i4_power ( int i, int j )

/******************************************************************************/
/*
  Purpose:

    I4_POWER returns the value of I^J.

  Licensing:

    This code is distributed under the GNU LGPL license.

  Modified:

    23 October 2007

  Author:

    John Burkardt

  Parameters:

    Input, int I, J, the base and the power.  J should be nonnegative.

    Output, int I4_POWER, the value of I^J.
*/
{
  int k;
  int value;

  if ( j < 0 )
  {
    if ( i == 1 )
    {
      value = 1;
    }
    else if ( i == 0 )
    {
      fprintf ( stderr, "\n" );
      fprintf ( stderr, "I4_POWER - Fatal error!\n" );
      fprintf ( stderr, "  I^J requested, with I = 0 and J negative.\n" );
      exit ( 1 );
    }
    else
    {
      value = 0;
    }
  }
  else if ( j == 0 )
  {
    if ( i == 0 )
    {
      fprintf ( stderr, "\n" );
      fprintf ( stderr, "I4_POWER - Fatal error!\n" );
      fprintf ( stderr, "  I^J requested, with I = 0 and J = 0.\n" );
      exit ( 1 );
    }
    else
    {
      value = 1;
    }
  }
  else if ( j == 1 )
  {
    value = i;
  }
  else
  {
    value = 1;
    for ( k = 1; k <= j; k++ )
    {
      value = value * i;
    }
  }
  return value;
}
/******************************************************************************/

int i4_uniform_ab ( int a, int b, int *seed )

/******************************************************************************/
/*
  Purpose:

    I4_UNIFORM_AB returns a scaled pseudorandom I4 between A and B.

  Discussion:

    The pseudorandom number should be uniformly distributed
    between A and B.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    24 May 2012

  Author:

    John Burkardt

  Reference:

    Paul Bratley, Bennett Fox, Linus Schrage,
    A Guide to Simulation,
    Second Edition,
    Springer, 1987,
    ISBN: 0387964673,
    LC: QA76.9.C65.B73.

    Bennett Fox,
    Algorithm 647:
    Implementation and Relative Efficiency of Quasirandom
    Sequence Generators,
    ACM Transactions on Mathematical Software,
    Volume 12, Number 4, December 1986, pages 362-376.

    Pierre L'Ecuyer,
    Random Number Generation,
    in Handbook of Simulation,
    edited by Jerry Banks,
    Wiley, 1998,
    ISBN: 0471134031,
    LC: T57.62.H37.

    Peter Lewis, Allen Goodman, James Miller,
    A Pseudo-Random Number Generator for the System/360,
    IBM Systems Journal,
    Volume 8, Number 2, 1969, pages 136-143.

  Parameters:

    Input, int A, B, the limits of the interval.

    Input/output, int *SEED, the "seed" value, which should NOT be 0.
    On output, SEED has been updated.

    Output, int I4_UNIFORM_AB, a number between A and B.
*/
{
  int c;
  int i4_huge = 2147483647;
  int k;
  float r;
  int value;

  if ( *seed == 0 )
  {
    fprintf ( stderr, "\n" );
    fprintf ( stderr, "I4_UNIFORM_AB - Fatal error!\n" );
    fprintf ( stderr, "  Input value of SEED = 0.\n" );
    exit ( 1 );
  }
/*
  Guaranteee A <= B.
*/
  if ( b < a )
  {
    c = a;
    a = b;
    b = c;
  }

  k = *seed / 127773;

  *seed = 16807 * ( *seed - k * 127773 ) - k * 2836;

  if ( *seed < 0 )
  {
    *seed = *seed + i4_huge;
  }

  r = ( float ) ( *seed ) * 4.656612875E-10;
/*
  Scale R to lie between A-0.5 and B+0.5.
*/
  r = ( 1.0 - r ) * ( ( float ) ( a ) - 0.5 ) 
    +         r   * ( ( float ) ( b ) + 0.5 );
/*
  Round R to the nearest integer.
*/
  value = round ( r );
/*
  Guarantee that A <= VALUE <= B.
*/
  if ( value < a )
  {
    value = a;
  }
  if ( b < value )
  {
    value = b;
  }

  return value;
}
/******************************************************************************/

void i4mat_print ( int m, int n, int a[], char *title )

/******************************************************************************/
/*
  Purpose:

    I4MAT_PRINT prints an I4MAT.

  Discussion:

    An I4MAT is an MxN array of I4's, stored by (I,J) -> [I+J*M].

  Licensing:

    This code is distributed under the GNU LGPL license.

  Modified:

    28 May 2008

  Author:

    John Burkardt

  Parameters:

    Input, int M, the number of rows in A.

    Input, int N, the number of columns in A.

    Input, int A[M*N], the M by N matrix.

    Input, char *TITLE, a title.
*/
{
  i4mat_print_some ( m, n, a, 1, 1, m, n, title );

  return;
}
/******************************************************************************/

void i4mat_print_some ( int m, int n, int a[], int ilo, int jlo, int ihi,
  int jhi, char *title )

/******************************************************************************/
/*
  Purpose:

    I4MAT_PRINT_SOME prints some of an I4MAT.

  Discussion:

    An I4MAT is an MxN array of I4's, stored by (I,J) -> [I+J*M].

  Licensing:

    This code is distributed under the GNU LGPL license.

  Modified:

    20 August 2010

  Author:

    John Burkardt

  Parameters:

    Input, int M, the number of rows of the matrix.
    M must be positive.

    Input, int N, the number of columns of the matrix.
    N must be positive.

    Input, int A[M*N], the matrix.

    Input, int ILO, JLO, IHI, JHI, designate the first row and
    column, and the last row and column to be printed.

    Input, char *TITLE, a title.
*/
{
# define INCX 10

  int i;
  int i2hi;
  int i2lo;
  int j;
  int j2hi;
  int j2lo;

  fprintf ( stdout, "\n" );
  fprintf ( stdout, "%s\n", title );

  if ( m <= 0 || n <= 0 )
  {
    fprintf ( stdout, "\n" );
    fprintf ( stdout, "  (None)\n" );
    return;
  }
/*
  Print the columns of the matrix, in strips of INCX.
*/
  for ( j2lo = jlo; j2lo <= jhi; j2lo = j2lo + INCX )
  {
    j2hi = j2lo + INCX - 1;
    j2hi = i4_min ( j2hi, n );
    j2hi = i4_min ( j2hi, jhi );

    fprintf ( stdout, "\n" );
/*
  For each column J in the current range...

  Write the header.
*/
    fprintf ( stdout, "  Col:" );
    for ( j = j2lo; j <= j2hi; j++ )
    {
      fprintf ( stdout, "  %6d", j - 1 );
    }
    fprintf ( stdout, "\n" );
    fprintf ( stdout, "  Row\n" );
    fprintf ( stdout, "\n" );
/*
  Determine the range of the rows in this strip.
*/
    i2lo = i4_max ( ilo, 1 );
    i2hi = i4_min ( ihi, m );

    for ( i = i2lo; i <= i2hi; i++ )
    {
/*
  Print out (up to INCX) entries in row I, that lie in the current strip.
*/
      fprintf ( stdout, "%5d:", i - 1 );
      for ( j = j2lo; j <= j2hi; j++ )
      {
        fprintf ( stdout, "  %6d", a[i-1+(j-1)*m] );
      }
      fprintf ( stdout, "\n" );
    }
  }

  return;
# undef INCX
}
/******************************************************************************/

void i4vec_heap_d ( int n, int a[] )

/******************************************************************************/
/*
  Purpose:

    I4VEC_HEAP_D reorders an I4VEC into a descending heap.

  Discussion:

    An I4VEC is a vector of I4's.

    A heap is an array A with the property that, for every index J,
    A[J] >= A[2*J+1] and A[J] >= A[2*J+2], (as long as the indices
    2*J+1 and 2*J+2 are legal).

  Diagram:

                  A(0)
                /      \
            A(1)         A(2)
          /     \        /  \
      A(3)       A(4)  A(5) A(6)
      /  \       /   \
    A(7) A(8)  A(9) A(10)

  Licensing:

    This code is distributed under the GNU LGPL license.

  Modified:

    30 April 1999

  Author:

    John Burkardt

  Reference:

    Albert Nijenhuis, Herbert Wilf,
    Combinatorial Algorithms,
    Academic Press, 1978, second edition,
    ISBN 0-12-519260-6.

  Parameters:

    Input, int N, the size of the input array.

    Input/output, int A[N].
    On input, an unsorted array.
    On output, the array has been reordered into a heap.
*/
{
  int i;
  int ifree;
  int key;
  int m;
/*
  Only nodes (N/2)-1 down to 0 can be "parent" nodes.
*/
  for ( i = ( n / 2 ) - 1; 0 <= i; i-- )
  {
/*
  Copy the value out of the parent node.
  Position IFREE is now "open".
*/
    key = a[i];
    ifree = i;

    for ( ; ; )
    {
/*
  Positions 2*IFREE + 1 and 2*IFREE + 2 are the descendants of position
  IFREE.  (One or both may not exist because they equal or exceed N.)
*/
      m = 2 * ifree + 1;
/*
  Does the first position exist?
*/
      if ( n <= m )
      {
        break;
      }
      else
      {
/*
  Does the second position exist?
*/
        if ( m + 1 < n )
        {
/*
  If both positions exist, take the larger of the two values,
  and update M if necessary.
*/
          if ( a[m] < a[m+1] )
          {
            m = m + 1;
          }
        }
/*
  If the large descendant is larger than KEY, move it up,
  and update IFREE, the location of the free position, and
  consider the descendants of THIS position.
*/
        if ( key < a[m] )
        {
          a[ifree] = a[m];
          ifree = m;
        }
        else
        {
          break;
        }
      }
    }
/*
  When you have stopped shifting items up, return the item you
  pulled out back to the heap.
*/
    a[ifree] = key;
  }

  return;
}
/******************************************************************************/

int *i4vec_indicator_new ( int n )

/******************************************************************************/
/*
  Purpose:

    I4VEC_INDICATOR_NEW sets an I4VEC to the indicator vector.

  Discussion:

    An I4VEC is a vector of I4's.

  Licensing:

    This code is distributed under the GNU LGPL license.

  Modified:

    26 August 2008

  Author:

    John Burkardt

  Parameters:

    Input, int N, the number of elements of A.

    Output, int I4VEC_INDICATOR_NEW[N], the array.
*/
{
  int *a;
  int i;

  a = ( int * ) malloc ( n * sizeof ( int ) );

  for ( i = 0; i < n; i++ )
  {
    a[i] = i + 1;
  }
  return a;
}
/******************************************************************************/

int i4vec_max ( int n, int a[] )

/******************************************************************************/
/*
  Purpose:

    I4VEC_MAX returns the value of the maximum element in an I4VEC.

  Discussion:

    An I4VEC is a vector of I4's.

  Licensing:

    This code is distributed under the GNU LGPL license.

  Modified:

    17 May 2003

  Author:

    John Burkardt

  Parameters:

    Input, int N, the number of entries in the array.

    Input, int A[N], the array to be checked.

    Output, int IVEC_MAX, the value of the maximum element.  This
    is set to 0 if N <= 0.
*/
{
  int i;
  int value;

  if ( n <= 0 )
  {
    return 0;
  }

  value = a[0];

  for ( i = 1; i < n; i++ )
  {
    if ( value < a[i] )
    {
      value = a[i];
    }
  }

  return value;
}
/******************************************************************************/

void i4vec_print ( int n, int a[], char *title )

/******************************************************************************/
/*
  Purpose:

    I4VEC_PRINT prints an I4VEC.

  Discussion:

    An I4VEC is a vector of I4's.

  Licensing:

    This code is distributed under the GNU LGPL license.

  Modified:

    14 November 2003

  Author:

    John Burkardt

  Parameters:

    Input, int N, the number of components of the vector.

    Input, int A[N], the vector to be printed.

    Input, char *TITLE, a title.
*/
{
  int i;

  fprintf ( stdout, "\n" );
  fprintf ( stdout, "%s\n", title );
  fprintf ( stdout, "\n" );

  for ( i = 0; i < n; i++ )
  {
    fprintf ( stdout, "  %6d: %8d\n", i, a[i] );
  }
  return;
}
/******************************************************************************/

void i4vec_sort_heap_a ( int n, int a[] )

/******************************************************************************/
/*
  Purpose:

    I4VEC_SORT_HEAP_A ascending sorts an I4VEC using heap sort.

  Discussion:

    An I4VEC is a vector of I4's.

  Licensing:

    This code is distributed under the GNU LGPL license.

  Modified:

    30 April 1999

  Author:

    John Burkardt

  Reference:

    Albert Nijenhuis, Herbert Wilf,
    Combinatorial Algorithms,
    Academic Press, 1978, second edition,
    ISBN 0-12-519260-6.

  Parameters:

    Input, int N, the number of entries in the array.

    Input/output, int A[N].
    On input, the array to be sorted;
    On output, the array has been sorted.
*/
{
  int n1;
  int temp;

  if ( n <= 1 )
  {
    return;
  }
/*
  1: Put A into descending heap form.
*/
  i4vec_heap_d ( n, a );
/*
  2: Sort A.

  The largest object in the heap is in A[0].
  Move it to position A[N-1].
*/
  temp = a[0];
  a[0] = a[n-1];
  a[n-1] = temp;
/*
  Consider the diminished heap of size N1.
*/
  for ( n1 = n-1; 2 <= n1; n1-- )
  {
/*
  Restore the heap structure of the initial N1 entries of A.
*/
    i4vec_heap_d ( n1, a );
/*
  Take the largest object from A[0] and move it to A[N1-1].
*/
    temp = a[0];
    a[0] = a[n1-1];
    a[n1-1] = temp;
  }

  return;
}
/******************************************************************************/

int i4vec_sorted_unique_count ( int n, int a[] )

/******************************************************************************/
/*
  Purpose:

    I4VEC_SORTED_UNIQUE_COUNT counts unique elements in a sorted I4VEC.

  Discussion:

    An I4VEC is a vector of I4's.

    Because the array is sorted, this algorithm is O(N).

  Licensing:

    This code is distributed under the GNU LGPL license.

  Modified:

    08 June 2010

  Author:

    John Burkardt

  Parameters:

    Input, int N, the number of elements of A.

    Input, int A[N], the sorted array to examine.

    Output, int I4VEC_SORTED_UNIQUE_COUNT, the number of unique elements of A.
*/
{
  int i;
  int unique_num;

  unique_num = 0;

  if ( n < 1 )
  {
    return unique_num;
  }

  unique_num = 1;

  for ( i = 1; i < n; i++ )
  {
    if ( a[i-1] != a[i] )
    {
      unique_num = unique_num + 1;
    }
  }

  return unique_num;
}
/******************************************************************************/

void pruefer_to_tree_arc ( int nnode, int iarray[], int inode[], int jnode[] )

/******************************************************************************/
/*
  Purpose:

    PRUEFER_TO_TREE_ARC is given a Pruefer code, and computes the tree.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    15 July 2013

  Reference:

    Dennis Stanton, Dennis White,
    Constructive Combinatorics,
    Springer Verlag, New York, 1986.

  Parameters:

    Input, int NNODE, the number of nodes.

    Input, int IARRAY[NNODE-2], the Pruefer code of the tree.

    Output, int INODE[NNODE-1], JNODE[NNODE-1], the edge
    array of the tree.  The I-th edge joins nodes INODE(I) and JNODE(I).
*/
{
  int i;
  int ii;
  int *iwork;
  int j;
/*
  Initialize IWORK(I) to count the number of neighbors of node I.
  The Pruefer code uses each node one less time than its total
  number of neighbors.
*/
  iwork = ( int * ) malloc ( nnode * sizeof ( int ) );

  for ( i = 0; i < nnode; i++ )
  {
    iwork[i] = 1;
  }
  for ( i = 0; i < nnode - 2; i++ )
  {
    iwork[iarray[i]-1] = iwork[iarray[i]-1] + 1;
  }

  for ( i = 0; i < nnode - 1; i++ )
  {
    inode[i] = -1;
    jnode[i] = -1;
  }
/*
  Now process each entry in the Pruefer code.
*/
  for ( i = 0; i < nnode - 2; i++ )
  {
    ii = -1;
    for ( j = 0; j < nnode; j++ )
    {
      if ( iwork[j] == 1 )
      {
        ii = j;
      }
    }
    inode[i] = ii + 1;
    jnode[i] = iarray[i];
    iwork[ii] = 0;
    iwork[iarray[i]-1] = iwork[iarray[i]-1] - 1;
  }

  inode[nnode-2] = iarray[nnode-3];
 
  if ( iarray[nnode-3] != 1 )
  {
    jnode[nnode-2] = 1;
  }
  else
  {
    jnode[nnode-2] = 2;
  }

  free ( iwork );

  return;
}
/******************************************************************************/

void pruefer_to_tree_2 ( int nnode, int iarray[], int itree[] )

/******************************************************************************/
/*
  Purpose:

    PRUEFER_TO_TREE_2 produces the edge list of a tree from its Pruefer code.

  Discussion:

    One can thus exhibit all trees on N nodes, produce
    one at random, find the M-th one on the list, etc, by
    manipulating the Pruefer codes.

    For every labeled tree on N nodes, there is a unique N-2 tuple
    of integers A1 through AN-2, with each A between 1 and N.  There
    are N^(N-2) such sequences, and each one is associated with exactly
    one tree.

    This routine apparently assumes that the Pruefer code is
    generated by taking the LOWEST labeled terminal node each time.
    This is not consistent with PRUEFER_TO_TREE and TREE_TO_PRUEFER.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    27 June 2013

  Author:

    Original FORTRAN77 version by Albert Nijenhuis, Herbert Wilf.
    C version by John Burkardt.

  Reference:

    Albert Nijenhuis. Herbert Wilf,
    Combinatorial Algorithms,
    Academic Press, 1978, second edition,
    ISBN 0-12-519260-6.

  Parameters:

    Input, int NNODE, number of nodes in desired tree.

    Input, int IARRAY[NNODE].  IARRAY(I), I = 1, NNODE-2 
    is the Pruefer code for the tree.

    Output, int ITREE[NNODE]; the I-th edge of the tree
    joins nodes I and ITREE(I).
*/
{
  int i;
  int ir;
  int j;
  int k;
  int kp;
  int l;

  for ( i = 0; i < nnode; i++ )
  {
    itree[i] = 0;
  }
 
  for ( i = nnode - 2; 1 <= i; i-- )
  {
    l = iarray[i-1];
 
    if ( itree[l-1] == 0 )
    {
      iarray[i-1] = - l;
      itree[l-1] = - 1;
    }
  }
  iarray[nnode-2] = nnode;
/*
  Find next index K so that ITREE(K) is 0.
*/
  k = 1;

 
  while ( itree[k-1] != 0 )
  {
    k = k + 1;
  }
 
  j = 0;
  kp = k;
 
  for ( ; ; )
   {
    j = j + 1;
    ir = abs ( iarray[j-1] );
    itree[kp-1] = ir;
 
    if ( j == nnode - 1 )
    {
      break;
    }
 
    if ( 0 < iarray[j-1] )
    {
      while ( itree[k-1] != 0 )
      {
        k = k + 1;
      }
      kp = k;
      continue;
    }
 
    if ( k < ir )
    {
      itree[ir-1] = 0;
      while ( itree[k-1] != 0 )
      {
        k = k + 1;
      }
      kp = k;
      continue;
    } 
    kp = ir;
  }
/*
  Restore the signs of IARRAY.
*/
  for ( i = 0; i < nnode - 2; i++ )
  {
    iarray[i] = abs ( iarray[i] );
  }
 
  return;
}
/******************************************************************************/

int *pruefer_to_tree_2_new ( int nnode, int iarray[] )

/******************************************************************************/
/*
  Purpose:

    PRUEFER_TO_TREE_2_NEW produces the edge list of a tree from its Pruefer code.

  Discussion:

    One can thus exhibit all trees on N nodes, produce
    one at random, find the M-th one on the list, etc, by
    manipulating the Pruefer codes.

    For every labeled tree on N nodes, there is a unique N-2 tuple
    of integers A1 through AN-2, with each A between 1 and N.  There
    are N^(N-2) such sequences, and each one is associated with exactly
    one tree.

    This routine apparently assumes that the Pruefer code is
    generated by taking the LOWEST labeled terminal node each time.
    This is not consistent with PRUEFER_TO_TREE and TREE_TO_PRUEFER.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    27 June 2013

  Author:

    Original FORTRAN77 version by Albert Nijenhuis, Herbert Wilf.
    C version by John Burkardt.

  Reference:

    Albert Nijenhuis. Herbert Wilf,
    Combinatorial Algorithms,
    Academic Press, 1978, second edition,
    ISBN 0-12-519260-6.

  Parameters:

    Input, int NNODE, number of nodes in desired tree.

    Input, int IARRAY[NNODE].  IARRAY(I), I = 1, NNODE-2 
    is the Pruefer code for the tree.

    Output, int PRUEFER_TO_TREE_2_NEW[NNODE]; the I-th edge of the tree
    joins nodes I and ITREE(I).
*/
{
  int *itree;

  itree = ( int * ) malloc ( nnode * sizeof ( int ) );

  pruefer_to_tree_2 ( nnode, iarray, itree );
 
  return itree;
}
/******************************************************************************/

double r8_uniform_01 ( int *seed )

/******************************************************************************/
/*
  Purpose:

    R8_UNIFORM_01 returns a pseudorandom R8 scaled to [0,1].

  Discussion:

    This routine implements the recursion

      seed = 16807 * seed mod ( 2^31 - 1 )
      r8_uniform_01 = seed / ( 2^31 - 1 )

    The integer arithmetic never requires more than 32 bits,
    including a sign bit.

    If the initial seed is 12345, then the first three computations are

      Input     Output      R8_UNIFORM_01
      SEED      SEED

         12345   207482415  0.096616
     207482415  1790989824  0.833995
    1790989824  2035175616  0.947702

  Licensing:

    This code is distributed under the GNU LGPL license.

  Modified:

    11 August 2004

  Author:

    John Burkardt

  Reference:

    Paul Bratley, Bennett Fox, Linus Schrage,
    A Guide to Simulation,
    Springer Verlag, pages 201-202, 1983.

    Pierre L'Ecuyer,
    Random Number Generation,
    in Handbook of Simulation
    edited by Jerry Banks,
    Wiley Interscience, page 95, 1998.

    Bennett Fox,
    Algorithm 647:
    Implementation and Relative Efficiency of Quasirandom
    Sequence Generators,
    ACM Transactions on Mathematical Software,
    Volume 12, Number 4, pages 362-376, 1986.

    P A Lewis, A S Goodman, J M Miller,
    A Pseudo-Random Number Generator for the System/360,
    IBM Systems Journal,
    Volume 8, pages 136-143, 1969.

  Parameters:

    Input/output, int *SEED, the "seed" value.  Normally, this
    value should not be 0.  On output, SEED has been updated.

    Output, double R8_UNIFORM_01, a new pseudorandom variate, strictly between
    0 and 1.
*/
{
  int i4_huge = 2147483647;
  int k;
  double r;

  k = *seed / 127773;

  *seed = 16807 * ( *seed - k * 127773 ) - k * 2836;

  if ( *seed < 0 )
  {
    *seed = *seed + i4_huge;
  }

  r = ( ( double ) ( *seed ) ) * 4.656612875E-10;

  return r;
}
/******************************************************************************/

void timestamp ( void )

/******************************************************************************/
/*
  Purpose:

    TIMESTAMP prints the current YMDHMS date as a time stamp.

  Example:

    31 May 2001 09:45:54 AM

  Licensing:

    This code is distributed under the GNU LGPL license.

  Modified:

    24 September 2003

  Author:

    John Burkardt

  Parameters:

    None
*/
{
# define TIME_SIZE 40

  static char time_buffer[TIME_SIZE];
  const struct tm *tm;
  time_t now;

  now = time ( NULL );
  tm = localtime ( &now );

  strftime ( time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm );

  fprintf ( stdout, "%s\n", time_buffer );

  return;
# undef TIME_SIZE
}
/******************************************************************************/

void tree_arc_center ( int nnode, int inode[], int jnode[], int center[], 
  int *eccent, int *parity )

/******************************************************************************/
/*
  Purpose:

    TREE_ARC_CENTER computes the center, eccentricity, and parity of a tree.

  Discussion:

    A tree is an undirected graph of N nodes, which uses N-1 edges,
    and is connected.  

    A graph with N-1 edges is not guaranteed to be a tree, and so this
    routine must first check that condition before proceeding.

    The edge distance between two nodes I and J is the minimum number of
    edges that must be traversed in a path from I and J.

    The eccentricity of a node I is the maximum edge distance between
    node I and the other nodes J in the graph.

    The radius of a graph is the minimum eccentricity over all nodes
    in the graph.

    The diameter of a graph is the maximum eccentricity over all nodes
    in the graph.

    The center of a graph is the set of nodes whose eccentricity is 
    equal to the radius, that is, the set of nodes of minimum eccentricity.

    For a tree, the center is either a single node, or a pair of
    neighbor nodes.

    The parity of the tree is 1 if the center is a single node, or 2 if
    the center is 2 nodes.

    The center of a tree can be found by removing all "leaves", that is,
    nodes of degree 1.  This step is repeated until only 1 or 2 nodes
    are left.

    Thanks to Alexander Sax for pointing out that a previous version of the
    code was failing when the tree had an odd parity, that is, a single
    center node, 15 April 2013.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    16 July 2013

  Author:

    John Burkardt

  Parameters:

    Input, int NNODE, the number of nodes.

    Input, int INODE[NNODE-1], JNODE[NNODE-1], the edges of
    the tree.  Edge I connects nodes INODE(I) and JNODE(I).

    Output, int CENTER[2].  CENTER(1) is the index of the
    first node in the center.  CENTER(2) is 0 if there is only one node
    in the center, or else the index of the second node.

    Output, int *ECCENT, the eccentricity of the nodes in 
    the center, and the radius of the the tree.

    Output, int *PARITY, the parity of the tree, which is
    normally 1 or 2.
*/
{
  int *degree;
  int i;
  int iedge;
  int ileaf;
  int j;
  int *list;
  int nedge;
  int nleaf;
  int nnode2;
  int result;

  *eccent = 0;
  center[0] = 0;
  center[1] = 0;
  *parity = 0;

  if ( nnode <= 0 )
  {
    fprintf ( stderr, "\n" );
    fprintf ( stderr, "TREE_ARC_CENTER - Fatal error!\n" );
    fprintf ( stderr, "  NNODE <= 0.\n" );
    exit ( 1 );
  }
  else if ( nnode == 1 )
  {
    *eccent = 0;
    center[0] = 1;
    center[1] = 0;
    *parity = 1;
    return;
  }
  else if ( nnode == 2 )
  {
    *eccent = 1;
    center[0] = 1;
    center[1] = 2;
    *parity = 2;
    return;
  }
/*
  Is this graph really a tree?
*/
  nedge = nnode - 1;
  result = graph_arc_is_tree ( nedge, inode, jnode, nnode );

  if ( result == 0 )
  {
    fprintf ( stderr, "\n" );
    fprintf ( stderr, "TREE_ARC_CENTER - Fatal error!\n" );
    fprintf ( stderr, "  This graph is NOT a tree.\n" );
    exit ( 1 );
  }
/*
  Compute the degrees.
*/
  degree = graph_arc_degree ( nnode, nedge, inode, jnode );
/*
  Defoliate the tree.
*/
  nnode2 = nnode;
  list = ( int * ) malloc ( nnode * sizeof ( int ) );

  for ( ; ; )
  {
    *eccent = *eccent + 1;
/*
  Find and mark the leaves.
*/
    nleaf = 0;

    for ( i = 1; i <= nnode; i++ )
    {
      if ( degree[i-1] == 1 )
      {
        nleaf = nleaf + 1;
        list[nleaf-1] = i;
      }
    }
/*
  Delete the leaves.
*/
    for ( ileaf = 1; ileaf <= nleaf; ileaf++ )
    {
      i = list[ileaf-1];

      iedge = 0;
      j = 0;

      for ( ; ; )
      {
        iedge = iedge + 1;

        if ( nedge < iedge )
        {
          fprintf ( stderr, "\n" );
          fprintf ( stderr, "TREE_ARC_CENTER - Fatal error!\n" );
          fprintf ( stderr, "  Data or algorithm failure.\n" );
          exit ( 1 );
        }

        if ( inode[iedge-1] == i )
        {
          j = jnode[iedge-1];
          inode[iedge-1] = - inode[iedge-1];
          jnode[iedge-1] = - jnode[iedge-1];
        }
        else if ( jnode[iedge-1] == i )
        {
          j = inode[iedge-1];
          inode[iedge-1] = - inode[iedge-1];
          jnode[iedge-1] = - jnode[iedge-1];
        }

        if ( j != 0 )
        {
          break;
        }
      }

      degree[i-1] = -1;
      nnode2 = nnode2 - 1;
      degree[j-1] = degree[j-1] - 1;
/*
  If the other node has degree 0, we must have just finished
  stripping all leaves from the tree, leaving a single node.
  Don't kill it here.  It is our odd center.

     if ( degree(j) == 0 )
     {
       nnode2 = nnode2 - 1
     }
*/
    }
/*
  Find the remaining nodes.
*/
    nnode2 = 0;

    for ( i = 1; i <= nnode; i++ )
    {
      if ( 0 <= degree[i-1] )
      {
        nnode2 = nnode2 + 1;
        list[nnode2-1] = i;
      }
    }
/*
  If at least 3, more pruning is required.
*/
    if ( nnode2 < 3 )
    {
      break;
    }
  }
/*
  If only one or two nodes left, we are done.
*/
  *parity = nnode2;

  for ( i = 0; i < nnode2; i++ )
  {
    center[i] = list[i];
  }
  for ( i = 0; i < nedge; i++ )
  {
    inode[i] = abs ( inode[i] );
    jnode[i] = abs ( jnode[i] );
  }

  free ( list );

  return;
}
/******************************************************************************/

void tree_arc_diam ( int nnode, int inode[], int jnode[], int *diam, 
  int label[], int *n1, int *n2 )

/******************************************************************************/
/*
  Purpose:

    TREE_ARC_DIAM computes the "diameter" of a tree.

  Discussion:

    A tree is an undirected graph of N nodes, which uses N-1 edges,
    and is connected.  

    A graph with N-1 edges is not guaranteed to be a tree, and so this
    routine must first check that condition before proceeding.

    The diameter of a graph is the length of the longest possible
    path that never repeats an edge.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    19 July 2013

  Author:

    John Burkardt

  Parameters:

    Input, int NNODE, the number of nodes.

    Input, int INODE[NNODE-1], JNODE[NNODE-1], the edges 
    of the tree.  Edge I connects nodes INODE(I) and JNODE(I).

    Output, int *DIAM, the length of the longest path 
    in the tree.

    Output, int LABEL[NNODE], marks the path between 
    nodes N1 and N2.  Node I is in this path if LABEL(I) is 1.

    Output, int *N1, *N2, the indices of two nodes in the 
    tree which are separated by DIAM edges.
*/
{
  int *degree;
  int i;
  int invals;
  int j;
  int k;
  int kstep;
  int nabe;
  int nedge;
  int result;

  if ( nnode <= 0 )
  {
    *diam = 0;
    fprintf ( stderr, "\n" );
    fprintf ( stderr, "TREE_ARC_DIAM - Fatal error!\n" );
    fprintf ( stderr, "  NNODE <= 0.\n" );
    exit ( 1 );
  }

  if ( nnode == 1 )
  {
    *diam = 0;
    *n1 = 1;
    *n2 = 1;
    return;
  }
/*
  Is this graph really a tree?
*/
  nedge = nnode - 1;
  result = graph_arc_is_tree ( nedge, inode, jnode, nnode );

  if ( result == 0 )
  {
    fprintf ( stderr, "\n" );
    fprintf ( stderr, "TREE_ARC_DIAM - Fatal error!\n" );
    fprintf ( stderr, "  This graph is NOT a tree.\n" );
    exit ( 1 );
  }

  label = i4vec_indicator_new ( nnode );
/*
  On step K:

    Identify the terminal and interior nodes.

    If there are no interior nodes left, 

      then there are just two nodes left at all.  The diameter is 2*K-1, 
      and a maximal path extends between the nodes whose labels are 
      contained in the two remaining terminal nodes.

    Else

      The label of each terminal node is passed to its interior neighbor.
      If more than one label arrives, take any one.

      The terminal nodes are removed.
*/
  kstep = 0;
  degree = ( int * ) malloc ( nnode * sizeof ( int ) );

  for ( ; ; )
  {
    kstep = kstep + 1;
/*
  Compute the degree of each node.
*/
    for ( j = 1; j <= nnode; j++ )
    {
      degree[j-1] = 0;   
    }
    for ( j = 1; j <= nedge; j++ )
    {
      k = inode[j-1];
      if ( 0 < k )
      {
        degree[k-1] = degree[k-1] + 1;
      }

      k = jnode[j-1];
      if ( 0 < k )
      {
        degree[k-1] = degree[k-1] + 1;
      }
    }
/*
  Count the number of interior nodes.
*/
    invals = 0;
    for ( i = 1; i <= nnode; i++ )
    {
      if ( 1 < degree[i-1] )
      {
        invals = invals + 1;
      }
    }
/*
  If there are 1 or 0 interior nodes, it's time to stop.
*/
    if ( invals == 1 )
    {
      *diam = 2 * kstep;
      break;
    }
    else if ( invals == 0 )
    {
      *diam = 2 * kstep - 1;
      break;
    }
/*
  If there are at least two interior nodes, then chop off the 
  terminal nodes and pass their labels inward.
*/
    for ( k = 1; k <= nnode; k++ )
    {
      if ( degree[k-1] == 1 )
      {
        for ( j = 1; j <= nedge; j++ )
        {
          if ( inode[j-1] == k )
          {
            nabe = jnode[j-1];
            label[nabe-1] = label[k-1];
            inode[j-1] = - inode[j-1];
            jnode[j-1] = - jnode[j-1];
          }
          else if ( jnode[j-1] == k )
          {
            nabe = inode[j-1];
            label[nabe-1] = label[k-1];
            inode[j-1] = - inode[j-1];
            jnode[j-1] = - jnode[j-1];
          }
        }
      }
    }
  }
/*
  Now get the labels from two of the remaining terminal nodes.
  The nodes represented by these labels will be a diameter apart.
*/
  *n1 = 0;
  *n2 = 0;

  for ( i = 1; i <= nnode; i++ );
  {
    if ( degree[i-1] == 1 )
    {
      if ( *n1 == 0 )
      {
        *n1 = label[i-1];
      }
      else if ( *n2 == 0 )
      {
        *n2 = label[i-1];
      }
    }
  }
/*
  Set the labels of the interior node (if any) and nodes marked
  N1 and N2 to 1, and all others to 0.  This will label the nodes on the path.
*/
  if ( invals == 1 )
  {
    for ( i = 1; i <= nnode; i++ )
    {
      if ( 1 < degree[i-1] )
      {
        label[i-1] = 1;
      }
    }
  }

  for ( i = 1; i <= nnode; i++ )
  {
    if ( label[i-1] == *n1 || label[i-1] == *n2 )
    {
      label[i-1] = 1;
    }
    else
    {
      label[i-1] = 0;
    }
  }
/*
  Clean up the arrays.
*/
  for ( j = 1; j <= nedge; j++ )
  {
    inode[j-1] = abs ( inode[j-1] );
    k = inode[j-1];
    jnode[j-1] = abs ( jnode[j-1] );
    k = jnode[j-1];
  }

  free ( degree );

  return;
}
/******************************************************************************/

void tree_arc_random ( int nnode, int *seed, int code[], int inode[], 
  int jnode[] )

/******************************************************************************/
/*
  Purpose:

    TREE_ARC_RANDOM selects a random labeled tree and its Pruefer code.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    16 July 2013

  Author:

    John Burkardt

  Parameters:

    Input, int NNODE, the number of nodes.

    Input/output, int *SEED, a seed for the random number 
    generator.

    Output, int CODE[NNODE-2], the Pruefer code for the 
    labeled tree.

    Output, int INODE[NNODE-1], JNODE[NNODE-1], the edge 
    array for the tree.
*/
{
  int i;

  if ( nnode <= 2 )
  {
    return;
  }

  vec_random ( nnode-2, nnode, seed, code );
 
  for ( i = 0; i < nnode - 2; i++ )
  {
    code[i] = code[i] + 1;
  }
  pruefer_to_tree_arc ( nnode, code, inode, jnode );
 
  return;
}
/******************************************************************************/

int *tree_arc_to_pruefer ( int nnode, int inode[], int jnode[] )

/******************************************************************************/
/*
  Purpose:

    TREE_ARC_TO_PRUEFER is given a labeled tree, and computes its Pruefer code.

  Discussion:

    A tree is an undirected graph of N nodes, which uses N-1 edges,
    and is connected.  

    A graph with N-1 edges is not guaranteed to be a tree, and so this
    routine must first check that condition before proceeding.

    The Pruefer code is a correspondence between all labeled trees of
    N nodes, and all list of N-2 integers between 1 and N (with repetition
    allowed).  The number of labeled trees on N nodes is therefore N^(N-2).

    The Pruefer code is constructed from the tree as follows:

    A terminal node on the tree is defined as a node with only one neighbor.

    Consider the set of all terminal nodes on the tree.  Take the one
    with the highest label, I.  Record the label of its neighbor, J.
    Delete node I and the edge between node I and J.

    J is the first entry in the Pruefer code for the tree.   Repeat
    the operation a total of N-2 times to get the complete code.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    18 July 2013

  Reference:

    Dennis Stanton, Dennis White,
    Constructive Combinatorics,
    Springer Verlage, New York, 1986.

  Parameters:

    Input, int NNODE, the number of nodes.

    Input, int INODE[NNODE-1], JNODE[NNODE-1], the edge array 
    of the tree.  The I-th edge joins nodes INODE(I) and JNODE(I).

    Output, int TREE_ARC_TO_PRUEFER[NNODE-2], the Pruefer code of the tree.
*/
{
  int *code;
  int *degree;
  int i;
  int i2;
  int iterm;
  int j;
  int jsave;
  int nedge;
  int result;
/*
  Is this graph really a tree?
*/
  nedge = nnode - 1;
  result = graph_arc_is_tree ( nedge, inode, jnode, nnode );

  if ( result == 0 )
  {
    fprintf ( stderr, "\n" );
    fprintf ( stderr, "TREE_ARC_TO_PRUEFER - Fatal error!\n" );
    fprintf ( stderr, "  This graph is NOT a tree.\n" );
    exit ( 1 );
  }

  code = ( int * ) malloc ( ( nnode - 2 ) * sizeof ( int ) );
/*
  Compute the degree of each node.
*/
  nedge = nnode - 1;
  degree =  graph_arc_degree ( nnode, nedge, inode, jnode );
/*
  Compute the next term of the Pruefer code.
*/
  for ( i = 1; i <= nnode - 2; i++ )
  {
/*
  Find the terminal node with the highest label.
*/
    iterm = 0;
 
    for ( j = 1; j <= nnode; j++ )
    {
      if ( degree[j-1] == 1 )
      {
        iterm = j;
      }
    }
/*
  Find the edge that includes this node, and note the
  index of the other node.
*/
    for ( j = 1; j < nnode - 1; j++ )
    {
      jsave = j;
 
      if ( inode[j-1] == iterm )
      {
        i2 = 2;
        break;
      }
      else if ( jnode[j-1] == iterm )
      {
        i2 = 1;
        break;
      }
    }
/*
  Delete the edge from the tree.
*/
    degree[inode[jsave-1]-1] = degree[inode[jsave-1]-1] - 1;
    degree[jnode[jsave-1]-1] = degree[jnode[jsave-1]-1] - 1;
/*
  Add the neighbor of the node to the Pruefer code.
*/
    if ( i2 == 1 )
    {
      code[i-1] = inode[jsave-1];
    }
    else
    {
      code[i-1] = jnode[jsave-1];
    }
/*
  Negate the nodes in the edge list to mark them as used.
*/
    inode[jsave-1] = - inode[jsave-1];
    jnode[jsave-1] = - jnode[jsave-1];
  }
/*
  Before returning, restore the original form of the edge list.
*/
  for ( i = 1; i <= nnode - 1; i++ )
  {
    inode[i-1] = abs ( inode[i-1] );
    jnode[i-1] = abs ( jnode[i-1] );
  }

  free ( degree );

  return code;
}
/******************************************************************************/

int tree_enum ( int nnode )

/******************************************************************************/
/*
  Purpose:

    TREE_ENUM enumerates the labeled trees on NNODE nodes.

  Discussion:

    The formula is due to Cauchy.

  Example:

    NNODE      NTREE

    0              1
    1              1
    2              1
    3              3
    4             16
    5            125
    6           1296
    7          16807
    8         262144
    9        4782969
   10      100000000

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    14 July 2013

  Author:

    John Burkardt

  Parameters:

    Input, int NNODE, the number of nodes in each tree.
    NNODE must normally be at least 3, but for this routine,
    any value of NNODE is allowed.  Values of NNODE greater than 10
    will probably overflow.

    Output, int TREE_ENUM, the number of distinct labeled trees.
*/
{
  int ntree;

  if ( nnode < 0 )
  {
    ntree = 0;
  }
  else if ( nnode == 0 )
  {
    ntree = 1;
  }
  else if ( nnode == 1 )
  {
    ntree = 1;
  }
  else if ( nnode == 2 )
  {
    ntree = 1;
  }
  else
  {
    ntree = i4_power ( nnode, nnode - 2 );
  }

  return ntree;
}
/******************************************************************************/

void tree_parent_next ( int nnode, int code[], int itree[], int *more )

/******************************************************************************/
/*
  Purpose:

    TREE_PARENT_NEXT generates, one at a time, all labeled trees.

  Discussion:

    The routine also returns the corresponding Pruefer codes.

  Formula:

    There are N^(N-2) labeled trees on N nodes (Cayley's formula).

    The number of trees in which node I has degree D(I) is the
    multinomial coefficient: ( N-2; D(1)-1, D(2)-1, ..., D(N)-1 ).

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    18 July 2013

  Parameters:

    Input, int NNODE, the number of nodes to be used in 
    the trees.

    Input/output, int CODE[NNODE].  The first NNODE-2 entries 
    of CODE contain the Pruefer code for the given labeled tree.

    Output, int ITREE[NNODE].  The first NNODE-1 entries 
    of ITREE describe the edges that go between the nodes.  Each pair
    (I, ITREE(I)) represents an edge.  Thus if ITREE(5) = 3,
    there is an edge from node 3 to node 5.

    Input/output, int *MORE.  On the first call only, the
    user is required to set MORE = .FALSE.  Then call the routine, and
    it will return information about the first tree
    as well as setting MORE to the value .TRUE.
    Keep calling to get another tree until MORE is .FALSE.
    on return, at which point there are no more trees.
*/
{
  int i;

  if ( *more )
  {
    for ( i = 0; i < nnode - 2; i++ )
    {
      code[i] = code[i] - 1;
     }
  }

  vec_next ( nnode-2, nnode, code, more );
 
  for ( i = 0; i < nnode - 2; i++ )
  {
    code[i] = code[i] + 1;
  }

  pruefer_to_tree_2 ( nnode, code, itree );
 
  return;
}
/******************************************************************************/

void tree_parent_to_arc ( int nnode, int parent[], int *nedge, int inode[], 
  int jnode[] )

/******************************************************************************/
/*
  Purpose:

    TREE_PARENT_TO_ARC converts a tree from parent to arc representation.

  Discussion:

    Parent representation lists the parent node of each node.  For a
    tree of N nodes, one node has a parent of 0, representing a null link.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    18 July 2013

  Author:

    John Burkardt

  Parameters:

    Input, int NNODE, the number of nodes in the tree.

    Input, int PARENT[NNODE], the parent node representation 
    of the tree.

    Output, int *NEDGE, the number of edges, normally NNODE-1.

    Output, int INODE[NEDGE], JNODE[NEDGE], pairs of nodes
    that define the links.
*/
{
  int i;

  *nedge = 0;

  for ( i = 1; i <= nnode; i++ )
  {
    if ( parent[i-1] != 0 )
    {
      *nedge = *nedge + 1;
      inode[*nedge-1] = i;
      jnode[*nedge-1] = parent[i-1];
    }
  }

  return;
}
/******************************************************************************/

int tree_rb_enum ( int n )

/******************************************************************************/
/*
  Purpose:

    TREE_RB_ENUM returns the number of rooted binary trees with N nodes.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    14 July 2013

  Author:

    John Burkardt

  Parameters:

    Input, int N, the number of nodes in the rooted 
    binary tree.  N should be odd.

    Output, int TREE_RB_ENUM, the number of rooted binary trees 
    with N nodes.
*/
{
  int *c;
  int m;
  int num;

  if ( n < 0 )
  {
    num = 0;
  }
  else if ( n == 0 )
  {
    num = 1;
  }
  else if ( ( n % 2 ) == 0 )
  {
    num = 0;
  }
  else
  {
    m = ( n - 1 ) / 2;
    c = catalan ( m );
    num = c[m];
    free ( c );
  }

  return num;
}
/******************************************************************************/

void tree_rb_lex_next ( int n, int a[], int *more )

/******************************************************************************/
/*
  Purpose:

    TREE_RB_LEX_NEXT generates rooted binary trees in lexicographic order.

  Discussion:

    The information definining the tree of N nodes is stored in a vector 
    of 0's and 1's, in preorder traversal form.  Essentially, the
    shape of the tree is traced out with a pencil that starts at the root,
    and stops at the very last null leaf.  The first time that a (non-null) 
    node is encountered, a 1 is added to the vector, and the left 
    descendant of the node is visited next.  When the path returns from
    the first descendant, the second descendant is visited.  When then path
    returns again, the path passes back up from the node to its parent.
    A null leaf is encountered only once, and causes a zero to be added to 
    the vector, and the path goes back up to the parent node.  

    The lexicographic order is used to order the vectors of 1's and 0's.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    05 August 2013

  Reference:

    Frank Ruskey,
    Combinatorial Generation,
    To appear.

  Parameters:

    Input, int N, the number of nodes in the rooted binary
    tree.  N should be odd.

    Input/output, int A[N], the preorder traversal form for
    the previous/next rooted binary tree.

    Output, logical *MORE, is TRUE if the next rooted binary tree was
    returned on this call, or FALSE if there are no more rooted binary
    trees, and the output of the previous call was the last one.
*/
{
  int i;
  int k;
  int p;
  int q;

  if ( ! *more )
  {
    for ( i = 1; i <= n - 2; i = i + 2 )
    {
      a[i-1] = 1;
    }
    for ( i = 2; i <= n - 1; i = i + 2 )
    {
      a[i-1] = 0;
    }
    a[n-1] = 0;
    *more = 1;
    return;
  }
/*
  Find the last 1 in A.
*/
  k = n;
  while ( a[k-1] == 0 )
  {
    k = k - 1;
  }
  q = n - k - 1;
/*
  Find the last 0 preceding the last 1 in A.
  If there is none, then we are done, because 11...1100..00 
  is the final element.
*/
  for ( ; ; )
  {
    if ( k == 1 )
    {
      *more = 0;
      return;
    }

    if ( a[k-1] == 0 )
    {
      break;
    }
    k = k - 1;
  }

  p = n - k - q - 1;
  a[k-1] = 1;
  for ( i = k + 1; i <= n - 2 * p + 1; i++ )
  {
    a[i-1] = 0;
  }
  for ( i = n - 2 * p + 2; i <= n - 2; i = i + 2 )
  {
    a[i-1] = 1;
  }
  for ( i = n - 2 * p + 3; i <= n - 1; i = i + 2 )
  {
    a[i-1] = 0;
  }
  a[n-1] = 0;

  return;
}
/******************************************************************************/

int *tree_rb_to_parent ( int n, int a[] )

/******************************************************************************/
/*
  Purpose:

    TREE_RB_TO_PARENT converts rooted binary tree to parent node representation.

  Discussion:

    Parent node representation of a tree assigns to each node a "parent" node,
    which represents the first link of the path between the node and the 
    root node.  The root node itself is assigned a parent of 0.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    18 July 2013

  Author:

    John Burkardt

  Parameters:

    Input, int N, the number of nodes in the tree.

    Input, int A[N], the preorder traversal form for the
    rooted binary tree.

    Output, int TREE_RB_TO_PARENT[N], the parent node representation 
    of the tree.
*/
{
  int dad;
  int k;
  int node;
  int node_num;
  int *parent;
  int *use;

  parent = ( int * ) malloc ( n * sizeof ( int ) );
  use = ( int * ) malloc ( n * sizeof ( int ) );

  node = 0;
  node_num = 0;

  for ( k = 1; k <= n; k++ )
  {
    dad = node;
    node_num = node_num + 1;
    node = node_num;
    parent[node-1] = dad;

    if ( a[k-1] == 1 )
    {
      use[node-1] = 0;
    }
    else
    {
      use[node-1] = 2;

      while ( use[node-1] == 2 )
      {
        node = dad;
        if ( node == 0 )
        {
          break;
        }
        use[node-1] = use[node-1] + 1;
        dad = parent[node-1];
      }
    }
  }
  free ( use );

  return parent;
}
/******************************************************************************/

void tree_rb_yule ( int *n, int *seed, int a[] )

/******************************************************************************/
/*
  Purpose:

    TREE_RB_YULE adds two nodes to a rooted binary tree using the Yule model.

  Discussion:

    The Yule model is a simulation of how an evolutionary family tree
    develops.  We start with a root node.  The internal nodes of the tree 
    are inactive and never change.  Each pendant or leaf node of the
    tree represents a biological family that can spontaneously "fission",
    developing two new distinct sub families.  In graphical terms, the node
    becomes internal, with two new leaf nodes depending from it.

    The tree is stored in inorder traversal form.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    05 August 2013

  Author:

    John Burkardt

  Parameters:

    Input/output, int *N, the number of nodes in the input
    tree.  On output, this number has been increased, usually by 2.

    Input/output, int *SEED, a seed for the random number
    generator.

    Input/output, int A[*], the preorder traversal form 
    for the rooted binary tree.  The number of entries in A is N.
*/
{
  int i;
  int ileaf;
  int j;
  int jleaf;
  int nleaf;

  if ( *n <= 0 )
  {
    *n = 1;
    a[0] = 0;
    return;
  }
/*
  Count the expected number of leaves, which are the 0 values.
*/
  nleaf = ( *n + 1 ) / 2;
/*
  Choose a random number between 1 and NLEAF.
*/
  ileaf = i4_uniform_ab ( 1, nleaf, seed );
/*
  Locate leaf number ILEAF.
*/
  j = 0;
  jleaf = 0;
  for ( i = 1; i <= *n; i++ )
  {
    if ( a[i-1] == 0 )
    {
      jleaf = jleaf + 1;
    }
    if ( jleaf == ileaf )
    {
      j = i;
      break;
    }
  }
/*
  Replace '0' by '100'
*/
  for ( i = *n; j <= i; i-- )
  {
    a[i+1] = a[i-1];
  }
  a[j-1] = 1;
  a[j] = 0;

  *n = *n + 2;

  return;
}
/******************************************************************************/

int *tree_rooted_code ( int nnode, int parent[] )

/******************************************************************************/
/*
  Purpose:

    TREE_ROOTED_CODE returns the code of a rooted tree.

  Discussion:

    This code for a rooted tree depends on the node ordering, so it's actually
    the code for a labeled rooted tree.  To eliminate the effects of node
    labeling, one could choose as the code for a tree the maximum of all
    the codes associated with the different possible labelings of the tree.
    There are more effective ways of arriving at this code than simply
    generating all possible codes and comparing them.  

    For a tree with NNODES, the code is a list of 2*NNODE 0's and 1's,
    describing a traversal of the tree starting at an imaginary node 0,
    moving "down" to the root (a code entry of 1), and then moving
    "down" (1) or "up" (0) as the tree is traversed in a depth first
    manner.  The final move must be from the root up to the imaginary
    node 0.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    18 July 2013

  Author:

    John Burkardt

  Parameters:

    Input, int NNODE, the number of nodes.

    Input, int PARENT[NNODE], is the parent node of each node.
    The node with parent 0 is the root.

    Output, int TREE_ROOTED_CODE[2*NNODE], the code for the tree.
*/
{
  int *code;
  int father;
  int i;
  int k;
  int son;

  code = ( int * ) malloc ( 2 * nnode * sizeof ( int ) );
/*
  Find the root.
*/
  father = 0;
  for ( i = 1; i <= nnode; i++ )
  {
    if ( parent[i-1] == 0 )
    {
      k = 1;
      code[0] = 1;
      father = i;
      break;
    }
  }

  if ( father == 0 )
  {
    fprintf ( stderr, "\n" );
    fprintf ( stderr, "TREE_ROOTED_CODE - Fatal error!\n" );
    fprintf ( stderr, "  Could not find the root.\n" );
    exit ( 1 );
  }

  while ( father != 0 ) 
  {
    k = k + 1;
    code[k-1] = 0;
    for ( son = 1; son <= nnode; son++ )
    {
      if ( parent[son-1] == father )
      {
        code[k-1] = 1;
        father = son;
        break;
      }
    }

    if ( code[k-1] == 0 )
    {
      parent[father-1] = - parent[father-1];
      father = - parent[father-1];
    }
  }

  for ( i = 0; i < nnode; i++ )
  {
    parent[i] = - parent[i];
  }

  return code;
}
/******************************************************************************/

int tree_rooted_code_compare ( int nnode, int npart, int code1[], int code2[] )

/******************************************************************************/
/*
  Purpose:

    TREE_ROOTED_CODE_COMPARE compares a portion of the code for two rooted trees.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    18 July 2013

  Author:

    John Burkardt

  Parameters:

    Input, int NNODE, the number of nodes.

    Input, int NPART, the number of nodes for which the code
    has been determined.  This determines the portion of the codes to be
    compared.  We expect 0 <= NPART <= NNODE.

    Input, int CODE1[2*NNODE], CODE2[2*NNODE], the two 
    rooted tree codes to be compared.

    Output, int TREE_ROOTED_CODE_COMPARE, the result of the comparison.
    -1, CODE1 < CODE2,
     0, CODE1 = CODE2,
    +1, CODE1 > CODE2.
*/
{
  int i;
  int ihi;
  int result;

  result = 0;

  if ( npart <= 0 )
  {
    return result;
  }

  ihi = 2 * nnode;
  if ( npart < nnode )
  {
    ihi = 2 * npart;
  }

  for ( i = 0; i < ihi; i++ )
  {
    if ( code1[i] < code2[i] )
    {
      result = -1;
      return result;
    }
    else if ( code2[i] < code1[i] )
    {
      result = +1;
      return result;
    }
  }

  return result;
}
/******************************************************************************/

void tree_rooted_depth ( int nnode, int parent[], int *depth, int depth_node[] )

/******************************************************************************/
/*
  Purpose:

    TREE_ROOTED_DEPTH returns the depth of a rooted tree.

  Discussion:

    The depth of any node of a rooted tree is the number of edges in 
    the shortest path from the root to the node.

    The depth of the rooted tree is the maximum of the depths
    of all the nodes.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    17 July 2013

  Author:

    John Burkardt

  Parameters:

    Input, int NNODE, the number of nodes.

    Input, int PARENT[NNODE], is the parent node of each node.
    The node with parent 0 is the root.

    Output, int *DEPTH, the depth of the tree.

    Output, int DEPTH_NODE[NNODE], the depth of each node.
*/
{
  int i;
  int j;
  int root;
/*
  Find the root.
*/
  root = -1;
  for ( i = 1; i <= nnode; i++ )
  {
    if ( parent[i-1] == 0 )
    {
      root = i;
      break;
    }
  }

  if ( root == -1 )
  {
    fprintf ( stderr, "\n" );
    fprintf ( stderr, "TREE_ROOTED_DEPTH - Fatal error!\n" );
    fprintf ( stderr, "  Could not find the root.\n" );
    exit ( 1 );
  }
/*
  Determine the depth of each node by moving towards the node.
  If you reach a node whose depth is already known, stop early.
*/
  for ( i = 0; i < nnode; i++ )
  {
    depth_node[i] = 0;
  }

  for ( i = 1; i <= nnode; i++ )
  {
    j = i;

    while ( j != root )
    {
      depth_node[i-1] = depth_node[i-1] + 1;
      j = parent[j-1];

      if ( 0 < depth_node[j-1] )
      {
        depth_node[i-1] = depth_node[i-1] + depth_node[j-1];
        break;
      }
    }
  }
/*
  Determine the maximum depth.
*/
  *depth = i4vec_max ( nnode, depth_node );

  return;
}
/******************************************************************************/

int *tree_rooted_enum ( int nnode )

/******************************************************************************/
/*
  Purpose:

    TREE_ROOTED_ENUM counts the number of unlabeled rooted trees.

  Example:

    Input    Output

      1         1
      2         1
      3         2
      4         4
      5         9
      6        20
      7        48
      8       115
      9       286
     10       719

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    14 July 2013

  Author:

    Original FORTRAN77 version by Albert Nijenhuis, Herbert Wilf.
    C version by John Burkardt.

  Reference:

    Albert Nijenhuis, Herbert Wilf,
    Combinatorial Algorithms,
    Academic Press, 1978, second edition,
    ISBN 0-12-519260-6.

  Parameters:

    Input, int NNODE, the number of nodes.

    Output, int TREE_ROOTED_ENUM[NNODE].  NTREE(I) is the number of 
    rooted, unlabeled trees on I nodes, for I = 1, 2, ... NNODE.
*/
{
  int i;
  int id;
  int isum;
  int itd;
  int j;
  int nlast;
  int *ntree;

  ntree = ( int * ) malloc ( nnode * sizeof ( int ) );

  ntree[0] = 1;
 
  for ( nlast = 2; nlast <= nnode; nlast++ )
  {
    isum = 0;
 
    for ( id = 1; id <= nlast - 1; id++ )
    {
      i = nlast;
      itd = ntree[id-1] * id;
 
      for ( j = 1; j <= nlast - 1; j++ )
      {
        i = i - id;

        if ( i <= 0 )
        {
          break;
        }
        isum = isum + ntree[i-1] * itd;
      }
    }
    ntree[nlast-1] = isum / ( nlast - 1 );
  }
  return ntree;
}
/******************************************************************************/

int *tree_rooted_random ( int nnode, int *seed )

/******************************************************************************/
/*
  Purpose:

    TREE_ROOTED_RANDOM selects a random unlabeled rooted tree.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    17 July 2013

  Author:

    Original FORTRAN77 version by Albert Nijenhuis, Herbert Wilf.
    C version by John Burkardt.

  Reference:

    Albert Nijenhuis, Herbert Wilf,
    Combinatorial Algorithms,
    Academic Press, 1978, second edition,
    ISBN 0-12-519260-6.

  Parameters:

    Input, int NNODE, the number of nodes.

    Input/output, int *SEED, a seed for the random 
    number generator.

    Output, int TREE_ROOTED_RANDOM[NNODE].  (I,ITREE(I)) is the I-th edge
    of the output tree for I = 2,NNODE.  ITREE(1)=0.
*/
{
  int i;
  int id;
  int is1;
  int is2;
  int itd;
  int *itree;
  int iz;
  int j;
  int l;
  int ll;
  int ls;
  int m;
  int *ntree;
  int nval;
  double r;
  int *stack;

  if ( nnode <= 0  )
  {
    fprintf ( stderr, "\n" );
    fprintf ( stderr, "TREE_ROOTED_RANDOM - Fatal error!\n" );
    fprintf ( stderr, "  NNODE = %d\n", nnode );
    fprintf ( stderr, "  but NNODE must be at least 1.\n" );
    exit ( 1 );
  }

  itree = ( int * ) malloc ( nnode * sizeof ( int ) );
  stack = ( int * ) malloc ( 2 * nnode * sizeof ( int ) );
/*
  Compute a table of the number of such trees for a given number of nodes.
*/
  ntree = tree_rooted_enum ( nnode );
/*
  Now select one such tree at random.
*/
  l = 0;

  nval = nnode;
  is1 = 0;
  is2 = 0;
  
  for ( ; ; )
  {
    while ( 2 < nval )
    {
      r = r8_uniform_01 ( seed );

      iz = ( int ) ( ( nval - 1 ) * ntree[nval-1] * r );

      id = 0;
  
      id = id + 1;
      itd = id * ntree[id-1];
      m = nval;
      j = 0;
 
      for ( ; ; )
      {
        j = j + 1;
        m = m - id;

        if ( m < 1 )
        {
          id = id + 1;
          itd = id * ntree[id-1];
          m = nval;
          j = 0;
          continue;
        }

        iz = iz - ntree[m-1] * itd;
        if ( iz < 0 )
        {
          break;
        }
      }
      is1 = is1 + 1;
      stack[0+(is1-1)*2] = j;
      stack[1+(is1-1)*2] = id;
      nval = m;
    }
 
    itree[is2] = l;
    l = is2 + 1;
    is2 = is2 + nval;

    if ( 1 < nval )
    {
      itree[is2-1] = is2 - 1;
    }
 
    for ( ; ; )
    {
      nval = stack[1+(is1-1)*2];
 
      if ( nval != 0 )
      {
        stack[1+(is1-1)*2] = 0;
        break;
      }
 
      j = stack[0+(is1-1)*2];
      is1 = is1 - 1;
      m = is2 - l + 1;
      ll = itree[l-1];
      ls = l + ( j - 1 ) * m - 1;
 
      if ( j != 1 )
      {
        for ( i = l; i <= ls; i++ )
        {
          itree[i+m-1] = itree[i-1] + m;
          if ( ( ( i - l ) % m ) == 0 )
          {
            itree[i+m-1] = ll;
          }
        }
      }
      is2 = ls + m;
 
      if ( is2 == nnode )
      {
        free ( ntree );
        free ( stack );
        return itree;
      }
      l = ll;
    }
  }
 
  free ( ntree );
  free ( stack );

  return itree;
}
/******************************************************************************/

void vec_next ( int n, int ibase, int iarray[], int *more )

/******************************************************************************/
/*
  Purpose:

    VEC_NEXT generates all N-vectors of integers modulo a given base.

  Discussion:

    The items are produced one at a time.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    14 July 2013

  Parameters:

    Input, int N, the size of the vectors to be used.

    Input, int IBASE, the base to be used.  IBASE = 2 will
    give vectors of 0's and 1's, for instance.

    Input/output, int IARRAY[N].  On each return,
    IARRAY will contain entries in the range 0 to IBASE-1.

    Input/output, int *MORE.  Set this variable 0 before
    the first call.  Normally, MORE will be returned 1 but
    once all the vectors have been generated, MORE will be
    reset 0 and you should stop calling the program.
*/
{
  int i;
  static int kount = 0;
  static int last = 0;
  int nn;

  if ( ! ( *more ) )
  {
    kount = 1;
    last = i4_power ( ibase, n );
    *more = 1;
    for ( i = 0; i < n; i++ )
    {
      iarray[i] = 0;
    }
  }
  else
  {
    kount = kount + 1;

    if ( kount == last )
    {
      *more = 0;
    }

    iarray[n-1] = iarray[n-1] + 1;
 
    for ( i = 1; i <= n; i++ )
    {
      nn = n - i;

      if ( iarray[nn] < ibase )
      {
        return;
      }

      iarray[nn] = 0;

      if ( nn != 0 )
      {
        iarray[nn-1] = iarray[nn-1] + 1;
      }
    }
  }
  return;
}
/******************************************************************************/

void vec_random ( int n, int base, int *seed, int a[] )

/******************************************************************************/
/*
  Purpose:

    VEC_RANDOM selects a random N-vector of integers modulo a given base.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    21 March 2009

  Author:

    John Burkardt

  Parameters:

    Input, int N, the size of the vector to be generated.

    Input, int BASE, the base to be used.

    Input/output, int *SEED, a random number seed.

    Output, int A[N], a list of N random values between
    0 and BASE-1.
*/
{
  int i;

  for ( i = 0; i < n; i++ )
  {
    a[i] = i4_uniform_ab ( 0, base-1, seed );
  }

  return;
}
/******************************************************************************/

int *vec_random_new ( int n, int base, int *seed )

/******************************************************************************/
/*
  Purpose:

    VEC_RANDOM_NEW selects a random N-vector of integers modulo a given base.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    21 March 2009

  Author:

    John Burkardt

  Parameters:

    Input, int N, the size of the vector to be generated.

    Input, int BASE, the base to be used.

    Input/output, int *SEED, a random number seed.

    Output, int VEC_RANDOM_NEW[N], a list of N random values between
    0 and BASE-1.
*/
{
  int *a;
  int i;

  a = ( int * ) malloc ( n * sizeof ( int ) );

  for ( i = 0; i < n; i++ )
  {
    a[i] = i4_uniform_ab ( 0, base-1, seed );
  }

  return a;
}