# include # include # include # include # include using namespace std; # include "toms515.hpp" //****************************************************************************80 int binom ( int n, int k ) //****************************************************************************80 // // Purpose: // // BINOM computes the binomial coefficient. // // Discussion: // // This is ACM algorithm 160 translated to Fortran. // // It calculates the number of combinations of N things taken K at a time. // // Modified: // // 01 April 2016 // // Author: // // Bill Buckles, Matthew Lybanon // // Reference: // // Bill Buckles, Matthew Lybanon, // Algorithm 515: Generation of a Vector from the Lexicographical Index, // ACM Transactions on Mathematical Software, // Volume 3, Number 2, June 1977, pages 180-182. // // Parameters: // // Input, int N, K, the parameters for the binomial // coefficient. // // Output, int BINOM, the binomial coefficient. // { int i; int k1; int p; int r; k1 = k; p = n - k1; if ( k1 < p ) { p = k1; k1 = n - p; } if ( p == 0 ) { r = 1; } else { r = k1 + 1; } for ( i = 2; i <= p; i++ ) { r = ( r * ( k1 + i ) ) / i; } return r; } //****************************************************************************80 int *comb ( int n, int p, int l ) //****************************************************************************80 // // Purpose: // // COMB selects a subset of order P from a set of order N. // // Discussion: // // This subroutine finds the combination set of N things taken // P at a time for a given lexicographic index. // // Modified: // // 01 April 2016 // // Author: // // Bill Buckles, Matthew Lybanon // // Reference: // // Bill Buckles, Matthew Lybanon, // Algorithm 515: Generation of a Vector from the Lexicographical Index, // ACM Transactions on Mathematical Software, // Volume 3, Number 2, June 1977, pages 180-182. // // Parameters: // // Input, int N, the number of things in the set. // // Input, int P, the number of things in each combination. // 0 < P < N. // // Input, int L, the lexicographic index of the // desired combination. 1 <= L <= choose(N,P). // // Output, int COMB[P], the combination set. // { int *c; int i; int k; int p1; int r; c = new int[p]; // // Special case: P = 1 // if ( p == 1 ) { c[0] = l; return c; } // // Initialize lower bound index. // k = 0; // // Select elements in ascending order. // p1 = p - 1; c[0] = 0; for ( i = 1; i <= p1; i++ ) { // // Update lower bound as the previously selected element. // if ( 1 < i ) { c[i-1] = c[i-2]; } // // Check validity of each entry. // for ( ; ; ) { c[i-1] = c[i-1] + 1; r = binom ( n - c[i-1], p - i ); k = k + r; if ( l <= k ) { break; } } k = k - r; } c[p-1] = c[p1-1] + l - k; return c; } //****************************************************************************80 bool i4_choose_check ( int n, int k ) //****************************************************************************80 // // Purpose: // // I4_CHOOSE_CHECK reports whether the binomial coefficient can be computed. // // Licensing: // // This code is distributed under the MIT license. // // Modified: // // 28 March 2016 // // Author: // // John Burkardt // // Parameters: // // Input, int N, K, the binomial parameters. // // Output, bool I4_CHOOSE_CHECK is: // TRUE, if C(N,K) < maximum integer. // FALSE, otherwise. // { bool check; double choose_nk_log; const int i4_huge = 2147483647; double i4_huge_log; i4_huge_log = log ( ( double ) ( i4_huge ) ); choose_nk_log = lgamma ( ( double ) ( n + 1 ) ) - lgamma ( ( double ) ( k + 1 ) ) - lgamma ( ( double ) ( n - k + 1 ) ); check = ( choose_nk_log < i4_huge_log ); return check; } //****************************************************************************80 int i4_uniform_ab ( int a, int b, int &seed ) //****************************************************************************80 // // 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 MIT license. // // Modified: // // 02 October 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, a number between A and B. // { int c; const int i4_huge = 2147483647; int k; float r; int value; if ( seed == 0 ) { cerr << "\n"; cerr << "I4_UNIFORM_AB - Fatal error!\n"; cerr << " Input value of SEED = 0.\n"; exit ( 1 ); } // // Guarantee 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 ); // // 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; } return value; } //****************************************************************************80 void timestamp ( ) //****************************************************************************80 // // Purpose: // // TIMESTAMP prints the current YMDHMS date as a time stamp. // // Example: // // 31 May 2001 09:45:54 AM // // Licensing: // // This code is distributed under the MIT license. // // Modified: // // 08 July 2009 // // Author: // // John Burkardt // // Parameters: // // None // { # define TIME_SIZE 40 static char time_buffer[TIME_SIZE]; const struct std::tm *tm_ptr; std::time_t now; now = std::time ( NULL ); tm_ptr = std::localtime ( &now ); std::strftime ( time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm_ptr ); std::cout << time_buffer << "\n"; return; # undef TIME_SIZE }