# include # include # include # include # include # include # include "ubvec.h" /******************************************************************************/ 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 GNU LGPL 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_max ( int i1, int i2 ) /******************************************************************************/ /* Purpose: I4_MAX returns the maximum of two I4's. Licensing: This code is distributed under the GNU LGPL license. Modified: 29 August 2006 Author: John Burkardt Parameters: Input, int I1, I2, are two integers to be compared. Output, int I4_MAX, the larger of I1 and I2. */ { int value; if ( i2 < i1 ) { value = i1; } else { value = i2; } return value; } /******************************************************************************/ int i4_min ( int i1, int i2 ) /******************************************************************************/ /* Purpose: I4_MIN returns the smaller of two I4's. Licensing: This code is distributed under the GNU LGPL license. Modified: 29 August 2006 Author: John Burkardt Parameters: Input, int I1, I2, two integers to be compared. Output, int I4_MIN, the smaller of I1 and I2. */ { int value; if ( i1 < i2 ) { value = i1; } else { value = i2; } return value; } /******************************************************************************/ int i4_power ( int i, int j ) /******************************************************************************/ /* Purpose: I4_POWER returns the value of I^J. Licensing: This code is distributed under the GNU LGPL license. Modified: 23 October 2007 Author: John Burkardt Parameters: Input, int I, J, the base and the power. J should be nonnegative. Output, int I4_POWER, the value of I^J. */ { int k; int value; if ( j < 0 ) { if ( i == 1 ) { value = 1; } else if ( i == 0 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "I4_POWER - Fatal error!\n" ); fprintf ( stderr, " I^J requested, with I = 0 and J negative.\n" ); exit ( 1 ); } else { value = 0; } } else if ( j == 0 ) { if ( i == 0 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "I4_POWER - Fatal error!\n" ); fprintf ( stderr, " I^J requested, with I = 0 and J = 0.\n" ); exit ( 1 ); } else { value = 1; } } else if ( j == 1 ) { value = i; } else { value = 1; for ( k = 1; k <= j; k++ ) { value = value * i; } } return value; } /******************************************************************************/ int i4_uniform_ab ( int a, int b, int *seed ) /******************************************************************************/ /* Purpose: I4_UNIFORM_AB returns a scaled pseudorandom I4 between A and B. Discussion: The pseudorandom number should be uniformly distributed between A and B. Licensing: This code is distributed under the GNU LGPL license. Modified: 24 May 2012 Author: John Burkardt Reference: Paul Bratley, Bennett Fox, Linus Schrage, A Guide to Simulation, Second Edition, Springer, 1987, ISBN: 0387964673, LC: QA76.9.C65.B73. Bennett Fox, Algorithm 647: Implementation and Relative Efficiency of Quasirandom Sequence Generators, ACM Transactions on Mathematical Software, Volume 12, Number 4, December 1986, pages 362-376. Pierre L'Ecuyer, Random Number Generation, in Handbook of Simulation, edited by Jerry Banks, Wiley, 1998, ISBN: 0471134031, LC: T57.62.H37. Peter Lewis, Allen Goodman, James Miller, A Pseudo-Random Number Generator for the System/360, IBM Systems Journal, Volume 8, Number 2, 1969, pages 136-143. Parameters: Input, int A, B, the limits of the interval. Input/output, int *SEED, the "seed" value, which should NOT be 0. On output, SEED has been updated. Output, int I4_UNIFORM_AB, a number between A and B. */ { int c; const int i4_huge = 2147483647; int k; float r; int value; if ( *seed == 0 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "I4_UNIFORM_AB - Fatal error!\n" ); fprintf ( stderr, " Input value of SEED = 0.\n" ); exit ( 1 ); } /* Guaranteee A <= B. */ if ( b < a ) { c = a; a = b; b = c; } k = *seed / 127773; *seed = 16807 * ( *seed - k * 127773 ) - k * 2836; if ( *seed < 0 ) { *seed = *seed + i4_huge; } r = ( float ) ( *seed ) * 4.656612875E-10; /* Scale R to lie between A-0.5 and B+0.5. */ r = ( 1.0 - r ) * ( ( float ) ( a ) - 0.5 ) + r * ( ( float ) ( b ) + 0.5 ); /* Round R to the nearest integer. */ value = round ( r ); /* Guarantee that A <= VALUE <= B. */ if ( value < a ) { value = a; } if ( b < value ) { value = b; } return value; } /******************************************************************************/ void i4vec_print ( int n, int a[], char *title ) /******************************************************************************/ /* Purpose: I4VEC_PRINT prints an I4VEC. Discussion: An I4VEC is a vector of I4's. Licensing: This code is distributed under the GNU LGPL license. Modified: 14 November 2003 Author: John Burkardt Parameters: Input, int N, the number of components of the vector. Input, int A[N], the vector to be printed. Input, char *TITLE, a title. */ { int i; fprintf ( stdout, "\n" ); fprintf ( stdout, "%s\n", title ); fprintf ( stdout, "\n" ); for ( i = 0; i < n; i++ ) { fprintf ( stdout, " %6d: %8d\n", i, a[i] ); } return; } /******************************************************************************/ 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 GNU LGPL 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 GNU LGPL license. Modified: 02 June 2015 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 ); if ( 0 < title_len ) { printf ( "\n" ); printf ( "%s\n", title ); } if ( 0 < n ) { for ( ilo = 1; ilo <= n; ilo = ilo + 5 ) { ihi = i4_min ( ilo + 5 - 1, n ); for ( i = ilo; i <= ihi; i++ ) { printf ( "%12d", a[i-1] ); } printf ( "\n" ); } } else { printf ( " (empty vector)\n" ); } 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 GNU LGPL 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 GNU LGPL 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 morse_thue ( unsigned int ui ) /******************************************************************************/ /* Purpose: MORSE_THUE generates a Morse_Thue number. Discussion: The Morse_Thue sequence can be defined in a number of ways. A) Start with the string containing the single letter '0'; then repeatedly apply the replacement rules '0' -> '01' and '1' -> '10' to the letters of the string. The Morse_Thue sequence is the resulting letter sequence. B) Starting with the string containing the single letter '0', repeatedly append the binary complement of the string to itself. Thus, '0' becomes '0' + '1' or '01', then '01' becomes '01' + '10', which becomes '0110' + '1001', and so on. C) Starting with I = 0, the I-th Morse-Thue number is determined by taking the binary representation of I, adding the digits, and computing the remainder modulo 2. Example: I binary S -- ------ -- 0 0 0 1 1 1 2 10 1 3 11 0 4 100 1 5 101 0 6 110 0 7 111 1 8 1000 1 Licensing: This code is distributed under the GNU LGPL license. Modified: 30 December 2014 Author: John Burkardt Parameters: Input, unsigned int UI, the index of the Morse-Thue number. Output, int MORSE_THUE, the Morse-Thue number of index I. */ { int i; const int nbits = 32; int s; int *ubvec; /* Expand UI into binary form. */ ubvec = ui4_to_ubvec ( ui, nbits ); /* Sum the 1's in the binary representation. */ s = 0; for ( i = 0; i < nbits; i++ ) { s = s + ubvec[i]; } /* Take the value modulo 2. */ s = s % 2; free ( ubvec ); return s; } /******************************************************************************/ unsigned int nim_sum ( unsigned int ui, unsigned int uj ) /******************************************************************************/ /* Purpose: NIM_SUM computes the Nim sum of two integers. Discussion: If K is the Nim sum of I and J, then each bit of K is the exclusive OR of the corresponding bits of I and J. Example: I J K I_2 J_2 K_2 ---- ---- ---- ---------- ---------- ---------- 0 0 0 0 0 0 1 0 1 1 0 1 1 1 0 1 1 0 2 7 5 10 111 101 11 28 23 1011 11100 10111 Licensing: This code is distributed under the GNU LGPL license. Modified: 01 December 2015 Author: John Burkardt Parameters: Input, unsigned int UI, UJ, the integers to be Nim-summed. Output, unsigned int NIM_SUM, the Nim sum of I and J. */ { const int nbits = 32; int *ubvec1; int *ubvec2; int *ubvec3; unsigned int value; ubvec1 = ui4_to_ubvec ( ui, nbits ); ubvec2 = ui4_to_ubvec ( uj, nbits ); ubvec3 = ubvec_xor ( nbits, ubvec1, ubvec2 ); value = ubvec_to_ui4 ( nbits, ubvec3 ); free ( ubvec1 ); free ( ubvec2 ); free ( ubvec3 ); return value; } /******************************************************************************/ void timestamp ( ) /******************************************************************************/ /* Purpose: TIMESTAMP prints the current YMDHMS date as a time stamp. Example: 17 June 2014 09:45:54 AM Licensing: This code is distributed under the GNU LGPL license. Modified: 17 June 2014 Author: John Burkardt Parameters: None */ { # define TIME_SIZE 40 static char time_buffer[TIME_SIZE]; const struct tm *tm; time_t now; now = time ( NULL ); tm = localtime ( &now ); strftime ( time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm ); printf ( "%s\n", time_buffer ); return; # undef TIME_SIZE } /******************************************************************************/ int *ubvec_add ( int n, int ubvec1[], int ubvec2[] ) /******************************************************************************/ /* Purpose: UBVEC_ADD adds two unsigned binary vectors. Discussion: A UBVEC is a vector of binary digits representing an unsigned integer. UBVEC[N-1] contains the units digit, UBVEC[N-2] the coefficient of 2, UBVEC[N-3] the coefficient of 4 and so on, so that printing the digits in order gives the binary form of the number. Example: N = 4 UBVEC1 + UBVEC2 = UBVEC3 ( 0 0 0 1 ) + ( 0 0 1 1 ) = ( 0 1 0 0 ) 1 + 3 = 4 Licensing: This code is distributed under the GNU LGPL license. Modified: 21 March 2009 Author: John Burkardt Parameters: Input, int N, the length of the vectors. Input, int UBVEC1[N], UBVEC2[N], the vectors to be added. Output, int UBVEC_ADD[N], the sum of the two input vectors. */ { int i; int *ubvec3; ubvec3 = ( int * ) malloc ( n * sizeof ( int ) ); for ( i = 0; i < n; i++ ) { ubvec3[i] = ubvec1[i] + ubvec2[i]; } for ( i = n - 1; 0 <= i; i-- ) { while ( 2 <= ubvec3[i] ) { ubvec3[i] = ubvec3[i] - 2; if ( 0 < i ) { ubvec3[i-1] = ubvec3[i-1] + 1; } else { fprintf ( stderr, "\n" ); fprintf ( stderr, "UBVEC_ADD - Fatal error!\n" ); fprintf ( stderr, " Addition overflow.\n" ); exit ( 1 ); } } } return ubvec3; } /******************************************************************************/ int *ubvec_and ( int n, int ubvec1[], int ubvec2[] ) /******************************************************************************/ /* Purpose: UBVEC_AND computes the AND of two unsigned binary vectors. Discussion: A UBVEC is a vector of N binary digits. A UBVEC can be interpreted as a binary representation of an unsigned integer, with the first entry being the coefficient of 2^(N-1) and the last entry the coefficient of 1. UBVEC # ----- -- 00000 0 00001 1 00010 2 10000 16 Licensing: This code is distributed under the GNU LGPL license. Modified: 01 December 2015 Author: John Burkardt Parameters: Input, int N, the length of the vectors. Input, int UBVEC1[N], UBVEC2[N], the vectors. Input, int UBVEC_AND[N], the AND of the two vectors. */ { int i; int *ubvec3; ubvec3 = ( int * ) malloc ( n * sizeof ( int ) ); for ( i = 0; i < n; i++ ) { ubvec3[i] = i4_min ( ubvec1[i], ubvec2[i] ); } return ubvec3; } /******************************************************************************/ bool ubvec_check ( int n, int ubvec[] ) /******************************************************************************/ /* Purpose: UBVEC_CHECK checks an unsigned binary vector. Discussion: The only check made is that the entries are all 0 or 1. A UBVEC is a vector of N binary digits. A UBVEC can be interpreted as a binary representation of an unsigned integer, with the first entry being the coefficient of 2^(N-1) and the last entry the coefficient of 1. UBVEC # ----- -- 00000 0 00001 1 00010 2 10000 16 Licensing: This code is distributed under the GNU LGPL license. Modified: 01 December 2015 Author: John Burkardt Parameters: Input, int N, the length of the vectors. Input, int UBVEC[N], the vector to be checked. Output, bool UBVEC_CHECK. TRUE, the data is legal. FALSE, the data is illegal. */ { bool check; int i; check = true; for ( i = 0; i < n; i++ ) { if ( ubvec[i] < 0 || 1 < ubvec[i] ) { check = false; break; } } return check; } /******************************************************************************/ int *ubvec_complement1 ( int n, int ubvec1[] ) /******************************************************************************/ /* Purpose: UBVEC_COMPLEMENT1 computes the one's complement of an unsigned binary vector. Discussion: A UBVEC is a vector of N binary digits. A UBVEC can be interpreted as a binary representation of an unsigned integer, with the first entry being the coefficient of 2^(N-1) and the last entry the coefficient of 1. UBVEC # ----- -- 00000 0 00001 1 00010 2 10000 16 Licensing: This code is distributed under the GNU LGPL license. Modified: 01 December 2015 Author: John Burkardt Parameters: Input, int N, the length of the vectors. Input, int UBVEC1[N], the vector to be complemented. Output, int UBVEC_COMPLEMENT1[N], the complemented vector. */ { int i; int *ubvec2; ubvec2 = ( int * ) malloc ( n * sizeof ( int ) ); for ( i = 0; i < n; i++ ) { ubvec2[i] = 1 - ubvec1[i]; } return ubvec2; } /******************************************************************************/ int ubvec_enum ( int n ) /******************************************************************************/ /* Purpose: UBVEC_ENUM enumerates the unsigned binary vectors of length N. Discussion: A UBVEC is a vector of N binary digits. A UBVEC can be interpreted as a binary representation of an unsigned integer, with the first entry being the coefficient of 2^(N-1) and the last entry the coefficient of 1. UBVEC # ----- -- 00000 0 00001 1 00010 2 10000 16 Licensing: This code is distributed under the GNU LGPL license. Modified: 01 December 2015 Author: John Burkardt Parameters: Input, int N, the length of the vectors. Output, int UBVEC_ENUM, the number of binary vectors. */ { int value; value = i4_power ( 2, n ); return value; } /******************************************************************************/ void ubvec_next ( int n, int ubvec[] ) /******************************************************************************/ /* Purpose: UBVEC_NEXT generates the next UBVEC. Discussion: The vectors are produced in the order: (0,0,...,0), (0,0,...,1), ... (1,1,...,1) and the "next" vector after (1,1,...,1) is (0,0,...,0). That is, we allow wrap around. A UBVEC is a vector of N binary digits. A UBVEC can be interpreted as a binary representation of an unsigned integer, with the first entry being the coefficient of 2^(N-1) and the last entry the coefficient of 1. UBVEC # ----- -- 00000 0 00001 1 00010 2 10000 16 Example: N = 3 Input Output ----- ------ 0 0 0 => 0 0 1 0 0 1 => 0 1 0 0 1 0 => 0 1 1 0 1 1 => 1 0 0 1 0 0 => 1 0 1 1 0 1 => 1 1 0 1 1 0 => 1 1 1 1 1 1 => 0 0 0 Licensing: This code is distributed under the GNU LGPL license. Modified: 01 December 2015 Author: John Burkardt Parameters: Input, int N, the dimension of the vectors. Input/output, int UBVEC[N], on output, the successor to the input vector. */ { int i; for ( i = n - 1; 0 <= i; i-- ) { if ( ubvec[i] == 0 ) { ubvec[i] = 1; return; } ubvec[i] = 0; } return; } /******************************************************************************/ void ubvec_next_gray ( int n, int t[] ) /******************************************************************************/ /* Purpose: UBVEC_NEXT_GRAY computes the next UBVEC in the Gray code. Licensing: This code is distributed under the GNU LGPL license. Modified: 01 December 201 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. */ { int i; int weight; weight = i4vec_sum ( n, t ); if ( ( weight % 2 ) == 0 ) { if ( t[n-1] == 0 ) { t[n-1] = 1; } else { t[n-1] = 0; } 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; } return; } } /* The final element was input. Return the first element. */ for ( i = 0; i < n; i++ ) { t[i] = 0; } } return; } /******************************************************************************/ void ubvec_next_grlex ( int n, int ubvec[] ) /******************************************************************************/ /* Purpose: UBVEC_NEXT_GRLEX generates the next UBVEC in GRLEX order. Discussion: N = 3 Input Output ----- ------ 0 0 0 => 0 0 1 0 0 1 => 0 1 0 0 1 0 => 1 0 0 1 0 0 => 0 1 1 0 1 1 => 1 0 1 1 0 1 => 1 1 0 1 1 0 => 1 1 1 1 1 1 => 0 0 0 A UBVEC is a vector of N binary digits. A UBVEC can be interpreted as a binary representation of an unsigned integer, with the first entry being the coefficient of 2^(N-1) and the last entry the coefficient of 1. UBVEC # ----- -- 00000 0 00001 1 00010 2 10000 16 Licensing: This code is distributed under the GNU LGPL license. Modified: 01 December 2015 Author: John Burkardt Parameters: Input, int N, the dimension. Input/output, int UBVEC[N], the binary vector whose successor is desired, and on output, the successor. */ { int i; int o; int s; int z; /* Initialize locations of 0 and 1. */ if ( ubvec[0] == 0 ) { z = 0; o = -1; } else { z = -1; o = 0; } /* Moving from right to left, search for a "1", preceded by a "0". */ for ( i = n - 1; 1 <= i; i-- ) { if ( ubvec[i] == 1 ) { o = i; if ( ubvec[i-1] == 0 ) { z = i - 1; break; } } } /* UBVEC = 0 */ if ( o == -1 ) { ubvec[n-1] = 1; } /* 01 never occurs. So for sure, B(1) = 1. */ else if ( z == -1 ) { s = i4vec_sum ( n, ubvec ); if ( s == n ) { for ( i = 0; i < n; i++ ) { ubvec[i] = 0; } } else { for ( i = 0; i < n - s - 1; i++ ) { ubvec[i] = 0; } for ( i = n - s - 1; i < n; i++ ) { ubvec[i] = 1; } } } /* Found the rightmost "01" string. Replace it by "10". Shift following 1's to the right. */ else { ubvec[z] = 1; ubvec[o] = 0; s = 0; for ( i = o + 1; i < n; i++ ) { s = s + ubvec[i]; } for ( i = o + 1; i < n - s; i++ ) { ubvec[i] = 0; } for ( i = n - s; i < n; i++ ) { ubvec[i] = 1; } } return; } /******************************************************************************/ int *ubvec_or ( int n, int ubvec1[], int ubvec2[] ) /******************************************************************************/ /* Purpose: UBVEC_OR computes the OR of two unsigned binary vectors. Discussion: A UBVEC is a vector of N binary digits. A UBVEC can be interpreted as a binary representation of an unsigned integer, with the first entry being the coefficient of 2^(N-1) and the last entry the coefficient of 1. UBVEC # ----- -- 00000 0 00001 1 00010 2 10000 16 Licensing: This code is distributed under the GNU LGPL license. Modified: 01 December 2015 Author: John Burkardt Parameters: Input, int N, the length of the vectors. Input, int UBVEC1[N], UBVEC2[N], the vectors. Input, int UBVEC_OR[N], the OR of the two vectors. */ { int i; int *ubvec3; ubvec3 = ( int * ) malloc ( n * sizeof ( int ) ); for ( i = 0; i < n; i++ ) { ubvec3[i] = i4_max ( ubvec1[i], ubvec2[i] ); } return ubvec3; } /******************************************************************************/ void ubvec_print ( int n, int bvec[], char *title ) /******************************************************************************/ /* Purpose: UBVEC_PRINT prints a UBVEC, with an optional title. Discussion: A UBVEC is a vector of binary digits representing an unsigned integer. UBVEC[N-1] contains the units digit, UBVEC[N-2] the coefficient of 2, UBVEC[N-3] the coefficient of 4 and so on, so that printing the digits in order gives the binary form of the number. Licensing: This code is distributed under the GNU LGPL license. Modified: 01 December 2006 Author: John Burkardt Parameters: Input, int N, the number of components of the vector. Input, int BVEC[N], the vector to be printed. Input, char *TITLE, a title to be printed first. TITLE may be blank. */ { int i; int ihi; int ilo; if ( 0 < strlen ( title ) ) { printf ( "\n" ); printf ( "%s\n", title ); } for ( ilo = 0; ilo < n; ilo = ilo + 70 ) { ihi = i4_min ( ilo + 70 - 1, n - 1 ); printf ( " " ); for ( i = ilo; i <= ihi; i++ ) { printf ( "%d", bvec[i] ); } printf ( "\n" ); } return; } /******************************************************************************/ int *ubvec_random ( int n, int *seed ) /******************************************************************************/ /* Purpose: UBVEC_RANDOM returns a pseudorandom UBVEC. Discussion: A UBVEC is a vector of N binary digits. A UBVEC can be interpreted as a binary representation of an unsigned integer, with the first entry being the coefficient of 2^(N-1) and the last entry the coefficient of 1. UBVEC # ----- -- 00000 0 00001 1 00010 2 10000 16 Licensing: This code is distributed under the GNU LGPL license. Modified: 01 December 2015 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 N, the order of the vector. Input/output, int *SEED, the "seed" value, which should NOT be 0. On output, SEED has been updated. Output, int UBVEC_RANDOM[N], a pseudorandom binary vector. */ { const int i4_huge = 2147483647; const int i4_huge_half = 1073741823; int i; int k; int *ubvec; if ( seed == 0 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "UBVEC_RANDOM - Fatal error!\n" ); fprintf ( stderr, " Input value of SEED = 0.\n" ); exit ( 1 ); } ubvec = ( 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; } if ( i4_huge_half < *seed ) { ubvec[i] = 0; } else { ubvec[i] = 1; } } return ubvec; } /******************************************************************************/ int ubvec_rank_gray ( int n, int ubvec[] ) /******************************************************************************/ /* Purpose: UBVEC_RANK_GRAY ranks a UBVEC according to the Gray ordering. Licensing: This code is distributed under the GNU LGPL license. Modified: 01 December 2015 Author: John Burkardt Parameters: Input, int N, the number of components of the vector. Input, int UBVEC[N], the vector to be printed. Output, int UBVEC_RANK, the rank of the UBVEC. */ { int rank; int ui4; ui4 = ubvec_to_ui4 ( n, ubvec ); rank = ui4_rank_gray ( ui4 ); return rank; } /******************************************************************************/ int *ubvec_reverse ( int n, int ubvec1[] ) /******************************************************************************/ /* Purpose: UBVEC_REVERSE reverses a UBVEC. Discussion: A UBVEC is a vector of N binary digits. A UBVEC can be interpreted as a binary representation of an unsigned integer, with the first entry being the coefficient of 2^(N-1) and the last entry the coefficient of 1. UBVEC # ----- -- 00000 0 00001 1 00010 2 10000 16 Licensing: This code is distributed under the GNU LGPL license. Modified: 19 September 2015 Author: John Burkardt Parameters: Input, int N, the length of the vectors. Input, int UBVEC1[N], the vector to be reversed. Output, int UBVEC_REVERSE[N], the reversed vector. */ { int i; int *ubvec2; ubvec2 = ( int * ) malloc ( n * sizeof ( int ) ); for ( i = 0; i < n; i++ ) { ubvec2[i] = ubvec1[n-1-i]; } return ubvec2; } /******************************************************************************/ int *ubvec_unrank_gray ( int rank, int n ) /******************************************************************************/ /* Purpose: UBVEC_UNRANK_GRAY unranks a UBVEC. Licensing: This code is distributed under the GNU LGPL license. Modified: 01 December 2015 Author: John Burkardt Parameters: Input, int RANK, the rank of the UBVEC. 0 <= RANK < 2^N. Input, int N, the size of the UBVEC. Output, int UBVEC_UNRANK_GRAY[N], the UBVEC of given rank. */ { int *ubvec; int ui4; ui4 = ui4_unrank_gray ( rank ); ubvec = ui4_to_ubvec ( ui4, n ); return ubvec; } /******************************************************************************/ int *ubvec_unrank_grlex ( int rank, int n ) /******************************************************************************/ /* Purpose: UBVEC_UNRANK_GRLEX unranks a UBVEC using the GRLEX ordering. Discussion: Licensing: This code is distributed under the GNU LGPL license. Modified: 01 December 2015 Author: John Burkardt Parameters: Input, int RANK, the rank. 0 <= RANK < 2^N. Input, int N, the size of the UBVEC. Output, int UBVEC_UNRANK_GRLEX[N], the UBVEC of the given rank. */ { int *b; int i; int k; int mk; int mk_old; int mk_plus; int rank_k; int *t; b = ( int * ) malloc ( n * sizeof ( n ) ); mk = 0; for ( k = 0; k <= n; k++ ) { mk_old = mk; mk_plus = i4_choose ( n, k ); mk = mk_old + mk_plus; if ( rank < mk ) { rank_k = rank - mk_old; t = ksubset_colex_unrank ( rank_k, k, n ); for ( i = 0; i < n; i++ ) { b[i] = 0; } for ( i = 0; i < k; i++ ) { b[n-t[i]] = 1; } free ( t ); return b; } } /* If we got here, the rank is too large. */ fprintf ( stderr, "\n" ); fprintf ( stderr, "UBVEC_UNRANK_GRLEX - Fatal error!\n" ); fprintf ( stderr, " Input value of rank is too high.\n" ); exit ( 1 ); } /******************************************************************************/ unsigned int ubvec_to_ui4 ( int n, int bvec[] ) /******************************************************************************/ /* Purpose: UBVEC_TO_UI4 makes an unsigned integer from an unsigned binary vector. Discussion: A UBVEC is a vector of binary digits representing an unsigned integer. UBVEC[N-1] contains the units digit, UBVEC[N-2] the coefficient of 2, UBVEC[N-3] the coefficient of 4 and so on, so that printing the digits in order gives the binary form of the number. Example: N = 4 BVEC binary I ---------- ----- -- 1 2 3 4 ---------- 0 0 0 1 1 1 0 0 1 0 10 2 0 0 1 1 11 3 0 1 0 0 100 4 1 0 0 1 1001 9 1 1 1 1 1111 15 Licensing: This code is distributed under the GNU LGPL license. Modified: 21 March 2009 Author: John Burkardt Parameters: Input, int N, the dimension of the vector. Input, int BVEC[N], the binary representation. Input, unsigned int UBVEC_TO_UI4, the integer value. */ { const unsigned int base = 2; int i; unsigned int ui4; ui4 = 0; for ( i = 0; i < n; i++ ) { ui4 = base * ui4 + ( unsigned int ) bvec[i]; } return ui4; } /******************************************************************************/ int *ubvec_xor ( int n, int ubvec1[], int ubvec2[] ) /******************************************************************************/ /* Purpose: UBVEC_XOR computes the exclusive OR of two unsigned binary vectors. Discussion: A UBVEC is a vector of N binary digits. A UBVEC can be interpreted as a binary representation of an unsigned integer, with the first entry being the coefficient of 2^(N-1) and the last entry the coefficient of 1. UBVEC # ----- -- 00000 0 00001 1 00010 2 10000 16 Licensing: This code is distributed under the GNU LGPL license. Modified: 01 December 2015 Author: John Burkardt Parameters: Input, int N, the length of the vectors. Input, int UBVEC1[N], UBVEC2[N], the vectors to be XOR'ed. Input, int UBVEC3[N], the exclusive OR of the two vectors. */ { int i; int *ubvec3; ubvec3 = ( int * ) malloc ( n * sizeof ( int ) ); for ( i = 0; i < n; i++ ) { ubvec3[i] = ( ( ubvec1[i] + ubvec2[i] ) % 2 ); } return ubvec3; } /******************************************************************************/ int ui4_rank_gray ( unsigned int ui4 ) /******************************************************************************/ /* Purpose: UI4_RANK_GRAY ranks a Gray code. Discussion: This routine is entirely arithmetical, and does not require access to bit testing and setting routines. Given the number GRAY, its ranking is the order in which it would be visited in the Gray code ordering. The Gray code ordering begins Rank Gray Gray (Dec) (Bin) 0 0 0000 1 1 0001 2 3 0011 3 2 0010 4 6 0110 5 7 0111 6 5 0101 7 4 0100 8 12 0110 etc This routine is given a Gray code, and has to return the rank. Example: Gray Gray Rank (Dec) (Bin) 0 0 0 1 1 1 2 10 3 3 11 2 4 100 7 5 101 6 6 110 4 7 111 5 8 1000 15 9 1001 14 10 1010 12 11 1011 13 12 1100 8 13 1101 9 14 1110 11 15 1111 10 16 10000 31 Licensing: This code is distributed under the GNU LGPL license. Modified: 01 December 2015 Author: John Burkardt Parameters: Input, unsigned int UI4, the Gray code to be ranked. Output, int UI4_RANK_GRAY, the rank of UI4. */ { int k; bool last; bool next; int rank; int two_k; unsigned int ui4_copy; ui4_copy = ui4; if ( ui4_copy == 0 ) { rank = 0; return rank; } /* Find TWO_K, the largest power of 2 less than or equal to GRAY. */ k = 0; two_k = 1; while ( 2 * two_k <= ui4_copy ) { two_k = two_k * 2; k = k + 1; } rank = two_k; last = true; ui4_copy = ui4_copy - two_k; while ( 0 < k ) { two_k = two_k / 2; k = k - 1; next = ( two_k <= ui4_copy && ui4_copy < two_k * 2 ); if ( next ) { ui4_copy = ui4_copy - two_k; } if ( next != last ) { rank = rank + two_k; last = true; } else { last = false; } } return rank; } /******************************************************************************/ int *ui4_to_ubvec ( unsigned int ui4, int n ) /******************************************************************************/ /* Purpose: UI4_TO_UBVEC makes a unsigned binary vector from an integer. Discussion: A UBVEC is a vector of binary digits representing an unsigned integer. UBVEC[N-1] contains the units digit, UBVEC[N-2] the coefficient of 2, UBVEC[N-3] the coefficient of 4 and so on, so that printing the digits in order gives the binary form of the number. To guarantee that there will be enough space for any value of I, it would be necessary to set N = 32. Example: UI4 BVEC binary 0 1 2 3 4 5 1 2 4 8 16 32 -- ---------------- ------ 1 1 0 0 0 0 1 1 2 0 1 0 0 1 0 10 3 1 1 0 0 1 1 11 4 0 0 0 1 0 0 100 9 0 0 1 0 0 1 1001 57 1 1 1 0 1 1 110111 Licensing: This code is distributed under the GNU LGPL license. Modified: 01 December 2015 Author: John Burkardt Parameters: Input, unsigned int UI4, an integer to be represented. Input, int N, the dimension of the vector. Output, int UI4_TO_UBVEC[N], the binary representation. */ { int i; int *ubvec; ubvec = ( int * ) malloc ( n * sizeof ( int ) ); for ( i = n - 1; 0 <= i; i-- ) { ubvec[i] = ui4 % 2; ui4 = ui4 / 2; } return ubvec; } /******************************************************************************/ int ui4_unrank_gray ( int rank ) /******************************************************************************/ /* Purpose: UI4_UNRANK_GRAY unranks a Gray code. Discussion: This routine is entirely arithmetical, and does not require access to bit testing and setting routines. The binary values of the Gray codes of successive integers differ in just one bit. The sequence of Gray codes for 0 to (2^N)-1 can be interpreted as a Hamiltonian cycle on a graph of the cube in N dimensions. Example: Rank Gray Gray (Dec) (Bin) 0 0 0 1 1 1 2 3 11 3 2 10 4 6 110 5 7 111 6 5 101 7 4 100 8 12 1100 9 14 1001 10 12 1010 11 13 1011 12 8 1100 13 9 1101 14 11 1110 15 10 1111 16 31 10000 Licensing: This code is distributed under the GNU LGPL license. Modified: 01 December 2015 Author: John Burkardt Parameters: Input, int RANK, the integer whose Gray code is desired. Output, int UI4_UNRANK_GRAY, the Gray code of the given rank. */ { int k; bool last; bool next; int rank_copy; int two_k; unsigned int ui4; if ( rank <= 0 ) { ui4 = 0; return ui4; } rank_copy = rank; k = 0; two_k = 1; while ( 2 * two_k <= rank_copy ) { two_k = two_k * 2; k = k + 1; } ui4 = two_k; rank_copy = rank_copy - two_k; next = true; while ( 0 < k ) { two_k = two_k / 2; k = k - 1; last = next; next = ( two_k <= rank_copy && rank_copy <= two_k * 2 ); if ( next != last ) { ui4 = ui4 + two_k; } if ( next ) { rank_copy = rank_copy - two_k; } } return ui4; }