# include # include # include # include # include # include "toms515.h" /******************************************************************************/ int binom ( int n, int k ) /******************************************************************************/ /* 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; } /******************************************************************************/ int *comb ( int n, int p, int l ) /******************************************************************************/ /* 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 = ( int * ) malloc ( p * sizeof ( int ) ); /* 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; } /******************************************************************************/ int i4_choose_check ( int n, int k ) /******************************************************************************/ /* 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, int I4_CHOOSE_CHECK is: TRUE, if C(N,K) < maximum integer. FALSE, otherwise. */ { int 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; } /******************************************************************************/ 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 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; 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 timestamp ( ) /******************************************************************************/ /* 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: 24 September 2003 Author: John Burkardt Parameters: None */ { # define TIME_SIZE 40 static char time_buffer[TIME_SIZE]; const struct tm *tm; time_t now; now = time ( NULL ); tm = localtime ( &now ); strftime ( time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm ); fprintf ( stdout, "%s\n", time_buffer ); return; # undef TIME_SIZE }