# include # include # include # include # include # include "uniform.h" /******************************************************************************/ void bvec_print ( int n, int bvec[], char *title ) /******************************************************************************/ /* Purpose: bvec_print() prints a BVEC, with an optional title. Discussion: A BVEC is a vector of binary digits representing an integer. BVEC[0] is 0 for positive values and 1 for negative values, which are stored in 2's complement form. For positive values, BVEC[N-1] contains the units digit, BVEC[N-2] the coefficient of 2, BVEC[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 MIT license. Modified: 01 December 2006 Author: John Burkardt Input: int N, the number of components of the vector. int BVEC[N], the vector to be printed. 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 *bvec_uniform_new ( int n ) /******************************************************************************/ /* Purpose: bvec_uniform_new() returns a random binary vector. Discussion: A BVEC is a vector of binary digits representing an integer. BVEC[0] is 0 for positive values and 1 for negative values, which are stored in 2's complement form. For positive values, BVEC[N-1] contains the units digit, BVEC[N-2] the coefficient of 2, BVEC[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 MIT license. Modified: 25 December 2014 Author: John Burkardt Input: int N, the length of the vectors. Output: int BVEC_UNIFORM_NEW[N], the randomly selected vector. */ { int i; int *bvec; bvec = ( int * ) malloc ( n * sizeof ( int ) ); for ( i = 0; i < n; i++ ) { bvec[i] = ( lrand48 ( ) % 2 ); } return bvec; } /******************************************************************************/ float complex c4_uniform_01 ( ) /******************************************************************************/ /* Purpose: c4_uniform_01() returns a unit pseudorandom C4. Discussion: The angle should be uniformly distributed between 0 and 2 * PI, the square root of the radius uniformly distributed between 0 and 1. This results in a uniform distribution of values in the unit circle. Licensing: This code is distributed under the MIT license. Modified: 11 January 2007 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. Output: float complex C4_UNIFORM_01, a pseudorandom complex value. */ { float r; const float r4_pi = 3.1415926E+00; float theta; float complex value; r = sqrt ( drand48 ( ) ); theta = 2.0 * r4_pi * drand48 ( ); value = r * ( cos ( theta ) + I * sin ( theta ) ); return value; } /******************************************************************************/ void c4mat_print ( int m, int n, float complex a[], char *title ) /******************************************************************************/ /* Purpose: c4mat_print() prints a C4MAT. Discussion: A C4MAT is a matrix of single precision complex values. Licensing: This code is distributed under the MIT license. Modified: 09 March 2014 Author: John Burkardt Input: int M, N, the number of rows and columns in the matrix. float complex A[M*N], the matrix. char *TITLE, a title. */ { c4mat_print_some ( m, n, a, 1, 1, m, n, title ); return; } /******************************************************************************/ void c4mat_print_some ( int m, int n, float complex a[], int ilo, int jlo, int ihi, int jhi, char *title ) /******************************************************************************/ /* Purpose: c4mat_print_some() prints some of a C4MAT. Discussion: A C4MAT is a matrix of float complex values. Licensing: This code is distributed under the MIT license. Modified: 09 March 2014 Author: John Burkardt Input: int M, N, the number of rows and columns in the matrix. float complex A[M*N], the matrix. int ILO, JLO, IHI, JHI, the first row and column, and the last row and column to be printed. char *TITLE, a title. */ { float complex c; int i; int i2hi; int i2lo; int inc; int incx = 4; int j; int j2; int j2hi; int j2lo; printf ( "\n" ); printf ( "%s\n", title ); /* Print the columns of the matrix, in strips of INCX. */ for ( j2lo = jlo; j2lo <= jhi; j2lo = j2lo + incx ) { j2hi = j2lo + incx - 1; if ( n < j2hi ) { j2hi = n; } if ( jhi < j2hi ) { j2hi = jhi; } inc = j2hi + 1 - j2lo; printf ( "\n" ); printf ( " Col: " ); for ( j = j2lo; j <= j2hi; j++ ) { j2 = j + 1 - j2lo; printf ( " %10d", j ); } printf ( "\n" ); printf ( " Row\n" ); printf ( " ---\n" ); /* Determine the range of the rows in this strip. */ i2lo = 1; if ( i2lo < ilo ) { i2lo = ilo; } i2hi = m; if ( ihi < i2hi ) { i2hi = ihi; } for ( i = i2lo; i <= i2hi; i++ ) { /* Print out (up to) INCX entries in row I, that lie in the current strip. */ for ( j2 = 1; j2 <= inc; j2++ ) { j = j2lo - 1 + j2; c = a[i-1+(j-1)*m]; printf ( " %8g %8g", creal ( c ), cimag ( c ) ); } printf ( "\n" ); } } return; } /******************************************************************************/ void c4mat_uniform_01 ( int m, int n, float complex c[] ) /******************************************************************************/ /* Purpose: c4mat_uniform_01() returns a unit pseudorandom C4MAT. Discussion: The angles should be uniformly distributed between 0 and 2 * PI, the square roots of the radius uniformly distributed between 0 and 1. This results in a uniform distribution of values in the unit circle. Licensing: This code is distributed under the MIT license. Modified: 12 January 2007 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. Input: int M, N, the number of rows and columns in the matrix. Output: float complex C[M*N], the pseudorandom complex matrix. */ { int i; int j; float r; const float r4_pi = 3.1415926; float theta; for ( j = 0; j < n; j++ ) { for ( i = 0; i < m; i++ ) { r = sqrt ( drand48 ( ) ); theta = 2.0 * r4_pi * ( drand48 ( ) ); c[i+j*m] = r * ( cos ( theta ) + I * sin ( theta ) ); } } return; } /******************************************************************************/ float complex *c4mat_uniform_01_new ( int m, int n ) /******************************************************************************/ /* Purpose: c4mat_uniform_01_new() returns a unit pseudorandom C4MAT. Discussion: The angles should be uniformly distributed between 0 and 2 * PI, the square roots of the radius uniformly distributed between 0 and 1. This results in a uniform distribution of values in the unit circle. Licensing: This code is distributed under the MIT license. Modified: 12 January 2007 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. Input: int M, N, the number of rows and columns in the matrix. Output: float complex C4MAT_UNIFORM_01_NEW[M*N], the pseudorandom complex matrix. */ { float complex *c; int i; int j; float r; const float r4_pi = 3.1415926; float theta; c = ( float complex * ) malloc ( m * n * sizeof ( float complex ) ); for ( j = 0; j < n; j++ ) { for ( i = 0; i < m; i++ ) { r = sqrt ( drand48 ( ) ); theta = 2.0 * r4_pi * ( drand48 ( ) ); c[i+j*m] = r * ( cos ( theta ) + I * sin ( theta ) ); } } return c; } /******************************************************************************/ void c4vec_print ( int n, float complex a[], char *title ) /******************************************************************************/ /* Purpose: c4vec_print() prints a C4VEC. Discussion: A C4VEC is a vector of float complex values. Licensing: This code is distributed under the MIT license. Modified: 09 March 2014 Author: John Burkardt Input: int N, the number of components of the vector. float complex A[N], the vector to be printed. 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, " %8d: %14f %14f\n", i, creal( a[i] ), cimag ( a[i] ) ); } return; } /******************************************************************************/ void c4vec_uniform_01 ( int n, float complex c[] ) /******************************************************************************/ /* Purpose: c4vec_uniform_01() returns a unit pseudorandom C4VEC. Discussion: The angles should be uniformly distributed between 0 and 2 * PI, the square roots of the radius uniformly distributed between 0 and 1. This results in a uniform distribution of values in the unit circle. Licensing: This code is distributed under the MIT license. Modified: 12 January 2007 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. Input: int N, the number of values to compute. Output: float complex C[N], the pseudorandom complex vector. */ { int i; float r; const float r4_pi = 3.1415926; float theta; for ( i = 0; i < n; i++ ) { r = sqrt ( drand48 ( ) ); theta = 2.0 * r4_pi * ( drand48 ( ) ); c[i] = r * ( cos ( theta ) + I * sin ( theta ) ); } return; } /******************************************************************************/ float complex *c4vec_uniform_01_new ( int n ) /******************************************************************************/ /* Purpose: c4vec_uniform_01_new() returns a unit pseudorandom C4VEC. Discussion: The angles should be uniformly distributed between 0 and 2 * PI, the square roots of the radius uniformly distributed between 0 and 1. This results in a uniform distribution of values in the unit circle. Licensing: This code is distributed under the MIT license. Modified: 12 January 2007 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. Input: int N, the number of values to compute. Output: float complex C4VEC_UNIFORM_01_NEW[N], the pseudorandom complex vector. */ { float complex *c; int i; float r; const float r4_pi = 3.1415926; float theta; c = ( float complex * ) malloc ( n * sizeof ( float complex ) ); for ( i = 0; i < n; i++ ) { r = sqrt ( drand48 ( ) ); theta = 2.0 * r4_pi * ( drand48 ( ) ); c[i] = r * ( cos ( theta ) + I * sin ( theta ) ); } return c; } /******************************************************************************/ double complex c8_uniform_01 ( ) /******************************************************************************/ /* Purpose: c8_uniform_01() returns a unit pseudorandom C8. Discussion: The angle should be uniformly distributed between 0 and 2 * PI, the square root of the radius uniformly distributed between 0 and 1. This results in a uniform distribution of values in the unit circle. Licensing: This code is distributed under the MIT license. Modified: 11 January 2007 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. Output: double complex C8_UNIFORM_01, a pseudorandom complex value. */ { double r; const double r8_pi = 3.141592653589793; double theta; double complex value; r = sqrt ( drand48 ( ) ); theta = 2.0 * r8_pi * ( drand48 ( ) ); value = r * ( cos ( theta ) + I * sin ( theta ) ); return value; } /******************************************************************************/ void c8mat_print ( int m, int n, double complex a[], char *title ) /******************************************************************************/ /* Purpose: c8mat_print() prints a C8MAT. Discussion: A C8MAT is a matrix of double precision complex values. Licensing: This code is distributed under the MIT license. Modified: 08 July 2011 Author: John Burkardt Input: int M, N, the number of rows and columns in the matrix. double complex A[M*N], the matrix. char *TITLE, a title. */ { c8mat_print_some ( m, n, a, 1, 1, m, n, title ); return; } /******************************************************************************/ void c8mat_print_some ( int m, int n, double complex a[], int ilo, int jlo, int ihi, int jhi, char *title ) /******************************************************************************/ /* Purpose: c8mat_print_some() prints some of a C8MAT. Discussion: A C8MAT is a matrix of double precision complex values. Licensing: This code is distributed under the MIT license. Modified: 08 July 2011 Author: John Burkardt Input: int M, N, the number of rows and columns in the matrix. double complex A[M*N], the matrix. int ILO, JLO, IHI, JHI, the first row and column, and the last row and column to be printed. char *TITLE, a title. */ { double complex c; int i; int i2hi; int i2lo; int inc; int incx = 4; int j; int j2; int j2hi; int j2lo; printf ( "\n" ); printf ( "%s\n", title ); /* Print the columns of the matrix, in strips of INCX. */ for ( j2lo = jlo; j2lo <= jhi; j2lo = j2lo + incx ) { j2hi = j2lo + incx - 1; if ( n < j2hi ) { j2hi = n; } if ( jhi < j2hi ) { j2hi = jhi; } inc = j2hi + 1 - j2lo; printf ( "\n" ); printf ( " Col: " ); for ( j = j2lo; j <= j2hi; j++ ) { j2 = j + 1 - j2lo; printf ( " %10d", j ); } printf ( "\n" ); printf ( " Row\n" ); printf ( " ---\n" ); /* Determine the range of the rows in this strip. */ i2lo = 1; if ( i2lo < ilo ) { i2lo = ilo; } i2hi = m; if ( ihi < i2hi ) { i2hi = ihi; } for ( i = i2lo; i <= i2hi; i++ ) { /* Print out (up to) INCX entries in row I, that lie in the current strip. */ for ( j2 = 1; j2 <= inc; j2++ ) { j = j2lo - 1 + j2; c = a[i-1+(j-1)*m]; printf ( " %8g %8g", creal ( c ), cimag ( c ) ); } printf ( "\n" ); } } return; } /******************************************************************************/ void c8mat_uniform_01 ( int m, int n, double complex c[] ) /******************************************************************************/ /* Purpose: c8mat_uniform_01() returns a unit pseudorandom C8MAT. Discussion: The angles should be uniformly distributed between 0 and 2 * PI, the square roots of the radius uniformly distributed between 0 and 1. This results in a uniform distribution of values in the unit circle. Licensing: This code is distributed under the MIT license. Modified: 12 January 2007 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. Input: int M, N, the number of rows and columns in the matrix. Output: double complex C[M*N], the pseudorandom complex matrix. */ { int i; int j; double r; const double r8_pi = 3.141592653589793; double theta; for ( j = 0; j < n; j++ ) { for ( i = 0; i < m; i++ ) { r = sqrt ( drand48 ( ) ); theta = 2.0 * r8_pi * ( drand48 ( ) ); c[i+j*m] = r * ( cos ( theta )+ I * sin ( theta ) ); } } return; } /******************************************************************************/ double complex *c8mat_uniform_01_new ( int m, int n ) /******************************************************************************/ /* Purpose: c8mat_uniform_01_new() returns a unit pseudorandom C8MAT. Discussion: The angles should be uniformly distributed between 0 and 2 * PI, the square roots of the radius uniformly distributed between 0 and 1. This results in a uniform distribution of values in the unit circle. Licensing: This code is distributed under the MIT license. Modified: 12 January 2007 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. Input: int M, N, the number of rows and columns in the matrix. Output: double complex C8MAT_UNIFORM_01_NEW[M*N], the pseudorandom complex matrix. */ { double complex *c; int i; int j; double r; const double r8_pi = 3.141592653589793; double theta; c = ( double complex * ) malloc ( m * n * sizeof ( double complex ) ); for ( j = 0; j < n; j++ ) { for ( i = 0; i < m; i++ ) { r = sqrt ( drand48 ( ) ); theta = 2.0 * r8_pi * ( drand48 ( ) ); c[i+j*m] = r * ( cos ( theta )+ I * sin ( theta ) ); } } return c; } /******************************************************************************/ void c8vec_print ( int n, double complex a[], char *title ) /******************************************************************************/ /* Purpose: c8vec_print() prints a C8VEC. Discussion: A C8VEC is a vector of double complex values. Licensing: This code is distributed under the MIT license. Modified: 27 January 2013 Author: John Burkardt Input: int N, the number of components of the vector. double complex A[N], the vector to be printed. 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, " %8d: %14f %14f\n", i, creal( a[i] ), cimag ( a[i] ) ); } return; } /******************************************************************************/ void c8vec_uniform_01 ( int n, double complex c[] ) /******************************************************************************/ /* Purpose: c8vec_uniform_01() returns a unit pseudorandom C8VEC. Discussion: The angles should be uniformly distributed between 0 and 2 * PI, the square roots of the radius uniformly distributed between 0 and 1. This results in a uniform distribution of values in the unit circle. Licensing: This code is distributed under the MIT license. Modified: 12 January 2007 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. Input: int N, the number of values to compute. Output: double complex C[N], the pseudorandom vector. */ { int i; double r; const double r8_pi = 3.141592653589793; double theta; for ( i = 0; i < n; i++ ) { r = sqrt ( drand48 ( ) ); theta = 2.0 * r8_pi * ( drand48 ( ) ); c[i] = r * ( cos ( theta ) + I * sin ( theta ) ); } return; } /******************************************************************************/ double complex *c8vec_uniform_01_new ( int n ) /******************************************************************************/ /* Purpose: c8vec_uniform_01_new() returns a unit pseudorandom C8VEC. Discussion: The angles should be uniformly distributed between 0 and 2 * PI, the square roots of the radius uniformly distributed between 0 and 1. This results in a uniform distribution of values in the unit circle. Licensing: This code is distributed under the MIT license. Modified: 12 January 2007 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. Input: int N, the number of values to compute. Output: double complex C8VEC_UNIFORM_01_NEW[N], the pseudorandom vector. */ { double complex *c; int i; double r; const double r8_pi = 3.141592653589793; double theta; c = ( double complex * ) malloc ( n * sizeof ( double complex ) ); for ( i = 0; i < n; i++ ) { r = sqrt ( drand48 ( ) ); theta = 2.0 * r8_pi * ( drand48 ( ) ); c[i] = r * ( cos ( theta ) + I * sin ( theta ) ); } return c; } /*****************************************************************************/ char ch_uniform_ab ( char a, char b ) /******************************************************************************/ /* Purpose: ch_uniform_ab() returns a scaled CH. Licensing: This code is distributed under the MIT license. Modified: 22 May 2008 Author: John Burkardt Input: char A, B, the minimum and maximum acceptable characters. Output: char CH_UNIFORM, the randomly chosen character. */ { char c; float r; r = r4_uniform_01 ( ); c = a + ( char ) ( r * ( float ) ( b + 1 - a ) ); return c; } /******************************************************************************/ int congruence ( int a, int b, int c, int *error ) /******************************************************************************/ /* Purpose: congruence() solves a congruence of the form A * X = C ( mod B ). Discussion: A, B and C are given integers. The equation is solvable if and only if the greatest common divisor of A and B also divides C. Licensing: This code is distributed under the MIT license. Modified: 15 November 2004 Author: John Burkardt Reference: Eric Weisstein, editor, CRC Concise Encylopedia of Mathematics, CRC Press, 1998, page 446. Input: int A, B, C, the coefficients of the Diophantine equation. int *ERROR, error flag, is 1 if an error occurred.. int CONGRUENCE, the solution of the Diophantine equation. X will be between 0 and B-1. */ { # define N_MAX 100 int a_copy; int a_mag; int a_sign; int b_copy; int b_mag; int c_copy; int g; int k; int n; int q[N_MAX]; int swap; int x; int y; int z; /* Defaults for output parameters. */ *error = 0; x = 0; y = 0; /* Special cases. */ if ( a == 0 && b == 0 && c == 0 ) { x = 0; return x; } else if ( a == 0 && b == 0 && c != 0 ) { *error = 1; x = 0; return x; } else if ( a == 0 && b != 0 && c == 0 ) { x = 0; return x; } else if ( a == 0 && b != 0 && c != 0 ) { x = 0; if ( ( c % b ) != 0 ) { *error = 2; } return x; } else if ( a != 0 && b == 0 && c == 0 ) { x = 0; return x; } else if ( a != 0 && b == 0 && c != 0 ) { x = c / a; if ( ( c % a ) != 0 ) { *error = 3; return x; } return x; } else if ( a != 0 && b != 0 && c == 0 ) { /* g = i4_gcd ( a, b ); */ /* x = b / g; */ x = 0; return x; } /* Now handle the "general" case: A, B and C are nonzero. Step 1: Compute the GCD of A and B, which must also divide C. */ g = i4_gcd ( a, b ); if ( ( c % g ) != 0 ) { *error = 4; return x; } a_copy = a / g; b_copy = b / g; c_copy = c / g; /* Step 2: Split A and B into sign and magnitude. */ a_mag = abs ( a_copy ); a_sign = i4_sign ( a_copy ); b_mag = abs ( b_copy ); /* Another special case, A_MAG = 1 or B_MAG = 1. */ if ( a_mag == 1 ) { x = a_sign * c_copy; return x; } else if ( b_mag == 1 ) { x = 0; return x; } /* Step 3: Produce the Euclidean remainder sequence. */ if ( b_mag <= a_mag ) { swap = 0; q[0] = a_mag; q[1] = b_mag; } else { swap = 1; q[0] = b_mag; q[1] = a_mag; } n = 3; for ( ; ; ) { q[n-1] = ( q[n-3] % q[n-2] ); if ( q[n-1] == 1 ) { break; } n = n + 1; if ( N_MAX < n ) { *error = 1; fprintf ( stderr, "\n" ); fprintf ( stderr, "CONGRUENCE - Fatal error!\n" ); fprintf ( stderr, " Exceeded number of iterations.\n" ); exit ( 1 ); } } /* Step 4: Now go backwards to solve X * A_MAG + Y * B_MAG = 1. */ y = 0; for ( k = n; 2 <= k; k-- ) { x = y; y = ( 1 - x * q[k-2] ) / q[k-1]; } /* Step 5: Undo the swapping. */ if ( swap == 1 ) { z = x; x = y; y = z; } /* Step 6: Now apply signs to X and Y so that X * A + Y * B = 1. */ x = x * a_sign; /* Step 7: Multiply by C, so that X * A + Y * B = C. */ x = x * c_copy; /* Step 8: Now force 0 <= X < B. */ x = x % b; /* Step 9: Force positivity. */ if ( x < 0 ) { x = x + b; } return x; # undef N_MAX } /******************************************************************************/ char digit_to_ch ( int i ) /******************************************************************************/ /* Purpose: digit_to_ch() returns the base 10 digit character corresponding to a digit. Example: I C ----- --- 0 '0' 1 '1' ... ... 9 '9' 10 '*' -83 '*' Licensing: This code is distributed under the MIT license. Modified: 16 June 2003 Author: John Burkardt Input: int I, the digit, which should be between 0 and 9. char DIGIT_TO_CH, the appropriate character '0' through '9' or '*'. */ { char c; if ( 0 <= i && i <= 9 ) { c = '0' + i; } else { c = '*'; } return c; } /******************************************************************************/ int i4_gcd ( int i, int j ) /******************************************************************************/ /* Purpose: i4_gcd() finds the greatest common divisor of I and J. Discussion: Only the absolute values of I and J are considered, so that the result is always nonnegative. If I or J is 0, I4_GCD is returned as max ( 1, abs ( I ), abs ( J ) ). If I and J have no common factor, I4_GCD is returned as 1. Otherwise, using the Euclidean algorithm, I4_GCD is the largest common factor of I and J. Licensing: This code is distributed under the MIT license. Modified: 07 May 2003 Author: John Burkardt Input: int I, J, two numbers whose greatest common divisor is desired. int I4_GCD, the greatest common divisor of I and J. */ { int ip; int iq; int ir; /* Return immediately if either I or J is zero. */ if ( i == 0 ) { return i4_max ( 1, abs ( j ) ); } else if ( j == 0 ) { return i4_max ( 1, abs ( i ) ); } /* Set IP to the larger of I and J, IQ to the smaller. This way, we can alter IP and IQ as we go. */ ip = i4_max ( abs ( i ), abs ( j ) ); iq = i4_min ( abs ( i ), abs ( j ) ); /* Carry out the Euclidean algorithm. */ for ( ; ; ) { ir = ip % iq; if ( ir == 0 ) { break; } ip = iq; iq = ir; } return iq; } /******************************************************************************/ 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 Input: int I4_HUGE, a "huge" integer. */ { return 2147483647; } /******************************************************************************/ int i4_log_10 ( int i ) /******************************************************************************/ /* Purpose: i4_log_10() returns the whole part of the logarithm base 10 of an I4. Discussion: It should be the case that 10^I4_LOG_10(I) <= |I| < 10^(I4_LOG_10(I)+1). (except for I = 0). The number of decimal digits in I is I4_LOG_10(I) + 1. Example: I I4_LOG_10(I) 0 0 1 0 2 0 9 0 10 1 11 1 99 1 100 2 101 2 999 2 1000 3 1001 3 9999 3 10000 4 10001 4 Licensing: This code is distributed under the MIT license. Modified: 17 June 2003 Author: John Burkardt Input: int I, the integer. int I4_LOG_10, the whole part of the logarithm of abs ( I ). */ { int ten_pow; int value; i = abs ( i ); ten_pow = 10; value = 0; while ( ten_pow <= i ) { ten_pow = ten_pow * 10; value = value + 1; } 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 Input: int I1, I2, are two integers to be compared. 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 Input: int I1, I2, two integers to be compared. int I4_MIN, the smaller of I1 and I2. */ { int value; if ( i1 < i2 ) { value = i1; } else { value = i2; } return value; } /******************************************************************************/ int i4_sign ( int i ) /******************************************************************************/ /* Purpose: i4_sign() returns the sign of an I4. Discussion: The sign of 0 and all positive integers is taken to be +1. The sign of all negative integers is -1. Licensing: This code is distributed under the MIT license. Modified: 06 May 2003 Author: John Burkardt Input: int I, the integer whose sign is desired. int I4_SIGN, the sign of I. */ { int value; if ( i < 0 ) { value = - 1; } else { value = 1; } return value; } /******************************************************************************/ void i4_swap ( int *i, int *j ) /******************************************************************************/ /* Purpose: i4_swap() switches two I4's. Licensing: This code is distributed under the MIT license. Modified: 07 January 2002 Author: John Burkardt Input: Iint *I, *J: two values Output: int *I, *J: the values have been interchanged. */ { int k; k = *i; *i = *j; *j = k; return; } /******************************************************************************/ char *i4_to_s ( int i ) /******************************************************************************/ /* Purpose: i4_to_s() converts an I4 to a string. Example: INTVAL S 1 1 -1 -1 0 0 1952 1952 123456 123456 1234567 1234567 Licensing: This code is distributed under the MIT license. Modified: 13 March 2004 Author: John Burkardt Input: int I, an integer to be converted. char *I4_TO_S, the representation of the integer. */ { int digit; int j; int length; int ten_power; char *s; length = i4_log_10 ( i ); ten_power = ( int ) pow ( ( double ) 10, ( double ) length ); if ( i < 0 ) { length = length + 1; } /* Add one position for the trailing null. */ length = length + 1; s = ( char * ) malloc ( length * sizeof ( char ) ); if ( i == 0 ) { s[0] = '0'; s[1] = '\0'; return s; } /* Now take care of the sign. */ j = 0; if ( i < 0 ) { s[j] = '-'; j = j + 1; i = abs ( i ); } /* Find the leading digit of I, strip it off, and stick it into the string. */ while ( 0 < ten_power ) { digit = i / ten_power; s[j] = digit_to_ch ( digit ); j = j + 1; i = i - digit * ten_power; ten_power = ten_power / 10; } /* Tack on the trailing NULL. */ s[j] = '\0'; j = j + 1; return s; } /******************************************************************************/ int i4_uniform_ab ( int a, int b ) /******************************************************************************/ /* 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. Input: int A, B, the limits of the interval. Output: int I4_UNIFORM_AB, a number between A and B. */ { int c; float r; int value; /* Guarantee A <= B. */ if ( b < a ) { c = a; a = b; b = c; } r = drand48 ( ); /* Scale R to lie between A-0.5 and B+0.5. */ r = ( 1.0 - r ) * ( ( float ) ( a ) - 0.5 ) + r * ( ( float ) ( b ) + 0.5 ); /* Round R to the nearest integer. */ value = round ( r ); /* Guarantee that A <= VALUE <= B. */ if ( value < a ) { value = a; } if ( b < value ) { value = b; } return value; } /******************************************************************************/ void i4mat_print ( int m, int n, int a[], char *title ) /******************************************************************************/ /* Purpose: i4mat_print() prints an I4MAT. Discussion: An I4MAT is an MxN array of I4's, stored by (I,J) -> [I+J*M]. Licensing: This code is distributed under the MIT license. Modified: 28 May 2008 Author: John Burkardt Input: int M, the number of rows in A. int N, the number of columns in A. int A[M*N], the M by N matrix. 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 Input: int M, the number of rows of the matrix. M must be positive. int N, the number of columns of the matrix. N must be positive. int A[M*N], the matrix. int ILO, JLO, IHI, JHI, designate the first row and column, and the last row and column to be printed. 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; if ( n < j2hi ) { j2hi = n; } if ( jhi < j2hi ) { 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. */ if ( 1 < ilo ) { i2lo = ilo; } else { i2lo = 1; } if ( m < ihi ) { i2hi = m; } else { i2hi = ihi; } 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 i4mat_uniform_ab ( int m, int n, int a, int b, int x[] ) /******************************************************************************/ /* Purpose: i4mat_uniform_ab() returns a scaled pseudorandom 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: 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. Input: int M, N, the number of rows and columns. int A, B, the limits of the pseudorandom values. Output: int X[M*N], a matrix of pseudorandom values. */ { int c; int i; int j; float r; int value; /* Guaranteee A <= B. */ if ( b < a ) { c = a; a = b; b = c; } for ( j = 0; j < n; j++ ) { for ( i = 0; i < m; i++ ) { r = drand48 ( ); /* 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+j*m] = value; } } return; } /******************************************************************************/ int *i4mat_uniform_ab_new ( int m, int n, int a, int b ) /******************************************************************************/ /* Purpose: i4mat_uniform_ab_new() returns a scaled pseudorandom 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: 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. Input: int M, N, the number of rows and columns. int A, B, the limits of the pseudorandom values. Output: int I4MAT_UNIFORM_AB_NEW[M*N], a matrix of pseudorandom values. */ { int c; int i; int j; float r; int value; int *x; /* Guaranteee A <= B. */ if ( b < a ) { c = a; a = b; b = c; } x = ( int * ) malloc ( m * n * sizeof ( int ) ); for ( j = 0; j < n; j++ ) { for ( i = 0; i < m; i++ ) { r = drand48 ( ); /* 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+j*m] = value; } } return x; } /******************************************************************************/ int i4vec_max ( int n, int a[] ) /******************************************************************************/ /* Purpose: i4vec_max() returns the value of the maximum element in an I4VEC. Licensing: This code is distributed under the MIT license. Modified: 17 May 2003 Author: John Burkardt Input: int N, the number of entries in the array. int A[N], the array to be checked. 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; } /******************************************************************************/ float i4vec_mean ( int n, int x[] ) /******************************************************************************/ /* Purpose: i4vec_mean() returns the mean of an I4VEC. Licensing: This code is distributed under the MIT license. Modified: 01 May 1999 Author: John Burkardt Input: int N, the number of entries in the vector. int X[N], the vector whose mean is desired. float I4VEC_MEAN, the mean, or average, of the vector entries. */ { int i; float mean; mean = 0.0; for ( i = 0; i < n; i++ ) { mean = mean + ( float ) x[i]; } mean = mean / ( float ) n; return mean; } /******************************************************************************/ int i4vec_min ( int n, int a[] ) /******************************************************************************/ /* Purpose: i4vec_min() returns the minimum element in an I4VEC. Licensing: This code is distributed under the MIT license. Modified: 17 May 2003 Author: John Burkardt Input: int N, the number of entries in the array. int A[N], the array to be checked. int I4VEC_MIN, the value of the minimum 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 ( a[i] < value ) { value = a[i]; } } return value; } /******************************************************************************/ void i4vec_print ( int n, int a[], char *title ) /******************************************************************************/ /* Purpose: i4vec_print() prints an I4VEC. Discussion: An I4VEC is a vector of I4's. Licensing: This code is distributed under the MIT license. Modified: 14 November 2003 Author: John Burkardt Input: int N, the number of components of the vector. int A[N], the vector to be printed. 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_uniform_ab ( int n, int a, int b, int x[] ) /******************************************************************************/ /* Purpose: i4vec_uniform_ab() 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: 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. Input: integer N, the dimension of the vector. int A, B, the limits of the interval. Output: int X[N], a vector of random values between A and B. */ { int c; int i; float r; int value; /* Guaranteee A <= B. */ if ( b < a ) { c = a; a = b; b = c; } for ( i = 0; i < n; i++ ) { r = drand48 ( ); /* 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; } /******************************************************************************/ int *i4vec_uniform_ab_new ( int n, int a, int b ) /******************************************************************************/ /* 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: 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. Input: integer N, the dimension of the vector. int A, B, the limits of the interval. Output: int I4VEC_UNIFORM_AB_NEW[N], a vector of random values between A and B. */ { int c; int i; float r; int value; int *x; /* Guarantee A <= B. */ if ( b < a ) { c = a; a = b; b = c; } x = ( int * ) malloc ( n * sizeof ( int ) ); for ( i = 0; i < n; i++ ) { r = drand48 ( ); /* 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; } /******************************************************************************/ float i4vec_variance ( int n, int x[] ) /******************************************************************************/ /* Purpose: i4vec_variance() returns the variance of an I4VEC. Licensing: This code is distributed under the MIT license. Modified: 01 May 1999 Author: John Burkardt Input: int N, the number of entries in the vector. int X[N], the vector whose variance is desired. float I4VEC_VARIANCE, the variance of the vector entries. */ { int i; float mean; float variance; mean = i4vec_mean ( n, x ); variance = 0.0; for ( i = 0; i < n; i++ ) { variance = variance + ( ( float ) x[i] - mean ) * ( ( float ) x[i] - mean ); } if ( 1 < n ) { variance = variance / ( float ) ( n - 1 ); } else { variance = 0.0; } return variance; } /******************************************************************************/ int l4_uniform ( ) /******************************************************************************/ /* Purpose: l4_uniform() returns a pseudorandom L4. Discussion: An L4 is a LOGICAL value. Licensing: This code is distributed under the MIT license. Modified: 23 April 2008 Author: John Burkardt Output: int l4_uniform: a pseudorandom logical value. */ { double r; int value; r = drand48 ( ); value = ( r < 0.5 ); return value; } /******************************************************************************/ void l4mat_print ( int m, int n, int a[], char *title ) /******************************************************************************/ /* Purpose: l4mat_print() prints an L4MAT. Discussion: An L4MAT is an array of L4 values. Licensing: This code is distributed under the MIT license. Modified: 03 November 2011 Author: John Burkardt Input: int M, the number of rows in A. int N, the number of columns in A. int A[M*N], the matrix. char *TITLE, a title. */ { l4mat_print_some ( m, n, a, 0, 0, m - 1, n - 1, title ); return; } /******************************************************************************/ void l4mat_print_some ( int m, int n, int a[], int ilo, int jlo, int ihi, int jhi, char *title ) /******************************************************************************/ /* Purpose: l4mat_print_some() prints some of an L4MAT. Discussion: An L4MAT is an array of L4 values. Licensing: This code is distributed under the MIT license. Modified: 03 November 2011 Author: John Burkardt Input: int M, N, the number of rows and columns. int A[M*N], an M by N matrix to be printed. int ILO, JLO, the first row and column to print. int IHI, JHI, the last row and column to print. char *TITLE, a title. */ { int i; int i2hi; int i2lo; int incx = 35; int j; int j2hi; int j2lo; printf ( "\n" ); printf ( "%s\n", title ); for ( j2lo = i4_max ( jlo, 0 ); j2lo <= i4_min ( jhi, n - 1 ); j2lo = j2lo + incx ) { j2hi = j2lo + incx - 1; if ( n - 1 < j2hi ) { j2hi = n - 1; } if ( jhi < j2hi ) { j2hi = jhi; } printf ( "\n" ); if ( 100 <= j2hi ) { printf ( " " ); for ( j = j2lo; j <= j2hi; j++ ) { printf ( " %1d", j / 100 ); } printf ( "\n" ); } if ( 10 <= j2hi ) { printf ( " " ); for ( j = j2lo; j <= j2hi; j++ ) { printf ( " %1d", ( ( j / 10 ) % 10 ) ); } printf ( "\n" ); } printf ( " Col " ); for ( j = j2lo; j <= j2hi; j++ ) { printf ( " %1d", ( j % 10 ) ); } printf ( "\n" ); printf ( " Row\n" ); printf ( "\n" ); i2lo = 0; if ( i2lo < ilo ) { i2lo = ilo; } i2hi = m - 1; if ( ihi < i2hi ) { i2hi = ihi; } for ( i = i2lo; i <= i2hi; i++ ) { printf ( "%5d:", i ); for ( j = j2lo; j <= j2hi; j++ ) { printf ( " %1d", a[i+j*m] ); } printf ( "\n" ); } } return; } /******************************************************************************/ int *l4mat_uniform_new ( int m, int n ) /******************************************************************************/ /* Purpose: l4mat_uniform_new() returns a pseudorandom L4MAT. Discussion: An LMAT is a two dimensional array of LOGICAL values. Licensing: This code is distributed under the MIT license. Modified: 23 April 2008 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. Input: int M, N, the order of the matrix. Output: int L4MAT_UNIFORM_NEW[M*N], a pseudorandom logical matrix. */ { int i; int j; int *l4mat; double r; l4mat = ( int * ) malloc ( m * n * sizeof ( int ) ); for ( j = 0; j < n; j++ ) { for ( i = 0; i < m; i++ ) { r = drand48 ( ); l4mat[i+j*m] = ( r < 0.5 ); } } return l4mat; } /******************************************************************************/ void l4vec_print ( int n, int a[], char *title ) /******************************************************************************/ /* Purpose: l4vec_print() prints an L4VEC. Licensing: This code is distributed under the MIT license. Modified: 06 May 2014 Author: John Burkardt Input: int N, the number of components of the vector. int A[N], the (logical) vector to be printed. char *TITLE, a title. */ { int i; printf ( "\n" ); printf ( "%s\n", title ); printf ( "\n" ); for ( i = 0; i < n; i++ ) { if ( a[i] == 0 ) { printf ( " %8d: F\n", i ); } else { printf ( " %8d: T\n", i ); } } return; } /******************************************************************************/ int *l4vec_uniform_new ( int n ) /******************************************************************************/ /* Purpose: l4vec_uniform_new() returns a pseudorandom L4VEC. Discussion: An L4VEC is a vector of LOGICAL values. Licensing: This code is distributed under the MIT license. Modified: 23 April 2008 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 LEcuyer, 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. Input: int N, the order of the vector. Output: int L4VEC_UNIFORM_NEW[N], a pseudorandom logical vector. /*/ { int i; int *l4vec; double r; l4vec = ( int * ) malloc ( n * sizeof ( int ) ); for ( i = 0; i < n; i++ ) { r = drand48 ( ); l4vec[i] = ( r < 0.5 ); } return l4vec; } /******************************************************************************/ int power_mod ( int a, int n, int m ) /******************************************************************************/ /* Purpose: power_mod() computes mod ( A^N, M ). Discussion: Some programming tricks are used to speed up the computation, and to allow computations in which A^N is much too large to store in a real word. First, for efficiency, the power A^N is computed by determining the binary expansion of N, then computing A, A^2, A^4, and so on by repeated squaring, and multiplying only those factors that contribute to A^N. Secondly, the intermediate products are immediately "mod'ed", which keeps them small. For instance, to compute mod ( A^13, 11 ), we essentially compute 13 = 1 + 4 + 8 A^13 = A * A^4 * A^8 mod ( A^13, 11 ) = mod ( A, 11 ) * mod ( A^4, 11 ) * mod ( A^8, 11 ). Fermat's little theorem says that if P is prime, and A is not divisible by P, then ( A^(P-1) - 1 ) is divisible by P. Licensing: This code is distributed under the MIT license. Modified: 16 November 2004 Author: John Burkardt Input: int A, the base of the expression to be tested. A should be nonnegative. int N, the power to which the base is raised. N should be nonnegative. int M, the divisor against which the expression is tested. M should be positive. int POWER_MOD, the remainder when A**N is divided by M. */ { long long int a_square2; int d; long long int m2; int x; long long int x2; if ( a < 0 ) { return -1; } if ( m <= 0 ) { return -1; } if ( n < 0 ) { return -1; } /* A_SQUARE contains the successive squares of A. */ a_square2 = ( long long int ) a; m2 = ( long long int ) m; x2 = ( long long int ) 1; while ( 0 < n ) { d = n % 2; if ( d == 1 ) { x2 = ( x2 * a_square2 ) % m2; } a_square2 = ( a_square2 * a_square2 ) % m2; n = ( n - d ) / 2; } /* Ensure that 0 <= X. */ while ( x2 < 0 ) { x2 = x2 + m2; } x = ( int ) x2; return x; } /******************************************************************************/ int r4_nint ( float x ) /******************************************************************************/ /* Purpose: r4_nint() returns the nearest integer to an R4. Example: X R4_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 Input: float X, the value. int R4_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; } /******************************************************************************/ float r4_uniform_ab ( float a, float b ) /******************************************************************************/ /* Purpose: r4_uniform_ab() returns a scaled pseudorandom R4. Discussion: The pseudorandom number should be uniformly distributed between A and B. Licensing: This code is distributed under the MIT license. Modified: 19 April 2011 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. Input: float A, B, the limits of the interval. Output: float R4_UNIFORM_AB, a number strictly between A and B. */ { float value; value = a + ( b - a ) * drand48 ( ); return value; } /******************************************************************************/ float r4_uniform_01 ( ) /******************************************************************************/ /* Purpose: r4_uniform_01() returns a unit pseudorandom R4. Licensing: This code is distributed under the MIT license. Modified: 16 November 2004 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. Output: float R4_UNIFORM_01, a new pseudorandom variate, strictly between 0 and 1. */ { float value; value = drand48 ( ); return value; } /******************************************************************************/ void r4mat_print ( int m, int n, float a[], char *title ) /******************************************************************************/ /* Purpose: r4mat_print() prints an R4MAT. Discussion: An R4MAT is a doubly dimensioned array of R4's, which may be stored as a vector in column-major order. Entry A(I,J) is stored as A[I+J*M] Licensing: This code is distributed under the MIT license. Modified: 28 May 2008 Author: John Burkardt Input: int M, the number of rows in A. int N, the number of columns in A. float A[M*N], the M by N matrix. char *TITLE, a title. */ { r4mat_print_some ( m, n, a, 1, 1, m, n, title ); return; } /******************************************************************************/ void r4mat_print_some ( int m, int n, float a[], int ilo, int jlo, int ihi, int jhi, char *title ) /******************************************************************************/ /* Purpose: r4mat_print_some() prints some of an R4MAT. Discussion: An R4MAT is a doubly dimensioned array of R4's, which may be stored as a vector in column-major order. Licensing: This code is distributed under the MIT license. Modified: 20 August 2010 Author: John Burkardt Input: int M, the number of rows of the matrix. M must be positive. int N, the number of columns of the matrix. N must be positive. float A[M*N], the matrix. int ILO, JLO, IHI, JHI, designate the first row and column, and the last row and column to be printed. char *TITLE, a title. */ { # define INCX 5 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 5. */ for ( j2lo = jlo; j2lo <= jhi; j2lo = j2lo + INCX ) { j2hi = j2lo + INCX - 1; if ( n < j2hi ) { j2hi = n; } if ( jhi < j2hi ) { 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, " %7d ", j - 1 ); } fprintf ( stdout, "\n" ); fprintf ( stdout, " Row\n" ); fprintf ( stdout, "\n" ); /* Determine the range of the rows in this strip. */ if ( 1 < ilo ) { i2lo = ilo; } else { i2lo = 1; } if ( m < ihi ) { i2hi = m; } else { i2hi = ihi; } for ( i = i2lo; i <= i2hi; i++ ) { /* Print out (up to) 5 entries in row I, that lie in the current strip. */ fprintf ( stdout, "%5d:", i - 1 ); for ( j = j2lo; j <= j2hi; j++ ) { fprintf ( stdout, " %14f", a[i-1+(j-1)*m] ); } fprintf ( stdout, "\n" ); } } return; # undef INCX } /******************************************************************************/ void r4mat_uniform_01 ( int m, int n, float r[] ) /******************************************************************************/ /* Purpose: r4mat_uniform_01() returns a unit pseudorandom R4MAT. Licensing: This code is distributed under the MIT license. Modified: 12 January 2007 Author: John Burkardt Input: int M, N, the number of rows and columns. Output: float R[M*N], a matrix of pseudorandom values. */ { int i; int j; for ( j = 0; j < n; j++ ) { for ( i = 0; i < m; i++ ) { r[i+j*m] = drand48 ( ); } } return; } /******************************************************************************/ float *r4mat_uniform_01_new ( int m, int n ) /******************************************************************************/ /* Purpose: r4mat_uniform_01_new() returns a unit pseudorandom R4MAT. Licensing: This code is distributed under the MIT license. Modified: 12 January 2007 Author: John Burkardt Input: int M, N, the number of rows and columns. Output: float R4MAT_UNIFORM_01_NEW[M*N], a matrix of pseudorandom values. */ { int i; int j; float *r; r = ( float * ) malloc ( m * n * sizeof ( float ) ); for ( j = 0; j < n; j++ ) { for ( i = 0; i < m; i++ ) { r[i+j*m] = drand48 ( ); } } return r; } /******************************************************************************/ void r4mat_uniform_ab ( int m, int n, float b, float c, float r[] ) /******************************************************************************/ /* Purpose: r4mat_uniform_ab() returns a scaled pseudorandom R4MAT. Licensing: This code is distributed under the MIT license. Modified: 23 April 2008 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. Input: int M, N, the number of rows and columns. float B, C, the limits of the pseudorandom values. Output: float R[M*N], a matrix of pseudorandom values. */ { int i; int j; for ( j = 0; j < n; j++ ) { for ( i = 0; i < m; i++ ) { r[i+j*m] = b + ( c - b ) * drand48 ( ); } } return; } /******************************************************************************/ float *r4mat_uniform_ab_new ( int m, int n, float b, float c ) /******************************************************************************/ /* Purpose: r4mat_uniform_ab_new() returns a scaled pseudorandom R4MAT. Licensing: This code is distributed under the MIT license. Modified: 23 April 2008 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. Input: int M, N, the number of rows and columns. float B, C, the limits of the pseudorandom values. Output: float R4MAT_UNIFORM_AB_NEW[M*N], a matrix of pseudorandom values. */ { int i; int j; float *r; r = ( float * ) malloc ( m * n * sizeof ( float ) ); for ( j = 0; j < n; j++ ) { for ( i = 0; i < m; i++ ) { r[i+j*m] = b + ( c - b ) * drand48 ( ); } } return r; } /******************************************************************************/ void r4vec_print ( int n, float a[], char *title ) /******************************************************************************/ /* Purpose: r4vec_print() prints an R4VEC. Discussion: An R4VEC is a vector of R4's. Licensing: This code is distributed under the MIT license. Modified: 08 April 2009 Author: John Burkardt Input: int N, the number of components of the vector. float A[N], the vector to be printed. 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, " %8d: %14f\n", i, a[i] ); } return; } /******************************************************************************/ void r4vec_uniform_ab ( int n, float b, float c, float r[] ) /******************************************************************************/ /* Purpose: r4vec_uniform_ab() returns a scaled pseudorandom R4VEC. Licensing: This code is distributed under the MIT license. Modified: 23 April 2008 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. Input: int N, the number of entries in the vector. float B, C, the lower and upper limits of the pseudorandom values. Output: float R[N], the vector of pseudorandom values. */ { int i; for ( i = 0; i < n; i++ ) { r[i] = b + ( c - b ) * drand48 ( ); } return; } /******************************************************************************/ float *r4vec_uniform_ab_new ( int n, float b, float c ) /******************************************************************************/ /* Purpose: r4vec_uniform_ab_new() returns a scaled pseudorandom R4VEC. Licensing: This code is distributed under the MIT license. Modified: 23 April 2008 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. Input: int N, the number of entries in the vector. float B, C, the lower and upper limits of the pseudorandom values. Output: float R4VEC_UNIFORM_AB_NEW[N], the vector of pseudorandom values. */ { int i; float *r; r = ( float * ) malloc ( n * sizeof ( float ) ); for ( i = 0; i < n; i++ ) { r[i] = b + ( c - b ) * drand48 ( ); } return r; } /******************************************************************************/ void r4vec_uniform_01 ( int n, float r[] ) /******************************************************************************/ /* Purpose: r4vec_uniform_01() returns a unit pseudorandom R4VEC. Licensing: This code is distributed under the MIT license. Modified: 12 January 2007 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. Input: int N, the number of entries in the vector. Output: float R[N], the vector of pseudorandom values. */ { int i; for ( i = 0; i < n; i++ ) { r[i] = drand48 ( ); } return; } /******************************************************************************/ float *r4vec_uniform_01_new ( int n ) /******************************************************************************/ /* Purpose: r4vec_uniform_01_new() returns a unit pseudorandom R4VEC. Licensing: This code is distributed under the MIT license. Modified: 12 January 2007 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. Input: int N, the number of entries in the vector. Output: float R4VEC_UNIFORM_01_NEW[N], the vector of pseudorandom values. */ { int i; float *r; r = ( float * ) malloc ( n * sizeof ( float ) ); for ( i = 0; i < n; i++ ) { r[i] = drand48 ( ); } return r; } /******************************************************************************/ int r8_nint ( double x ) /******************************************************************************/ /* Purpose: r8_nint() returns the nearest I4 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: 26 August 2004 Author: John Burkardt Input: double X, the value. int R8_NINT, the nearest integer to X. */ { int s; if ( x < 0.0 ) { s = -1; } else { s = 1; } return ( s * ( int ) ( fabs ( x ) + 0.5 ) ); } /******************************************************************************/ double r8_uniform_ab ( double a, double b ) /******************************************************************************/ /* Purpose: r8_uniform_ab() returns a scaled pseudorandom R8. Discussion: The pseudorandom number should be uniformly distributed between A and B. Licensing: This code is distributed under the MIT license. Modified: 19 April 2011 Author: John Burkardt Input: double A, B, the limits of the interval. Output: double R8_UNIFORM_AB, a number strictly between A and B. */ { double value; value = a + ( b - a ) * drand48 ( ); return value; } /******************************************************************************/ double r8_uniform_01 ( ) /******************************************************************************/ /* Purpose: r8_uniform_01() returns a unit pseudorandom R8. Licensing: This code is distributed under the MIT license. Modified: 11 August 2004 Author: John Burkardt Output: double R8_UNIFORM_01, a new pseudorandom variate, strictly between 0 and 1. */ { double r; r = drand48 ( ); return r; } /******************************************************************************/ double *r8col_uniform_abvec_new ( int m, int n, double a[], double b[] ) /******************************************************************************/ /* Purpose: r8col_uniform_abvec_new() fills an R8COL with scaled pseudorandom numbers. Discussion: An R8COL is an array of R8 values, regarded as a set of column vectors. The user specifies a minimum and maximum value for each row. Licensing: This code is distributed under the MIT license. Modified: 28 December 2014 Author: John Burkardt Reference: Paul Bratley, Bennett Fox, Linus Schrage, A Guide to Simulation, Springer Verlag, pages 201-202, 1983. Bennett Fox, Algorithm 647: Implementation and Relative Efficiency of Quasirandom Sequence Generators, ACM Transactions on Mathematical Software, Volume 12, Number 4, pages 362-376, 1986. Philip Lewis, Allen Goodman, James Miller, A Pseudo-Random Number Generator for the System/360, IBM Systems Journal, Volume 8, pages 136-143, 1969. Input: int M, N, the number of rows and columns. double A[M], B[M], the upper and lower limits. Output: double R8COL_UNIFORM_ABVEC_NEW[M*N], a matrix of pseudorandom values. */ { int i; int j; double *r; r = ( double * ) malloc ( m * n * sizeof ( double ) ); for ( j = 0; j < n; j++ ) { for ( i = 0; i < m; i++ ) { r[i+j*m] = a[i] + ( b[i] - a[i] ) * drand48 ( ); } } return r; } /******************************************************************************/ void r8mat_print ( int m, int n, double a[], char *title ) /******************************************************************************/ /* Purpose: r8mat_print() prints an R8MAT. Discussion: An R8MAT is a doubly dimensioned array of R8 values, stored as a vector in column-major order. Entry A(I,J) is stored as A[I+J*M] Licensing: This code is distributed under the MIT license. Modified: 28 May 2008 Author: John Burkardt Input: int M, the number of rows in A. int N, the number of columns in A. double A[M*N], the M by N matrix. char *TITLE, a title. */ { r8mat_print_some ( m, n, a, 1, 1, m, n, title ); return; } /******************************************************************************/ void r8mat_print_some ( int m, int n, double a[], int ilo, int jlo, int ihi, int jhi, char *title ) /******************************************************************************/ /* Purpose: r8mat_print_some() prints some of an R8MAT. Discussion: An R8MAT is a doubly dimensioned array of R8 values, stored as a vector in column-major order. Licensing: This code is distributed under the MIT license. Modified: 26 June 2013 Author: John Burkardt Input: int M, the number of rows of the matrix. M must be positive. int N, the number of columns of the matrix. N must be positive. double A[M*N], the matrix. int ILO, JLO, IHI, JHI, designate the first row and column, and the last row and column to be printed. char *TITLE, a title. */ { # define INCX 5 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 5. */ for ( j2lo = jlo; j2lo <= jhi; j2lo = j2lo + INCX ) { j2hi = j2lo + INCX - 1; if ( n < j2hi ) { j2hi = n; } if ( jhi < j2hi ) { 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, " %7d ", j - 1 ); } fprintf ( stdout, "\n" ); fprintf ( stdout, " Row\n" ); fprintf ( stdout, "\n" ); /* Determine the range of the rows in this strip. */ if ( 1 < ilo ) { i2lo = ilo; } else { i2lo = 1; } if ( m < ihi ) { i2hi = m; } else { i2hi = ihi; } for ( i = i2lo; i <= i2hi; i++ ) { /* Print out (up to) 5 entries in row I, that lie in the current strip. */ fprintf ( stdout, "%5d:", i - 1 ); for ( j = j2lo; j <= j2hi; j++ ) { fprintf ( stdout, " %14g", a[i-1+(j-1)*m] ); } fprintf ( stdout, "\n" ); } } return; # undef INCX } /******************************************************************************/ void r8mat_uniform_01 ( int m, int n, double r[] ) /******************************************************************************/ /* Purpose: r8mat_uniform_01() returns a unit pseudorandom R8MAT. Licensing: This code is distributed under the MIT license. Modified: 03 October 2005 Author: John Burkardt Input: int M, N, the number of rows and columns. Output: double R[M*N], a matrix of pseudorandom values. */ { int i; int j; for ( j = 0; j < n; j++ ) { for ( i = 0; i < m; i++ ) { r[i+j*m] = drand48 ( ); } } return; } /******************************************************************************/ double *r8mat_uniform_01_new ( int m, int n ) /******************************************************************************/ /* Purpose: r8mat_uniform_01_new() returns a unit pseudorandom R8MAT. Licensing: This code is distributed under the MIT license. Modified: 03 October 2005 Author: John Burkardt Input: int M, N, the number of rows and columns. Output: double R8MAT_UNIFORM_01_NEW[M*N], a matrix of pseudorandom values. */ { int i; int j; double *r; r = ( double * ) malloc ( m * n * sizeof ( double ) ); for ( j = 0; j < n; j++ ) { for ( i = 0; i < m; i++ ) { r[i+j*m] = drand48 ( ); } } return r; } /******************************************************************************/ void r8mat_uniform_ab ( int m, int n, double a, double b, double r[] ) /******************************************************************************/ /* Purpose: r8mat_uniform_ab() returns a scaled pseudorandom R8MAT. Licensing: This code is distributed under the MIT license. Modified: 03 October 2005 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. Input: int M, N, the number of rows and columns. double A, B, the limits of the pseudorandom values. Output: double R[M*N], a matrix of pseudorandom values. */ { int i; int j; for ( j = 0; j < n; j++ ) { for ( i = 0; i < m; i++ ) { r[i+j*m] = a + ( b - a ) * drand48 ( ); } } return; } /******************************************************************************/ double *r8mat_uniform_ab_new ( int m, int n, double a, double b ) /******************************************************************************/ /* Purpose: r8mat_uniform_ab_new() returns a scaled pseudorandom R8MAT. Licensing: This code is distributed under the MIT license. Modified: 03 October 2005 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. Input: int M, N, the number of rows and columns. double A, B, the limits of the pseudorandom values. Output: double R8MAT_UNIFORM_AB_NEW[M*N], a matrix of pseudorandom values. */ { int i; int j; double *r; r = ( double * ) malloc ( m * n * sizeof ( double ) ); for ( j = 0; j < n; j++ ) { for ( i = 0; i < m; i++ ) { r[i+j*m] = a + ( b - a ) * drand48 ( ); } } return r; } /******************************************************************************/ double *r8row_uniform_abvec_new ( int m, int n, double a[], double b[] ) /******************************************************************************/ /* Purpose: r8row_uniform_abvec_new() fills an R8ROW with scaled pseudorandom numbers. Discussion: An R8ROW is an array of R8 values, regarded as a set of row vectors. The user specifies a minimum and maximum value for each column. Licensing: This code is distributed under the MIT license. Modified: 15 January 2012 Author: John Burkardt Reference: Paul Bratley, Bennett Fox, Linus Schrage, A Guide to Simulation, Springer Verlag, pages 201-202, 1983. Bennett Fox, Algorithm 647: Implementation and Relative Efficiency of Quasirandom Sequence Generators, ACM Transactions on Mathematical Software, Volume 12, Number 4, pages 362-376, 1986. Philip Lewis, Allen Goodman, James Miller, A Pseudo-Random Number Generator for the System/360, IBM Systems Journal, Volume 8, pages 136-143, 1969. Input: int M, N, the number of rows and columns. double A[N], B[N], the upper and lower limits. Output: double R8ROW_UNIFORM_ABVEC_NEW[M*N], a matrix of pseudorandom values. */ { int i; int j; double *r; r = ( double * ) malloc ( m * n * sizeof ( double ) ); for ( i = 0; i < m; i++ ) { for ( j = 0; j < n; j++ ) { r[i+j*m] = a[j] + ( b[j] - a[j] ) * drand48 ( ); } } return r; } /******************************************************************************/ void r8vec_copy ( int n, double a1[], double a2[] ) /******************************************************************************/ /* Purpose: r8vec_copy() copies an R8VEC. Licensing: This code is distributed under the MIT license. Modified: 03 July 2005 Author: John Burkardt Input: int N, the number of entries in the vectors. double A1[N], the vector to be copied. double A2[N], the copy of A1. */ { int i; for ( i = 0; i < n; i++ ) { a2[i] = a1[i]; } return; } /******************************************************************************/ double *r8vec_normal_01_new ( int n ) /******************************************************************************/ /* Purpose: r8vec_normal_01_new() returns a unit pseudonormal R8VEC. Discussion: The standard normal probability distribution function (PDF) has mean 0 and standard deviation 1. This routine can generate a vector of values on one call. It has the feature that it should provide the same results in the same order no matter how we break up the task. The Box-Muller method is used, which is efficient, but generates an even number of values each time. On any call to this routine, an even number of new values are generated. Depending on the situation, one value may be left over. In that case, it is saved for the next call. Licensing: This code is distributed under the MIT license. Modified: 18 October 2004 Author: John Burkardt Input: int N, the number of values desired. If N is negative, then the code will flush its internal memory; in particular, if there is a saved value to be used on the next call, it is instead discarded. Output: double R8VEC_NORMAL_01_NEW[N], a sample of the standard normal PDF. Local: int MADE, records the number of values that have been computed. On input with negative N, this value overwrites the return value of N, so the user can get an accounting of how much work has been done. double R[N+1], is used to store some uniform random values. Its dimension is N+1, but really it is only needed to be the smallest even number greater than or equal to N. int SAVED, is 0 or 1 depending on whether there is a single saved value left over from the previous call. int X_LO, X_HI, records the range of entries of X that we need to compute. This starts off as 1:N, but is adjusted if we have a saved value that can be immediately stored in X(1), and so on. float Y, the value saved from the previous call, if SAVED is 1. */ { # define R8_PI 3.141592653589793 int i; int m; static int made = 0; double *r; static int saved = 0; double *x; int x_hi; int x_lo; static double y = 0.0; x = ( double * ) malloc ( n * sizeof ( double ) ); /* I'd like to allow the user to reset the internal data. But this won't work properly if we have a saved value Y. I'm making a crock option that allows the user to signal explicitly that any internal memory should be flushed, by passing in a negative value for N. */ if ( n < 0 ) { made = 0; saved = 0; y = 0.0; return NULL; } else if ( n == 0 ) { return NULL; } /* Record the range of X we need to fill in. */ x_lo = 1; x_hi = n; /* Use up the old value, if we have it. */ if ( saved == 1 ) { x[0] = y; saved = 0; x_lo = 2; } /* Maybe we don't need any more values. */ if ( x_hi - x_lo + 1 == 0 ) { } /* If we need just one new value, do that here to avoid null arrays. */ else if ( x_hi - x_lo + 1 == 1 ) { r = r8vec_uniform_01_new ( 2 ); x[x_hi-1] = sqrt ( - 2.0 * log ( r[0] ) ) * cos ( 2.0 * R8_PI * r[1] ); y = sqrt ( - 2.0 * log ( r[0] ) ) * sin ( 2.0 * R8_PI * r[1] ); saved = 1; made = made + 2; free ( r ); } /* If we require an even number of values, that's easy. */ else if ( ( x_hi - x_lo + 1 ) % 2 == 0 ) { m = ( x_hi - x_lo + 1 ) / 2; r = r8vec_uniform_01_new ( 2 * m ); for ( i = 0; i <= 2*m-2; i = i + 2 ) { x[x_lo+i-1] = sqrt ( - 2.0 * log ( r[i] ) ) * cos ( 2.0 * R8_PI * r[i+1] ); x[x_lo+i ] = sqrt ( - 2.0 * log ( r[i] ) ) * sin ( 2.0 * R8_PI * r[i+1] ); } made = made + x_hi - x_lo + 1; free ( r ); } /* If we require an odd number of values, we generate an even number, and handle the last pair specially, storing one in X(N), and saving the other for later. */ else { x_hi = x_hi - 1; m = ( x_hi - x_lo + 1 ) / 2 + 1; r = r8vec_uniform_01_new ( 2 * m ); for ( i = 0; i <= 2*m-4; i = i + 2 ) { x[x_lo+i-1] = sqrt ( - 2.0 * log ( r[i] ) ) * cos ( 2.0 * R8_PI * r[i+1] ); x[x_lo+i ] = sqrt ( - 2.0 * log ( r[i] ) ) * sin ( 2.0 * R8_PI * r[i+1] ); } i = 2*m - 2; x[x_lo+i-1] = sqrt ( - 2.0 * log ( r[i] ) ) * cos ( 2.0 * R8_PI * r[i+1] ); y = sqrt ( - 2.0 * log ( r[i] ) ) * sin ( 2.0 * R8_PI * r[i+1] ); saved = 1; made = made + x_hi - x_lo + 2; free ( r ); } return x; # undef R8_PI } /******************************************************************************/ void r8vec_print ( int n, double a[], char *title ) /******************************************************************************/ /* Purpose: r8vec_print() prints an R8VEC. Discussion: An R8VEC is a vector of R8's. Licensing: This code is distributed under the MIT license. Modified: 08 April 2009 Author: John Burkardt Input: int N, the number of components of the vector. double A[N], the vector to be printed. 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, " %8d: %14g\n", i, a[i] ); } return; } /******************************************************************************/ void r8vec_uniform_01 ( int n, double r[] ) /******************************************************************************/ /* Purpose: r8vec_uniform_01() returns a unit pseudorandom R8VEC. Licensing: This code is distributed under the MIT license. Modified: 19 August 2004 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. Input: int N, the number of entries in the vector. Output: double R[N], the vector of pseudorandom values. */ { int i; for ( i = 0; i < n; i++ ) { r[i] = drand48 ( ); } return; } /******************************************************************************/ double *r8vec_uniform_01_new ( int n ) /******************************************************************************/ /* Purpose: r8vec_uniform_01_new() returns a unit pseudorandom R8VEC. Licensing: This code is distributed under the MIT license. Modified: 19 August 2004 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. Input: int N, the number of entries in the vector. Output: double R8VEC_UNIFORM_01_NEW[N], the vector of pseudorandom values. */ { int i; double *r; r = ( double * ) malloc ( n * sizeof ( double ) ); for ( i = 0; i < n; i++ ) { r[i] = drand48 ( ); } return r; } /******************************************************************************/ void r8vec_uniform_ab ( int n, double a, double b, double r[] ) /******************************************************************************/ /* Purpose: r8vec_uniform_ab() returns a scaled pseudorandom R8VEC. Discussion: Each dimension ranges from A to B. Licensing: This code is distributed under the MIT license. Modified: 30 January 2005 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. Input: int N, the number of entries in the vector. double A, B, the lower and upper limits of the pseudorandom values. Output: double R[N], the vector of pseudorandom values. */ { int i; for ( i = 0; i < n; i++ ) { r[i] = a + ( b - a ) * drand48 ( ); } return; } /******************************************************************************/ double *r8vec_uniform_ab_new ( int n, double a, double b ) /******************************************************************************/ /* Purpose: r8vec_uniform_ab_new() returns a scaled pseudorandom R8VEC. Discussion: Each dimension ranges from A to B. Licensing: This code is distributed under the MIT license. Modified: 30 January 2005 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. Input: int N, the number of entries in the vector. double A, B, the lower and upper limits of the pseudorandom values. Output: double R8VEC_UNIFORM_AB_NEW[N], the vector of pseudorandom values. */ { int i; double *r; r = ( double * ) malloc ( n * sizeof ( double ) ); for ( i = 0; i < n; i++ ) { r[i] = a + ( b - a ) * drand48 ( ); } return r; } /******************************************************************************/ void r8vec_uniform_abvec ( int n, double a[], double b[], double r[] ) /******************************************************************************/ /* Purpose: r8vec_uniform_abvec() returns a scaled pseudorandom R8VEC. Discussion: Dimension I ranges from A[I] to B[I]. Licensing: This code is distributed under the MIT license. Modified: 30 January 2005 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. Input: int N, the number of entries in the vector. double A[N], B[N], the lower and upper limits of the pseudorandom values. Output: double R[N], the vector of pseudorandom values. */ { int i; for ( i = 0; i < n; i++ ) { r[i] = a[i] + ( b[i] - a[i] ) * drand48 ( ); } return; } /******************************************************************************/ double *r8vec_uniform_abvec_new ( int n, double a[], double b[] ) /******************************************************************************/ /* Purpose: r8vec_uniform_abvec_new() returns a scaled pseudorandom R8VEC. Discussion: Dimension I ranges from A[I] to B[I]. Licensing: This code is distributed under the MIT license. Modified: 30 January 2005 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. Input: int N, the number of entries in the vector. double A[N], B[N], the lower and upper limits of the pseudorandom values. Output: double R8VEC_UNIFORM_ABVEC_NEW[N], the vector of pseudorandom values. */ { int i; double *r; r = ( double * ) malloc ( n * sizeof ( double ) ); for ( i = 0; i < n; i++ ) { r[i] = a[i] + ( b[i] - a[i] ) * drand48 ( ); } return r; } /******************************************************************************/ double *r8vec_uniform_unit_new ( int m ) /******************************************************************************/ /* Purpose: r8vec_uniform_unit_new() generates a random unit vector. Licensing: This code is distributed under the MIT license. Modified: 04 October 2012 Author: John Burkardt Input: int M, the dimension of the space. Output: double R8VEC_UNIFORM_UNIT_NEW[M], a random direction vector, with unit norm. */ { double *a; int i; double norm; /* Take M random samples from the normal distribution. */ a = r8vec_normal_01_new ( m ); /* Compute the norm. */ norm = 0.0; for ( i = 0; i < m; i++ ) { norm = norm + a[i] * a[i]; } norm = sqrt ( norm ); /* Normalize. */ for ( i = 0; i < m; i++ ) { a[i] = a[i] / norm; } return a; }