# include # include # include # include # include # include "combo.h" /******************************************************************************/ void backtrack ( int l, int iarray[], int *indx, int *k, int *nstack, int stack[], int maxstack ) /******************************************************************************/ /* Purpose: backtrack() supervises a backtrack search. Discussion: The routine builds a vector, one element at a time, which is required to satisfy some condition. At any time, the partial vector may be discovered to be unsatisfactory, but the routine records information about where the last arbitrary choice was made, so that the search can be carried out efficiently, rather than starting out all over again. Licensing: This code is distributed under the MIT license. Modified: 25 July 2011 Author: Original FORTRAN77 version by Albert Nijenhuis, Herbert Wilf. C 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, int L, the length of the completed candidate vector. Input/output, int IARRAY[L], the candidate vector. Input/output, int *INDX. On input, set INDX = 0 to start a search. On output: 1, a complete output vector has been determined. 2, candidates are needed. 3, no more possible vectors exist. Input/output, int *K, the current length of the candidate vector. Input/output, int *NSTACK, the current length of the stack. Input/output, int STACK[MAXSTACK], a list of candidates for positions 1 through K. Input, int MAXSTACK, the maximum length of the stack. */ { /* If this is the first call, request a candidate for position 1. */ if ( *indx == 0 ) { *k = 1; *nstack = 0; *indx = 2; return; } /* Examine the stack. */ for ( ; ; ) { *nstack = *nstack - 1; /* If there are candidates for position K, take the first available one off the stack, and increment K. This may cause K to reach the desired value of L, in which case we need to signal the user that a complete set of candidates is being returned. */ if ( stack[*nstack] != 0 ) { iarray[*k-1] = stack[*nstack-1]; stack[*nstack-1] = stack[*nstack] - 1; if ( *k != l ) { *k = *k + 1; *indx = 2; } else { *indx = 1; } break; } /* If there are no candidates for position K, then decrement K. If K is still positive, repeat the examination of the stack. */ else { *k = *k - 1; if ( *k <= 0 ) { *indx = 3; break; } } } return; } /******************************************************************************/ int bal_seq_check ( int n, int t[] ) /******************************************************************************/ /* Purpose: bal_seq_check() checks a balanced sequence. Licensing: This code is distributed under the MIT license. Modified: 25 November 2015 Author: John Burkardt Reference: Donald Kreher, Douglas Simpson, Combinatorial Algorithms, CRC Press, 1998, ISBN: 0-8493-3988-X, LC: QA164.K73. Parameters: Input, int N, the number of 0's (and 1's) in the sequence. N must be positive. Input, int T[2*N], a balanced sequence. Output, int BAL_SEQ_CHECK. 1, the data is legal. 0, the data is not legal. */ { int check; int i; int one_count; int zero_count; check = 1; if ( n < 1 ) { check = 0; return check; } one_count = 0; zero_count = 0; for ( i = 0; i < 2 * n; i++ ) { if ( t[i] == 0 ) { zero_count = zero_count + 1; } else if ( t[i] == 1 ) { one_count = one_count + 1; } else { check = 0; return check; } if ( zero_count < one_count ) { check = 0; return check; } } if ( one_count != zero_count ) { check = 0; } return check; } /******************************************************************************/ int bal_seq_enum ( int n ) /******************************************************************************/ /* Purpose: bal_seq_enum() enumerates the balanced sequences. Licensing: This code is distributed under the MIT license. Modified: 24 July 2011 Author: John Burkardt Reference: Donald Kreher, Douglas Simpson, Combinatorial Algorithms, CRC Press, 1998, ISBN: 0-8493-3988-X, LC: QA164.K73. Parameters: Input, int N, the number of 0's (and 1's) in the sequence. N must be nonnegative. Output, int BAL_SEQ_ENUM, the number of balanced sequences. */ { int value; value = i4_choose ( 2 * n, n ) / ( n + 1 ); return value; } /******************************************************************************/ int bal_seq_rank ( int n, int t[] ) /******************************************************************************/ /* Purpose: bal_seq_rank() ranks a balanced sequence. Licensing: This code is distributed under the MIT license. Modified: 24 July 2011 Author: John Burkardt Reference: Donald Kreher, Douglas Simpson, Combinatorial Algorithms, CRC Press, 1998, ISBN: 0-8493-3988-X, LC: QA164.K73. Parameters: Input, int N, the number of 0's (and 1's) in the sequence. N must be positive. Input, int T[2*N], a balanced sequence. Output, int BAL_SEQ_RANK, the rank of the balanced sequence. */ { int check; int mxy; int rank; int x; int y; /* Check. */ check = bal_seq_check ( n, t ); if ( ! check ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "BAL_SEQ_RANK(): Fatal error!\n" ); fprintf ( stderr, " The input array is illegal.\n" ); exit ( 1 ); } y = 0; rank = 0; for ( x = 1; x <= 2 * n - 1; x++ ) { if ( t[x-1] == 0 ) { y = y + 1; } else { mxy = mountain ( n, x, y + 1 ); rank = rank + mxy; y = y - 1; } } return rank; } /******************************************************************************/ void bal_seq_successor ( int n, int t[], int *rank ) /******************************************************************************/ /* Purpose: bal_seq_successor() computes the lexical balanced sequence successor. Licensing: This code is distributed under the MIT license. Modified: 26 July 2011 Author: John Burkardt Reference: Donald Kreher, Douglas Simpson, Combinatorial Algorithms, CRC Press, 1998, ISBN: 0-8493-3988-X, LC: QA164.K73. Parameters: Input, int N, the number of 0's (and 1's) in the sequence. N must be positive. Input/output, int T[2*N], on input, a balanced sequence, and on output, its lexical successor. Input/output, int *RANK, the rank. If RANK = -1 on input, then the routine understands that this is the first call, and that the user wishes the routine to supply the first element in the ordering, which has RANK = 0. In general, the input value of RANK is increased by 1 for output, unless the very last element of the ordering was input, in which case the output value of RANK is 0. */ { int check; int i; int j; int open; int open_index; int slot; int slot_index; int slot_ones; /* Return the first element. */ if ( *rank == -1 ) { for ( i = 0; i < n; i++ ) { t[i] = 0; } for ( i = n; i < 2 * n; i++ ) { t[i] = 1; } *rank = 0; return; } /* Check. */ check = bal_seq_check ( n, t ); if ( ! check ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "BAL_SEQ_SUCCESSOR(): Fatal error!\n" ); fprintf ( stderr, " The input array is illegal.\n" ); exit ( 1 ); } /* After the I-th 0 there is a 'slot' with the capacity to hold between 0 and I ones. The first element of the sequence has all the 1''s cowering behind the N-th 0. We seek to move a 1 to the left, and to do it lexically, we will move a 1 to the rightmost slot that is under capacity. Find the slot. */ slot = 0; slot_index = 0; slot_ones = 0; open = 0; open_index = 0; for ( i = 1; i <= 2 * n; i++ ) { if ( t[i-1] == 0 ) { if ( 0 < slot ) { if ( slot_ones < slot ) { open = slot; open_index = slot_index; } } slot = slot + 1; slot_index = i; } else { slot_ones = slot_ones + 1; } } /* If OPEN is not 0, then preserve the string up to the OPEN-th 0, preserve the 1''s that follow, but then write a 1, then all the remaining 0's and all the remaining 1's. */ if ( open != 0 ) { j = open_index + 1; while ( t[j-1] == 1 ) { j = j + 1; } t[j-1] = 1; for ( i = open + 1; i <= n; i++ ) { j = j + 1; t[j-1] = 0; } for ( i = j + 1; i <= 2 * n; i++ ) { t[i-1] = 1; } } /* If OPEN is 0, the last element was input. Return the first one. */ else { for ( i = 0; i < n; i++ ) { t[i] = 0; } for ( i = n; i < 2 * n; i++ ) { t[i] = 1; } *rank = 0; return; } *rank = *rank + 1; return; } /******************************************************************************/ int *bal_seq_to_tableau ( int n, int t[] ) /******************************************************************************/ /* Purpose: bal_seq_to_tableau() converts a balanced sequence to a 2 by N tableau. Licensing: This code is distributed under the MIT license. Modified: 26 July 2011 Author: John Burkardt Reference: Donald Kreher, Douglas Simpson, Combinatorial Algorithms, CRC Press, 1998, ISBN: 0-8493-3988-X, LC: QA164.K73. Parameters: Input, int N, the number of 0's (and 1's) in the sequence. N must be positive. Input, int T[2*N], a balanced sequence. Output, int BAL_SEQ_TO_TABLEAU[2*N], a 2 by N tableau. */ { int c[2]; int check; int i; int r; int *tab; /* Check. */ check = bal_seq_check ( n, t ); if ( ! check ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "BAL_SEQ_TO_TABLEAU(): Fatal error!\n" ); fprintf ( stderr, " The input array is illegal.\n" ); exit ( 1 ); } tab = ( int * ) malloc ( 2 * n * sizeof ( int ) ); c[0] = 0; c[1] = 0; for ( i = 1; i <= 2 * n; i++ ) { r = t[i-1] + 1; c[r-1] = c[r-1] + 1; tab[r-1+(c[r-1]-1)*2] = i; } return tab; } /******************************************************************************/ int *bal_seq_unrank ( int rank, int n ) /******************************************************************************/ /* Purpose: bal_seq_unrank() unranks a balanced sequence. Licensing: This code is distributed under the MIT license. Modified: 26 July 2011 Author: John Burkardt Reference: Donald Kreher, Douglas Simpson, Combinatorial Algorithms, CRC Press, 1998, ISBN: 0-8493-3988-X, LC: QA164.K73. Parameters: Input, int RANK, the rank of the balanced sequence. Input, int N, the number of 0's (and 1's) in the sequence. N must be positive. Output, int BAL_SEQ_UNRANK[2*N], a balanced sequence. */ { int low; int m; int nseq; int *t; int x; int y; /* Check. */ if ( n < 1 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "BAL_SEQ_UNRANK(): Fatal error!\n" ); fprintf ( stderr, " Input N is illegal.\n" ); exit ( 1 ); } nseq = bal_seq_enum ( n ); if ( rank < 0 || nseq < rank ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "BAL_SEQ_UNRANK(): Fatal error!\n" ); fprintf ( stderr, " The input rank is illegal.\n" ); fprintf ( stderr, " Input rank = %d\n", rank ); fprintf ( stderr, " But 0 <= rank < %d is required.\n", nseq ); exit ( 1 ); } t = ( int * ) malloc ( 2 * n * sizeof ( int ) ); y = 0; low = 0; for ( x = 0; x < 2 * n; x++ ) { m = mountain ( n, x + 1, y + 1 ); if ( rank <= low + m - 1 ) { y = y + 1; t[x] = 0; } else { low = low + m; y = y - 1; t[x] = 1; } } return t; } /******************************************************************************/ int *bell_numbers ( int m ) /******************************************************************************/ /* Purpose: bell_numbers() computes the Bell numbers. Discussion: There are B(M) restricted growth functions of length M. There are B(M) partitions of a set of M objects. B(M) is the sum of the Stirling numbers of the second kind, S(M,N), for N = 0 to M. Licensing: This code is distributed under the MIT license. Modified: 24 July 2011 Author: John Burkardt Reference: Donald Kreher, Douglas Simpson, Combinatorial Algorithms, CRC Press, 1998, ISBN: 0-8493-3988-X, LC: QA164.K73. Parameters: Input, int M, indicates how many Bell numbers are to compute. M must be nonnegative. Output, int bell_numbers[M+1], the first M+1 Bell numbers. */ { int *b; int i; int j; b = ( int * ) malloc ( ( m + 1 ) * sizeof ( int ) ); b[0] = 1; for ( j = 1; j <= m; j++ ) { b[j] = 0; for ( i = 0; i < j; i++ ) { b[j] = b[j] + i4_choose ( j - 1, i ) * b[i]; } } return b; } /******************************************************************************/ void bell_values ( int *n_data, int *n, int *c ) /******************************************************************************/ /* Purpose: bell_values() returns some values of the Bell numbers. Discussion: 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. The Bell numbers were named for Eric Temple Bell. In Mathematica, the function can be evaluated by Sum[StirlingS2[n,m],{m,1,n}] Definition: 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. 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 MIT 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. 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 Bell number. Output, int *C, the value of the Bell number. */ { # define N_MAX 11 static int c_vec[N_MAX] = { 1, 1, 2, 5, 15, 52, 203, 877, 4140, 21147, 115975 }; static 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 } /******************************************************************************/ int cycle_check ( int n, int ncycle, int t[], int index[] ) /******************************************************************************/ /* Purpose: cycle_check() checks a permutation in cycle form. Licensing: This code is distributed under the MIT license. Modified: 31 November 2015 Author: John Burkardt Reference: Donald Kreher, Douglas Simpson, Combinatorial Algorithms, CRC Press, 1998, ISBN: 0-8493-3988-X, LC: QA164.K73. Parameters: Input, int N, the number of items permuted. N must be positive. Input, int NCYCLE, the number of cycles. 1 <= NCYCLE <= N. Input, int T[N], INDEX[NCYCLE], describes the permutation as a collection of NCYCLE cycles. The first cycle is T(1) -> T(2) -> ... -> T(INDEX(1)) -> T(1). Output, int CYCLE_CHECK. 1, the data is legal. 0, the data is not legal. */ { int check; int i; int ifind; int iseek; check = 1; /* N must be at least 1. */ if ( n < 1 ) { check = 0; return check; } /* 1 <= NCYCLE <= N. */ if ( ncycle < 1 || n < ncycle ) { check = 0; return check; } /* 1 <= INDEX(I) <= N. */ for ( i = 0; i < ncycle; i++ ) { if ( index[i] < 1 || n < index[i] ) { check = 0; return check; } } /* The INDEX values sum to N. */ if ( i4vec_sum ( ncycle, index ) != n ) { check = 0; return check; } /* 1 <= T(I) <= N. */ for ( i = 0; i < n; i++ ) { if ( t[i] < 1 || n < t[i] ) { check = 0; return check; } } /* Verify that every value from 1 to N occurs in T. */ for ( iseek = 1; iseek <= n; iseek++ ) { ifind = -1; for ( i = 0; i < n; i++ ) { if ( t[i] == iseek ) { ifind = i + 1; break; } } if ( ifind == -1 ) { check = 0; return check; } } return check; } /******************************************************************************/ int *cycle_to_perm ( int n, int ncycle, int t[], int index[] ) /******************************************************************************/ /* Purpose: cycle_to_perm() converts a permutation from cycle to array form. Licensing: This code is distributed under the MIT license. Modified: 26 November 2015 Author: John Burkardt Reference: Donald Kreher, Douglas Simpson, Combinatorial Algorithms, CRC Press, 1998, ISBN: 0-8493-3988-X, LC: QA164.K73. Parameters: Input, int N, the number of items permuted. N must be positive. Input, int NCYCLE, the number of cycles. 1 <= NCYCLE <= N. Input, int T[N], INDEX[NCYCLE], describes the permutation as a collection of NCYCLE cycles. The first cycle is T(1) -> T(2) -> ... -> T(INDEX(1)) -> T(1). Output, int CYCLE_TO_PERM[N], describes the permutation using a single array. For each index I, I -> P(I). */ { int check; int i; int j; int jhi; int jlo; int *p; /* Check. */ check = cycle_check ( n, ncycle, t, index ); if ( ! check ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "CYCLE_TO_PERM(): Fatal error!\n" ); fprintf ( stderr, " CYCLE_CHECK rejects the input data.\n" ); exit ( 1 ); } p = ( int * ) malloc ( n * sizeof ( int ) ); jhi = 0; for ( i = 1; i <= ncycle; i++ ) { jlo = jhi + 1; jhi = jhi + index[i-1]; for ( j = jlo; j <= jhi; j++ ) { if ( j < jhi ) { p[t[j-1]-1] = t[j]; } else { p[t[j-1]-1] = t[jlo-1]; } } } return p; } /******************************************************************************/ int dist_enum ( int k, int m ) /******************************************************************************/ /* Purpose: dist_enum() returns the number of distributions of indistinguishable objects. Licensing: This code is distributed under the MIT license. Modified: 24 July 2011 Author: John Burkardt Parameters: Input, int K, the number of distinguishable "slots". Input, int M, the number of indistinguishable objects. Output, int DIST_ENUM, the number of distributions of M indistinguishable objects about K distinguishable slots. */ { int value; value = i4_choose ( m + k - 1, m ); return value; } /******************************************************************************/ void dist_next ( int k, int m, int q[], int *leftmost, int *more ) /******************************************************************************/ /* Purpose: dist_next() returns the next distribution of indistinguishable objects. Discussion: A distribution of M objects into K parts is an ordered sequence of K nonnegative integers which sum to M. This is similar to a partition of a set into K subsets, except that here the order matters. That is, (1,1,2) and (1,2,1) are considered to be different distributions. On the first call to this routine, the user should set MORE = FALSE, to signal that this is a startup for the given computation. The routine will return the first distribution, and set MORE = TRUE. If the user calls again, with MORE = TRUE, the next distribution is being requested. If the routine returns with MORE = TRUE, then that distribution was found and returned. However, if the routine returns with MORE = FALSE, then no more distributions were found; the enumeration of distributions has terminated. A "distribution of M indistinguishable objects into K slots" is sometimes called a "composition of the integer M into K parts". Example: K = 3, M = 5 0 0 5 0 1 4 0 2 3 0 3 2 0 4 1 0 5 0 1 0 4 1 1 3 1 2 2 1 3 1 1 4 0 2 0 3 2 1 2 2 2 1 2 3 0 3 0 2 3 1 1 3 2 0 4 0 1 4 1 0 5 0 0 Licensing: This code is distributed under the MIT license. Modified: 26 July 2011 Author: John Burkardt Reference: Robert Fenichel, Algorithm 329: Distribution of Indistinguishable Objects into Distinguishable Slots, Communications of the ACM, Volume 11, Number 6, June 1968, page 430. Parameters: Input, int K, the number of distinguishable "slots". Input, int M, the number of indistinguishable objects. Input/output, int Q[K], the number of objects in each slot. Input/output, int *LEFTMOST, set this to 0 before the first call. Input/output, int *MORE, used by the user to start the computation, and by the routine to stop the computation. */ { int i; /* The startup call. */ if ( ! ( *more ) ) { *more = 1; for ( i = 0; i < k - 1; i++ ) { q[i] = 0; } q[k-1] = m; *leftmost = k + 1; } /* There are no more distributions. Reset Q to the first distribution in the sequence. */ else if ( q[0] == m ) { *more = 0; for ( i = 0; i < k - 1; i++ ) { q[i] = 0; } q[k-1] = m; *leftmost = k + 1; } else if ( *leftmost < k + 1 ) { *leftmost = *leftmost - 1; q[k-1] = q[*leftmost-1] - 1; q[*leftmost-1] = 0; q[*leftmost-2] = q[*leftmost-2] + 1; if ( q[k-1] != 0 ) { *leftmost = k + 1; } } else { if ( q[k-1] == 1 ) { *leftmost = k; } q[k-1] = q[k-1] - 1; q[k-2] = q[k-2] + 1; } return; } /******************************************************************************/ int edge_check ( int n_node, int n_edge, int t[] ) /******************************************************************************/ /* Purpose: edge_check() checks a graph stored by edges. Licensing: This code is distributed under the MIT license. Modified: 26 November 2015 Author: John Burkardt Parameters: Input, int N_NODE, the number of nodes in the graph. N_NODE must be positive. Input, int N_EDGE, the number of edges in the graph. N_EDGE must be positive. Input, int T(2,N_EDGE), describes the edges of the tree as pairs of nodes. Output, int EDGE_CHECK. 1, the data is legal. 0, the data is not legal. */ { int check; int i; int j; int j2; check = 1; if ( n_node < 1 ) { check = 0; return check; } if ( n_edge < 1 ) { check = 0; return check; } /* Every edge must join two legal nodes. */ for ( i = 0; i < 2; i++ ) { for ( j = 0; j < n_edge; j++ ) { if ( t[i+j*2] < 1 || n_node < t[i+j*2] ) { check = 0; return check; } } } /* Every edge must join distinct nodes. */ for ( j = 0; j < n_edge; j++ ) { if ( t[0+j*2] == t[1+j*2] ) { check = 0; return check; } } /* Every edge must be distinct. */ for ( j = 0; j < n_edge - 1; j++ ) { for ( j2 = j + 1; j2 < n_edge; j2++ ) { if ( t[0+j*2] == t[0+j2*2] && t[1+j*2] == t[1+j2*2] ) { check = 0; return check; } else if ( t[0+j*2] == t[1+j2*2] && t[1+j*2] == t[0+j2*2] ) { check = 0; return check; } } } return check; } /******************************************************************************/ int *edge_degree ( int n_node, int n_edge, int t[] ) /******************************************************************************/ /* Purpose: edge_degree() returns the degree of the nodes of a graph stored by edges. Licensing: This code is distributed under the MIT license. Modified: 26 November 2015 Author: John Burkardt Reference: Donald Kreher, Douglas Simpson, Combinatorial Algorithms, CRC Press, 1998, ISBN: 0-8493-3988-X, LC: QA164.K73. Parameters: Input, int N_NODE, the number of nodes in the graph. N_NODE must be positive. Input, int N_EDGE, the number of edges in the graph. N_EDGE must be positive. Input, int T[2*N_EDGE], describes the edges of the tree as pairs of nodes. Output, int EDGE_DEGREE[N_NODE], the degree of each node. */ { int check; int *d; int i; int j; /* Check. */ check = edge_check ( n_node, n_edge, t ); if ( ! check ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "EDGE_DEGREE(): Fatal error!\n" ); fprintf ( stderr, " The input array is illegal.\n" ); exit ( 1 ); } /* Compute the degree of each node. */ d = ( int * ) malloc ( n_node * sizeof ( int ) ); for ( i = 0; i < n_node; i++ ) { d[i] = 0; } for ( j = 0; j < n_edge; j++ ) { d[t[0+j*2]-1] = d[t[0+j*2]-1] + 1; d[t[1+j*2]-1] = d[t[1+j*2]-1] + 1; } return d; } /******************************************************************************/ int edge_enum ( int n_node ) /******************************************************************************/ /* Purpose: edge_enum() enumerates the maximum number of edges in a graph on N_NODE nodes. Licensing: This code is distributed under the MIT license. Modified: 24 July 2011 Author: John Burkardt Parameters: Input, int N_NODE, the number of nodes in the graph. N_NODE must be positive. Output, int EDGE_ENUM, the maximum number of edges in a graph on N_NODE nodes. */ { int value; value = ( n_node * ( n_node - 1 ) ) / 2; return value; } /******************************************************************************/ void gamma_log_values ( int *n_data, double *x, double *fx ) /******************************************************************************/ /* Purpose: gamma_log_values() returns some values of the Log Gamma function. Discussion: In Mathematica, the function can be evaluated by: Log[Gamma[x]] Licensing: This code is distributed under the MIT license. Modified: 14 August 2004 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, double *X, the argument of the function. Output, double *FX, the value of the function. */ { # define N_MAX 20 static double fx_vec[N_MAX] = { 0.1524063822430784E+01, 0.7966778177017837E+00, 0.3982338580692348E+00, 0.1520596783998375E+00, 0.0000000000000000E+00, -0.4987244125983972E-01, -0.8537409000331584E-01, -0.1081748095078604E+00, -0.1196129141723712E+00, -0.1207822376352452E+00, -0.1125917656967557E+00, -0.9580769740706586E-01, -0.7108387291437216E-01, -0.3898427592308333E-01, 0.00000000000000000E+00, 0.69314718055994530E+00, 0.17917594692280550E+01, 0.12801827480081469E+02, 0.39339884187199494E+02, 0.71257038967168009E+02 }; static double x_vec[N_MAX] = { 0.20E+00, 0.40E+00, 0.60E+00, 0.80E+00, 1.00E+00, 1.10E+00, 1.20E+00, 1.30E+00, 1.40E+00, 1.50E+00, 1.60E+00, 1.70E+00, 1.80E+00, 1.90E+00, 2.00E+00, 3.00E+00, 4.00E+00, 10.00E+00, 20.00E+00, 30.00E+00 }; if ( *n_data < 0 ) { *n_data = 0; } *n_data = *n_data + 1; if ( N_MAX < *n_data ) { *n_data = 0; *x = 0.0; *fx = 0.0; } else { *x = x_vec[*n_data-1]; *fx = fx_vec[*n_data-1]; } return; # undef N_MAX } /******************************************************************************/ int gray_code_check ( int n, int t[] ) /******************************************************************************/ /* Purpose: gray_code_check() checks a Gray code element. Licensing: This code is distributed under the MIT license. Modified: 25 November 2015 Author: John Burkardt Input: int N, the number of digits in each element. N must be positive. int T[N], an element of the Gray code. Each entry T(I) is either 0 or 1. Output: int GRAY_CODE_CHECK. 1, the data is legal. 0, the data is not legal. */ { int check; int i; check = 1; if ( n < 1 ) { check = 0; return check; } for ( i = 0; i < n; i++ ) { if ( t[i] != 0 && t[i] != 1 ) { check = 0; return check; } } return check; } /******************************************************************************/ int gray_code_enum ( int n ) /******************************************************************************/ /* Purpose: gray_code_enum() enumerates the Gray codes on N digits. Licensing: This code is distributed under the MIT license. Modified: 24 July 2011 Author: John Burkardt Input: int N, the number of digits in each element. N must be nonnegative. Output: int GRAY_CODE_ENUM, the number of distinct elements. */ { int value; value = i4_power ( 2, n ); return value; } /******************************************************************************/ int *gray_code_random ( int n ) /******************************************************************************/ /* Purpose: gray_code_random() randomly selects a Gray code of N digits. Licensing: This code is distributed under the MIT license. Modified: 11 September 2022 Author: John Burkardt Reference: Donald Kreher, Douglas Simpson, Combinatorial Algorithms, CRC Press, 1998, ISBN: 0-8493-3988-X, LC: QA164.K73. Input: int N, the number of digits in each element. N must be positive. Output: int GRAY_CODE_RANDOM[N], the random Gray code. */ { int gray_num; int rank; int *t; /* Check. */ if ( n < 1 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "gray_code_random(): Fatal error!\n" ); fprintf ( stderr, " Input N is illegal.\n" ); exit ( 1 ); } /* Compute GRAY_NUM, the number of Gray codes on N digits. */ gray_num = gray_code_enum ( n ); /* Choose RANK between 1 and GRAY_NUM. */ rank = 1 + ( rand ( ) % gray_num ); /* Compute the Gray code of given RANK. */ t = gray_code_unrank ( rank, n ); return t; } /******************************************************************************/ int gray_code_rank ( int n, int t[] ) /******************************************************************************/ /* Purpose: gray_code_rank() computes the rank of a Gray code element. Licensing: This code is distributed under the MIT license. Modified: 24 July 2011 Author: John Burkardt Reference: Donald Kreher, Douglas Simpson, Combinatorial Algorithms, CRC Press, 1998, ISBN: 0-8493-3988-X, LC: QA164.K73. Input: int N, the number of digits in each element. N must be positive. int T[N], an element of the Gray code. Each entry is either 0 or 1. Output: int GRAY_CODE_RANK, the rank of the element. */ { int b; int check; int i; int rank; /* Check. */ check = gray_code_check ( n, t ); if ( ! check ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "GRAY_CODE_RANK(): Fatal error!\n" ); fprintf ( stderr, " The input array is illegal.\n" ); exit ( 1 ); } rank = 0; b = 0; for ( i = n - 1; 0 <= i; i-- ) { if ( t[n-i-1] != 0 ) { b = 1 - b; } if ( b == 1 ) { rank = rank + i4_power ( 2, i ); } } return rank; } /******************************************************************************/ void gray_code_successor ( int n, int t[], int *rank ) /******************************************************************************/ /* Purpose: gray_code_successor() computes the binary reflected Gray code successor. Example: 000, 001, 011, 010, 110, 111, 101, 100, after which the sequence repeats. Discussion: In the original code, the successor of the element that has an initial 1 followed by N-1 zeroes is undefined. In this version, the successor is the element with N zeroes. Licensing: This code is distributed under the MIT license. Modified: 26 July 2011 Author: John Burkardt Reference: Donald Kreher, Douglas Simpson, Combinatorial Algorithms, CRC Press, 1998, ISBN: 0-8493-3988-X, LC: QA164.K73. Parameters: Input, int N, the number of digits in each element. N must be positive. Input/output, int T[N]. On input, T contains an element of the Gray code, that is, each entry T(I) is either 0 or 1. On output, T contains the successor to the input value; this is an element of the Gray code, which differs from the input value in a single position. Input/output, int *RANK, the rank. If RANK = -1 on input, then the routine understands that this is the first call, and that the user wishes the routine to supply the first element in the ordering, which has RANK = 0. In general, the input value of RANK is increased by 1 for output, unless the very last element of the ordering was input, in which case the output value of RANK is 0. */ { int check; int i; int weight; /* Return the first element. */ if ( *rank == -1 ) { for ( i = 0; i < n; i++ ) { t[i] = 0; } *rank = 0; return; } /* Check. */ check = gray_code_check ( n, t ); if ( ! check ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "GRAY_CODE_SUCCESSOR(): Fatal error!\n" ); fprintf ( stderr, " The input array is illegal.\n" ); exit ( 1 ); } weight = i4vec_sum ( n, t ); if ( ( weight % 2 ) == 0 ) { if ( t[n-1] == 0 ) { t[n-1] = 1; } else { t[n-1] = 0; } *rank = *rank + 1; return; } else { for ( i = n - 1; 1 <= i; i-- ) { if ( t[i] == 1 ) { if ( t[i-1] == 0 ) { t[i-1] = 1; } else { t[i-1] = 0; } *rank = *rank + 1; return; } } /* The final element was input. Return the first element. */ for ( i = 0; i < n; i++ ) { t[i] = 0; } *rank = 0; } return; } /******************************************************************************/ int *gray_code_unrank ( int rank, int n ) /******************************************************************************/ /* Purpose: gray_code_unrank() computes the Gray code element of given rank. Licensing: This code is distributed under the MIT license. Modified: 25 July 2011 Author: John Burkardt Reference: Donald Kreher, Douglas Simpson, Combinatorial Algorithms, CRC Press, 1998, ISBN: 0-8493-3988-X, LC: QA164.K73. Input: int RANK, the rank of the element. 0 <= RANK <= 2^N. int N, the number of digits in each element. N must be positive. Output: int GRAY_CODE_UNRANK[N], the element of the Gray code which has the given rank. */ { int b; int bprime; int i; int ngray; int rank_copy; int *t; /* Check. */ if ( n < 1 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "gray_code_unrank(): Fatal error!\n" ); fprintf ( stderr, " Input N is illegal.\n" ); exit ( 1 ); } ngray = gray_code_enum ( n ); if ( rank < 0 || ngray < rank ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "gray_code_unrank(): Fatal error!\n" ); fprintf ( stderr, " The input rank is illegal.\n" ); exit ( 1 ); } t = ( int * ) malloc ( n * sizeof ( int ) ); rank_copy = rank; for ( i = 0; i < n; i++ ) { t[i] = 0; } bprime = 0; for ( i = n - 1; 0 <= i; i-- ) { b = rank_copy / i4_power ( 2, i ); if ( b != bprime ) { t[n-i-1] = 1; } bprime = b; rank_copy = rank_copy - b * i4_power ( 2, i ); } return t; } /******************************************************************************/ int i4_choose ( int n, int k ) /******************************************************************************/ /* Purpose: i4_choose() computes the binomial coefficient C(N,K). Discussion: The value is calculated in such a way as to avoid overflow and roundoff. The calculation is done in integer arithmetic. The formula used is: C(N,K) = N! / ( K! * (N-K)! ) Licensing: This code is distributed under the MIT license. Modified: 09 December 2013 Author: John Burkardt Reference: ML Wolfson, HV Wright, Algorithm 160: Combinatorial of M Things Taken N at a Time, Communications of the ACM, Volume 6, Number 4, April 1963, page 161. Parameters: Input, int N, K, are the values of N and K. Output, int I4_CHOOSE, the number of combinations of N things taken K at a time. */ { int i; int mn; int mx; int value; if ( k < n - k ) { mn = k; mx = n - k; } else { mn = n - k; mx = k; } if ( mn < 0 ) { value = 0; } else if ( mn == 0 ) { value = 1; } else { value = mx + 1; for ( i = 2; i <= mn; i++ ) { value = ( value * ( mx + i ) ) / i; } } return value; } /******************************************************************************/ int i4_factorial ( int n ) /******************************************************************************/ /* Purpose: i4_factorial() computes the factorial of N. Discussion: factorial ( N ) = product ( 1 <= I <= N ) I Licensing: This code is distributed under the MIT license. Modified: 26 June 2008 Author: John Burkardt Parameters: Input, int N, the argument of the factorial function. If N is less than 1, the function value is returned as 1. 0 <= N <= 13 is required. Output, int I4_FACTORIAL, the factorial of N. */ { int i; int value; value = 1; if ( 13 < n ) { fprintf ( stderr, "I4_FACTORIAL(): Fatal error!\n" ); fprintf ( stderr, " I4_FACTORIAL(N) cannot be computed as an integer\n" ); fprintf ( stderr, " for 13 < N.\n" ); fprintf ( stderr, " Input value N = %d\n", n ); exit ( 1 ); } for ( i = 1; i <= n; i++ ) { value = value * i; } return value; } /******************************************************************************/ void i4_factorial_values ( int *n_data, int *n, int *fn ) /******************************************************************************/ /* Purpose: i4_factorial_values() returns values of the factorial function. Discussion: 0! = 1 I! = Product ( 1 <= J <= I ) I In Mathematica, the function can be evaluated by: n! Licensing: This code is distributed under the MIT license. Modified: 27 July 2011 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 argument of the function. Output, int *FN, the value of the function. */ { # define N_MAX 13 static int fn_vec[N_MAX] = { 1, 1, 2, 6, 24, 120, 720, 5040, 40320, 362880, 3628800, 39916800, 479001600 }; static int n_vec[N_MAX] = { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12 }; if ( *n_data < 0 ) { *n_data = 0; } *n_data = *n_data + 1; if ( N_MAX < *n_data ) { *n_data = 0; *n = 0; *fn = 0; } else { *n = n_vec[*n_data-1]; *fn = fn_vec[*n_data-1]; } return; # undef N_MAX } /******************************************************************************/ int i4_fall ( int x, int n ) /******************************************************************************/ /* Purpose: i4_fall() computes the falling factorial function [X]_N. Discussion: Note that the number of "injections" or 1-to-1 mappings from a set of N elements to a set of M elements is [M]_N. The number of permutations of N objects out of M is [M]_N. Moreover, the Stirling numbers of the first kind can be used to convert a falling factorial into a polynomial, as follows: [X]_N = S^0_N + S^1_N * X + S^2_N * X^2 + ... + S^N_N X^N. Formula: [X]_N = X * ( X - 1 ) * ( X - 2 ) * ... * ( X - N + 1 ). Licensing: This code is distributed under the MIT license. Modified: 08 May 2003 Author: John Burkardt Parameters: Input, int X, the argument of the falling factorial function. Input, int N, the order of the falling factorial function. If N = 0, FALL = 1, if N = 1, FALL = X. Note that if N is negative, a "rising" factorial will be computed. Output, int I4_FALL, the value of the falling factorial function. */ { int i; int value; value = 1; if ( 0 < n ) { for ( i = 1; i <= n; i++ ) { value = value * x; x = x - 1; } } else if ( n < 0 ) { for ( i = -1; n <= i; i-- ) { value = value * x; x = x + 1; } } return value; } /******************************************************************************/ void i4_fall_values ( int *n_data, int *m, int *n, int *fmn ) /******************************************************************************/ /* Purpose: i4_fall_values() returns values of the integer falling factorial function. Discussion: The definition of the falling factorial function is (m)_n = (m)! / (m-n)! = ( m ) * ( m - 1 ) * ( m - 2 ) ... * ( m - n + 1 ) = Gamma ( m + 1 ) / Gamma ( m - n + 1 ) We assume 0 <= N <= M. In Mathematica, the function can be evaluated by: FactorialPower[m,n] Licensing: This code is distributed under the MIT license. Modified: 14 December 2014 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 *M, N, the arguments of the function. Output, int *FMN, the value of the function. */ { # define N_MAX 15 static int fmn_vec[N_MAX] = { 1, 5, 20, 60, 120, 120, 0, 1, 10, 4000, 90, 4896, 24, 912576, 0 }; static int m_vec[N_MAX] = { 5, 5, 5, 5, 5, 5, 5, 50, 10, 4000, 10, 18, 4, 98, 1 }; static int n_vec[N_MAX] = { 0, 1, 2, 3, 4, 5, 6, 0, 1, 1, 2, 3, 4, 3, 7 }; if ( *n_data < 0 ) { *n_data = 0; } *n_data = *n_data + 1; if ( N_MAX < *n_data ) { *n_data = 0; *m = 0; *n = 0; *fmn = 0; } else { *m = m_vec[*n_data-1]; *n = n_vec[*n_data-1]; *fmn = fmn_vec[*n_data-1]; } return; # undef N_MAX } /******************************************************************************/ int i4_huge ( ) /******************************************************************************/ /* Purpose: i4_huge() returns a "huge" I4. Licensing: This code is distributed under the MIT license. Modified: 29 August 2006 Author: John Burkardt Parameters: Output, int I4_HUGE, a "huge" integer. */ { static int value = 2147483647; return value; } /******************************************************************************/ int i4_max ( int i1, int i2 ) /******************************************************************************/ /* Purpose: i4_max() returns the maximum of two I4's. Licensing: This code is distributed under the MIT 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 MIT 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 MIT 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. Discussion: The pseudorandom number should be uniformly distributed between A and B. Licensing: This code is distributed under the MIT 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; } /******************************************************************************/ int *i4mat_copy_new ( int m, int n, int a1[] ) /******************************************************************************/ /* Purpose: i4mat_copy_new() copies an I4MAT to a "new" 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 MIT license. Modified: 27 August 2008 Author: John Burkardt Parameters: Input, int M, N, the number of rows and columns. Input, int A1[M*N], the matrix to be copied. Output, int I4MAT_COPY_NEW[M*N], the copy of A1. */ { int *a2; int i; int j; a2 = ( int * ) malloc ( m * n * sizeof ( int ) ); for ( j = 0; j < n; j++ ) { for ( i = 0; i < m; i++ ) { a2[i+j*m] = a1[i+j*m]; } } return a2; } /******************************************************************************/ 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 MIT 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 MIT 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_backtrack ( int n, int maxstack, int stack[], int x[], int *indx, int *k, int *nstack, int ncan[] ) /******************************************************************************/ /* Purpose: I4VEC_BACKTRACK supervises a backtrack search for an I4VEC. Discussion: The routine tries to construct an integer vector one index at a time, using possible candidates as supplied by the user. At any time, the partially constructed vector may be discovered to be unsatisfactory, but the routine records information about where the last arbitrary choice was made, so that the search can be carried out efficiently, rather than starting out all over again. First, call the routine with INDX = 0 so it can initialize itself. Now, on each return from the routine, if INDX is: 1, you've just been handed a complete candidate vector; Admire it, analyze it, do what you like. 2, please determine suitable candidates for position X(K). Return the number of candidates in NCAN(K), adding each candidate to the end of STACK, and increasing NSTACK. 3, you're done. Stop calling the routine; Licensing: This code is distributed under the MIT license. Modified: 13 July 2004 Author: Original FORTRAN77 version by Albert Nijenhuis, Herbert Wilf. C++ 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, int N, the number of positions to be filled in the vector. Input, int MAXSTACK, the maximum length of the stack. Input, int STACK[MAXSTACK], a list of all current candidates for all positions 1 through K. Input/output, int X[N], the partial or complete candidate vector. Input/output, int *INDX, a communication flag. On input, 0 to start a search. On output: 1, a complete output vector has been determined and returned in X(1:N); 2, candidates are needed for position X(K); 3, no more possible vectors exist. Input/output, int *K, if INDX=2, the current vector index being considered. Input/output, int *NSTACK, the current length of the stack. Input/output, int NCAN[N], lists the current number of candidates for positions 1 through K. */ { /* If this is the first call, request a candidate for position 1. */ if ( *indx == 0 ) { *k = 1; *nstack = 0; *indx = 2; return; } /* Examine the stack. */ for ( ; ; ) { /* If there are candidates for position K, take the first available one off the stack, and increment K. This may cause K to reach the desired value of N, in which case we need to signal the user that a complete set of candidates is being returned. */ if ( 0 < ncan[(*k)-1] ) { x[(*k)-1] = stack[(*nstack)-1]; *nstack = *nstack - 1; ncan[(*k)-1] = ncan[(*k)-1] - 1; if ( *k != n ) { *k = *k + 1; *indx = 2; } else { *indx = 1; } break; } /* If there are no candidates for position K, then decrement K. If K is still positive, repeat the examination of the stack. */ else { *k = *k - 1; if ( *k <= 0 ) { *indx = 3; break; } } } return; } /******************************************************************************/ int *i4vec_copy_new ( int n, int a1[] ) /******************************************************************************/ /* Purpose: I4VEC_COPY_NEW copies an I4VEC. Discussion: An I4VEC is a vector of I4's. Licensing: This code is distributed under the MIT license. Modified: 04 July 2008 Author: John Burkardt Parameters: Input, int N, the number of entries in the vectors. Input, int A1[N], the vector to be copied. Output, int I4VEC_COPY_NEW[N], the copy of A1. */ { int *a2; int i; a2 = ( int * ) malloc ( n * sizeof ( int ) ); for ( i = 0; i < n; i++ ) { a2[i] = a1[i]; } return a2; } /******************************************************************************/ int i4vec_dot_product ( int n, int x[], int y[] ) /******************************************************************************/ /* Purpose: I4VEC_DOT_PRODUCT computes the dot product of two I4VEC's. Discussion: An I4VEC is a vector of I4's. Licensing: This code is distributed under the MIT license. Modified: 19 December 2011 Author: John Burkardt Parameters: Input, int N, the size of the array. Input, int X[N], Y[N], the arrays. Output, int I4VEC_DOT_PRODUCT, the dot product of X and Y. */ { int i; int value; value = 0; for ( i = 0; i < n; i++ ) { value = value + x[i] * y[i]; } return value; } /******************************************************************************/ void i4vec_indicator1 ( int n, int a[] ) /******************************************************************************/ /* Purpose: I4VEC_INDICATOR1 sets an I4VEC to the indicator1 vector. Discussion: An I4VEC is a vector of I4's. Licensing: This code is distributed under the MIT license. Modified: 24 March 2009 Author: John Burkardt Parameters: Input, int N, the number of elements of A. Output, int A[N], the initialized array. */ { int i; for ( i = 0; i < n; i++ ) { a[i] = i + 1; } return; } /******************************************************************************/ 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 MIT 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; } /******************************************************************************/ int *i4vec_part1_new ( int n, int npart ) /******************************************************************************/ /* Purpose: I4VEC_PART1_NEW partitions an integer into NPART parts. Example: Input: N = 17, NPART = 5 Output: X = ( 13, 1, 1, 1, 1 ). Licensing: This code is distributed under the MIT license. Modified: 24 July 2011 Author: John Burkardt Parameters: Input, int N, the integer to be partitioned. N may be positive, zero, or negative. Input, int NPART, the number of entries in the array. 1 <= NPART <= N. Output, int I4VEC_PART1_NEW[NPART], the partition of N. The entries of X add up to N. X(1) = N + 1 - NPART, and all other entries are equal to 1. */ { int i; int *x; if ( npart < 1 || n < npart ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "I4VEC_PART1_NEW(): Fatal error!\n" ); fprintf ( stderr, " The input value of NPART is illegal.\n" ); exit ( 1 ); } x = ( int * ) malloc ( npart * sizeof ( int ) ); x[0] = n + 1 - npart; for ( i = 1; i < npart; i++ ) { x[i] = 1; } return x; } /******************************************************************************/ void i4vec_part2 ( int n, int npart, int x[] ) /******************************************************************************/ /* Purpose: I4VEC_PART2 partitions an integer into NPART nearly equal parts. Example: Input: N = 17, NPART = 5 Output: X = ( 4, 4, 3, 3, 3 ). Licensing: This code is distributed under the MIT license. Modified: 13 April 2013 i4vec Author: John Burkardt Parameters: Input, int N, the integer to be partitioned. N may be positive, zero, or negative. Input, int NPART, the number of entries in the array. 1 <= NPART Output, int X[NPART], the partition of N. The entries of X add up to N. The entries of X are either all equal, or differ by at most 1. The entries of X all have the same sign as N, and the "largest" entries occur first. */ { int i; int j; if ( npart < 1 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "I4VEC_PART2(): Fatal error!\n" ); fprintf ( stderr, " The input value of NPART is illegal.\n" ); exit ( 1 ); } for ( i = 0; i < npart; i++ ) { x[i] = 0; } if ( 0 < n ) { j = 1; for ( i = 1; i <= n; i++ ) { x[j-1] = x[j-1] + 1; j = j + 1; if ( npart < j ) { j = 1; } } } else if ( n < 0 ) { j = 1; for ( i = n; i <= -1; i++ ) { x[j-1] = x[j-1] - 1; j = j + 1; if ( npart < j ) { j = 1; } } } return; } /******************************************************************************/ int *i4vec_part2_new ( int n, int npart ) /******************************************************************************/ /* Purpose: I4VEC_PART2_NEW partitions an integer N into NPART nearly equal parts. Example: Input: N = 17, NPART = 5 Output: X = ( 4, 4, 3, 3, 3 ). Licensing: This code is distributed under the MIT license. Modified: 25 July 2011 Author: John Burkardt Parameters: Input, int N, the integer to be partitioned. N may be positive, zero, or negative. Input, int NPART, the number of entries in the array. 1 <= NPART Output, int I4VEC_PART2[NPART], the partition of N. The entries of X add up to N. The entries of X are either all equal, or differ by at most 1. The entries of X all have the same sign as N, and the "largest" entries occur first. */ { int *x; x = ( int * ) malloc ( npart * sizeof ( int ) ); i4vec_part2 ( n, npart, x ); return x; } /******************************************************************************/ 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 MIT 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_reverse ( int n, int a[] ) /******************************************************************************/ /* Purpose: I4VEC_REVERSE reverses the elements of an I4VEC. Discussion: An I4VEC is a vector of I4's. Example: Input: N = 5, A = ( 11, 12, 13, 14, 15 ). Output: A = ( 15, 14, 13, 12, 11 ). Licensing: This code is distributed under the MIT license. Modified: 30 June 2009 Author: John Burkardt Parameters: Input, int N, the number of entries in the array. Input/output, int A[N], the array to be reversed. */ { int i; int j; for ( i = 0; i < n / 2; i++ ) { j = a[i]; a[i] = a[n-1-i]; a[n-1-i] = j; } return; } /******************************************************************************/ int i4vec_search_binary_a ( int n, int a[], int b ) /******************************************************************************/ /* Purpose: I4VEC_SEARCH_BINARY_A searches an ascending sorted I4VEC for a value. Discussion: An I4VEC is a vector of I4's. Binary search is used. Licensing: This code is distributed under the MIT license. Modified: 06 September 2008 Author: John Burkardt Reference: Donald Kreher, Douglas Simpson, Algorithm 1.9, Combinatorial Algorithms, CRC Press, 1998, page 26. Parameters: Input, int N, the number of elements in the vector. Input, int A[N], the array to be searched. A must be sorted in ascending order. Input, int B, the value to be searched for. Output, int I4VEC_SEARCH_BINARY_A, the result of the search. -1, B does not occur in A. I, A[I] = B. */ { int high; int index; int low; int mid; /* Check. */ if ( n <= 0 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "I4VEC_SEARCH_BINARY_A(): Fatal error!\n" ); fprintf ( stderr, " The array dimension N is less than 1.\n" ); exit ( 1 ); } index = -1; low = 1; high = n; while ( low <= high ) { mid = ( low + high ) / 2; if ( a[mid-1] == b ) { index = mid; break; } else if ( a[mid-1] < b ) { low = mid + 1; } else if ( b < a[mid-1] ) { high = mid - 1; } } return index; } /******************************************************************************/ int i4vec_search_binary_d ( int n, int a[], int b ) /******************************************************************************/ /* Purpose: I4VEC_SEARCH_BINARY_D searches a descending sorted I4VEC for a value. Discussion: An I4VEC is a vector of I4's. Binary search is used. Licensing: This code is distributed under the MIT license. Modified: 06 June 2010 Author: John Burkardt Reference: Donald Kreher, Douglas Simpson, Algorithm 1.9, Combinatorial Algorithms, CRC Press, 1998, page 26. Parameters: Input, int N, the number of elements in the vector. Input, int A[N], the array to be searched. A must be sorted in descending order. Input, int B, the value to be searched for. Output, int I4VEC_SEARCH_BINARY_D, the result of the search. -1, B does not occur in A. I, A[I] = B. */ { int high; int index; int low; int mid; /* Check. */ if ( n <= 0 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "I4VEC_SEARCH_BINARY_D(): Fatal error!\n" ); fprintf ( stderr, " The array dimension N is less than 1.\n" ); exit ( 1 ); } index = -1; low = 1; high = n; while ( low <= high ) { mid = ( low + high ) / 2; if ( a[mid-1] == b ) { index = mid; break; } else if ( b < a[mid-1] ) { low = mid + 1; } else if ( a[mid-1] < b ) { high = mid - 1; } } return index; } /******************************************************************************/ void i4vec_sort_insert_a ( int n, int a[] ) /******************************************************************************/ /* Purpose: I4VEC_SORT_INSERT_A uses an ascending insertion sort on an I4VEC. Discussion: An I4VEC is a vector of I4's. Licensing: This code is distributed under the MIT license. Modified: 25 July 2010 Author: John Burkardt Reference: Donald Kreher, Douglas Simpson, Algorithm 1.1, Combinatorial Algorithms, CRC Press, 1998, page 11. Parameters: Input, int N, the number of items in the vector. N must be positive. Input/output, int A[N]. On input, A contains data to be sorted. On output, the entries of A have been sorted in ascending order. */ { int i; int j; int x; for ( i = 1; i < n; i++ ) { x = a[i]; j = i; while ( 1 <= j && x < a[j-1] ) { a[j] = a[j-1]; j = j - 1; } a[j] = x; } return; } /******************************************************************************/ void i4vec_sort_insert_d ( int n, int a[] ) /******************************************************************************/ /* Purpose: I4VEC_SORT_INSERT_D uses a descending insertion sort on an I4VEC. Discussion: An I4VEC is a vector of I4's. Licensing: This code is distributed under the MIT license. Modified: 25 July 2010 Author: John Burkardt Reference: Donald Kreher, Douglas Simpson, Algorithm 1.1, Combinatorial Algorithms, CRC Press, 1998, page 11. Parameters: Input, int N, the number of items in the vector. N must be positive. Input/output, int A[N]. On input, A contains data to be sorted. On output, the entries of A have been sorted in ascending order. */ { int i; int j; int x; for ( i = 1; i < n; i++ ) { x = a[i]; j = i; while ( 1 <= j && a[j-1] < x ) { a[j] = a[j-1]; j = j - 1; } a[j] = x; } return; } /******************************************************************************/ int i4vec_sum ( int n, int a[] ) /******************************************************************************/ /* Purpose: I4VEC_SUM sums the entries of an I4VEC. Discussion: An I4VEC is a vector of I4's. Example: Input: A = ( 1, 2, 3, 4 ) Output: I4VEC_SUM = 10 Licensing: This code is distributed under the MIT license. Modified: 29 May 2003 Author: John Burkardt Parameters: Input, int N, the number of entries in the vector. Input, int A[N], the vector to be summed. Output, int I4VEC_SUM, the sum of the entries of A. */ { int i; int sum; sum = 0; for ( i = 0; i < n; i++ ) { sum = sum + a[i]; } return sum; } /******************************************************************************/ void i4vec_transpose_print ( int n, int a[], char *title ) /******************************************************************************/ /* Purpose: I4VEC_TRANSPOSE_PRINT prints an I4VEC "transposed". Discussion: An I4VEC is a vector of I4's. Example: A = { 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11 } TITLE = "My vector: " My vector: 1 2 3 4 5 6 7 8 9 10 11 Licensing: This code is distributed under the MIT license. Modified: 25 July 2010 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; int ihi; int ilo; int title_len; title_len = strlen ( title ); for ( ilo = 1; ilo <= n; ilo = ilo + 10 ) { if ( 0 < title_len ) { if ( ilo == 1 ) { printf ( "%s", title ); } else { for ( i = 1; i <= title_len; i++ ) { printf ( " " ); } } } ihi = i4_min ( ilo + 10 - 1, n ); for ( i = ilo; i <= ihi; i++ ) { printf ( " %6d", a[i-1] ); } printf ( "\n" ); } return; } /******************************************************************************/ int *i4vec_uniform_ab_new ( int n, int a, int b, int *seed ) /******************************************************************************/ /* Purpose: I4VEC_UNIFORM_AB_NEW returns a scaled pseudorandom I4VEC. Discussion: The pseudorandom numbers should be uniformly distributed between A and B. Licensing: This code is distributed under the MIT license. Modified: 13 January 2014 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, integer N, the dimension of the vector. 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 I4VEC_UNIFORM_AB_NEW[N], a vector of random values between A and B. */ { int c; int i; const int i4_huge = 2147483647; int k; float r; int value; int *x; if ( *seed == 0 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "I4VEC_UNIFORM_AB_NEW(): 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; } x = ( int * ) malloc ( n * sizeof ( int ) ); for ( i = 0; i < n; i++ ) { 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 ); /* Use rounding to convert R to an integer between A and B. */ value = round ( r ); /* Guarantee A <= VALUE <= B. */ if ( value < a ) { value = a; } if ( b < value ) { value = b; } x[i] = value; } return x; } /******************************************************************************/ void knapsack_01 ( int n, double mass_limit, double p[], double w[], double x[], double *mass, double *profit ) /******************************************************************************/ /* Purpose: KNAPSACK_01 solves the 0/1 knapsack problem. Discussion: The 0/1 knapsack problem is as follows: Given: a set of N objects, a profit P(I) and weight W(I) associated with each object, and a weight limit MASS_LIMIT, Determine: a set of choices X(I) which are 0 or 1, that maximizes the profit P = Sum ( 1 <= I <= N ) P(I) * X(I) subject to the constraint Sum ( 1 <= I <= N ) W(I) * X(I) <= MASS_LIMIT. This routine assumes that the objects have already been sorted in order of decreasing "profit density", P(I)/W(I). Licensing: This code is distributed under the MIT license. Modified: 28 July 2011 Author: John Burkardt Reference: Donald Kreher, Douglas Simpson, Combinatorial Algorithms, CRC Press, 1998, ISBN: 0-8493-3988-X, LC: QA164.K73. Parameters: Input, int N, the number of objects. Input, double MASS_LIMIT, the weight limit of the chosen objects. Input/output, double P[N], the "profit" or value of each object. P is assumed to be nonnegative. Input/output, double W[N], the "weight" or cost of each object. W is assumed to be nonnegative. Output, double X[N], the choice function for the objects. 0, the object was not taken. 1, the object was taken. Output, double *MASS, the total mass of the objects taken. Output, double *PROFIT, the total profit of the objects taken. */ { int i; int indx; int k; double mass_1; double mass_2; double mass_best; double mass_remaining; int maxstack = 100; int *ncan; int nstack; double profit_1; double profit_2; double profit_best; double *stack; double *x_best; ncan = ( int * ) malloc ( n * sizeof ( int ) ); stack = ( double * ) malloc ( maxstack * sizeof ( double ) ); x_best = ( double * ) malloc ( n * sizeof ( double ) ); nstack = 0; /* Initialize the "best so far" data. */ for ( i = 0; i < n; i++ ) { x_best[i] = 0.0; } profit_best = 0.0; mass_best = 0; /* Begin the backtracking solution. */ indx = 0; for ( ; ; ) { r8vec_backtrack ( n, maxstack, stack, x, &indx, &k, &nstack, ncan ); /* Got a new candidate. Compare it to the best so far. */ if ( indx == 1 ) { *profit = r8vec_dot_product ( n, p, x ); *mass = r8vec_dot_product ( n, w, x ); if ( profit_best < *profit || ( *profit == profit_best && *mass < mass_best ) ) { profit_best = *profit; mass_best = *mass; for ( i = 0; i < n; i++ ) { x_best[i] = x[i]; } } } /* Need candidates for X(K). X(K) = 1 is possible if: * adding W(K) to our mass doesn''t put us over our mass limit; * and adding P(K) to our current profit, and taking the best we could get using rational X for the remainder would put us over our current best. X(K) = 0 is always possible. */ else if ( indx == 2 ) { ncan[k-1] = 0; mass_1 = w[k-1]; for ( i = 0; i < k - 1; i++ ) { mass_1 = mass_1 + w[i] * x[i]; } if ( mass_1 <= mass_limit ) { mass_remaining = mass_limit - mass_1; profit_1 = p[k-1]; for ( i = 0; i < k - 1; i++ ) { profit_1 = profit_1 + p[i] * x[i]; } if ( k < n ) { knapsack_rational ( n - k, mass_remaining, p+k, w+k, x+k, &mass_2, &profit_2 ); } else { profit_2 = 0.0; } if ( profit_best < profit_1 + profit_2 ) { if ( maxstack <= nstack ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "KNAPSACK_01(): Fatal error!\n" ); fprintf ( stderr, " Exceeded stack space.\n" ); exit ( 1 ); } ncan[k-1] = ncan[k-1] + 1; nstack = nstack + 1; stack[nstack-1] = 1.0; } } if ( maxstack <= nstack ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "KNAPSACK_01(): Fatal error!\n" ); fprintf ( stderr, " Exceeded stack space.\n" ); exit ( 1 ); } ncan[k-1] = ncan[k-1] + 1; nstack = nstack + 1; stack[nstack-1] = 0.0; } /* Done. Return the best solution. */ else { *profit = profit_best; *mass = mass_best; for ( i = 0; i < n; i++ ) { x[i] = x_best[i]; } break; } } free ( ncan ); free ( stack ); free ( x_best ); return; } /******************************************************************************/ void knapsack_rational ( int n, double mass_limit, double p[], double w[], double x[], double *mass, double *profit ) /******************************************************************************/ /* Purpose: KNAPSACK_RATIONAL solves the rational knapsack problem. Discussion: The rational knapsack problem is a generalization of the 0/1 knapsack problem. It is mainly used to derive a bounding function for the 0/1 knapsack problem. The 0/1 knapsack problem is as follows: Given: a set of N objects, a profit P(I) and weight W(I) associated with each object, and a weight limit MASS_LIMIT, Determine: a set of choices X(I) which are 0 or 1, that maximizes the profit P = Sum ( 1 <= I <= N ) P(I) * X(I) subject to the constraint Sum ( 1 <= I <= N ) W(I) * X(I) <= MASS_LIMIT. By contrast, the rational knapsack problem allows the values X(I) to be any value between 0 and 1. A solution for the rational knapsack problem is known. Arrange the objects in order of their "profit density" ratios P(I)/W(I), and then take in order as many of these as you can. If you still have "room" in the weight constraint, then you should take the maximal fraction of the very next object, which will complete your weight limit, and maximize your profit. If should be obvious that, given the same data, a solution for the rational knapsack problem will always have a profit that is at least as high as for the 0/1 problem. Since the rational knapsack maximum profit is easily computed, this makes it a useful bounding function. Note that this routine assumes that the objects have already been arranged in order of the "profit density". Licensing: This code is distributed under the MIT license. Modified: 25 July 2011 Author: John Burkardt Reference: Donald Kreher, Douglas Simpson, Combinatorial Algorithms, CRC Press, 1998, ISBN: 0-8493-3988-X, LC: QA164.K73. Parameters: Input, int N, the number of objects. Input, double MASS_LIMIT, the weight limit of the chosen objects. Input, double P[N], the "profit" or value of each object. The entries of P are assumed to be nonnegative. Input, double W[N], the "weight" or cost of each object. The entries of W are assumed to be nonnegative. Output, double X[N], the choice function for the objects. 0.0, the object was not taken. 1.0, the object was taken. R, where 0 < R < 1, a fractional amount of the object was taken. Output, double *MASS, the total mass of the objects taken. Output, double *PROFIT, the total profit of the objects taken. */ { int i; *mass = 0.0; *profit = 0.0; for ( i = 0; i < n; i++ ) { if ( mass_limit <= *mass ) { x[i] = 0.0; } else if ( *mass + w[i] <= mass_limit ) { x[i] = 1.0; *mass = *mass + w[i]; *profit = *profit + p[i]; } else { x[i] = ( mass_limit - *mass ) / w[i]; *mass = mass_limit; *profit = *profit + p[i] * x[i]; } } return; } /******************************************************************************/ void knapsack_reorder ( int n, double p[], double w[] ) /******************************************************************************/ /* Purpose: KNAPSACK_REORDER reorders the knapsack data by "profit density". Discussion: This routine must be called to rearrange the data before calling routines that handle a knapsack problem. The "profit density" for object I is defined as P(I)/W(I). Licensing: This code is distributed under the MIT license. Modified: 26 July 2011 Author: John Burkardt Reference: Donald Kreher, Douglas Simpson, Combinatorial Algorithms, CRC Press, 1998, ISBN: 0-8493-3988-X, LC: QA164.K73. Parameters: Input, int N, the number of objects. Input/output, double P[N], the "profit" or value of each object. Input/output, double W[N], the "weight" or cost of each object. */ { int i; int j; double t; /* Rearrange the objects in order of "profit density". */ for ( i = 0; i < n; i++ ) { for ( j = i + 1; j < n; j++ ) { if ( p[i] * w[j] < p[j] * w[i] ) { t = p[i]; p[i] = p[j]; p[j] = t; t = w[i]; w[i] = w[j]; w[j] = t; } } } return; } /******************************************************************************/ int ksubset_colex_check ( int k, int n, int t[] ) /******************************************************************************/ /* Purpose: KSUBSET_COLEX_CHECK checks a K subset in colex form. Licensing: This code is distributed under the MIT license. Modified: 24 July 2011 Author: John Burkardt Reference: Donald Kreher, Douglas Simpson, Combinatorial Algorithms, CRC Press, 1998, ISBN: 0-8493-3988-X, LC: QA164.K73. Parameters: Input, int K, the number of elements each K subset must have. 0 <= K <= N. Input, int N, the number of elements in the master set. 0 <= N. Input, int T[K], describes a K subset. T(I) is the I-th element of the K subset. The elements must be listed in DESCENDING order. Output, int KSUBSET_COLEX_CHECK. 1, the data is legal. 0, the data is not legal. */ { int check; int i; int tmax; check = 1; if ( n < 0 ) { check = 0; return check; } if ( k < 0 || n < k ) { check = 0; return check; } tmax = n + 1; for ( i = 0; i < k; i++ ) { if ( t[i] <= 0 || tmax <= t[i] ) { check = 0; return check; } tmax = t[i]; } return check; } /******************************************************************************/ int ksubset_colex_rank ( int k, int n, int t[] ) /******************************************************************************/ /* Purpose: KSUBSET_COLEX_RANK computes the colex rank of a K subset. Licensing: This code is distributed under the MIT license. Modified: 26 July 2011 Author: John Burkardt Reference: Donald Kreher, Douglas Simpson, Combinatorial Algorithms, CRC Press, 1998, ISBN: 0-8493-3988-X, LC: QA164.K73. Parameters: Input, int K, the number of elements each K subset must have. 1 <= K <= N. Input, int N, the number of elements in the master set. N must be positive. Input, int T[K[, describes a K subset. T(I) is the I-th element of the K subset. The elements must be listed in DESCENDING order. Output, int KSUBSET_COLEX_RANK, the rank of the subset. */ { int check; int i; int rank; /* Check. */ check = ksubset_colex_check ( k, n, t ); if ( ! check ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "KSUBSET_COLEX_RANK(): Fatal error!\n" ); fprintf ( stderr, " The input array is illegal.\n" ); exit ( 1 ); } rank = 0; for ( i = 0; i < k; i++ ) { rank = rank + i4_choose ( t[i] - 1, k - i ); } return rank; } /******************************************************************************/ void ksubset_colex_successor ( int k, int n, int t[], int *rank ) /******************************************************************************/ /* Purpose: KSUBSET_COLEX_SUCCESSOR computes the K subset colex successor. Discussion: In the original code, there is a last element with no successor. Licensing: This code is distributed under the MIT license. Modified: 26 July 2011 Author: John Burkardt Reference: Donald Kreher, Douglas Simpson, Combinatorial Algorithms, CRC Press, 1998, ISBN: 0-8493-3988-X, LC: QA164.K73. Parameters: Input, int K, the number of elements each K subset must have. 1 <= K <= N. Input, int N, the number of elements in the master set. N must be positive. Input/output, int T[K], describes a K subset. T(I) is the I-th element. The elements must be listed in DESCENDING order. On input, T describes a K subset. On output, T describes the next K subset in the ordering. If the input T was the last in the ordering, then the output T will be the first. Input/output, int *RANK, the rank. If RANK = -1 on input, then the routine understands that this is the first call, and that the user wishes the routine to supply the first element in the ordering, which has RANK = 0. In general, the input value of RANK is increased by 1 for output, unless the very last element of the ordering was input, in which case the output value of RANK is 0. */ { int check; int i; /* Return the first element. */ if ( *rank == -1 ) { for ( i = 1; i <= k; i++ ) { t[i-1] = k + 1 - i; } *rank = 0; return; } /* Check. */ check = ksubset_colex_check ( k, n, t ); if ( ! check ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "KSUBSET_COLEX_SUCCESSOR(): Fatal error!\n" ); fprintf ( stderr, " The input array is illegal.\n" ); exit ( 1 ); } for ( i = k - 1; 1 <= i; i-- ) { if ( t[k-i] + 1 < t[k-i-1] ) { t[k-i] = t[k-i] + 1; *rank = *rank + 1; return; } } if ( t[0] < n ) { t[0] = t[0] + 1; for ( i = 1; i <= k - 1; i++ ) { t[k-i] = i; } *rank = *rank + 1; return; } /* The last K subset was input. Return the first one. */ for ( i = 1; i <= k; i++ ) { t[i-1] = k + 1 - i; } *rank = 0; return; } /******************************************************************************/ int *ksubset_colex_unrank ( int rank, int k, int n ) /******************************************************************************/ /* Purpose: KSUBSET_COLEX_UNRANK computes the K subset of given colex rank. Licensing: This code is distributed under the MIT license. Modified: 25 July 2011 Author: John Burkardt Reference: Donald Kreher, Douglas Simpson, Combinatorial Algorithms, CRC Press, 1998, ISBN: 0-8493-3988-X, LC: QA164.K73. Parameters: Input, int RANK, the rank of the K subset. Input, int K, the number of elements each K subset must have. 0 <= K <= N. Input, int N, the number of elements in the master set. N must be positive. Output, int KSUBSET_COLEX_UNRANK[K], describes the K subset of the given rank. T(I) is the I-th element. The elements must be listed in DESCENDING order. */ { int i; int nksub; int rank_copy; int *t; int x; /* Check. */ if ( n < 1 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "KSUBSET_COLEX_UNRANK(): Fatal error!\n" ); fprintf ( stderr, " Input N is illegal.\n" ); exit ( 1 ); } if ( k == 0 ) { t = ( int * ) malloc ( k * sizeof ( int ) ); return t; } if ( k < 0 || n < k ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "KSUBSET_COLEX_UNRANK(): Fatal error!\n" ); fprintf ( stderr, " Input K is illegal.\n" ); exit ( 1 ); } nksub = ksubset_enum ( k, n ); if ( rank < 0 || nksub < rank ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "KSUBSET_COLEX_UNRANK(): Fatal error!\n" ); fprintf ( stderr, " The input rank is illegal.\n" ); exit ( 1 ); } rank_copy = rank; x = n; t = ( int * ) malloc ( k * sizeof ( int ) ); for ( i = 1; i <= k; i++ ) { while ( rank_copy < i4_choose ( x, k + 1 - i ) ) { x = x - 1; } t[i-1] = x + 1; rank_copy = rank_copy - i4_choose ( x, k + 1 - i ); } return t; } /******************************************************************************/ int ksubset_enum ( int k, int n ) /******************************************************************************/ /* Purpose: KSUBSET_ENUM enumerates the K element subsets of an N set. Licensing: This code is distributed under the MIT license. Modified: 24 July 2011 Author: John Burkardt Parameters: Input, int K, the number of elements each K subset must have. 0 <= K <= N. Input, int N, the number of elements in the master set. 0 <= N. Output, int KSUBSET_ENUM, the number of distinct elements. */ { int value; value = i4_choose ( n, k ); return value; } /******************************************************************************/ int ksubset_lex_check ( int k, int n, int t[] ) /******************************************************************************/ /* Purpose: KSUBSET_LEX_CHECK checks a K subset in lex form. Licensing: This code is distributed under the MIT license. Modified: 26 July 2011 Author: John Burkardt Reference: Donald Kreher, Douglas Simpson, Combinatorial Algorithms, CRC Press, 1998, ISBN: 0-8493-3988-X, LC: QA164.K73. Parameters: Input, int K, the number of elements each K subset must have. 0 <= K <= N. Input, int N, the number of elements in the master set. 0 <= N. Input, int T[K], describes a K subset. T(I) is the I-th element of the K subset. The elements must be listed in ASCENDING order. Output, int KSUBSET_LEX_CHECK. 1, the data is legal. 0, the data is not legal. */ { int check; int i; int tmin; check = 1; if ( n < 0 ) { check = 0; return check; } if ( k < 0 || n < k ) { check = 0; return check; } tmin = 0; for ( i = 0; i < k; i++ ) { if ( t[i] <= tmin || n < t[i] ) { check = 0; return check; } tmin = t[i]; } return check; } /******************************************************************************/ int ksubset_lex_rank ( int k, int n, int t[] ) /******************************************************************************/ /* Purpose: KSUBSET_LEX_RANK computes the lexicographic rank of a K subset. Licensing: This code is distributed under the MIT license. Modified: 25 July 2011 Author: John Burkardt Reference: Donald Kreher, Douglas Simpson, Combinatorial Algorithms, CRC Press, 1998, ISBN: 0-8493-3988-X, LC: QA164.K73. Parameters: Input, int K, the number of elements each K subset must have. 1 <= K <= N. Input, int N, the number of elements in the master set. N must be positive. Input, int T[K], describes a K subset. T(I) is the I-th element. The elements must be listed in ascending order. Output, int KSUBSET_LEX_RANK, the rank of the K subset. */ { int check; int i; int j; int rank; int tim1; /* Check. */ check = ksubset_lex_check ( k, n, t ); if ( ! check ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "KSUBSET_LEX_RANK(): Fatal error!\n" ); fprintf ( stderr, " The input array is illegal.\n" ); exit ( 1 ); } rank = 0; for ( i = 1; i <= k; i++ ) { if ( i == 1 ) { tim1 = 0; } else { tim1 = t[i-2]; } if ( tim1 + 1 <= t[i-1] - 1 ) { for ( j = tim1 + 1; j <= t[i-1] - 1; j++ ) { rank = rank + i4_choose ( n - j, k - i ); } } } return rank; } /******************************************************************************/ void ksubset_lex_successor ( int k, int n, int t[], int *rank ) /******************************************************************************/ /* Purpose: KSUBSET_LEX_SUCCESSOR computes the K subset lexicographic successor. Discussion: In the original code, there is a last element with no successor. Licensing: This code is distributed under the MIT license. Modified: 26 July 2011 Author: John Burkardt Reference: Donald Kreher, Douglas Simpson, Combinatorial Algorithms, CRC Press, 1998, ISBN: 0-8493-3988-X, LC: QA164.K73. Parameters: Input, int K, the number of elements each K subset must have. 1 <= K <= N. Input, int N, the number of elements in the master set. N must be positive. Input/output, int T[K], describes a K subset. T(I) is the I-th element. The elements must be listed in ascending order. On input, T describes a K subset. On output, T describes the next K subset in the ordering. If the input T was the last in the ordering, then the output T will be the first. Input/output, int *RANK, the rank. If RANK = -1 on input, then the routine understands that this is the first call, and that the user wishes the routine to supply the first element in the ordering, which has RANK = 0. In general, the input value of RANK is increased by 1 for output, unless the very last element of the ordering was input, in which case the output value of RANK is 0. */ { int check; int i; int isave; int j; /* Return the first element. */ if ( *rank == -1 ) { i4vec_indicator1 ( k, t ); *rank = 0; return; } /* Check. */ check = ksubset_lex_check ( k, n, t ); if ( ! check ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "KSUBSET_LEX_SUCCESSOR(): Fatal error!\n" ); fprintf ( stderr, " The input array is illegal.\n" ); exit ( 1 ); } isave = 0; for ( i = k; 1 <= i; i-- ) { if ( t[i-1] != n - k + i ) { isave = i; break; } } /* The last K subset was input. Return the first one. */ if ( isave == 0 ) { i4vec_indicator1 ( k, t ); *rank = 0; } else { for ( j = k; isave <= j; j-- ) { t[j-1] = t[isave-1] + 1 + j - isave; } *rank = *rank + 1; } return; } /******************************************************************************/ int *ksubset_lex_unrank ( int rank, int k, int n ) /******************************************************************************/ /* Purpose: KSUBSET_LEX_UNRANK computes the K subset of given lexicographic rank. Licensing: This code is distributed under the MIT license. Modified: 25 July 2011 Author: John Burkardt Reference: Donald Kreher, Douglas Simpson, Combinatorial Algorithms, CRC Press, 1998, ISBN: 0-8493-3988-X, LC: QA164.K73. Parameters: Input, int RANK, the rank of the K subset. Input, int K, the number of elements each K subset must have. 0 <= K <= N. Input, int N, the number of elements in the master set. N must be positive. Output, int KSUBSET_LEX_RANK[K], describes the K subset of the given rank. T(I) is the I-th element. The elements must be listed in ascending order. */ { int i; int nksub; int rank_copy; int *t; int x; /* Check. */ if ( n < 1 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "KSUBSET_LEX_UNRANK(): Fatal error!\n" ); fprintf ( stderr, " Input N is illegal.\n" ); exit ( 1 ); } if ( k == 0 ) { t = ( int * ) malloc ( k * sizeof ( int ) ); return t; } if ( k < 0 || n < k ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "KSUBSET_LEX_UNRANK(): Fatal error!\n" ); fprintf ( stderr, " Input K is illegal.\n" ); exit ( 1 ); } nksub = ksubset_enum ( k, n ); if ( rank < 0 || nksub < rank ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "KSUBSET_LEX_UNRANK(): Fatal error!\n" ); fprintf ( stderr, " Input rank is illegal.\n" ); exit ( 1 ); } t = ( int * ) malloc ( k * sizeof ( int ) ); rank_copy = rank; x = 1; for ( i = 1; i <= k; i++ ) { while ( i4_choose ( n - x, k - i ) <= rank_copy ) { rank_copy = rank_copy - i4_choose ( n - x, k - i ); x = x + 1; } t[i-1] = x; x = x + 1; } return t; } /******************************************************************************/ int ksubset_revdoor_rank ( int k, int n, int t[] ) /******************************************************************************/ /* Purpose: KSUBSET_REVDOOR_RANK computes the revolving door rank of a K subset. Licensing: This code is distributed under the MIT license. Modified: 25 July 2011 Author: John Burkardt Reference: Donald Kreher, Douglas Simpson, Combinatorial Algorithms, CRC Press, 1998, ISBN: 0-8493-3988-X, LC: QA164.K73. Parameters: Input, int K, the number of elements each K subset must have. 1 <= K <= N. Input, int N, the number of elements in the master set. N must be positive. Input, int T[K], describes a K subset. T(I) is the I-th element. The elements must be listed in ascending order. Output, int KSUBSET_REVDOOR_RANK, the rank of the K subset. */ { int check; int i; int rank; int s; /* Check. */ check = ksubset_lex_check ( k, n, t ); if ( ! check ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "KSUBSET_REVDOOR_RANK(): Fatal error!\n" ); fprintf ( stderr, " The input array is illegal.\n" ); exit ( 1 ); } if ( ( k % 2 ) == 0 ) { rank = 0; } else { rank = - 1; } s = 1; for ( i = k; 1 <= i; i-- ) { rank = rank + s * i4_choose ( t[i-1], i ); s = - s; } return rank; } /******************************************************************************/ void ksubset_revdoor_successor ( int k, int n, int t[], int *rank ) /******************************************************************************/ /* Purpose: KSUBSET_REVDOOR_SUCCESSOR computes the K subset revolving door successor. Discussion: After numerous attempts to implement the algorithm published in Kreher and Stinson, the Nijenhuis and Wilf version was implemented instead. The K and S algorithm is supposedly based on the N and W one. Licensing: This code is distributed under the MIT license. Modified: 26 July 2011 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. Donald Kreher, Douglas Simpson, Combinatorial Algorithms, CRC Press, 1998, ISBN: 0-8493-3988-X, LC: QA164.K73. Parameters: Input, int K, the number of elements each K subset must have. 1 <= K <= N. Input, int N, the number of elements in the master set. N must be positive. Input/output, int T[K], describes a K subset. T(I) is the I-th element. The elements must be listed in ascending order. On input, T describes a K subset. On output, T describes the next K subset in the ordering. If the input T was the last in the ordering, then the output T will be the first. Input/output, int *RANK, the rank. If RANK = -1 on input, then the routine understands that this is the first call, and that the user wishes the routine to supply the first element in the ordering, which has RANK = 0. In general, the input value of RANK is increased by 1 for output, unless the very last element of the ordering was input, in which case the output value of RANK is 0. */ { int check; int j; /* Return the first element. */ if ( *rank == - 1 ) { i4vec_indicator1 ( k, t ); *rank = 0; return; } /* Check. */ check = ksubset_lex_check ( k, n, t ); if ( ! check ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "KSUBSET_REVDOOR_SUCCESSOR(): Fatal error!\n" ); fprintf ( stderr, " The input array is illegal.\n" ); exit ( 1 ); } j = 0; for ( ; ; ) { if ( 0 < j || ( k % 2 ) == 0 ) { j = j + 1; if ( k < j ) { t[k-1] = k; *rank = 0; return; } if ( t[j-1] != j ) { t[j-1] = t[j-1] - 1; if ( j != 1 ) { t[j-2] = j - 1; } *rank = *rank + 1; return; } } j = j + 1; if ( j < k ) { if ( t[j-1] != t[j] - 1 ) { break; } } else { if ( t[j-1] != n ) { break; } } } t[j-1] = t[j-1] + 1; if ( j != 1 ) { t[j-2] = t[j-1] - 1; } *rank = *rank + 1; return; } /******************************************************************************/ int *ksubset_revdoor_unrank ( int rank, int k, int n ) /******************************************************************************/ /* Purpose: KSUBSET_REVDOOR_UNRANK computes the K subset of given revolving door rank. Licensing: This code is distributed under the MIT license. Modified: 25 July 2011 Author: John Burkardt Reference: Donald Kreher, Douglas Simpson, Combinatorial Algorithms, CRC Press, 1998, ISBN: 0-8493-3988-X, LC: QA164.K73. Parameters: Input, int RANK, the rank of the K subset. Input, int K, the number of elements each K subset must have. 1 <= K <= N. Input, int N, the number of elements in the master set. N must be positive. Output, int KSUBSET_REVDOOR_UNRANK[K], describes the K subset of the given rank. T(I) is the I-th element. The elements must be listed in ascending order. */ { int i; int nksub; int rank_copy; int *t; int x; /* Check. */ if ( n < 1 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "KSUBSET_REVDOOR_UNRANK(): Fatal error!\n" ); fprintf ( stderr, " Input N is illegal.\n" ); exit ( 1 ); } if ( k < 1 || n < k ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "KSUBSET_REVDOOR_UNRANK(): Fatal error!\n" ); fprintf ( stderr, " Input K is illegal.\n" ); exit ( 1 ); } nksub = ksubset_enum ( k, n ); if ( rank < 0 || nksub < rank ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "KSUBSET_REVDOOR_UNRANK(): Fatal error!\n" ); fprintf ( stderr, " The input rank is illegal.\n" ); exit ( 1 ); } rank_copy = rank; t = ( int * ) malloc ( k * sizeof ( int * ) ); x = n; for ( i = k; 1 <= i; i-- ) { while ( rank_copy < i4_choose ( x, i ) ) { x = x - 1; } t[i-1] = x + 1; rank_copy = i4_choose ( x + 1, i ) - rank_copy - 1; } return t; } /******************************************************************************/ void marriage ( int n, int prefer[], int rank[], int fiancee[], int next[] ) /******************************************************************************/ /* Purpose: MARRIAGE finds a stable set of marriages for given preferences. Discussion: Given a set of N men and N women who must be married in pairs, and information defining the relative rankings that each person assigns to the candidates of the opposite sex, this routine finds a stable set of marriages for them. A stable set of marriages is a pairing of the men and women with the stability property: if M1 marries W1 and M2 marries W2, then it is never the case that M1 and W2 would both prefer to be married to each other. An important application of stable marriage algorithms occurs in the annual matching of medical residents to hospitals. Licensing: This code is distributed under the MIT license. Modified: 28 July 2011 Author: John Burkardt Reference: Robert Sedgewick, Algorithms in C, Addison-Wesley, 1990, ISBN: 0-201-51425-7, LC: QA76.73.C15S43. Parameters: Input, int N, the number of pairs of men and women. Input, int PREFER[N*N]; for man I, the value of PREFER(I,J) represents his J-th preference for a wife. Input, int RANK[N*N]; for woman I, the value of RANK(I,J) represents her ranking of man number J. A value of 1 for RANK(I,J) means woman I ranks man J most preferable, while a value of N would mean she ranked him least preferable. Output, int FIANCEE[N]; for woman I, FIANCEE(I) is the man to whom she is now engaged. Output, int NEXT[N]; for man I, NEXT(I) is his preference ranking for the woman to whom he is now engaged. A value of 1 represents his first choice, a value of N his last. */ { int i; int m; int temp; int w; /* For man I, NEXT(I) is the woman I has most recently proposed to, and hence NEXT(I)+1 is the next one to try. */ for ( i = 0; i < n; i++ ) { next[i] = 0; } /* For woman I, FIANCEE(I) is the man she has agree to marry, or 0 if she has not agreed to any man yet. */ for ( i = 0; i < n; i++ ) { fiancee[i] = -1; } /* Start with an unengaged man, and end with an engaged woman. */ for ( i = 1; i <= n; i++ ) { m = i; for ( ; ; ) { next[m-1] = next[m-1] + 1; w = prefer[m-1+(next[m-1]-1)*n]; if ( fiancee[w-1] == -1 ) { fiancee[w-1] = m; break; } if ( rank[w-1+(m-1)*n] < rank[w-1+(fiancee[w-1]-1)*n] ) { temp = fiancee[w-1]; fiancee[w-1] = m; m = temp; } } } return; } /******************************************************************************/ int mountain ( int n, int x, int y ) /******************************************************************************/ /* Purpose: MOUNTAIN enumerates the mountains. Licensing: This code is distributed under the MIT license. Modified: 24 July 2011 Author: John Burkardt Reference: Donald Kreher, Douglas Simpson, Combinatorial Algorithms, CRC Press, 1998, ISBN: 0-8493-3988-X, LC: QA164.K73. Parameters: Input, int N, ... N must be positive. Input, int X, Y, ... 0 <= X <= 2 * N, Output, int MOUNTAIN, the value of the "mountain function" M ( N, X, Y ), which is the number of all mountain ranges from (X,Y) to (2*N,0) which do not drop below sea level. */ { int a; int b; int c; int value; /* Check. */ if ( n <= 0 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "MOUNTAIN(): Fatal error!\n" ); fprintf ( stderr, " N <= 0.\n" ); fprintf ( stderr, " N = %d\n", n ); exit ( 1 ); } else if ( x < 0 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "MOUNTAIN(): Fatal error!\n" ); fprintf ( stderr, " X < 0.\n" ); fprintf ( stderr, " X = %d\n", x ); exit ( 1 ); } else if ( 2 * n < x ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "MOUNTAIN(): Fatal error!\n" ); fprintf ( stderr, " 2 * N < X.\n" ); fprintf ( stderr, " X = %d\n", x ); fprintf ( stderr, " N = %d\n", n ); exit ( 1 ); } /* Special cases. */ if ( y < 0 ) { value = 0; } else if ( 2 * n < x + y ) { value = 0; } else if ( ( ( x + y ) % 2 ) == 1 ) { value = 0; } else { a = 2 * n - x; b = n - ( x + y ) / 2; c = n - 1 - ( x + y ) / 2; value = i4_choose ( a, b ) - i4_choose ( a, c ); } return value; } /******************************************************************************/ int npart_enum ( int n, int npart ) /******************************************************************************/ /* Purpose: NPART_ENUM enumerates the number of partitions of N with NPART parts. Licensing: This code is distributed under the MIT license. Modified: 24 July 2011 Author: John Burkardt Reference: Donald Kreher, Douglas Simpson, Combinatorial Algorithms, CRC Press, 1998, ISBN: 0-8493-3988-X, LC: QA164.K73. Parameters: Input, int N, the integer to be partitioned. Normally N must be positive, but for this routine any N is allowed. Input, int NPART, the number of parts of the partition. Normally, 1 <= NPART <= N is required, but for this routine any value of NPART is allowed. Output, int NPART_ENUM is the number of partitions of N with NPART parts. */ { int *p; int value; if ( n <= 0 ) { value = 0; } else if ( npart <= 0 || n < npart ) { value = 0; } else { p = npart_table ( n, npart ); value = p[n+npart*(n+1)]; free ( p ); } return value; } /******************************************************************************/ int *npart_rsf_lex_random ( int n, int npart, int *seed ) /******************************************************************************/ /* Purpose: npart_rsf_lex_random() returns a random RSF NPART partition. Licensing: This code is distributed under the MIT license. Modified: 25 July 2011 Author: John Burkardt Parameters: Input, int N, the integer to be partitioned. N must be positive. Input, int NPART, the number of parts of the partition. 1 <= NPART <= N. Input/output, int *SEED, a seed for the random number generator. Output, int NPART_RSF_LEX_RANDOM[NPART], contains the partition. A(1) through A(NPART) contain the nonzero integers which sum to N. */ { int *a; int npartitions; int rank; npartitions = npart_enum ( n, npart ); rank = i4_uniform_ab ( 1, npartitions, seed ); a = npart_rsf_lex_unrank ( rank, n, npart ); return a; } /******************************************************************************/ int npart_rsf_lex_rank ( int n, int npart, int a[] ) /******************************************************************************/ /* Purpose: NPART_RSF_LEX_RANK computes the lex rank of an RSF NPART partition. Licensing: This code is distributed under the MIT license. Modified: 25 July 2011 Author: John Burkardt Reference: Donald Kreher, Douglas Simpson, Combinatorial Algorithms, CRC Press, 1998, ISBN: 0-8493-3988-X, LC: QA164.K73. Parameters: Input, int N, the integer to be partitioned. N must be positive. Input, int NPART, the number of parts of the partition. 1 <= NPART <= N. Input, int A[NPART], contains the partition. A(1) through A(NPART) contain the nonzero integers which sum to N. Output, int NPART_RSF_LEX_RANK, the rank of the partition. */ { int *b; int check; int i; int ncopy; int npartcopy; int *p; int rank; /* Check. */ check = part_rsf_check ( n, npart, a ); if ( ! check ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "NPART_RSF_LEX_RANK(): Fatal error!\n" ); fprintf ( stderr, " The input array is illegal.\n" ); exit ( 1 ); } /* Get the table of partitions of N with NPART parts. */ p = npart_table ( n, npart ); /* Copy the partition "backwards". */ b = ( int * ) malloc ( npart * sizeof ( int ) ); for ( i = 1; i <= npart; i++ ) { b[i-1] = a[npart-i]; } rank = 0; ncopy = n; npartcopy = npart; while ( 0 < ncopy && 0 < npartcopy ) { if ( b[npartcopy-1] == 1 ) { ncopy = ncopy - 1; npartcopy = npartcopy - 1; } else { for ( i = 0; i < npartcopy; i++ ) { b[i] = b[i] - 1; } rank = rank + p[ncopy-1+(npartcopy-1)*(n+1)]; ncopy = ncopy - npartcopy; } } free ( b ); free ( p ); return rank; } /******************************************************************************/ void npart_rsf_lex_successor ( int n, int npart, int a[], int *rank ) /******************************************************************************/ /* Purpose: NPART_RSF_LEX_SUCCESSOR computes the RSF lex successor for NPART partitions. Licensing: This code is distributed under the MIT license. Modified: 26 July 2011 Author: John Burkardt Reference: Donald Kreher, Douglas Simpson, Combinatorial Algorithms, CRC Press, 1998, ISBN: 0-8493-3988-X, LC: QA164.K73. Parameters: Input, int N, the integer to be partitioned. N must be at least 1. Input, int NPART, the number of parts of the partition. 1 <= NPART <= N. Input/output, int A[NPART], contains the partition. A(1) through A(NPART) contain the nonzero integers which sum to N. Input/output, int *RANK, the rank. If RANK = -1 on input, then the routine understands that this is the first call, and that the user wishes the routine to supply the first element in the ordering, which has RANK = 0. In general, the input value of RANK is increased by 1 for output, unless the very last element of the ordering was input, in which case the output value of RANK is 0. */ { int check; int d; int i; int j; /* Return the first element. */ if ( *rank == -1 ) { if ( npart < 1 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "NPART_RSF_LEX_SUCCESSOR(): Fatal error!\n" ); fprintf ( stderr, " NPART < 1.\n" ); exit ( 1 ); } for ( i = 0; i < npart - 1; i++ ) { a[i] = 1; } a[npart-1] = n - ( npart - 1 ); *rank = 0; return; } /* Check. */ check = part_rsf_check ( n, npart, a ); if ( ! check ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "NPART_RSF_LEX_SUCCESSOR(): Fatal error!\n" ); fprintf ( stderr, " The input array is illegal.\n" ); exit ( 1 ); } /* Find the first index I for which A(NPART+1-I) + 1 < A(NPART). */ i = 2; for ( ; ; ) { if ( npart < i ) { break; } if ( a[npart-i] + 1 < a[npart-1] ) { break; } i = i + 1; } /* If no such index, we've reached the end of the line. */ if ( i == npart + 1 ) { for ( i = 0; i < npart - 1; i++ ) { a[i] = 1; } a[npart-1] = n - ( npart - 1 ); *rank = 0; return; } /* Otherwise, increment A(NPART+1-I), and adjust other entries. */ else { a[npart-i] = a[npart-i] + 1; d = - 1; for ( j = i - 1; 2 <= j; j-- ) { d = d + a[npart-j] - a[npart-i]; a[npart-j] = a[npart-i]; } a[npart-1] = a[npart-1] + d; } *rank = *rank + 1; return; } /******************************************************************************/ int *npart_rsf_lex_unrank ( int rank, int n, int npart ) /******************************************************************************/ /* Purpose: NPART_RSF_LEX_UNRANK unranks an RSF NPART partition in the lex ordering. Licensing: This code is distributed under the MIT license. Modified: 24 July 2011 Author: John Burkardt Reference: Donald Kreher, Douglas Simpson, Combinatorial Algorithms, CRC Press, 1998, ISBN: 0-8493-3988-X, LC: QA164.K73. Parameters: Input, int RANK, the rank of the partition. Input, int N, the integer to be partitioned. N must be positive. Input, int NPART, the number of parts of the partition. 1 <= NPART <= N. Output, int NPART_RSF_LEX_UNRANK[NPART], contains the partition. A(1) through A(NPART) contain the nonzero integers which sum to N. */ { int *a; int i; int ncopy; int npartcopy; int npartitions; int *p; int rank_copy; /* Check. */ if ( n <= 0 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "NPART_RSF_LEX_UNRANK(): Fatal error!\n" ); fprintf ( stderr, " The input N is illegal.\n" ); exit ( 1 ); } if ( npart < 1 || n < npart ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "NPART_RSF_LEX_UNRANK(): Fatal error!\n" ); fprintf ( stderr, " The input NPART is illegal.\n" ); exit ( 1 ); } npartitions = npart_enum ( n, npart ); if ( rank < 0 || npartitions < rank ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "NPART_RSF_LEX_UNRANK(): Fatal error!\n" ); fprintf ( stderr, " The input rank is illegal.\n" ); exit ( 1 ); } /* Get the table of partitions of N with NPART parts. */ p = npart_table ( n, npart ); a = ( int * ) malloc ( npart * sizeof ( int ) ); for ( i = 0; i < npart; i++ ) { a[i] = 0; } rank_copy = rank; ncopy = n; npartcopy = npart; while ( 0 < ncopy ) { if ( rank_copy < p[ncopy-1+(npartcopy-1)*(n+1)] ) { a[npart-npartcopy] = a[npart-npartcopy] + 1; ncopy = ncopy - 1; npartcopy = npartcopy - 1; } else { for ( i = 1; i <= npartcopy; i++ ) { a[npart-i] = a[npart-i] + 1; } rank_copy = rank_copy - p[ncopy-1+(npartcopy-1)*(n+1)]; ncopy = ncopy - npartcopy; } } return a; } /******************************************************************************/ void npart_sf_lex_successor ( int n, int npart, int a[], int *rank ) /******************************************************************************/ /* Purpose: NPART_SF_LEX_SUCCESSOR computes SF NPART partition. Licensing: This code is distributed under the MIT license. Modified: 26 July 2011 Author: John Burkardt Reference: Donald Kreher, Douglas Simpson, Combinatorial Algorithms, CRC Press, 1998, ISBN: 0-8493-3988-X, LC: QA164.K73. Parameters: Input, int N, the integer to be partitioned. N must be positive. Input, int NPART, the number of parts of the partition. 1 <= NPART <= N. Input/output, int A[NPART], contains the partition. A(1) through A(NPART) contain the nonzero integers which sum to N. The values in A must be in DESCENDING order. Input/output, int *RANK, the rank. If RANK = -1 on input, then the routine understands that this is the first call, and that the user wishes the routine to supply the first element in the ordering, which has RANK = 0. In general, the input value of RANK is increased by 1 for output, unless the very last element of the ordering was input, in which case the output value of RANK is 0. */ { int check; int i; int indx; int temp; /* Return the first element. */ if ( *rank == -1 ) { i4vec_part2 ( n, npart, a ); *rank = 0; return; } /* Check. */ check = part_sf_check ( n, npart, a ); if ( ! check ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "NPART_SF_LEX_SUCCESSOR(): Fatal error!\n" ); fprintf ( stderr, " The partition is illegal.\n" ); exit ( 1 ); } /* Find the last entry that is 2 or more. */ for ( i = npart; 1 <= i; i-- ) { if ( 1 < a[i-1] ) { indx = i; break; } } /* As long as the last nonunit occurs after the first position, have it donate 1 to the left. */ if ( 1 < indx ) { a[indx-1] = a[indx-1] - 1; a[indx-2] = a[indx-2] + 1; indx = indx - 1; for ( ; ; ) { if ( indx <= 1 ) { break; } if ( a[indx-1] <= a[indx-2] ) { break; } temp = a[indx-1]; a[indx-1] = a[indx-2]; a[indx-2] = temp; indx = indx - 1; } /* Sum the tail. */ temp = 0; for ( i = indx; i < npart; i++ ) { temp = temp + a[i]; } /* Partition the tail sum equally over the tail. */ i4vec_part2 ( temp, npart - indx, a+indx ); *rank = *rank + 1; } /* If A(2) through A(NPART) are 1, then this is the last element. Return the first one. */ else { i4vec_part2 ( n, npart, a ); *rank = 0; } return; } /******************************************************************************/ int *npart_table ( int n, int npart ) /******************************************************************************/ /* Purpose: NPART_TABLE tabulates the number of partitions of N having NPART parts. Licensing: This code is distributed under the MIT license. Modified: 24 July 2011 Author: John Burkardt Reference: Donald Kreher, Douglas Simpson, Combinatorial Algorithms, CRC Press, 1998, ISBN: 0-8493-3988-X, LC: QA164.K73. Parameters: Input, int N, the integer to be partitioned. N must be positive. Input, int NPART, the number of parts of the partition. 1 <= NPART <= N. Output, int NPART_TABLE[(N+1)*(NPART+1)], P(I,J) is the number of partitions of I having J parts. */ { int i; int j; int *p; p = ( int * ) malloc ( ( n + 1 ) * ( npart + 1 ) * sizeof ( int ) ); p[0+0*(n+1)] = 1; for ( i = 1; i <= n; i++ ) { p[i+0*(n+1)] = 0; } for ( i = 1; i <= n; i++ ) { for ( j = 1; j <= npart; j++ ) { if ( i < j ) { p[i+j*(n+1)] = 0; } else if ( i < 2 * j ) { p[i+j*(n+1)] = p[i-1+(j-1)*(n+1)]; } else { p[i+j*(n+1)] = p[i-1+(j-1)*(n+1)] + p[i-j+j*(n+1)]; } } } return p; } /******************************************************************************/ int part_enum ( int n ) /******************************************************************************/ /* Purpose: PART_ENUM enumerates the number of partitions of N. Licensing: This code is distributed under the MIT license. Modified: 24 July 2011 Author: John Burkardt Reference: Donald Kreher, Douglas Simpson, Combinatorial Algorithms, CRC Press, 1998, ISBN: 0-8493-3988-X, LC: QA164.K73. Parameters: Input, int N, the integer to be partitioned. Normally N must be positive, but for this routine any N is allowed. Output, int PART_ENUM is the number of partitions of N. */ { int *p; int value; if ( n < 0 ) { value = 0; } else { p = part_table ( n ); value = p[n]; free ( p ); } return value; } /******************************************************************************/ int part_rsf_check ( int n, int npart, int a[] ) /******************************************************************************/ /* Purpose: PART_RSF_CHECK checks a reverse standard form partition of an integer. Licensing: This code is distributed under the MIT license. Modified: 25 July 2011 Author: John Burkardt Reference: Donald Kreher, Douglas Simpson, Combinatorial Algorithms, CRC Press, 1998, ISBN: 0-8493-3988-X, LC: QA164.K73. Parameters: Input, int N, the integer to be partitioned. N must be positive. Input, int NPART, the number of parts of the partition. 1 <= NPART <= N. Input, int A[NPART], contains the partition. A(1) through A(NPART) contain the nonzero integers which sum to N. The entries must be in ASCENDING order. Output, int PART_RSF_CHECK. 1, the data is legal. 0, the data is not legal. */ { int check; int i; check = 1; if ( n < 1 ) { check = 0; return check; } if ( npart < 1 || n < npart ) { check = 0; return check; } /* Every entry must lie between 1 and N. */ for ( i = 0; i < npart; i++ ) { if ( a[i] < 1 || n < a[i] ) { check = 0; return check; } } /* The entries must be in ascending order. */ for ( i = 1; i < npart; i++ ) { if ( a[i] < a[i-1] ) { check = 0; return check; } } /* The entries must add up to N. */ if ( i4vec_sum ( npart, a ) != n ) { check = 0; return check; } return check; } /******************************************************************************/ int part_sf_check ( int n, int npart, int a[] ) /******************************************************************************/ /* Purpose: PART_SF_CHECK checks a standard form partition of an integer. Licensing: This code is distributed under the MIT license. Modified: 27 July 2011 Author: John Burkardt Reference: Donald Kreher, Douglas Simpson, Combinatorial Algorithms, CRC Press, 1998, ISBN: 0-8493-3988-X, LC: QA164.K73. Parameters: Input, int N, the integer to be partitioned. N must be positive. Input, int NPART, the number of parts of the partition. 1 <= NPART <= N. Input, int A[NPART], contains the partition. A(1) through A(NPART) contain the nonzero integers which sum to N. The entries must be in DESCENDING order. Output, int PART_SF_CHECK. 1, the data is legal. 0, the data is not legal. */ { int check; int i; check = 1; if ( n < 1 ) { check = 0; return check; } if ( npart < 1 || n < npart ) { check = 0; return check; } /* Every entry must lie between 1 and N. */ for ( i = 0; i < npart; i++ ) { if ( a[i] < 1 || n < a[i] ) { check = 0; return check; } } /* The entries must be in descending order. */ for ( i = 1; i < npart; i++ ) { if ( a[i-1] < a[i] ) { check = 0; return check; } } /* The entries must add up to N. */ if ( i4vec_sum ( npart, a ) != n ) { check = 0; return check; } return check; } /******************************************************************************/ int *part_sf_conjugate ( int n, int npart, int a[], int *npart2 ) /******************************************************************************/ /* Purpose: PART_SF_CONJUGATE computes the conjugate of a partition. Licensing: This code is distributed under the MIT license. Modified: 26 July 2011 Author: John Burkardt Reference: Donald Kreher, Douglas Simpson, Combinatorial Algorithms, CRC Press, 1998, ISBN: 0-8493-3988-X, LC: QA164.K73. Parameters: Input, int N, the integer to be partitioned. N must be positive. Input, int NPART, the number of parts of the partition. 1 <= NPART <= N. Input, int A[N], contains the partition. A(1) through A(NPART) contain the nonzero integers which sum to N. Output, int *NPART2, the number of parts of the conjugate partition. Output, int PART_SF_CONJUGATE[N], contains the conjugate partition. */ { int *b; int check; int i; int j; /* Check. */ check = part_sf_check ( n, npart, a ); if ( ! check ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "PART_SF_CONJUGATE(): Fatal error!\n" ); fprintf ( stderr, " The partition is illegal.\n" ); exit ( 1 ); } *npart2 = a[0]; b = ( int * ) malloc ( n * sizeof ( int ) ); for ( i = 0; i < *npart2; i++ ) { b[i] = 0; } for ( i = 0; i < npart; i++ ) { for ( j = 0; j < a[i]; j++ ) { b[j] = b[j] + 1; } } return b; } /******************************************************************************/ int part_sf_majorize ( int n, int nparta, int a[], int npartb, int b[] ) /******************************************************************************/ /* Purpose: PART_SF_MAJORIZE determines if partition A majorizes partition B. Discussion: The partitions must be in standard form. If A, with NPARTA parts, and B, with NPARTB parts, are both partitions of the same positive integer N, then we say that A majorizes B if, for every index K from 1 to N, it is true that sum ( 1 <= I <= K ) B(I) <= sum ( 1 <= I <= K ) A(I) where entries of A beyond index NPARTA, and of B beyond BPARTB are assumed to be 0. We say that A strictly majorizes B if A majorizes B, and for at least one index K the inequality is strict. For any two partitions of N, it is possible that A majorizes B, B majorizes A, both partitions majorize each other (in which case they are equal), or that neither majorizes the other. Licensing: This code is distributed under the MIT license. Modified: 26 July 2011 Author: John Burkardt Reference: Jack vanLint, Richard Wilson, A Course in Combinatorics, Cambridge, 1992, ISBN: 0-521-42260-4, LC: QA164.L56. Parameters: Input, int N, the integer to be partitioned. N must be positive. Input, int NPARTA, the number of parts in partition A. 1 <= NPARTA <= N. Input, int A[NPARTA], contains partition A in standard form. A(1) through A(NPARTA) contain nonzero integers which sum to N. Input, int NPARTB, the number of parts in partition B. 1 <= NPARTB <= N. Input, int B[NPARTB], contains partition B in standard form. B(1) through B(NPARTB) contain nonzero integers which sum to N. Output, int PART_SF_MAJORIZE, the result of the comparison. -2, A and B are incomparable, but would have been -1. -1, A < B, (A is strictly majorized by B), 0, A = B, (A and B are identical), +1, A > B, (A strictly majorizes B), +2, A and B are incomparable, but would have been +1. */ { int check; int i; int result; int suma; int sumb; /* Check. */ check = part_sf_check ( n, nparta, a ); if ( ! check ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "PART_SF_MAJORIZE(): Fatal error!\n" ); fprintf ( stderr, " The partition is illegal.\n" ); exit ( 1 ); } part_sf_check ( n, npartb, b ); if ( ! check ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "PART_SF_MAJORIZE(): Fatal error!\n" ); fprintf ( stderr, " The partition is illegal.\n" ); exit ( 1 ); } result = 0; suma = 0; sumb = 0; for ( i = 0; i < i4_min ( nparta, npartb ); i++ ) { if ( i < nparta ) { suma = suma + a[i]; } if ( i < npartb ) { sumb = sumb + b[i]; } if ( result == -1 ) { if ( sumb < suma ) { result = -2; return result; } } else if ( result == 0 ) { if ( suma < sumb ) { result = -1; } else if ( sumb < suma ) { result = +1; } } else if ( result == + 1 ) { if ( suma < sumb ) { result = +2; return result; } } } return result; } /******************************************************************************/ void part_successor ( int n, int *npart, int a[], int *rank ) /******************************************************************************/ /* Purpose: PART_SUCCESSOR computes the lexicographic partition successor. Discussion: part_successor() is "inspired by" the GenPartitions algorithm, but instead of relying on recursion, generates the partitions one at a time. Licensing: This code is distributed under the MIT license. Modified: 26 July 2011 Author: John Burkardt Reference: Donald Kreher, Douglas Simpson, Combinatorial Algorithms, CRC Press, 1998, ISBN: 0-8493-3988-X, LC: QA164.K73. Parameters: Input, int N, the integer to be partitioned. N must be positive. Input/output, int *NPART, the number of parts of the partition. 1 <= NPART <= N. Input/output, int A[N], contains the partition. A(1) through A(NPART) contain the nonzero integers which sum to N. Input/output, int *RANK, the rank. If RANK = -1 on input, then the routine understands that this is the first call, and that the user wishes the routine to supply the first element in the ordering, which has RANK = 0. In general, the input value of RANK is increased by 1 for output, unless the very last element of the ordering was input, in which case the output value of RANK is 0. */ { int asum; int check; int i; int ihi; int j; /* Return the first element. */ if ( *rank == -1 ) { for ( i = 0; i < n; i++ ) { a[i] = 1; } *npart = n; *rank = 0; return; } /* Check. */ check = part_sf_check ( n, *npart, a ); if ( ! check ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "PART_SUCCESSOR(): Fatal error!\n" ); fprintf ( stderr, " The partition is illegal.\n" ); exit ( 1 ); } /* If possible, increment the first intermediate position that is less than its left hand neighbor, and has at least one right hand neighbor. */ ihi = *npart - 1; for ( i = ihi; 2 <= i; i-- ) { if ( a[i-1] < a[i-2] ) { asum = - 1; for ( j = i + 1; j <= *npart; j++ ) { asum = asum + a[j-1]; } a[i-1] = a[i-1] + 1; for ( j = i + 1; j <= *npart; j++ ) { a[j-1] = 0; } *npart = i + asum; for ( j = i + 1; j <= *npart; j++ ) { a[j-1] = 1; } *rank = *rank + 1; return; } } /* A) there are two or more parts Increment the first, replace the rest by 1''s. */ if ( 2 <= *npart ) { a[0] = a[0] + 1; for ( j = 2; j <= *npart; j++ ) { a[j-1] = 0; } *npart = n - a[0] + 1; for ( j = 2; j <= *npart; j++ ) { a[j-1] = 1; } *rank = *rank + 1; } /* B) there is only one part. We have reached the last item. Return the first one. */ else if ( *npart == 1 ) { for ( i = 0; i < n; i++ ) { a[i] = 1; } *npart = n; *rank = 0; } return; } /******************************************************************************/ int *part_table ( int n ) /******************************************************************************/ /* Purpose: PART_TABLE tabulates the number of partitions of N. Licensing: This code is distributed under the MIT license. Modified: 26 July 2011 Author: John Burkardt Reference: Donald Kreher, Douglas Simpson, Combinatorial Algorithms, CRC Press, 1998, ISBN: 0-8493-3988-X, LC: QA164.K73. Parameters: Input, int N, the integer to be partitioned. N must be positive. Output, int P[N+1], P(I) is the number of partitions of I. */ { int i; int j; int *p; int psum; int sign; int w; int wprime; p = ( int * ) malloc ( ( n + 1 ) * sizeof ( int ) ); p[0] = 1; if ( n <= 0 ) { return p; } p[1] = 1; for ( i = 2; i <= n; i++ ) { sign = 1; psum = 0; w = 1; j = 1; wprime = w + j; while ( w < n ) { if ( 0 <= i - w ) { if ( sign == 1 ) { psum = psum + p[i-w]; } else { psum = psum - p[i-w]; } } if ( wprime <= i ) { if ( sign == 1 ) { psum = psum + p[i-wprime]; } else { psum = psum - p[i-wprime]; } } w = w + 3 * j + 1; j = j + 1; wprime = w + j; sign = - sign; } p[i] = psum; } return p; } /******************************************************************************/ int *partition_greedy ( int n, int a[] ) /******************************************************************************/ /* Purpose: PARTITION_GREEDY attacks the partition problem with a greedy algorithm. Discussion: Given a collection of N not necessarily distinct positive integers A(I), it is desired to partition the values into two groups, whose sums are as close as possible. Algorithm: Begin with sets 1 and 2 empty. Process the data in descending order of magnitude. The next item A(I) is added to set 1 or set 2, whichever has the smallest current sum. Stop as soon as all items have been allocated. Licensing: This code is distributed under the MIT license. Modified: 24 July 2011 Author: John Burkardt Reference: Brian Hayes, The Easiest Hard Problem, American Scientist, Volume 90, Number 2, March-April 2002, pages 113-117. Parameters: Input, int N, the number of values. N must be positive. Input/output, int A[N], a collection of positive values. On output, A has been sorted into descending order. Output, int PARTITION_GREEDY[N]; an entry is 0 if A[I] is part of set 0, and 1 if it is assigned to set 1. */ { int i; int *indx; int j; int sums[2]; sums[0] = 0; sums[1] = 0; i4vec_sort_insert_d ( n, a ); indx = ( int * ) malloc ( n * sizeof ( int ) ); for ( i = 0; i < n; i++ ) { if ( sums[0] < sums[1] ) { j = 0; } else { j = 1; } indx[i] = j; sums[j] = sums[j] + a[i]; } return indx; } /******************************************************************************/ int partn_enum ( int n, int nmax ) /******************************************************************************/ /* Purpose: PARTN_ENUM enumerates the partitions of N with maximum element NMAX. Licensing: This code is distributed under the MIT license. Modified: 24 July 2011 Author: John Burkardt Reference: Donald Kreher, Douglas Simpson, Combinatorial Algorithms, CRC Press, 1998, ISBN: 0-8493-3988-X, LC: QA164.K73. Parameters: Input, int N, the integer to be partitioned. Normally N must be positive, but for this routine any N is allowed. Input, int NMAX, the maximum element in the partition. Normally, 1 <= NMAX <= N is required, but for this routine any value of NMAX is allowed. Output, int PARTN_ENUM is the number of partitions of N with maximum element NMAX. */ { int *p; int value; if ( n <= 0 ) { value = 0; } else if ( nmax <= 0 || n < nmax ) { value = 0; } else { p = npart_table ( n, nmax ); value = p[n+nmax*(n+1)]; free ( p ); } return value; } /******************************************************************************/ int partn_sf_check ( int n, int nmax, int npart, int a[] ) /******************************************************************************/ /* Purpose: PARTN_SF_CHECK checks an SF partition of an integer with largest entry NMAX. Licensing: This code is distributed under the MIT license. Modified: 25 July 2011 Author: John Burkardt Reference: Donald Kreher, Douglas Simpson, Combinatorial Algorithms, CRC Press, 1998, ISBN: 0-8493-3988-X, LC: QA164.K73. Parameters: Input, int N, the integer to be partitioned. N must be positive. Input, int NMAX, the value of the largest entry. 1 <= NMAX <= N. Input, int NPART, the number of parts of the partition. 1 <= NPART <= N. Input, int A[NPART], contains the partition. A(1) through A(NPART) contain the nonzero integers which sum to N. The entries must be in DESCENDING order. Output, int PARTN_SF_CHECK. 1, the data is legal. 0, the data is not legal. */ { int asum; int check; int i; check = 1; if ( n < 1 ) { check = 0; return check; } if ( nmax < 1 || n < nmax ) { check = 0; return check; } if ( npart < 1 || n < npart ) { check = 0; return check; } /* Entry 1 must be NMAX. */ if ( a[0] != nmax ) { check = 0; return check; } /* Every entry must lie between 1 and N. */ for ( i = 0; i < npart; i++ ) { if ( a[i] < 1 || n < a[i] ) { check = 0; return check; } } /* The entries must be in descending order. */ for ( i = 1; i < npart; i++ ) { if ( a[i-1] < a[i] ) { check = 0; return check; } } /* The entries must add up to N. */ asum = i4vec_sum ( npart, a ); if ( asum != n ) { check = 0; return check; } return check; } /******************************************************************************/ void partn_successor ( int n, int nmax, int *npart, int a[], int *rank ) /******************************************************************************/ /* Purpose: PARTN_SUCCESSOR computes partitions whose largest part is NMAX. Licensing: This code is distributed under the MIT license. Modified: 28 July 2011 Author: John Burkardt Reference: Donald Kreher, Douglas Simpson, Combinatorial Algorithms, CRC Press, 1998, ISBN: 0-8493-3988-X, LC: QA164.K73. Parameters: Input, int N, the integer to be partitioned. N must be positive. Input, int NMAX, the maximum size of any part of the partition. 1 <= NMAX <= N. Input/output, int *NPART, the number of parts of the partition. 1 <= NPART <= N. Input/output, int A[N], contains the partition. A(1) through A(NPART) contain the nonzero integers which sum to N. Input/output, int *RANK, the rank. If RANK = -1 on input, then the routine understands that this is the first call, and that the user wishes the routine to supply the first element in the ordering, which has RANK = 0. In general, the input value of RANK is increased by 1 for output, unless the very last element of the ordering was input, in which case the output value of RANK is 0. */ { int check; int i; int index; int temp; /* Return the first element. */ if ( *rank == -1 ) { *npart = n + 1 - nmax; a[0] = nmax; for ( i = 1; i < *npart; i++ ) { a[i] = 1; } *rank = 0; return; } /* Check. */ check = partn_sf_check ( n, nmax, *npart, a ); if ( ! check ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "PARTN_SUCCESSOR(): Fatal error!\n" ); fprintf ( stderr, " The input array is illegal.\n" ); exit ( 1 ); } /* If there are at least two parts, and the next to last is not NMAX, then rob the last part and pay the next to the last part. Then, if the next to last part is too big, swap it leftwards. */ if ( 1 < *npart ) { if ( a[*npart-2] < nmax ) { a[*npart-1] = a[*npart-1] - 1; a[*npart-2] = a[*npart-2] + 1; index = *npart - 1; for ( ; ; ) { if ( index <= 1 ) { break; } if ( a[index-1] <= a[index-2] ) { break; } temp = a[index-2]; a[index-2] = a[index-1]; a[index-1] = temp; index = index - 1; } /* Sum the tail. */ temp = 0; for ( i = index; i < *npart; i++ ) { temp = temp + a[i]; } /* Spread the sum as 1''s. */ *npart = index + temp; for ( i = index; i < *npart; i++ ) { a[i] = 1; } *rank = *rank + 1; return; } } /* Otherwise, we have reached the last item. Return the first one. */ else { *npart = n + 1 - nmax; a[0] = nmax; for ( i = 1; i < *npart; i++ ) { a[i] = 1; } *rank = 0; return; } return; } /******************************************************************************/ int perm_check ( int n, int p[] ) /******************************************************************************/ /* Purpose: PERM_CHECK checks a representation of a permutation. Discussion: The routine is given N and P, a vector of length N. P is a legal represention of a permutation of the integers from 1 to N if and only if every integer from 1 to N occurs as a value of P(I) for some I between 1 and N. Licensing: This code is distributed under the MIT license. Modified: 25 July 2011 Author: John Burkardt Parameters: Input, int N, the number of values being permuted. N must be positive. Input, int P[N], the array to check. Output, int PERM_CHECK. 1, the data is legal. 0, the data is not legal. */ { int check; int i; int ifind; int iseek; check = 1; if ( n < 1 ) { check = 0; return check; } for ( i = 0; i < n; i++ ) { if ( p[i] < 1 || n < p[i] ) { check = 0; return check; } } for ( iseek = 1; iseek <= n; iseek++ ) { ifind = -1; for ( i = 0; i < n; i++ ) { if ( p[i] == iseek ) { ifind = i; break; } } if ( ifind == -1 ) { check = 0; return check; } } return check; } /******************************************************************************/ int perm_enum ( int n ) /******************************************************************************/ /* Purpose: PERM_ENUM enumerates the permutations on N digits. Licensing: This code is distributed under the MIT license. Modified: 24 July 2011 Author: John Burkardt Parameters: Input, int N, the number of values being permuted. N must be nonnegative. Output, int PERM_ENUM, the number of distinct elements. */ { int value; value = i4_factorial ( n ); return value; } /******************************************************************************/ int *perm_inv ( int n, int p[] ) /******************************************************************************/ /* Purpose: PERM_INV computes the inverse of a permutation. Licensing: This code is distributed under the MIT license. Modified: 25 July 2011 Author: John Burkardt Reference: Donald Kreher, Douglas Simpson, Combinatorial Algorithms, CRC Press, 1998, ISBN: 0-8493-3988-X, LC: QA164.K73. Parameters: Input, int N, the number of values being permuted. N must be positive. Input, int P[N], describes the permutation. P(I) is the item which is permuted into the I-th place by the permutation. Output, int PERM_INV[N], the inverse permutation. */ { int check; int i; int *pinv; /* Check. */ check = perm_check ( n, p ); if ( ! check ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "PERM_INV(): Fatal error!\n" ); fprintf ( stderr, " The permutation is not legal.\n" ); exit ( 1 ); } pinv = ( int * ) malloc ( n * sizeof ( int ) ); for ( i = 0; i < n; i++ ) { pinv[p[i]-1] = i + 1; } return pinv; } /******************************************************************************/ int perm_lex_rank ( int n, int p[] ) /******************************************************************************/ /* Purpose: PERM_LEX_RANK computes the lexicographic rank of a permutation. Discussion: The original code altered the input permutation. Licensing: This code is distributed under the MIT license. Modified: 25 July 2011 Author: John Burkardt Reference: Donald Kreher, Douglas Simpson, Combinatorial Algorithms, CRC Press, 1998, ISBN: 0-8493-3988-X, LC: QA164.K73. Parameters: Input, int N, the number of values being permuted. N must be positive. Input, int P[N], describes the permutation. P[I] is the item which is permuted into the I-th place by the permutation. Output, int PERM_LEX_RANK, the rank of the permutation. */ { int check; int i; int j; int *pcopy; int rank; /* Check. */ check = perm_check ( n, p ); if ( ! check ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "PERM_LEX_RANK(): Fatal error!\n" ); fprintf ( stderr, " The permutation is not legal.\n" ); exit ( 1 ); } rank = 0; pcopy = ( int * ) malloc ( n * sizeof ( int ) ); for ( i = 0; i < n; i++ ) { pcopy[i] = p[i]; } for ( j = 0; j < n; j++ ) { rank = rank + ( pcopy[j] - 1 ) * i4_factorial ( n - 1 - j ); for ( i = j + 1; i < n; i++ ) { if ( pcopy[j] < pcopy[i] ) { pcopy[i] = pcopy[i] - 1; } } } free ( pcopy ); return rank; } /******************************************************************************/ void perm_lex_successor ( int n, int p[], int *rank ) /******************************************************************************/ /* Purpose: PERM_LEX_SUCCESSOR computes the lexicographic permutation successor. Example: RANK Permutation 0 1 2 3 4 1 1 2 4 3 2 1 3 2 4 3 1 3 4 2 4 1 4 2 3 5 1 4 3 2 6 2 1 3 4 ... 23 4 3 2 1 Licensing: This code is distributed under the MIT license. Modified: 26 July 2011 Author: John Burkardt Reference: Donald Kreher, Douglas Simpson, Combinatorial Algorithms, CRC Press, 1998, ISBN: 0-8493-3988-X, LC: QA164.K73. Parameters: Input, int N, the number of values being permuted. N must be positive. Input/output, int P[N], describes the permutation. P(I) is the item which is permuted into the I-th place by the permutation. Input/output, int *RANK, the rank. If RANK = -1 on input, then the routine understands that this is the first call, and that the user wishes the routine to supply the first element in the ordering, which has RANK = 0. In general, the input value of RANK is increased by 1 for output, unless the very last element of the ordering was input, in which case the output value of RANK is 0. */ { int check; int i; int j; int temp; /* Return the first element. */ if ( *rank == -1 ) { i4vec_indicator1 ( n, p ); *rank = 0; return; } /* Check. */ check = perm_check ( n, p ); if ( ! check ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "PERM_LEX_SUCCESSOR(): Fatal error!\n" ); fprintf ( stderr, " The permutation is not legal.\n" ); exit ( 1 ); } /* Seek I, the highest index for which the next element is bigger. */ i = n - 1; for ( ; ; ) { if ( i <= 0 ) { break; } if ( p[i-1] <= p[i] ) { break; } i = i - 1; } /* If no I could be found, then we have reach the final permutation, N, N-1, ..., 2, 1. Time to start over again. */ if ( i == 0 ) { i4vec_indicator1 ( n, p ); *rank = 0; } else { /* Otherwise, look for the the highest index after I whose element is bigger than I''s. We know that I+1 is one such value, so the loop will never fail. */ j = n; while ( p[j-1] < p[i-1] ) { j = j - 1; } /* Interchange elements I and J. */ temp = p[i-1]; p[i-1] = p[j-1]; p[j-1] = temp; /* Reverse the elements from I+1 to N. */ i4vec_reverse ( n - i, p+i ); *rank = *rank + 1; } return; } /******************************************************************************/ int *perm_lex_unrank ( int rank, int n ) /******************************************************************************/ /* Purpose: PERM_LEX_UNRANK computes the permutation of given lexicographic rank. Licensing: This code is distributed under the MIT license. Modified: 26 July 2011 Author: John Burkardt Reference: Donald Kreher, Douglas Simpson, Combinatorial Algorithms, CRC Press, 1998, ISBN: 0-8493-3988-X, LC: QA164.K73. Parameters: Input, int RANK, the rank of the permutation. Input, int N, the number of values being permuted. N must be positive. Output, int PERM_LEX_UNRANK[N], describes the permutation. */ { int d; int i; int j; int nperm; int *p; int rank_copy; /* Check. */ if ( n < 1 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "PERM_LEX_UNRANK(): Fatal error!\n" ); fprintf ( stderr, " Input N is illegal.\n" ); exit ( 1 ); } nperm = perm_enum ( n ); if ( rank < 0 || nperm < rank ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "PERM_LEX_UNRANK(): Fatal error!\n" ); fprintf ( stderr, " The input rank is illegal.\n" ); exit ( 1 ); } rank_copy = rank; p = ( int * ) malloc ( n * sizeof ( int ) ); p[n-1] = 1; for ( j = 1; j <= n - 1; j++ ) { d = ( rank_copy % i4_factorial ( j + 1 ) ) / i4_factorial ( j ); rank_copy = rank_copy - d * i4_factorial ( j ); p[n-j-1] = d + 1; for ( i = n - j + 1; i <= n; i++ ) { if ( d < p[i-1] ) { p[i-1] = p[i-1] + 1; } } } return p; } /******************************************************************************/ int *perm_mul ( int n, int p[], int q[] ) /******************************************************************************/ /* Purpose: PERM_MUL computes the product of two permutations. Licensing: This code is distributed under the MIT license. Modified: 26 July 2011 Author: John Burkardt Reference: Donald Kreher, Douglas Simpson,inson, Combinatorial Algorithms, CRC Press, 1998, ISBN: 0-8493-3988-X, LC: QA164.K73. Parameters: Input, int N, the number of values being permuted. N must be positive. Input, int P[N], Q[N], describes the permutation factors. Output, int PERMN_MUL[N], the product permutation P * Q. R(I) = P(Q(I)). */ { int check; int i; int *r; /* Check. */ check = perm_check ( n, p ); if ( ! check ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "PERM_MUL(): Fatal error!\n" ); fprintf ( stderr, " The permutation is not legal.\n" ); exit ( 1 ); } check = perm_check ( n, q ); if ( ! check ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "PERM_MUL(): Fatal error!\n" ); fprintf ( stderr, " The permutation is not legal.\n" ); exit ( 1 ); } /* Use a temporary vector for the result, to avoid problems if some arguments are actually identified. */ r = ( int * ) malloc ( n * sizeof ( int ) ); for ( i = 0; i < n; i++ ) { r[i] = p[q[i]-1]; } return r; } /******************************************************************************/ int perm_parity ( int n, int p[] ) /******************************************************************************/ /* Purpose: PERM_PARITY computes the parity of a permutation. Discussion: The routine requires the use of a temporary array. A permutation is called "even" or "odd", depending on whether it is equivalent to an even or odd number of pairwise transpositions. This is known as the "parity" of the permutation. The "sign" of a permutation is +1 if it has even parity, and -1 if it has odd parity. Licensing: This code is distributed under the MIT license. Modified: 25 July 2011 Author: John Burkardt Reference: Donald Kreher, Douglas Simpson, Combinatorial Algorithms, CRC Press, 1998, ISBN: 0-8493-3988-X, LC: QA164.K73. Parameters: Input, int N, the number of values being permuted. N must be positive. Input, int P[N], describes the permutation. P(I) is the item which is permuted into the I-th place by the permutation. Output, int PERM_PARITY, the parity of the permutation. 0, the permutation has even parity. 1, the permutation has odd parity. */ { int *a; int c; int check; int i; int j; int parity; /* Check. */ check = perm_check ( n, p ); if ( ! check ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "PERM_PARITY(): Fatal error!\n" ); fprintf ( stderr, " The permutation is not legal.\n" ); exit ( 1 ); } a = ( int * ) malloc ( n * sizeof ( int ) ); for ( i = 0; i < n; i++ ) { a[i] = 0; } c = 0; for ( j = 1; j <= n; j++ ) { if ( a[j-1] == 0 ) { c = c + 1; a[j-1] = 1; i = j; while ( p[i-1] != j ) { i = p[i-1]; a[i-1] = 1; } } } parity = ( n - c ) % 2; free ( a ); return parity; } /******************************************************************************/ void perm_print ( int n, int p[], char *title ) /******************************************************************************/ /* Purpose: PERM_PRINT prints a permutation. Discussion: The permutation is assumed to be zero-based. Example: Input: P = 6 1 2 0 4 2 5 Printed output: "This is the permutation:" 0 1 2 3 4 5 6 6 1 2 0 4 2 5 Licensing: This code is distributed under the MIT license. Modified: 14 May 2011 Author: John Burkardt Parameters: Input, int N, the number of objects permuted. Input, int P[N], the permutation, in standard index form. Input, char *TITLE, a title. If no title is supplied, then only the permutation is printed. */ { int i; int ihi; int ilo; int inc = 20; if ( s_len_trim ( title ) != 0 ) { printf ( "\n" ); printf ( "%s\n", title ); for ( ilo = 0; ilo < n; ilo = ilo + inc ) { ihi = ilo + inc; if ( n < ihi ) { ihi = n; } printf ( "\n" ); printf ( " " ); for ( i = ilo; i < ihi; i++ ) { printf ( "%4d", i ); } printf ( "\n" ); printf ( " " ); for ( i = ilo; i < ihi; i++ ) { printf ( "%4d", p[i] ); } printf ( "\n" ); } } else { for ( ilo = 0; ilo < n; ilo = ilo + inc ) { ihi = ilo + inc; if ( n < ihi ) { ihi = n; } printf ( " " ); for ( i = ilo; i < ihi; i++ ) { printf ( "%4d", p[i] ); } printf ( "\n" ); } } return; } /******************************************************************************/ int *perm_random ( int n, int *seed ) /******************************************************************************/ /* Purpose: perm_random() selects a random permutation of 1, ..., N. Discussion: The algorithm is known as the Fisher-Yates or Knuth shuffle. Licensing: This code is distributed under the MIT license. Modified: 11 January 2016 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 objects to be permuted. Input/output, int *SEED, a seed for the random number generator. Output, int PERM_RANDOM[N], a permutation of ( 1, 2, ..., N ), in standard index form. */ { int i; int j; int *p; int t; p = ( int * ) malloc ( n * sizeof ( int ) ); for ( i = 0; i < n; i++ ) { p[i] = i + 1; } for ( i = 0; i < n - 1; i++ ) { j = i4_uniform_ab ( i, n - 1, seed ); t = p[i]; p[i] = p[j]; p[j] = t; } return p; } /******************************************************************************/ int perm_tj_rank ( int n, int p[] ) /******************************************************************************/ /* Purpose: PERM_TJ_RANK computes the Trotter-Johnson rank of a permutation. Licensing: This code is distributed under the MIT license. Modified: 25 July 2011 Author: John Burkardt Reference: Donald Kreher, Douglas Simpson, Combinatorial Algorithms, CRC Press, 1998, ISBN: 0-8493-3988-X, LC: QA164.K73. Parameters: Input, int N, the number of values being permuted. N must be positive. Input, int P[N], describes the permutation. P(I) is the item which is permuted into the I-th place by the permutation. Output, int PERM_TJ_RANK, the rank of the permutation. */ { int check; int i; int j; int k; int rank; /* Check. */ check = perm_check ( n, p ); if ( ! check ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "PERM_TJ_RANK(): Fatal error!\n" ); fprintf ( stderr, " The permutation is not legal.\n" ); exit ( 1 ); } rank = 0; for ( j = 2; j <= n; j++ ) { k = 1; i = 1; while ( p[i-1] != j ) { if ( p[i-1] < j ) { k = k + 1; } i = i + 1; } if ( ( rank % 2 ) == 0 ) { rank = j * rank + j - k; } else { rank = j * rank + k - 1; } } return rank; } /******************************************************************************/ void perm_tj_successor ( int n, int p[], int *rank ) /******************************************************************************/ /* Purpose: PERM_TJ_SUCCESSOR computes the Trotter-Johnson permutation successor. Licensing: This code is distributed under the MIT license. Modified: 26 July 2011 Author: John Burkardt Reference: Donald Kreher, Douglas Simpson, Combinatorial Algorithms, CRC Press, 1998, ISBN: 0-8493-3988-X, LC: QA164.K73. Parameters: Input, int N, the number of values being permuted. N must be positive. Input/output, int P[N], describes the permutation. P(I) is the item which is permuted into the I-th place by the permutation. Input/output, int *RANK, the rank. If RANK = -1 on input, then the routine understands that this is the first call, and that the user wishes the routine to supply the first element in the ordering, which has RANK = 0. In general, the input value of RANK is increased by 1 for output, unless the very last element of the ordering was input, in which case the output value of RANK is 0. */ { int check; int d; int done; int i; int m; int par; int *q; int st; int temp; /* Return the first element. */ if ( *rank == -1 ) { i4vec_indicator1 ( n, p ); *rank = 0; return; } /* Check. */ check = perm_check ( n, p ); if ( ! check ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "PERM_TJ_SUCCESSOR(): Fatal error!\n" ); fprintf ( stderr, " The permutation is not legal.\n" ); exit ( 1 ); } q = ( int * ) malloc ( n * sizeof ( int ) ); st = 0; for ( i = 0; i < n; i++ ) { q[i] = p[i]; } done = 0; m = n; while ( 1 < m && !done ) { d = 1; while ( q[d-1] != m ) { d = d + 1; } for ( i = d; i < m; i++ ) { q[i-1] = q[i]; } par = perm_parity ( m - 1, q ); if ( par == 1 ) { if ( d == m ) { m = m - 1; } else { temp = p[st+d-1]; p[st+d-1] = p[st+d]; p[st+d] = temp; done = 1; } } else { if ( d == 1 ) { m = m - 1; st = st + 1; } else { temp = p[st+d-1]; p[st+d-1] = p[st+d-2]; p[st+d-2] = temp; done = 1; } } } /* Last element was input. Return first one. */ if ( m == 1 ) { i4vec_indicator1 ( n, p ); *rank = 0; free ( q ); return; } *rank = *rank + 1; free ( q ); return; } /******************************************************************************/ int *perm_tj_unrank ( int rank, int n ) /******************************************************************************/ /* Purpose: PERM_TJ_UNRANK computes the permutation of given Trotter-Johnson rank. Licensing: This code is distributed under the MIT license. Modified: 25 July 2011 Author: John Burkardt Reference: Donald Kreher, Douglas Simpson, Combinatorial Algorithms, CRC Press, 1998, ISBN: 0-8493-3988-X, LC: QA164.K73. Parameters: Input, int RANK, the rank of the permutation. Input, int N, the number of values being permuted. N must be positive. Output, int PERM_TJ_UNRANK[N], describes the permutation. */ { int i; int j; int k; int jhi; int nperm; int *p; int r1; int r2; /* Check. */ if ( n < 1 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "PERM_TJ_UNRANK(): Fatal error!\n" ); fprintf ( stderr, " Input N is illegal.\n" ); exit ( 1 ); } nperm = perm_enum ( n ); if ( rank < 0 || nperm < rank ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "PERM_TJ_UNRANK(): Fatal error!\n" ); fprintf ( stderr, " The input rank is illegal.\n" ); exit ( 1 ); } p = ( int * ) malloc ( n * sizeof ( int ) ); p[0] = 1; r2 = 0; for ( j = 2; j <= n; j++ ) { /* Replace this ratio of factorials! */ r1 = ( rank * i4_factorial ( j ) ) / i4_factorial ( n ); k = r1 - j * r2; if ( ( r2 % 2 ) == 0 ) { jhi = j - k; } else { jhi = k + 1; } for ( i = j - 1; jhi <= i; i-- ) { p[i] = p[i-1]; } p[jhi-1] = j; r2 = r1; } return p; } /******************************************************************************/ void perm_to_cycle ( int n, int p[], int *ncycle, int t[], int index[] ) /******************************************************************************/ /* Purpose: PERM_TO_CYCLE converts a permutation from array to cycle form. Licensing: This code is distributed under the MIT license. Modified: 28 July 2011 Author: John Burkardt Reference: Donald Kreher, Douglas Simpson, Combinatorial Algorithms, CRC Press, 1998, ISBN: 0-8493-3988-X, LC: QA164.K73. Parameters: Input, int N, the number of values being permuted. N must be positive. Input, int P[N], describes the permutation using a single array. For each index I, I -> P(I). Output, int *NCYCLE, the number of cycles. 1 <= NCYCLE <= N. Output, int T[N], INDEX[N], describes the permutation as a collection of NCYCLE cycles. The first cycle is T(1) -> T(2) -> ... -> T(INDEX(1)) -> T(1). */ { int check; int i; int j; int nset; /* Check. */ check = perm_check ( n, p ); if ( ! check ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "PERM_TO_CYCLE(): Fatal error!\n" ); fprintf ( stderr, " The permutation is not legal.\n" ); exit ( 1 ); } /* Initialize. */ *ncycle = 0; for ( i = 0; i < n; i++ ) { index[i] = 0; } for ( i = 0; i < n; i++ ) { t[i] = 0; } nset = 0; /* Find the next unused entry. */ for ( i = 1; i <= n; i++ ) { if ( 0 < p[i-1] ) { *ncycle = *ncycle + 1; index[*ncycle-1] = 1; nset = nset + 1; t[nset-1] = p[i-1]; p[i-1] = - p[i-1]; for ( ; ; ) { j = t[nset-1]; if ( p[j-1] < 0 ) { break; } index[*ncycle-1] = index[*ncycle-1] + 1; nset = nset + 1; t[nset-1] = p[j-1]; p[j-1] = - p[j-1]; } } } /* If no unused entries remain, we are done. Restore the sign of the permutation and return. */ for ( i = 0; i < n; i++ ) { p[i] = - p[i]; } return; } /******************************************************************************/ int pruefer_check ( int n, int p[] ) /******************************************************************************/ /* Purpose: PRUEFER_CHECK checks a Pruefer code. Licensing: This code is distributed under the MIT license. Modified: 26 July 2011 Author: John Burkardt Reference: Donald Kreher, Douglas Simpson, Combinatorial Algorithms, CRC Press, 1998, ISBN: 0-8493-3988-X, LC: QA164.K73. Parameters: Input, int N, the number of nodes in the tree. N must be at least 3. Input, int P[N-2], the Pruefer code for the tree. Output, int PRUEFER_CHECK. 1, the data is legal. 0, the data is not legal. */ { int check; int i; check = 1; if ( n < 3 ) { check = 0; return check; } for ( i = 0; i < n - 2; i++ ) { if ( p[i] < 1 || n < p[i] ) { check = 0; return check; } } return check; } /******************************************************************************/ int pruefer_enum ( int n ) /******************************************************************************/ /* Purpose: PRUEFER_ENUM enumerates the Pruefer codes on N-2 digits. Licensing: This code is distributed under the MIT license. Modified: 29 August 2015 Author: John Burkardt Parameters: Input, int N, the number of digits in the code, plus 2. Output, int NCODE, the number of distinct elements. */ { int value; if ( n < 2 ) { value = 0; } else if ( n == 2 ) { value = 1; } else { value = i4_power ( n, n - 2 ); } return value; } /******************************************************************************/ int *pruefer_random ( int n ) /******************************************************************************/ /* Purpose: pruefer_random() returns a random Pruefer code of N-2 digits. Licensing: This code is distributed under the MIT license. Modified: 12 September 2022 Author: John Burkardt Reference: Donald Kreher, Douglas Simpson, Combinatorial Algorithms, CRC Press, 1998, ISBN: 0-8493-3988-X, LC: QA164.K73. Input: int N, the number of nodes in the tree. N must be at least 3. Output: int P[N-2], the random Pruefer code. */ { int i; int *p; p = ( int * ) malloc ( ( n - 2 ) * sizeof ( int ) ); for ( i = 0; i < n - 2; i++ ) { p[i] = 1 + ( rand ( ) % n ); } return p; } /******************************************************************************/ int pruefer_rank ( int n, int p[] ) /******************************************************************************/ /* Purpose: pruefer_rank() ranks a Pruefer code. Licensing: This code is distributed under the MIT license. Modified: 26 July 2011 Author: John Burkardt Reference: Donald Kreher, Douglas Simpson, Combinatorial Algorithms, CRC Press, 1998, ISBN: 0-8493-3988-X, LC: QA164.K73. Parameters: Input, int N, the number of nodes in the tree. N must be at least 3. Input, int P[N-2], the Pruefer code for the tree. Output, int PRUEFER_RANK, the rank of the Pruefer code. */ { int check; int i; int k; int rank; /* Check. */ check = pruefer_check ( n, p ); if ( ! check ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "PRUEFER_RANK(): Fatal error!\n" ); fprintf ( stderr, " Input array is illegal.\n" ); exit ( 1 ); } rank = 0; k = 1; for ( i = n - 3; 0 <= i; i-- ) { rank = rank + k * ( p[i] - 1 ); k = k * n; } return rank; } /******************************************************************************/ void pruefer_successor ( int n, int p[], int *rank ) /******************************************************************************/ /* Purpose: pruefer_successor() computes the lexical Pruefer sequence successor. Licensing: This code is distributed under the MIT license. Modified: 27 July 2011 Author: John Burkardt Reference: Donald Kreher, Douglas Simpson, Combinatorial Algorithms, CRC Press, 1998, ISBN: 0-8493-3988-X, LC: QA164.K73. Parameters: Input, int N, the number of nodes in the tree. N must be at least 3. Input/output, int P(N-2), on input, the Pruefer code for a tree, and on output, its lexical successor. Input/output, int *RANK, the rank. If RANK = -1 on input, then the routine understands that this is the first call, and that the user wishes the routine to supply the first element in the ordering, which has RANK = 0. In general, the input value of RANK is increased by 1 for output, unless the very last element of the ordering was input, in which case the output value of RANK is 0. */ { int check; int i; int j; /* Return the first element. */ if ( *rank == -1 ) { for ( i = 0; i < n - 2; i++ ) { p[i] = 1; } *rank = 0; return; } /* Check. */ check = pruefer_check ( n, p ); if ( ! check ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "PRUEFER_SUCCESSOR(): Fatal error!\n" ); fprintf ( stderr, " The input array is illegal.\n" ); exit ( 1 ); } j = n - 2; for ( ; ; ) { if ( p[j-1] != n ) { break; } j = j - 1; if ( j <= 0 ) { break; } } if ( j != 0 ) { p[j-1] = p[j-1] + 1; for ( i = j + 1; i <= n - 2; i++ ) { p[i-1] = 1; } *rank = *rank + 1; } else { for ( i = 0; i < n - 2; i++ ) { p[i] = 1; } *rank = 0; } return; } /******************************************************************************/ void pruefer_to_tree ( int n, int p[], int t[] ) /******************************************************************************/ /* Purpose: pruefer_to_tree() converts a Pruefer code to a tree. Discussion: The original code attempts to tack on an extra entry to P. Licensing: This code is distributed under the MIT license. Modified: 24 July 2011 Author: John Burkardt Reference: Donald Kreher, Douglas Simpson, Combinatorial Algorithms, CRC Press, 1998, ISBN: 0-8493-3988-X, LC: QA164.K73. Parameters: Input, int N, the number of nodes in the tree. N must be at least 3. Input, int P[N-2], the Pruefer code for the tree. Output, int T[2*(N-1)], describes the edges of the tree as pairs of nodes. */ { int check; int *d; int i; int j; int x; int y; /* Check. */ check = pruefer_check ( n, p ); if ( ! check ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "PRUEFER_TO_TREE(): Fatal error!\n" ); fprintf ( stderr, " The input array is illegal!\n" ); exit ( 1 ); } /* Initialize the tree to 0. */ for ( j = 0; j < n - 1; j++ ) { for ( i = 0; i < 2; i++ ) { t[i+j*2] = 0; } } d = ( int * ) malloc ( n * sizeof ( int ) ); for ( i = 0; i < n; i++ ) { d[i] = 1; } for ( i = 0; i < n - 2; i++ ) { d[p[i]-1] = d[p[i]-1] + 1; } for ( i = 1; i <= n - 1; i++ ) { x = n; while ( d[x-1] != 1 ) { x = x - 1; } if ( i == n - 1 ) { y = 1; } else { y = p[i-1]; } d[x-1] = d[x-1] - 1; d[y-1] = d[y-1] - 1; t[0+(i-1)*2] = x; t[1+(i-1)*2] = y; } free ( d ); return; } /******************************************************************************/ int *pruefer_to_tree_new ( int n, int p[] ) /******************************************************************************/ /* Purpose: pruefer_to_tree_new() converts a Pruefer code to a tree. Discussion: The original code attempts to tack on an extra entry to P. Licensing: This code is distributed under the MIT license. Modified: 26 July 2011 Author: John Burkardt Reference: Donald Kreher, Douglas Simpson, Combinatorial Algorithms, CRC Press, 1998, ISBN: 0-8493-3988-X, LC: QA164.K73. Parameters: Input, int N, the number of nodes in the tree. N must be at least 3. Input, int P[N-2], the Pruefer code for the tree. Output, int PRUEFER_TO_TREE_NEW[2*(N-1)], describes the edges of the tree as pairs of nodes. */ { int *t; t = ( int * ) malloc ( 2 * ( n - 1 ) * sizeof ( int ) ); pruefer_to_tree ( n, p, t ); return t; } /******************************************************************************/ int *pruefer_unrank ( int rank, int n ) /******************************************************************************/ /* Purpose: pruefer_unrank() unranks a Pruefer code. Licensing: This code is distributed under the MIT license. Modified: 24 July 2011 Author: John Burkardt Reference: Donald Kreher, Douglas Simpson, Combinatorial Algorithms, CRC Press, 1998, ISBN: 0-8493-3988-X, LC: QA164.K73. Input: int RANK, the rank of the Pruefer code. int N, the number of nodes in the tree. N must be at least 3. Output: int P[N-2], the Pruefer code for the tree. */ { int i; int ncode; int *p; int rank_copy; /* Check. */ if ( n < 1 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "PRUEFER_UNRANK(): Fatal error!\n" ); fprintf ( stderr, " Input N is illegal.\n" ); exit ( 1 ); } if ( n < 3 ) { p = NULL; return p; } ncode = pruefer_enum ( n ); if ( rank < 0 || ncode < rank ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "PRUEFER_UNRANK(): Fatal error!\n" ); fprintf ( stderr, " The input rank is illegal.\n" ); exit ( 1 ); } rank_copy = rank; p = ( int * ) malloc ( ( n - 2 ) * sizeof ( int ) ); for ( i = n - 3; 0 <= i; i-- ) { p[i] = ( rank_copy % n ) + 1; rank_copy = ( rank_copy - p[i] + 1 ) / n; } return p; } /******************************************************************************/ void queens ( int n, int iarray[], int k, int *nstack, int istack[], int maxstack ) /******************************************************************************/ /* Purpose: QUEENS finds possible positions for the K-th nonattacking queen. Discussion: The chessboard is N by N, and is being filled one column at a time, with a tentative solution to the nonattacking queen problem. So far, K-1 rows have been chosen, and we now need to provide a list of all possible rows that might be used in column K. Licensing: This code is distributed under the MIT license. Modified: 26 July 2011 Author: John Burkardt Parameters: Input, int N, the total number of queens to place, and the length of a side of the chessboard. Input, int IARRAY[N]. The first K-1 entries of IARRAY record the rows into which queens have already been placed. Input, int K, the column for which we need possible row positions for the next queen. Input/output, int *NSTACK, the current length of stack. On output, this has been updated. Input/output, int ISTACK[MAXSTACK]. On output, we have added the candidates, and the number of candidates, to the end of the stack. Input, int MAXSTACK, maximum dimension of ISTACK. */ { int diag; int irow; int jcol; int ncan; int row; ncan = 0; for ( irow = 1; irow <= n; irow++ ) { /* If row IROW has already been used, that is it. */ row = 0; for ( jcol = 1; jcol <= k - 1; jcol++ ) { if ( iarray[jcol-1] == irow ) { row = 1; } } if ( !row ) { diag = 0; for ( jcol = 1; jcol <= k - 1; jcol++ ) { if ( irow == iarray[jcol-1] + k - jcol || irow == iarray[jcol-1] - ( k - jcol ) ) { diag = 1; } } if ( !diag ) { ncan = ncan + 1; istack[*nstack] = irow; *nstack = *nstack + 1; } } } istack[*nstack] = ncan; *nstack = *nstack + 1; return; } /******************************************************************************/ double r8_choose ( int n, int k ) /******************************************************************************/ /* Purpose: R8_CHOOSE computes the combinatorial coefficient C(N,K). Discussion: Real arithmetic is used, and C(N,K) is computed directly, via Gamma functions, rather than recursively. C(N,K) is the number of distinct combinations of K objects chosen from a set of N distinct objects. A combination is like a set, in that order does not matter. C(N,K) = N! / ( (N-K)! * K! ) Example: The number of combinations of 2 things chosen from 5 is 10. C(5,2) = ( 5 * 4 * 3 * 2 * 1 ) / ( ( 3 * 2 * 1 ) * ( 2 * 1 ) ) = 10. The actual combinations may be represented as: (1,2), (1,3), (1,4), (1,5), (2,3), (2,4), (2,5), (3,4), (3,5), (4,5). Licensing: This code is distributed under the MIT license. Modified: 25 July 2011 Author: John Burkardt Parameters: Input, int N, the value of N. Input, int K, the value of K. Output, double R8_CHOOSE, the value of C(N,K) */ { double arg; double fack; double facn; double facnmk; double value; if ( n < 0 ) { value = 0.0; } else if ( k == 0 ) { value = 1.0; } else if ( k == 1 ) { value = ( double ) ( n ); } else if ( 1 < k && k < n - 1 ) { arg = ( double ) ( n + 1 ); facn = lgamma ( arg ); arg = ( double ) ( k + 1 ); fack = lgamma ( arg ); arg = ( double ) ( n - k + 1 ); facnmk = lgamma ( arg ); value = r8_nint ( exp ( facn - fack - facnmk ) ); } else if ( k == n - 1 ) { value = ( double ) ( n ); } else if ( k == n ) { value = 1.0; } else { value = 0.0; } return value; } /******************************************************************************/ int r8_nint ( double x ) /******************************************************************************/ /* Purpose: r8_nint() returns the nearest integer to an R8. Example: X R8_NINT 1.3 1 1.4 1 1.5 1 or 2 1.6 2 0.0 0 -0.7 -1 -1.1 -1 -1.6 -2 Licensing: This code is distributed under the MIT license. Modified: 05 May 2006 Author: John Burkardt Parameters: Input, double X, the value. Output, int R8_NINT, the nearest integer to X. */ { int s; int value; if ( x < 0.0 ) { s = - 1; } else { s = + 1; } value = s * ( int ) ( fabs ( x ) + 0.5 ); return value; } /******************************************************************************/ void r8vec_backtrack ( int n, int maxstack, double stack[], double x[], int *indx, int *k, int *nstack, int ncan[] ) /******************************************************************************/ /* Purpose: R8VEC_BACKTRACK supervises a backtrack search for a real vector. Discussion: The routine tries to construct a real vector one index at a time, using possible candidates as supplied by the user. At any time, the partially constructed vector may be discovered to be unsatisfactory, but the routine records information about where the last arbitrary choice was made, so that the search can be carried out efficiently, rather than starting out all over again. First, call the routine with INDX = 0 so it can initialize itself. Now, on each return from the routine, if INDX is: 1, you've just been handed a complete candidate vector; Admire it, analyze it, do what you like. 2, please determine suitable candidates for position X(K). Return the number of candidates in NCAN(K), adding each candidate to the end of STACK, and increasing NSTACK. 3, you're done. Stop calling the routine; Licensing: This code is distributed under the MIT license. Modified: 13 July 2004 Author: Original FORTRAN77 version by Albert Nijenhuis, Herbert Wilf. C 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, int N, the number of positions to be filled in the vector. Input, int MAXSTACK, the maximum length of the stack. Input, double STACK[MAXSTACK], a list of all current candidates for all positions 1 through K. Input/output, double X[N], the partial or complete candidate vector. Input/output, int *INDX, a communication flag. On input, 0 to start a search. On output: 1, a complete output vector has been determined and returned in X(1:N); 2, candidates are needed for position X(K); 3, no more possible vectors exist. Inout/output, int *K, if INDX=2, the current vector index being considered. Input/output, int *NSTACK, the current length of the stack. Input/output, int NCAN[N], lists the current number of candidates for positions 1 through K. */ { /* If this is the first call, request a candidate for position 1. */ if ( *indx == 0 ) { *k = 1; *nstack = 0; *indx = 2; return; } /* Examine the stack. */ for ( ; ; ) { /* If there are candidates for position K, take the first available one off the stack, and increment K. This may cause K to reach the desired value of N, in which case we need to signal the user that a complete set of candidates is being returned. */ if ( 0 < ncan[(*k)-1] ) { x[(*k)-1] = stack[(*nstack)-1]; *nstack = *nstack - 1; ncan[(*k)-1] = ncan[(*k)-1] - 1; if ( *k != n ) { *k = *k + 1; *indx = 2; } else { *indx = 1; } break; } /* If there are no candidates for position K, then decrement K. If K is still positive, repeat the examination of the stack. */ else { *k = *k - 1; if ( *k <= 0 ) { *indx = 3; break; } } } return; } /******************************************************************************/ double r8vec_dot_product ( int n, double a1[], double a2[] ) /******************************************************************************/ /* Purpose: R8VEC_DOT_PRODUCT computes the dot product of a pair of R8VEC's. Licensing: This code is distributed under the MIT license. Modified: 26 July 2007 Author: John Burkardt Parameters: Input, int N, the number of entries in the vectors. Input, double A1[N], A2[N], the two vectors to be considered. Output, double R8VEC_DOT_PRODUCT, the dot product of the vectors. */ { int i; double value; value = 0.0; for ( i = 0; i < n; i++ ) { value = value + a1[i] * a2[i]; } return value; } /******************************************************************************/ int rgf_check ( int m, int f[] ) /******************************************************************************/ /* Purpose: RGF_CHECK checks a restricted growth function. Licensing: This code is distributed under the MIT license. Modified: 25 July 2011 Author: John Burkardt Reference: Donald Kreher, Douglas Simpson, Combinatorial Algorithms, CRC Press, 1998, ISBN: 0-8493-3988-X, LC: QA164.K73. Parameters: Input, int M, the domain of the RGF is the integers from 1 to M. M must be positive. Input, int F(M), the restricted growth function. Output, int RGF_CHECK. 1, the data is legal. 0, the data is not legal. */ { int check; int fmax; int i; check = 1; if ( m <= 0 ) { check = 0; return check; } fmax = 0; for ( i = 0; i < m; i++ ) { if ( f[i] <= 0 || fmax + 1 < f[i] ) { check = 0; return check; } fmax = i4_max ( fmax, f[i] ); } return check; } /******************************************************************************/ int rgf_enum ( int m ) /******************************************************************************/ /* Purpose: RGF_ENUM enumerates the restricted growth functions on M. Licensing: This code is distributed under the MIT license. Modified: 28 July 2011 Author: John Burkardt Reference: Donald Kreher, Douglas Simpson, Combinatorial Algorithms, CRC Press, 1998, ISBN: 0-8493-3988-X, LC: QA164.K73. Parameters: Input, int M, the domain of the RGF is the integers from 1 to M. M must be positive. However, for the enumeration routine only, it is legal to call with any value of M. Output, int RGF_ENUM, the number of restricted growth functions. */ { int *b; int i; int j; int value; if ( m < 0 ) { value = 0; } else if ( m == 0 ) { value = 1; } else { b = ( int * ) malloc ( ( m + 1 ) * sizeof ( int ) ); b[0] = 1; for ( j = 1; j <= m; j++ ) { b[j] = 0; for ( i = 0; i < j; i++ ) { b[j] = b[j] + i4_choose ( j - 1, i ) * b[i]; } } value = b[m]; free ( b ); } return value; } /******************************************************************************/ int *rgf_g_table ( int m ) /******************************************************************************/ /* Purpose: RGF_G_TABLE tabulates the generalized restricted growth functions. Example: M = 6 D = 1 1 1 1 1 1 1 1 2 3 4 5 6 0 2 5 10 17 26 0 0 5 15 37 77 0 0 0 15 52 151 0 0 0 0 52 203 0 0 0 0 0 203 0 0 0 0 0 0 Licensing: This code is distributed under the MIT license. Modified: 25 July 2011 Author: John Burkardt Reference: Donald Kreher, Douglas Simpson, Combinatorial Algorithms, CRC Press, 1998, ISBN: 0-8493-3988-X, LC: QA164.K73. Parameters: Input, int M, indicates how many rows and columns are to be computed. M must be nonnegative. Output, int RGF_G_TABLE[(M+1)*(M+1)], the first M+1 rows and M+1 columns of the table of the number of generalized restricted growth functions. D(I,J) is the number of GRGF''s of length I with restriction parameter J. */ { int *d; int i; int j; d = ( int * ) malloc ( ( m + 1 ) * ( m + 1 ) * sizeof ( int ) ); for ( j = 0; j <= m; j++ ) { d[0+j*(m+1)] = 1; } for ( i = 1; i <= m; i++ ) { for ( j = 0; j <= m; j++ ) { if ( j <= m - i ) { d[i+j*(m+1)] = j * d[i-1+j*(m+1)] + d[i-1+(j+1)*(m+1)]; } else { d[i+j*(m+1)] = 0; } } } return d; } /******************************************************************************/ int rgf_rank ( int m, int f[] ) /******************************************************************************/ /* Purpose: RGF_RANK ranks a restricted growth function. Licensing: This code is distributed under the MIT license. Modified: 25 July 2011 Author: John Burkardt Reference: Donald Kreher, Douglas Simpson, Combinatorial Algorithms, CRC Press, 1998, ISBN: 0-8493-3988-X, LC: QA164.K73. Parameters: Input, int M, the domain of the RGF is the integers from 1 to M. M must be positive. Input, int F[M], the restricted growth function. Output, int RGF_RANK, the rank of the restricted growth function. */ { int check; int *d; int i; int j; int rank; /* Check. */ check = rgf_check ( m, f ); if ( ! check ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "RGF_RANK(): Fatal error!\n" ); fprintf ( stderr, " The input array is illegal!\n" ); exit ( 1 ); } /* Get the generalized restricted growth function table. */ d = rgf_g_table ( m ); rank = 0; j = 1; for ( i = 2; i <= m; i++ ) { rank = rank + ( f[i-1] - 1 ) * d[m-i+j*(m+1)]; j = i4_max ( j, f[i-1] ); } free ( d ); return rank; } /******************************************************************************/ void rgf_successor ( int m, int f[], int *rank ) /******************************************************************************/ /* Purpose: RGF_SUCCESSOR generates the next restricted growth function. Licensing: This code is distributed under the MIT license. Modified: 26 July 2011 Author: John Burkardt Reference: Donald Kreher, Douglas Simpson, Combinatorial Algorithms, CRC Press, 1998, ISBN: 0-8493-3988-X, LC: QA164.K73. Parameters: Input, int M, the domain of the RGF is the integers from 1 to M. M must be positive. Input/output, int F[M], the restricted growth function. Input/output, int *RANK, the rank. If RANK = -1 on input, then the routine understands that this is the first call, and that the user wishes the routine to supply the first element in the ordering, which has RANK = 0. In general, the input value of RANK is increased by 1 for output, unless the very last element of the ordering was input, in which case the output value of RANK is 0. */ { int check; int fmax; int i; int j; /* Return the first element. */ if ( *rank == -1 ) { for ( i = 0; i < m; i++ ) { f[i] = 1; } *rank = 0; return; } /* Check. */ check = rgf_check ( m, f ); if ( ! check ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "RGF_SUCCESSOR(): Fatal error!\n" ); fprintf ( stderr, " The input array is illegal!\n" ); exit ( 1 ); } /* Find the first position from the right which can be incremented. */ for ( i = m; 2 <= i; i-- ) { fmax = 1; for ( j = 2; j < i; j++ ) { fmax = i4_max ( fmax, f[j-1] ); } /* Increment the function at this position, and set later entries to 1. */ if ( f[i-1] != fmax + 1 ) { f[i-1] = f[i-1] + 1; for ( j = i + 1; j <= m; j++ ) { f[j-1] = 1; } *rank = *rank + 1; return; } } /* The final element was input. Return the first element. */ for ( i = 0; i < m; i++ ) { f[i] = 1; } *rank = 0; return; } /******************************************************************************/ void rgf_to_setpart ( int m, int f[], int *nsub, int s[], int index[] ) /******************************************************************************/ /* Purpose: RGF_TO_SETPART converts a restricted growth function to a set partition. Licensing: This code is distributed under the MIT license. Modified: 26 July 2011 Author: John Burkardt Reference: Donald Kreher, Douglas Simpson, Combinatorial Algorithms, CRC Press, 1998, ISBN: 0-8493-3988-X, LC: QA164.K73. Parameters: Input, int M, the domain of the RGF is the integers from 1 to M. M must be positive. Input, int F[M], the restricted growth function. Output, int *NSUB, the number of nonempty subsets into which the set is partitioned. Output, int S[M], describes the partition of a set of M objects into NSUB nonempty subsets. If element I of the superset belongs to subset J, then S(I) = J. Output, int INDEX[M], lists the location in S of the last element of each subset. Thus, the elements of subset 1 are S(1) through S(INDEX(1)), the elements of subset 2 are S(INDEX(1)+1) through S(INDEX(2)) and so on. */ { int check; int i; int j; int k; /* Check. */ check = rgf_check ( m, f ); if ( ! check ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "RGF_TO_SETPART(): Fatal error!\n" ); fprintf ( stderr, " The input array is illegal!\n" ); exit ( 1 ); } /* Determine the number of subsets. */ *nsub = i4vec_max ( m, f ); /* Initialize. */ for ( i = 0; i < m; i++ ) { s[i] = 0; } for ( i = 0; i < *nsub; i++ ) { index[i] = 0; } /* For each subset I, collect the indices of F which have value I. These are the elements of the I-th subset. */ k = 0; for ( i = 1; i <= *nsub; i++ ) { for ( j = 1; j <= m; j++ ) { if ( f[j-1] == i ) { k = k + 1; s[k-1] = j; } } index[i-1] = k; } return; } /******************************************************************************/ int *rgf_unrank ( int rank, int m ) /******************************************************************************/ /* Purpose: RGF_UNRANK returns the restricted growth function of a given rank. Licensing: This code is distributed under the MIT license. Modified: 25 July 2011 Author: John Burkardt Reference: Donald Kreher, Douglas Simpson, Combinatorial Algorithms, CRC Press, 1998, ISBN: 0-8493-3988-X, LC: QA164.K73. Parameters: Input, int RANK, the rank of the restricted growth function. Input, int M, the domain of the RGF is the integers from 1 to M. M must be positive. Output, int RGF_UNRANK[M], the restricted growth function. */ { int *d; int *f; int i; int j; int nrgf; int rank_copy; /* Check. */ if ( m < 1 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "RGF_UNRANK(): Fatal error!\n" ); fprintf ( stderr, " Input M is illegal.\n" ); exit ( 1 ); } nrgf = rgf_enum ( m ); if ( rank < 0 || nrgf < rank ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "RGF_UNRANK(): Fatal error!\n" ); fprintf ( stderr, " The input rank is illegal.\n" ); exit ( 1 ); } /* Get the generalized restricted growth function table. */ d = rgf_g_table ( m ); f = ( int * ) malloc ( m * sizeof ( int ) ); rank_copy = rank; j = 1; f[0] = 1; for ( i = 2; i <= m; i++ ) { if ( j * d[m-i+j*(m+1)] <= rank_copy ) { f[i-1] = j + 1; rank_copy = rank_copy - j * d[m-i+j*(m+1)]; j = j + 1; } else { f[i-1] = 1 + ( rank_copy / d[m-i+j*(m+1)] ); rank_copy = rank_copy % d[m-i+j*(m+1)]; } } free ( d ); return f; } /******************************************************************************/ int s_len_trim ( char *s ) /******************************************************************************/ /* Purpose: S_LEN_TRIM returns the length of a string to the last nonblank. Licensing: This code is distributed under the MIT license. Modified: 26 April 2003 Author: John Burkardt Parameters: Input, char *S, a pointer to a string. Output, int S_LEN_TRIM, the length of the string to the last nonblank. If S_LEN_TRIM is 0, then the string is entirely blank. */ { int n; char *t; n = strlen ( s ); t = s + strlen ( s ) - 1; while ( 0 < n ) { if ( *t != ' ' ) { return n; } t--; n--; } return n; } /******************************************************************************/ int setpart_check ( int m, int nsub, int s[], int index[] ) /******************************************************************************/ /* Purpose: SETPART_CHECK checks a set partition. Licensing: This code is distributed under the MIT license. Modified: 25 July 2011 Author: John Burkardt Reference: Donald Kreher, Douglas Simpson, Combinatorial Algorithms, CRC Press, 1998, ISBN: 0-8493-3988-X, LC: QA164.K73. Parameters: Input, int M, the number of elements of the set. M must be positive. Input, int NSUB, the number of nonempty subsets into which the set is partitioned. 1 <= NSUB <= M. Input, int S[M], contains the integers from 1 to M, grouped into subsets as described by INDEX. Input, int INDEX[NSUB], lists the location in S of the last element of each subset. Thus, the elements of subset 1 are S(1) through S(INDEX(1)), the elements of subset 2 are S(INDEX(1)+1) through S(INDEX(2)) and so on. Output, int SETPART_CHECK. 1, the data is legal. 0, the data is not legal. */ { int check; int i; int imin; int j; check = 1; /* Check M. */ if ( m < 1 ) { check = 0; return check; } /* Check NSUB. */ if ( nsub < 1 ) { check = 0; return check; } /* Check INDEX. */ imin = 0; for ( i = 0; i < nsub; i++ ) { if ( index[i] <= imin || m < index[i] ) { check = 0; return check; } imin = index[i]; } /* Check the elements of S. */ for ( i = 0; i < m; i++ ) { if ( s[i] <= 0 || m < s[i] ) { check = 0; return check; } for ( j = 0; j < i; j++ ) { if ( s[j] == s[i] ) { check = 0; return check; } } } return check; } /******************************************************************************/ int setpart_enum ( int m ) /******************************************************************************/ /* Purpose: SETPART_ENUM enumerates the partitions of a set of M elements. Licensing: This code is distributed under the MIT license. Modified: 27 July 2011 Author: John Burkardt Reference: Donald Kreher, Douglas Simpson, Combinatorial Algorithms, CRC Press, 1998, ISBN: 0-8493-3988-X, LC: QA164.K73. Parameters: Input, int M, the number of elements in the set. M must be positive. However, for the enumeration routine only, it is legal to call with any value of M. Output, int SETPART_ENUM, the number of partitions of the set. */ { int *b; int i; int j; int value; if ( m < 0 ) { value = 0; } else if ( m == 0 ) { value = 1; } else { b = ( int * ) malloc ( ( m + 1 ) * sizeof ( int ) ); b[0] = 1; for ( j = 1; j <= m; j++ ) { b[j] = 0; for ( i = 0; i < j; i++ ) { b[j] = b[j] + i4_choose ( j - 1, i ) * b[i]; } } value = b[m]; free ( b ); } return value; } /******************************************************************************/ int *setpart_to_rgf ( int m, int nsub, int s[], int index[] ) /******************************************************************************/ /* Purpose: SETPART_TO_RGF converts a set partition to a restricted growth function. Licensing: This code is distributed under the MIT license. Modified: 26 July 2011 Author: John Burkardt Reference: Donald Kreher, Douglas Simpson, Combinatorial Algorithms, CRC Press, 1998, ISBN: 0-8493-3988-X, LC: QA164.K73. Parameters: Input, int M, the number of elements of the set. M must be positive. Input, int NSUB, the number of nonempty subsets into which the set is partitioned. 1 <= NSUB <= M. Input, int S(M), contains the integers from 1 to M, grouped into subsets as described by INDEX. Input, int INDEX(NSUB), lists the location in S of the last element of each subset. Thus, the elements of subset 1 are S(1) through S(INDEX(1)), the elements of subset 2 are S(INDEX(1)+1) through S(INDEX(2)) and so on. Output, int SETPART_TO_RGF[M], the restricted growth function from M to NSUB. */ { int check; int *f; int i; int k; int khi; int klo; /* Check. */ check = setpart_check ( m, nsub, s, index ); if ( ! check ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "SETPART_TO_RGF(): Fatal error!\n" ); fprintf ( stderr, " The input array is illegal.\n" ); exit ( 1 ); } f = ( int * ) malloc ( m * sizeof ( int ) ); khi = 0; for ( i = 1; i <= nsub; i++ ) { klo = khi + 1; khi = index[i-1]; for ( k = klo; k <= khi; k++ ) { f[s[k-1]-1] = i; } } return f; } /******************************************************************************/ int *stirling_numbers1 ( int m, int n ) /******************************************************************************/ /* Purpose: stirling_numbers1() computes Stirling numbers of the first kind. Licensing: This code is distributed under the MIT license. Modified: 25 July 2011 Author: John Burkardt Reference: Donald Kreher, Douglas Simpson, Combinatorial Algorithms, CRC Press, 1998, ISBN: 0-8493-3988-X, LC: QA164.K73. Parameters: Input, int M, the maximum row to compute. M must be nonnegative. Input, int N, the maximum column to compute. N must be nonnegative. Output, int S[(M+1)*(N+1)], the first M+1 rows and N+1 columns of the table of Stirling numbers of the first kind. */ { int i; int j; int *s; s = ( int * ) malloc ( ( m + 1 ) * ( n + 1 ) * sizeof ( int ) ); s[0+0*(m+1)] = 1; for ( j = 1; j <= n; j++ ) { s[0+j*(m+1)] = 0; } for ( i = 1; i <= m; i++ ) { s[i+0*(m+1)] = 0; } /* This loop may be extraneous. */ for ( i = 0; i <= i4_min ( m, n - 1 ); i++ ) { s[i+(i+1)*(m+1)] = 0; } for ( i = 1; i <= m; i++ ) { for ( j = 1; j <= n; j++ ) { if ( j <= i ) { s[i+j*(m+1)] = s[i-1+(j-1)*(m+1)] - ( i - 1 ) * s[i-1+j*(m+1)]; } else { s[i+j*(m+1)] = 0; } } } return s; } /******************************************************************************/ int *stirling_numbers2 ( int m, int n ) /******************************************************************************/ /* Purpose: stirling_numbers2() computes Stirling numbers of the second kind. Discussion: The reference has a typographical error, referring to S(I-J,J-1) instead of S(I-1,J-1). Licensing: This code is distributed under the MIT license. Modified: 25 July 2011 Author: John Burkardt Reference: Donald Kreher, Douglas Simpson, Combinatorial Algorithms, CRC Press, 1998, ISBN: 0-8493-3988-X, LC: QA164.K73. Parameters: Input, int M, the maximum row to compute. M must be nonnegative. Input, int N, the maximum column to compute. N must be nonnegative. Output, int S[(M+1)*(N+1)], the first M+1 rows and N+1 columns of the table of Stirling numbers of the second kind. */ { int i; int j; int *s; s = ( int * ) malloc ( ( m + 1 ) * ( n + 1 ) * sizeof ( int ) ); s[0+0*(m+1)] = 1; for ( j = 1; j <= n; j++ ) { s[0+j*(m+1)] = 0; } for ( i = 1; i <= m; i++ ) { s[i+0*(m+1)] = 0; } /* This loop may be extraneous. */ for ( i = 0; i <= i4_min ( m, n - 1 ); i++ ) { s[i+(i+1)*(m+1)] = 0; } for ( i = 1; i <= m; i++ ) { for ( j = 1; j <= n; j++ ) { if ( j <= i ) { s[i+j*(m+1)] = j * s[i-1+j*(m+1)] + s[i-1+(j-1)*(m+1)]; } else { s[i+j*(m+1)] = 0; } } } return s; } /******************************************************************************/ int subset_check ( int n, int t[] ) /******************************************************************************/ /* Purpose: SUBSET_CHECK checks a subset. Licensing: This code is distributed under the MIT license. Modified: 28 July 2011 Author: John Burkardt Reference: Donald Kreher, Douglas Simpson, Combinatorial Algorithms, CRC Press, 1998, ISBN: 0-8493-3988-X, LC: QA164.K73. Parameters: Input, int N, the number of elements in the master set. N must be positive. Input, int T[N], the subset. If T(I) = 0, item I is not in the subset; if T(I) = 1, item I is in the subset. Output, int SUBSET_CHECK. 1, the data is legal. 0, the data is not legal. */ { int check; int i; check = 1; if ( n < 1 ) { check = 0; return check; } for ( i = 0; i < n; i++ ) { if ( t[i] != 0 && t[i] != 1 ) { check = 0; return check; } } return check; } /******************************************************************************/ int subset_colex_rank ( int n, int t[] ) /******************************************************************************/ /* Purpose: SUBSET_COLEX_RANK computes the colexicographic rank of a subset. Licensing: This code is distributed under the MIT license. Modified: 22 August 2011 Author: John Burkardt Reference: Donald Kreher, Douglas Simpson, Combinatorial Algorithms, CRC Press, 1998, ISBN: 0-8493-3988-X, LC: QA164.K73. Parameters: Input, int N, the number of items in the master set. N must be positive. Input, int T[N], the subset. If T(I) = 0, item I is not in the subset; if T(I) = 1, item I is in the subset. Output, int SUBSET_COLEX_RANK, the rank of the subset. */ { int check; int i; int rank; /* Check. */ check = subset_check ( n, t ); if ( ! check ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "SUBSET_COLEX_RANK(): Fatal error!\n" ); fprintf ( stderr, " The subset is illegal.\n" ); exit ( 1 ); } rank = 0; for ( i = 0; i < n; i++ ) { if ( t[i] == 1 ) { rank = rank + i4_power ( 2, i ); } } return rank; } /******************************************************************************/ void subset_colex_successor ( int n, int t[], int *rank ) /******************************************************************************/ /* Purpose: SUBSET_COLEX_SUCCESSOR computes the subset colexicographic successor. Discussion: In the original code, there is a last element with no successor. Licensing: This code is distributed under the MIT license. Modified: 26 July 2011 Author: John Burkardt Reference: Donald Kreher, Douglas Simpson, Combinatorial Algorithms, CRC Press, 1998, ISBN: 0-8493-3988-X, LC: QA164.K73. Parameters: Input, int N, the number of elements in the master set. N must be positive. Input/output, int T[N], describes a subset. T(I) is 0 if the I-th element of the master set is not in the subset, and is 1 if the I-th element is part of the subset. On input, T describes a subset. On output, T describes the next subset in the ordering. If the input T was the last in the ordering, then the output T will be the first. Input/output, int *RANK, the rank. If RANK = -1 on input, then the routine understands that this is the first call, and that the user wishes the routine to supply the first element in the ordering, which has RANK = 0. In general, the input value of RANK is increased by 1 for output, unless the very last element of the ordering was input, in which case the output value of RANK is 0. */ { int check; int i; /* Return the first element. */ if ( *rank == -1 ) { for ( i = 0; i < n; i++ ) { t[i] = 0; } *rank = 0; return; } /* Check. */ check = subset_check ( n, t ); if ( ! check ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "SUBSET_COLEX_RANK(): Fatal error!\n" ); fprintf ( stderr, " The subset is illegal.\n" ); exit ( 1 ); } for ( i = 0; i < n; i++ ) { if ( t[i] == 0 ) { t[i] = 1; *rank = *rank + 1; return; } else { t[i] = 0; } } *rank = 0; return; } /******************************************************************************/ int *subset_colex_unrank ( int rank, int n ) /******************************************************************************/ /* Purpose: SUBSET_COLEX_UNRANK computes the subset of given colexicographic rank. Licensing: This code is distributed under the MIT license. Modified: 25 July 2011 Author: John Burkardt Reference: Donald Kreher, Douglas Simpson, Combinatorial Algorithms, CRC Press, 1998, ISBN: 0-8493-3988-X, LC: QA164.K73. Parameters: Input, int RANK, the rank of the subset. Input, int N, the number of items in the master set. N must be positive. Output, int SUBSET_COLEX_UNRANK[N], the subsetof the given rank. If T(I) = 0, item I is not in the subset; if T(I) = 1, item I is in the subset. */ { int i; int nsub; int rank_copy; int *t; /* Check. */ if ( n < 1 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "SUBSET_COLEX_UNRANK(): Fatal error!\n" ); fprintf ( stderr, " Input N is illegal.\n" ); exit ( 1 ); } nsub = subset_enum ( n ); if ( rank < 0 || nsub < rank ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "SUBSET_COLEX_UNRANK(): Fatal error!\n" ); fprintf ( stderr, " The input rank is illegal.\n" ); exit ( 1 ); } rank_copy = rank; t = ( int * ) malloc ( n * sizeof ( int ) ); for ( i = 0; i < n; i++ ) { if ( ( rank_copy % 2 ) == 1 ) { t[i] = 1; } else { t[i] = 0; } rank_copy = rank_copy / 2; } return t; } /******************************************************************************/ int *subset_complement ( int n, int a[] ) /******************************************************************************/ /* Purpose: SUBSET_COMPLEMENT computes the complement of a set. Licensing: This code is distributed under the MIT license. Modified: 25 July 2011 Author: John Burkardt Reference: Donald Kreher, Douglas Simpson, Combinatorial Algorithms, CRC Press, 1998, ISBN: 0-8493-3988-X, LC: QA164.K73. Parameters: Input, int N, the order of the master set, of which A is a subset. N must be positive. Input, int A[N], a subset of the master set. A(I) = 0 if the I-th element is in the subset A, and is 1 otherwise. Output, int SUBSET_COMPLEMENT[N], the complement of A. */ { int *b; int check; int i; /* Check. */ check = subset_check ( n, a ); if ( ! check ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "SUBSET_COMPLEMENT(): Fatal error!\n" ); fprintf ( stderr, " The subset is illegal.\n" ); exit ( 1 ); } b = ( int * ) malloc ( n * sizeof ( int ) ); for ( i = 0; i < n; i++ ) { b[i] = 1 - a[i]; } return b; } /******************************************************************************/ int subset_distance ( int n, int t1[], int t2[] ) /******************************************************************************/ /* Purpose: SUBSET_DISTANCE computes the Hamming distance between two sets. Discussion: The sets T1 and T2 are assumed to be subsets of a set of N elements. Licensing: This code is distributed under the MIT license. Modified: 25 July 2011 Author: John Burkardt Reference: Donald Kreher, Douglas Simpson, Combinatorial Algorithms, CRC Press, 1998, ISBN: 0-8493-3988-X, LC: QA164.K73. Parameters: Input, int N, the order of the master set, of which T1 and T2 are subsets. N must be positive. Input, int T1[N], T2[N], two subsets of the master set. T1(I) = 0 if the I-th element is in the subset T1, and is 1 otherwise; T2 is defined similarly. Output, int SUBSET_DISTANCE, the Hamming distance between T1 and T2, defined as the number of elements of the master set which are in either T1 or T2 but not both. */ { int check; int dist; int i; /* Check. */ check = subset_check ( n, t1 ); if ( ! check ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "SUBSET_DISTANCE(): Fatal error!\n" ); fprintf ( stderr, " The subset is illegal.\n" ); exit ( 1 ); } check = subset_check ( n, t2 ); if ( ! check ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "SUBSET_DISTANCE(): Fatal error!\n" ); fprintf ( stderr, " The subset is illegal.\n" ); exit ( 1 ); } dist = 0; for ( i = 0; i < n; i++ ) { if ( ( t1[i] == 0 && t2[i] != 0 ) || ( t1[i] != 0 && t2[i] == 0 ) ) { dist = dist + 1; } } return dist; } /******************************************************************************/ int subset_enum ( int n ) /******************************************************************************/ /* Purpose: SUBSET_ENUM enumerates the subsets of a set with N elements. Licensing: This code is distributed under the MIT license. Modified: 24 July 2011 Author: John Burkardt Parameters: Input, int N, the number of elements in the set. N must be at least 0. Output, int SUBSET_ENUM, the number of distinct elements. */ { int value; value = i4_power ( 2, n ); return value; } /******************************************************************************/ int *subset_intersect ( int n, int a[], int b[] ) /******************************************************************************/ /* Purpose: SUBSET_INTERSECT computes the intersection of two sets. Licensing: This code is distributed under the MIT license. Modified: 25 July 2011 Author: John Burkardt Reference: Donald Kreher, Douglas Simpson, Combinatorial Algorithms, CRC Press, 1998, ISBN: 0-8493-3988-X, LC: QA164.K73. Parameters: Input, int N, the order of the master set, of which A and B are subsets. N must be positive. Input, int A[N], B[N], two subsets of the master set. A(I) = 0 if the I-th element is in the subset A, and is 1 otherwise; B is defined similarly. Output, int SUBSET_INTERSECT[N], the intersection of A and B. */ { int *c; int check; int i; /* Check. */ check = subset_check ( n, a ); if ( ! check ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "SUBSET_COLEX_RANK(): Fatal error!\n" ); fprintf ( stderr, " The subset is illegal.\n" ); exit ( 1 ); } check = subset_check ( n, b ); if ( ! check ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "SUBSET_COLEX_RANK(): Fatal error!\n" ); fprintf ( stderr, " The subset is illegal.\n" ); exit ( 1 ); } c = ( int * ) malloc ( n * sizeof ( int ) ); for ( i = 0; i < n; i++ ) { c[i] = i4_min ( a[i], b[i] ); } return c; } /******************************************************************************/ int subset_lex_rank ( int n, int t[] ) /******************************************************************************/ /* Purpose: SUBSET_LEX_RANK computes the lexicographic rank of a subset. Licensing: This code is distributed under the MIT license. Modified: 22 August 2011 Author: John Burkardt Reference: Donald Kreher, Douglas Simpson, Combinatorial Algorithms, CRC Press, 1998, ISBN: 0-8493-3988-X, LC: QA164.K73. Parameters: Input, int N, the number of items in the master set. N must be positive. Input, int T[N], the subset. If T(I) = 0, item I is not in the subset; if T(I) = 1, item I is in the subset. Output, int SUBSET_LEX_RANK, the rank of the subset. */ { int check; int i; int rank; /* Check. */ check = subset_check ( n, t ); if ( ! check ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "SUBSET_LEX_RANK(): Fatal error!\n" ); fprintf ( stderr, " The subset is illegal.\n" ); exit ( 1 ); } rank = 0; for ( i = 0; i < n; i++ ) { if ( t[i] == 1 ) { rank = rank + i4_power ( 2, n - i - 1 ); } } return rank; } /******************************************************************************/ void subset_lex_successor ( int n, int t[], int *rank ) /******************************************************************************/ /* Purpose: SUBSET_LEX_SUCCESSOR computes the subset lexicographic successor. Discussion: In the original code, there is a last element with no successor. Licensing: This code is distributed under the MIT license. Modified: 27 July 2011 Author: John Burkardt Reference: Donald Kreher, Douglas Simpson, Combinatorial Algorithms, CRC Press, 1998, ISBN: 0-8493-3988-X, LC: QA164.K73. Parameters: Input, int N, the number of elements in the master set. N must be positive. Input/output, int T[N], describes a subset. T(I) is 0 if the I-th element of the master set is not in the subset, and is 1 if the I-th element is part of the subset. On input, T describes a subset. On output, T describes the next subset in the ordering. If the input T was the last in the ordering, then the output T will be the first. Input/output, int *RANK, the rank. If RANK = -1 on input, then the routine understands that this is the first call, and that the user wishes the routine to supply the first element in the ordering, which has RANK = 0. In general, the input value of RANK is increased by 1 for output, unless the very last element of the ordering was input, in which case the output value of RANK is 0. */ { int check; int i; /* Return the first element. */ if ( *rank == -1 ) { for ( i = 0; i < n; i++ ) { t[i] = 0; } *rank = 0; return; } /* Check. */ check = subset_check ( n, t ); if ( ! check ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "SUBSET_LEX_SUCCESSOR(): Fatal error!\n" ); fprintf ( stderr, " The subset is illegal.\n" ); exit ( 1 ); } for ( i = n - 1; 0 <= i; i-- ) { if ( t[i] == 0 ) { t[i] = 1; *rank = *rank + 1; return; } else { t[i] = 0; } } *rank = 0; return; } /******************************************************************************/ int *subset_lex_unrank ( int rank, int n ) /******************************************************************************/ /* Purpose: SUBSET_LEX_UNRANK computes the subset of given lexicographic rank. Licensing: This code is distributed under the MIT license. Modified: 25 July 2011 Author: John Burkardt Reference: Donald Kreher, Douglas Simpson, Combinatorial Algorithms, CRC Press, 1998, ISBN: 0-8493-3988-X, LC: QA164.K73. Parameters: Input, int RANK, the rank of the subset. Input, int N, the number of items in the master set. N must be positive. Output, int SUBSET_LEX_UNRANK[N], the subset of the given rank. If T(I) = 0, item I is not in the subset; if T(I) = 1, item I is in the subset. */ { int i; int nsub; int rank_copy; int *t; /* Check. */ if ( n < 1 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "SUBSET_LEX_UNRANK(): Fatal error!\n" ); fprintf ( stderr, " Input N is illegal.\n" ); exit ( 1 ); } nsub = subset_enum ( n ); if ( rank < 0 || nsub < rank ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "SUBSET_LEX_UNRANK(): Fatal error!\n" ); fprintf ( stderr, " The input rank is illegal.\n" ); exit ( 1 ); } rank_copy = rank; t = ( int * ) malloc ( n * sizeof ( int ) ); for ( i = n - 1; 0 <= i; i-- ) { if ( ( rank_copy % 2 ) == 1 ) { t[i] = 1; } else { t[i] = 0; } rank_copy = rank_copy / 2; } return t; } /******************************************************************************/ int *subset_random ( int n, int *seed ) /******************************************************************************/ /* Purpose: subset_random() returns a random subset. Licensing: This code is distributed under the MIT license. Modified: 24 December 2015 Author: John Burkardt Parameters: Input, int N, the size of the set. Input/output, int *SEED, a seed for the random number generator. Output, int SUBSET_RANDOM[N], defines the subset using 0 and 1 values. */ { int *s; s = i4vec_uniform_ab_new ( n, 0, 1, seed ); return s; } /******************************************************************************/ int *subset_union ( int n, int a[], int b[] ) /******************************************************************************/ /* Purpose: SUBSET_UNION computes the union of two sets. Licensing: This code is distributed under the MIT license. Modified: 25 July 2011 Author: John Burkardt Reference: Donald Kreher, Douglas Simpson, Combinatorial Algorithms, CRC Press, 1998, ISBN: 0-8493-3988-X, LC: QA164.K73. Parameters: Input, int N, the order of the master set, of which A and B are subsets. N must be positive. Input, int A[N], B[N], two subsets of the master set. A(I) = 0 if the I-th element is in the subset A, and is 1 otherwise; B is defined similarly. Output, int SUBSET_UNION[N], the union of A and B. */ { int *c; int check; int i; /* Check. */ check = subset_check ( n, a ); if ( ! check ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "SUBSET_UNION(): Fatal error!\n" ); fprintf ( stderr, " The subset is illegal.\n" ); exit ( 1 ); } check = subset_check ( n, b ); if ( ! check ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "SUBSET_UNION(): Fatal error!\n" ); fprintf ( stderr, " The subset is illegal.\n" ); exit ( 1 ); } c = ( int * ) malloc ( n * sizeof ( int ) ); for ( i = 0; i < n; i++ ) { c[i] = i4_max ( a[i], b[i] ); } return c; } /******************************************************************************/ int subset_weight ( int n, int t[] ) /******************************************************************************/ /* Purpose: SUBSET_WEIGHT computes the Hamming weight of a set. Discussion: The Hamming weight is simply the number of elements in the set. Licensing: This code is distributed under the MIT license. Modified: 24 July 2011 Author: John Burkardt Reference: Donald Kreher, Douglas Simpson, Combinatorial Algorithms, CRC Press, 1998, ISBN: 0-8493-3988-X, LC: QA164.K73. Parameters: Input, int N, the order of the master set, of which T is a subset. N must be positive. Input, int T[N], defines the subset T. T(I) is 1 if I is an element of T, and 0 otherwise. Output, int SUBSET_WEIGHT, the Hamming weight of the subset T. */ { int check; int weight; /* Check. */ check = subset_check ( n, t ); if ( ! check ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "SUBSET_WEIGHT(): Fatal error!\n" ); fprintf ( stderr, " The subset is illegal.\n" ); exit ( 1 ); } weight = i4vec_sum ( n, t ); return weight; } /******************************************************************************/ int *subset_xor ( int n, int a[], int b[] ) /******************************************************************************/ /* Purpose: SUBSET_XOR computes the symmetric difference of two sets. Licensing: This code is distributed under the MIT license. Modified: 25 July 2011 Author: John Burkardt Reference: Donald Kreher, Douglas Simpson, Combinatorial Algorithms, CRC Press, 1998, ISBN: 0-8493-3988-X, LC: QA164.K73. Parameters: Input, int N, the order of the master set, of which A and B are subsets. N must be positive. Input, int A[N], B[N], two subsets of the master set. A(I) = 0 if the I-th element is in the subset A, and is 1 otherwise; B is defined similarly. Output, int SUBSET_XOR[N], the symmetric difference of A and B. */ { int *c; int check; int i; /* Check. */ check = subset_check ( n, a ); if ( ! check ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "SUBSET_XOR(): Fatal error!\n" ); fprintf ( stderr, " The subset is illegal.\n" ); exit ( 1 ); } check = subset_check ( n, b ); if ( ! check ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "SUBSET_XOR(): Fatal error!\n" ); fprintf ( stderr, " The subset is illegal.\n" ); exit ( 1 ); } c = ( int * ) malloc ( n * sizeof ( int ) ); for ( i = 0; i < n; i++ ) { c[i] = i4_max ( a[i], b[i] ) - i4_min ( a[i], b[i] ); } return c; } /******************************************************************************/ int subsetsum_swap ( int n, int a[], int sum_desired, int index[] ) /******************************************************************************/ /* Purpose: SUBSETSUM_SWAP seeks a solution of the subset sum problem by swapping. Discussion: Given a collection of N not necessarily distinct positive integers A(I), and a positive integer SUM_DESIRED, select a subset of the values so that their sum is as close as possible to SUM_DESIRED without exceeding it. Algorithm: Start with no values selected, and SUM_ACHIEVED = 0. Consider each element A(I): If A(I) is not selected and SUM_ACHIEVED + A(I) <= SUM_DESIRED, select A(I). If A(I) is still not selected, and there is a selected A(J) such that SUM_GOT < SUM_ACHIEVED + A(I) - A(J), select A(I) and deselect A(J). If no items were selected on this sweep, exit. Otherwise, repeat the search. Licensing: This code is distributed under the MIT license. Modified: 25 July 2011 Author: John Burkardt Reference: Donald Kreher, Douglas Simpson, Combinatorial Algorithms, CRC Press, 1998, ISBN: 0-8493-3988-X, LC: QA164.K73. Parameters: Input, int N, the number of values. N must be positive. Input/output, int A[N], a collection of positive values. On output, A has been sorted into descending order. Input, int SUM_DESIRED, the desired sum. Output, int INDEX[N]; INDEX(I) is 1 if A(I) is part of the sum, and 0 otherwise. Output, int SUBSETSUM_SWAP, the sum of the selected elements. */ { int i; int j; int nmove; int sum_achieved; /* Initialize. */ sum_achieved = 0; for ( i = 0; i < n; i++ ) { index[i] = 0; } /* Sort into descending order. */ i4vec_sort_insert_d ( n, a ); for ( ; ; ) { nmove = 0; for ( i = 0; i < n; i++ ) { if ( index[i] == 0 ) { if ( sum_achieved + a[i] <= sum_desired ) { index[i] = 1; sum_achieved = sum_achieved + a[i]; nmove = nmove + 1; continue; } } if ( index[i] == 0 ) { for ( j = 0; j < n; j++ ) { if ( index[j] == 1 ) { if ( sum_achieved < sum_achieved + a[i] - a[j] && sum_achieved + a[i] - a[j] <= sum_desired ) { index[j] = 0; index[i] = 1; nmove = nmove + 2; sum_achieved = sum_achieved + a[i] - a[j]; break; } } } } } if ( nmove <= 0 ) { break; } } return sum_achieved; } /******************************************************************************/ int tableau_check ( int n, int tab[] ) /******************************************************************************/ /* Purpose: TABLEAU_CHECK checks a 2 by N tableau. Licensing: This code is distributed under the MIT license. Modified: 25 December 2015 Author: John Burkardt Reference: Donald Kreher, Douglas Simpson, Combinatorial Algorithms, CRC Press, 1998, ISBN: 0-8493-3988-X, LC: QA164.K73. Parameters: Input, int N, the number of columns in the tableau. N must be positive. Input, int TAB[2*N], a 2 by N tableau. Output, int TABLEAU_CHECK. 1, the data is legal. 0, the data is not legal. */ { int check; int i; int j; check = 1; if ( n < 1 ) { check = 0; return check; } /* The entries must be between 1 and 2*N. */ for ( i = 0; i < 2; i++ ) { for ( j = 0; j < n; j++ ) { if ( tab[i+j*2] < 1 || 2 * n < tab[i+j*2] ) { check = 0; return check; } } } /* The entries must be increasing to the right. */ for ( i = 0; i < 2; i++ ) { for ( j = 1; j < n; j++ ) { if ( tab[i+j*2] <= tab[i+(j-1)*2] ) { check = 0; return check; } } } /* The entries must be increasing down. */ i = 1; for ( j = 0; j < n; j++ ) { if ( tab[i+j*2] <= tab[i-1+j*2] ) { check = 0; return check; } } return check; } /******************************************************************************/ int tableau_enum ( int n ) /******************************************************************************/ /* Purpose: TABLEAU_ENUM enumerates the 2 by N standard tableaus. Licensing: This code is distributed under the MIT license. Modified: 24 July 2011 Author: John Burkardt Reference: Donald Kreher, Douglas Simpson, Combinatorial Algorithms, CRC Press, 1998, ISBN: 0-8493-3988-X, LC: QA164.K73. Parameters: Input, int N, the number of columns in the tableau. N must be nonnegative. Output, int TABLEAU_ENUM, the number of 2 by N standard tableaus. */ { int value; value = i4_choose ( 2 * n, n ) / ( n + 1 ); return value; } /******************************************************************************/ int *tableau_to_bal_seq ( int n, int tab[] ) /******************************************************************************/ /* Purpose: tableau_to_bal_seq() converts a 2 by N tableau to a balanced sequence. Licensing: This code is distributed under the MIT license. Modified: 25 July 2011 Author: John Burkardt Reference: Donald Kreher, Douglas Simpson, Combinatorial Algorithms, CRC Press, 1998, ISBN: 0-8493-3988-X, LC: QA164.K73. Input: int N, the number of 0's (and 1's) in the sequence. N must be positive. int TAB[2*N], a 2 by N tableau. Output: int TABLEAU_TO_BAL_SEQ[2*N], a balanced sequence. */ { int check; int i; int j; int *t; /* Check. */ check = tableau_check ( n, tab ); if ( ! check ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "tableau_to_seq(): Fatal error!\n" ); fprintf ( stderr, " The tableau is illegal.\n" ); exit ( 1 ); } t = ( int * ) malloc ( 2 * n * sizeof ( int ) ); for ( i = 0; i < 2; i++ ) { for ( j = 0; j < n; j++ ) { t[tab[i+j*2]-1] = i; } } return t; } /******************************************************************************/ int tree_check ( int n, int t[] ) /******************************************************************************/ /* Purpose: tree_check() checks a tree. Licensing: This code is distributed under the MIT license. Modified: 29 July 2011 Author: John Burkardt Reference: Donald Kreher, Douglas Simpson, Combinatorial Algorithms, CRC Press, 1998, ISBN: 0-8493-3988-X, LC: QA164.K73. Input: int N, the number of nodes in the tree. N must be positive. int T[2*(N-1)], describes the edges of the tree as pairs of nodes. Output: int TREE_CHECK. 1, the data is legal. 0, the data is not legal. */ { int check; int *d; int i; int j; int k; int x; int y; check = 1; if ( n < 1 ) { check = 0; return check; } for ( i = 0; i < 2; i++ ) { for ( j = 0; j < n - 1; j++ ) { if ( t[i+j*2] < 1 || n < t[i+j*2] ) { check = 0; return check; } } } /* Compute the degree of each node. */ d = edge_degree ( n, n - 1, t ); /* Delete a node of degree 1, N-1 times. */ for ( k = 1; k <= n - 1; k++ ) { x = 1; while ( d[x-1] != 1 ) { x = x + 1; if ( n < x ) { check = 0; return check; } } /* Find its neighbor. */ j = 1; for ( ; ; ) { if ( t[0+(j-1)*2] == x ) { y = t[1+(j-1)*2]; break; } if ( t[1+(j-1)*2] == x ) { y = t[0+(j-1)*2]; break; } j = j + 1; if ( n - 1 < j ) { check = 0; return check; } } /* Delete the edge. */ t[0+(j-1)*2] = - t[0+(j-1)*2]; t[1+(j-1)*2] = - t[1+(j-1)*2]; d[x-1] = d[x-1] - 1; d[y-1] = d[y-1] - 1; } for ( j = 0; j < n - 1; j++ ) { for ( i = 0; i < 2; i++ ) { t[i+j*2] = - t[i+j*2]; } } free ( d ); return check; } /******************************************************************************/ int tree_enum ( int n ) /******************************************************************************/ /* Purpose: tree_enum() enumerates the trees on N nodes. Licensing: This code is distributed under the MIT license. Modified: 24 July 2011 Author: John Burkardt Input: int N, the number of nodes in each tree. N must normally be at least 3, but for this routine, any value of N is allowed. Output: int TREE_ENUM, the number of distinct elements. */ { int value; if ( n < 1 ) { value = 0; } else if ( n == 1 ) { value = 1; } else if ( n == 2 ) { value = 1; } else { value = i4_power ( n, n - 2 ); } return value; } /******************************************************************************/ int *tree_random ( int n ) /******************************************************************************/ /* Purpose: tree_random() randomly selects a tree on N vertices. Licensing: This code is distributed under the MIT license. Modified: 11 September 2022 Author: John Burkardt Reference: Donald Kreher, Douglas Simpson, Combinatorial Algorithms, CRC Press, 1998, ISBN: 0-8493-3988-X, LC: QA164.K73. Input: int N, the number of nodes in the tree. N must be at least 3. Output: int T[2*(N-1)], describes the edges of the tree as pairs of nodes. */ { int *p; int rank; int *t; int tree_num; /* Check. */ if ( n < 1 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "tree_random(): Fatal error!\n" ); fprintf ( stderr, " Input N is illegal.\n" ); exit ( 1 ); } /* Compute the number of trees. */ tree_num = tree_enum ( n ); /* Choose RANK betweeen 1 and TREE_NUM. */ rank = 1 + ( rand ( ) % tree_num ); /* Compute the Pruefer code P of the given RANK. */ p = pruefer_unrank ( rank, n ); /* Convert the Pruefer code P to a tree T. */ t = pruefer_to_tree_new ( n, p ); free ( p ); return t; } /******************************************************************************/ int tree_rank ( int n, int t[] ) /******************************************************************************/ /* Purpose: tree_rank() ranks a tree. Licensing: This code is distributed under the MIT license. Modified: 25 July 2011 Author: John Burkardt Reference: Donald Kreher, Douglas Simpson, Combinatorial Algorithms, CRC Press, 1998, ISBN: 0-8493-3988-X, LC: QA164.K73. Input: int N, the number of nodes in the tree. N must be at least 3. int T[2*(N-1)], describes the edges of the tree as pairs of nodes. Output: int RANK, the rank of the tree. */ { int check; int *p; int rank; /* Check the tree. */ check = tree_check ( n, t ); if ( ! check ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "tree_rank(): Fatal error!\n" ); fprintf ( stderr, " The tree is illegal.\n" ); exit ( 1 ); } /* Convert the tree to a Pruefer code. */ p = tree_to_pruefer ( n, t ); /* Find the rank of the Pruefer code. */ rank = pruefer_rank ( n, p ); free ( p ); return rank; } /******************************************************************************/ void tree_successor ( int n, int t[], int *rank ) /******************************************************************************/ /* Purpose: tree_successor() returns the successor of a tree. Licensing: This code is distributed under the MIT license. Modified: 26 July 2011 Author: John Burkardt Reference: Donald Kreher, Douglas Simpson, Combinatorial Algorithms, CRC Press, 1998, ISBN: 0-8493-3988-X, LC: QA164.K73. Input: int N, the number of nodes in the tree. N must be at least 3. int T[2*(N-1)], describes the edges of the tree as pairs of nodes. int *RANK, the rank of the tree. Output: int T[2*(N-1)], describes the edges of the successor tree as pairs of nodes. int *RANK, the rank of the successor tree. */ { int check; int i; int *p; /* Return the first element. */ if ( *rank == -1 ) { p = ( int * ) malloc ( ( n - 2 ) * sizeof ( int ) ); for ( i = 0; i < n - 2; i++ ) { p[i] = 1; } pruefer_to_tree ( n, p, t ); *rank = 0; free ( p ); return; } /* Check the tree. */ check = tree_check ( n, t ); if ( ! check ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "tree_successor(): Fatal error!\n" ); fprintf ( stderr, " The tree is illegal.\n" ); exit ( 1 ); } /* Convert the tree to a Pruefer code. */ p = tree_to_pruefer ( n, t ); /* Find the successor of the Pruefer code. */ pruefer_successor ( n, p, rank ); /* Convert the Pruefer code to the tree. */ pruefer_to_tree ( n, p, t ); free ( p ); return; } /******************************************************************************/ int *tree_to_pruefer ( int n, int t[] ) /******************************************************************************/ /* Purpose: tree_to_pruefer() converts a tree to a Pruefer code. Licensing: This code is distributed under the MIT license. Modified: 25 July 2011 Author: John Burkardt Reference: Donald Kreher, Douglas Simpson, Combinatorial Algorithms, CRC Press, 1998, ISBN: 0-8493-3988-X, LC: QA164.K73. Input: int N, the number of nodes in the tree. N must be positive. int T[2*(N-1)], describes the edges of the tree as pairs of nodes. Output: int TREE_TO_PRUEFER[N-2], the Pruefer code for the tree. */ { int check; int *d; int i; int j; int k; int *p; int x; int y; /* Check. */ check = tree_check ( n, t ); if ( ! check ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "tree_to_pruefer(): Fatal error!\n" ); fprintf ( stderr, " The tree is illegal.\n" ); exit ( 1 ); } /* Compute the degree of each node. */ d = edge_degree ( n, n - 1, t ); p = ( int * ) malloc ( ( n - 2 ) * sizeof ( int ) ); for ( j = 1; j <= n - 2; j++ ) { /* Find a node of degree 1. */ x = n; while ( d[x-1] != 1 ) { x = x - 1; } /* Find its neighbor. */ k = 1; for ( ; ; ) { if ( t[0+(k-1)*2] == x ) { y = t[1+(k-1)*2]; break; } if ( t[1+(k-1)*2] == x ) { y = t[0+(k-1)*2]; break; } k = k + 1; } /* Store the neighbor. */ p[j-1] = y; /* Delete the edge from the tree. */ d[x-1] = d[x-1] - 1; d[y-1] = d[y-1] - 1; t[0+(k-1)*2] = - t[0+(k-1)*2]; t[1+(k-1)*2] = - t[1+(k-1)*2]; } /* Remove the negative signs from the first N-2 columns of the tree. */ for ( j = 0; j < n - 2; j++ ) { for ( i = 0; i < 2; i++ ) { t[i+j*2] = - t[i+j*2]; } } free ( d ); return p; } /******************************************************************************/ int *tree_unrank ( int rank, int n ) /******************************************************************************/ /* Purpose: tree_unrank() unranks a tree. Licensing: This code is distributed under the MIT license. Modified: 27 July 2011 Author: John Burkardt Reference: Donald Kreher, Douglas Simpson, Combinatorial Algorithms, CRC Press, 1998, ISBN: 0-8493-3988-X, LC: QA164.K73. Input: int RANK, the rank of the tree. int N, the number of nodes in the tree. N must be at least 3. Output: int T[2*(N-1)], describes the edges of the tree as pairs of nodes. */ { int *p; int *t; int tree_num; /* Check. */ if ( n < 1 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "tree_unrank(): Fatal error!\n" ); fprintf ( stderr, " Input N is illegal.\n" ); exit ( 1 ); } tree_num = tree_enum ( n ); if ( rank < 0 || tree_num < rank ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "tree_unrank(): Fatal error!\n" ); fprintf ( stderr, " The input rank is illegal.\n" ); exit ( 1 ); } /* Unrank the Pruefer code. */ p = pruefer_unrank ( rank, n ); /* Convert the Pruefer code to a tree. */ t = pruefer_to_tree_new ( n, p ); free ( p ); return t; }