# include # include # include # include # include # include "normal.h" /******************************************************************************/ float complex c4_normal_01 ( int *seed ) /******************************************************************************/ /* Purpose: C4_NORMAL_01 returns a unit pseudonormal C4. Licensing: This code is distributed under the GNU LGPL license. Modified: 03 March 2015 Author: John Burkardt Parameters: Input/output, int *SEED, a seed for the random number generator. Output, float complex C4_NORMAL_01, a unit pseudonormal value. */ { const float r4_pi = 3.141592653589793; float complex value; float v1; float v2; float x_c; float x_r; v1 = r4_uniform_01 ( seed ); v2 = r4_uniform_01 ( seed ); x_r = sqrt ( - 2.0 * log ( v1 ) ) * cos ( 2.0 * r4_pi * v2 ); x_c = sqrt ( - 2.0 * log ( v1 ) ) * sin ( 2.0 * r4_pi * v2 ); value = x_r + x_c * I; return value; } /******************************************************************************/ double complex c8_normal_01 ( int *seed ) /******************************************************************************/ /* Purpose: C8_NORMAL_01 returns a unit pseudonormal C8. Licensing: This code is distributed under the GNU LGPL license. Modified: 03 March 2015 Author: John Burkardt Parameters: Input/output, int *SEED, a seed for the random number generator. Output, double complex C8_NORMAL_01, a unit pseudonormal value. */ { const double r4_pi = 3.141592653589793; double complex value; double v1; double v2; double x_c; double x_r; v1 = r8_uniform_01 ( seed ); v2 = r8_uniform_01 ( seed ); x_r = sqrt ( - 2.0 * log ( v1 ) ) * cos ( 2.0 * r4_pi * v2 ); x_c = sqrt ( - 2.0 * log ( v1 ) ) * sin ( 2.0 * r4_pi * v2 ); value = x_r + x_c * I; return value; } /******************************************************************************/ int i4_normal_ab ( double a, double b, int *seed ) /******************************************************************************/ /* Purpose: I4_NORMAL_AB returns a scaled pseudonormal I4. Discussion: The normal probability distribution function (PDF) is sampled, with mean A and standard deviation B. Licensing: This code is distributed under the GNU LGPL license. Modified: 02 July 2006 Author: John Burkardt Parameters: Input, double A, the mean of the PDF. Input, double B, the standard deviation of the PDF. Input/output, int *SEED, a seed for the random number generator. Output, int I4_NORMAL_AB, a sample of the normal PDF. */ { double r1; double r2; const double r8_pi = 3.141592653589793; int value; double x; r1 = r8_uniform_01 ( seed ); r2 = r8_uniform_01 ( seed ); x = sqrt ( - 2.0 * log ( r1 ) ) * cos ( 2.0 * r8_pi * r2 ); value = ( int ) round ( a + b * x ); return value; } /**********************************************************************/ void i4_normal_ab_test ( ) /**********************************************************************/ /* Purpose: I4_NORMAL_AB_TEST tests I4_NORMAL_AB. Licensing: This code is distributed under the GNU LGPL license. Modified: 02 July 2006 Author: John Burkardt */ { int i; double mu; int r; int seed; double sigma; printf ( "\n" ); printf ( "I4_NORMAL_AB_TEST\n" ); printf ( " I4_NORMAL_AB computes pseudonormal integers\n" ); printf ( " with mean MU and standard deviation SIGMA.\n" ); mu = 70.0; sigma = 10.0; seed = 123456789; printf ( "\n" ); printf ( " The mean = %f\n", mu ); printf ( " The standard deviation = %f\n", sigma ); printf ( " SEED = %d\n", seed ); printf ( "\n" ); for ( i = 1; i <= 10; i++ ) { r = i4_normal_ab ( mu, sigma, &seed ); printf ( " %8d %8d\n", i, r ); } return; } /******************************************************************************/ long long int i8_normal_ab ( double a, double b, long long int *seed ) /******************************************************************************/ /* Purpose: I8_NORMAL_AB returns a scaled pseudonormal I8. Discussion: The normal probability distribution function (PDF) is sampled, with mean A and standard deviation B. Licensing: This code is distributed under the GNU LGPL license. Modified: 02 July 2006 Author: John Burkardt Parameters: Input, double A, the mean of the PDF. Input, double B, the standard deviation of the PDF. Input/output, long long int *SEED, a seed for the random number generator. Output, long long int I8_NORMAL_AB, a sample of the normal PDF. */ { int seed_int; double value_double; long long int value_long_int; seed_int = ( int ) *seed; value_double = a + b * r8_normal_01 ( &seed_int ); if ( value_double < 0.0 ) { value_long_int = ( long long int ) ( value_double - 0.5 ); } else { value_long_int = ( long long int ) ( value_double + 0.5 ); } *seed = ( long long int ) seed_int; return value_long_int; } /******************************************************************************/ float r4_normal_01 ( int *seed ) /******************************************************************************/ /* Purpose: R4_NORMAL_01 returns a unit pseudonormal R4. Discussion: The standard normal probability distribution function (PDF) has mean 0 and standard deviation 1. The Box-Muller method is used, which is efficient, but generates two values at a time. Licensing: This code is distributed under the GNU LGPL license. Modified: 05 June 2013 Author: John Burkardt Parameters: Input/output, int *SEED, a seed for the random number generator. Output, float R4_NORMAL_01, a normally distributed random value. */ { float r1; float r2; const double r4_pi = 3.141592653589793; float x; r1 = r4_uniform_01 ( seed ); r2 = r4_uniform_01 ( seed ); x = sqrt ( - 2.0 * log ( r1 ) ) * cos ( 2.0 * r4_pi * r2 ); return x; } /******************************************************************************/ float r4_normal_ab ( float a, float b, int *seed ) /******************************************************************************/ /* Purpose: R4_NORMAL_AB returns a scaled pseudonormal R4. Discussion: The normal probability distribution function (PDF) is sampled, with mean A and standard deviation B. Licensing: This code is distributed under the GNU LGPL license. Modified: 29 June 2006 Author: John Burkardt Parameters: Input, float A, the mean of the PDF. Input, float B, the standard deviation of the PDF. Input/output, int *SEED, a seed for the random number generator. Output, float R4_NORMAL_AB, a sample of the normal PDF. */ { float value; value = a + b * r4_normal_01 ( seed ); return value; } /******************************************************************************/ float r4_uniform_01 ( int *seed ) /******************************************************************************/ /* Purpose: R4_UNIFORM_01 returns a unit pseudorandom R4. Discussion: This routine implements the recursion seed = 16807 * seed mod ( 2^31 - 1 ) r4_uniform_01 = seed / ( 2^31 - 1 ) The integer arithmetic never requires more than 32 bits, including a sign bit. If the initial seed is 12345, then the first three computations are Input Output R4_UNIFORM_01 SEED SEED 12345 207482415 0.096616 207482415 1790989824 0.833995 1790989824 2035175616 0.947702 Licensing: This code is distributed under the GNU LGPL license. Modified: 16 November 2004 Author: John Burkardt Reference: Paul Bratley, Bennett Fox, Linus Schrage, A Guide to Simulation, Springer Verlag, pages 201-202, 1983. Pierre L'Ecuyer, Random Number Generation, in Handbook of Simulation edited by Jerry Banks, Wiley Interscience, page 95, 1998. 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. Peter Lewis, Allen Goodman, James Miller, A Pseudo-Random Number Generator for the System/360, IBM Systems Journal, Volume 8, pages 136-143, 1969. Parameters: Input/output, int *SEED, the "seed" value. Normally, this value should not be 0. On output, SEED has been updated. Output, float R4_UNIFORM_01, a new pseudorandom variate, strictly between 0 and 1. */ { const int i4_huge = 2147483647; int k; float r; k = *seed / 127773; *seed = 16807 * ( *seed - k * 127773 ) - k * 2836; if ( *seed < 0 ) { *seed = *seed + i4_huge; } /* Although SEED can be represented exactly as a 32 bit integer, it generally cannot be represented exactly as a 32 bit real number! */ r = ( float ) ( *seed ) * 4.656612875E-10; return r; } /******************************************************************************/ 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 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 GNU LGPL license. Modified: 03 March 2015 Author: John Burkardt Parameters: Input, int M, the number of rows in A. Input, int N, the number of columns in A. Input, float A[M*N], the M by N matrix. Input, 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 values, stored as a vector in column-major order. Licensing: This code is distributed under the GNU LGPL license. Modified: 03 March 2015 Author: John Burkardt Parameters: Input, int M, the number of rows of the matrix. M must be positive. Input, int N, the number of columns of the matrix. N must be positive. Input, float A[M*N], the matrix. Input, int ILO, JLO, IHI, JHI, designate the first row and column, and the last row and column to be printed. Input, char *TITLE, a title. */ { # define INCX 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 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 GNU LGPL license. Modified: 03 March 2015 Author: John Burkardt Parameters: Input, int N, the number of components of the vector. Input, float A[N], the vector to be printed. Input, char *TITLE, a title. */ { int i; fprintf ( stdout, "\n" ); fprintf ( stdout, "%s\n", title ); fprintf ( stdout, "\n" ); for ( i = 0; i < n; i++ ) { fprintf ( stdout, " %8d: %14g\n", i, a[i] ); } return; } /******************************************************************************/ float *r4vec_uniform_01_new ( int n, int *seed ) /******************************************************************************/ /* Purpose: R4VEC_UNIFORM_01_NEW returns a unit pseudorandom R4VEC. Discussion: This routine implements the recursion seed = 16807 * seed mod ( 2^31 - 1 ) unif = seed / ( 2^31 - 1 ) The integer arithmetic never requires more than 32 bits, including a sign bit. Licensing: This code is distributed under the GNU LGPL license. Modified: 03 March 2015 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. Peter Lewis, Allen Goodman, James Miller, A Pseudo-Random Number Generator for the System/360, IBM Systems Journal, Volume 8, pages 136-143, 1969. Parameters: Input, int N, the number of entries in the vector. Input/output, int *SEED, a seed for the random number generator. Output, float R4VEC_UNIFORM_01_NEW[N], the vector of pseudorandom values. */ { int i; const int i4_huge = 2147483647; int k; float *r; r = ( float * ) malloc ( n * sizeof ( float ) ); for ( i = 0; i < n; i++ ) { k = *seed / 127773; *seed = 16807 * ( *seed - k * 127773 ) - k * 2836; if ( *seed < 0 ) { *seed = *seed + i4_huge; } r[i] = ( float ) ( *seed ) * 4.656612875E-10; } return r; } /******************************************************************************/ double r8_normal_01 ( int *seed ) /******************************************************************************/ /* Purpose: R8_NORMAL_01 returns a unit pseudonormal R8. Discussion: The standard normal probability distribution function (PDF) has mean 0 and standard deviation 1. Because this routine uses the Box Muller method, it requires pairs of uniform random values to generate a pair of normal random values. This means that on every other call, the code can use the second value that it calculated. However, if the user has changed the SEED value between calls, the routine automatically resets itself and discards the saved data. Licensing: This code is distributed under the GNU LGPL license. Modified: 06 August 2013 Author: John Burkardt Parameters: Input/output, int *SEED, a seed for the random number generator. Output, double R8_NORMAL_01, a normally distributed random value. */ { double r1; double r2; const double r8_pi = 3.141592653589793; double x; r1 = r8_uniform_01 ( seed ); r2 = r8_uniform_01 ( seed ); x = sqrt ( - 2.0 * log ( r1 ) ) * cos ( 2.0 * r8_pi * r2 ); return x; } /******************************************************************************/ double r8_normal_ab ( double a, double b, int *seed ) /******************************************************************************/ /* Purpose: R8_NORMAL_AB returns a scaled pseudonormal R8. Discussion: The normal probability distribution function (PDF) is sampled, with mean A and standard deviation B. Licensing: This code is distributed under the GNU LGPL license. Modified: 21 November 2004 Author: John Burkardt Parameters: Input, double A, the mean of the PDF. Input, double B, the standard deviation of the PDF. Input/output, int *SEED, a seed for the random number generator. Output, double R8_NORMAL_AB, a sample of the normal PDF. */ { double value; value = a + b * r8_normal_01 ( seed ); return value; } /******************************************************************************/ double r8_uniform_01 ( int *seed ) /******************************************************************************/ /* Purpose: R8_UNIFORM_01 returns a unit pseudorandom R8. Discussion: This routine implements the recursion seed = 16807 * seed mod ( 2^31 - 1 ) r8_uniform_01 = seed / ( 2^31 - 1 ) The integer arithmetic never requires more than 32 bits, including a sign bit. If the initial seed is 12345, then the first three computations are Input Output R8_UNIFORM_01 SEED SEED 12345 207482415 0.096616 207482415 1790989824 0.833995 1790989824 2035175616 0.947702 Licensing: This code is distributed under the GNU LGPL license. Modified: 11 August 2004 Author: John Burkardt Reference: Paul Bratley, Bennett Fox, Linus Schrage, A Guide to Simulation, Springer Verlag, pages 201-202, 1983. Pierre L'Ecuyer, Random Number Generation, in Handbook of Simulation edited by Jerry Banks, Wiley Interscience, page 95, 1998. 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. Peter Lewis, Allen Goodman, James Miller, A Pseudo-Random Number Generator for the System/360, IBM Systems Journal, Volume 8, pages 136-143, 1969. Parameters: Input/output, int *SEED, the "seed" value. Normally, this value should not be 0. On output, SEED has been updated. Output, double R8_UNIFORM_01, a new pseudorandom variate, strictly between 0 and 1. */ { const int i4_huge = 2147483647; int k; double r; k = *seed / 127773; *seed = 16807 * ( *seed - k * 127773 ) - k * 2836; if ( *seed < 0 ) { *seed = *seed + i4_huge; } /* Although SEED can be represented exactly as a 32 bit integer, it generally cannot be represented exactly as a 32 bit real number! */ r = ( double ) ( *seed ) * 4.656612875E-10; return r; } /******************************************************************************/ void r8mat_normal_01 ( int m, int n, int *seed, double r[] ) /******************************************************************************/ /* Purpose: R8MAT_NORMAL_01 returns a unit pseudonormal R8MAT. Licensing: This code is distributed under the GNU LGPL license. Modified: 03 October 2005 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. Peter Lewis, Allen Goodman, James Miller, A Pseudo-Random Number Generator for the System/360, IBM Systems Journal, Volume 8, pages 136-143, 1969. Parameters: Input, int M, N, the number of rows and columns in the array. Input/output, int *SEED, the "seed" value, which should NOT be 0. On output, SEED has been updated. Output, double R[M*N], the array of pseudonormal values. */ { r8vec_normal_01 ( m * n, seed, r ); return; } /******************************************************************************/ double *r8mat_normal_01_new ( int m, int n, int *seed ) /******************************************************************************/ /* Purpose: R8MAT_NORMAL_01_NEW returns a unit pseudonormal R8MAT. Licensing: This code is distributed under the GNU LGPL license. Modified: 03 October 2005 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. Peter Lewis, Allen Goodman, James Miller, A Pseudo-Random Number Generator for the System/360, IBM Systems Journal, Volume 8, pages 136-143, 1969. Parameters: Input, int M, N, the number of rows and columns in the array. Input/output, int *SEED, the "seed" value, which should NOT be 0. On output, SEED has been updated. Output, double R8MAT_NORMAL_01_NEW[M*N], the array of pseudonormal values. */ { double *r; r = r8vec_normal_01_new ( m * n, seed ); return r; } /******************************************************************************/ void r8mat_normal_ab ( int m, int n, double a, double b, int *seed, double r[] ) /******************************************************************************/ /* Purpose: R8MAT_NORMAL_AB returns a scaled pseudonormal R8MAT. Licensing: This code is distributed under the GNU LGPL license. Modified: 03 October 2005 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. Peter Lewis, Allen Goodman, James Miller, A Pseudo-Random Number Generator for the System/360, IBM Systems Journal, Volume 8, pages 136-143, 1969. Parameters: Input, int M, N, the number of rows and columns in the array. Input, double A, B, the mean and standard deviation. Input/output, int *SEED, the "seed" value, which should NOT be 0. On output, SEED has been updated. Output, double R[M*N], the array of pseudonormal values. */ { r8vec_normal_ab ( m * n, a, b, seed, r ); return; } /******************************************************************************/ double *r8mat_normal_ab_new ( int m, int n, double a, double b, int *seed ) /******************************************************************************/ /* Purpose: R8MAT_NORMAL_AB_NEW returns a scaled pseudonormal R8MAT. Licensing: This code is distributed under the GNU LGPL license. Modified: 03 October 2005 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. Peter Lewis, Allen Goodman, James Miller, A Pseudo-Random Number Generator for the System/360, IBM Systems Journal, Volume 8, pages 136-143, 1969. Parameters: Input, int M, N, the number of rows and columns in the array. Input, double A, B, the mean and standard deviation. Input/output, int *SEED, the "seed" value, which should NOT be 0. On output, SEED has been updated. Output, double R8MAT_NORMAL_AB_NEW[M*N], the array of pseudonormal values. */ { double *r; r = r8vec_normal_ab_new ( m * n, a, b, seed ); 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 GNU LGPL license. Modified: 28 May 2008 Author: John Burkardt Parameters: Input, int M, the number of rows in A. Input, int N, the number of columns in A. Input, double A[M*N], the M by N matrix. Input, 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 GNU LGPL license. Modified: 26 June 2013 Author: John Burkardt Parameters: Input, int M, the number of rows of the matrix. M must be positive. Input, int N, the number of columns of the matrix. N must be positive. Input, double A[M*N], the matrix. Input, int ILO, JLO, IHI, JHI, designate the first row and column, and the last row and column to be printed. Input, char *TITLE, a title. */ { # define INCX 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 r8vec_normal_01 ( int n, int *seed, double x[] ) /******************************************************************************/ /* Purpose: R8VEC_NORMAL_01 returns a unit pseudonormal R8VEC. Discussion: The standard normal probability distribution function (PDF) has mean 0 and standard deviation 1. Licensing: This code is distributed under the GNU LGPL license. Modified: 06 August 2013 Author: John Burkardt Parameters: Input, int N, the number of values desired. Input/output, int *SEED, a seed for the random number generator. Output, double X[N], a sample of the standard normal PDF. Local parameters: Local, 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. Local, int X_LO, X_HI, records the range of entries of X that we need to compute. */ { int i; int m; double *r; const double r8_pi = 3.141592653589793; int x_hi; int x_lo; /* Record the range of X we need to fill in. */ x_lo = 1; x_hi = n; /* If we need just one new value, do that here to avoid null arrays. */ if ( x_hi - x_lo + 1 == 1 ) { r = r8vec_uniform_01_new ( 2, seed ); x[x_hi-1] = sqrt ( - 2.0 * log ( r[0] ) ) * cos ( 2.0 * r8_pi * r[1] ); 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, seed ); 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] ); } 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, seed ); 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] ); free ( r ); } return; } /******************************************************************************/ double *r8vec_normal_01_new ( int n, int *seed ) /******************************************************************************/ /* 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. Licensing: This code is distributed under the GNU LGPL license. Modified: 06 August 2013 Author: John Burkardt Parameters: Input, int N, the number of values desired. Input/output, int *SEED, a seed for the random number generator. Output, double R8VEC_NORMAL_01_NEW[N], a sample of the standard normal PDF. Local parameters: Local, 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. Local, int X_LO, X_HI, records the range of entries of X that we need to compute. */ { int i; int m; double *r; const double r8_pi = 3.141592653589793; double *x; int x_hi; int x_lo; x = ( double * ) malloc ( n * sizeof ( double ) ); /* Record the range of X we need to fill in. */ x_lo = 1; x_hi = n; /* If we need just one new value, do that here to avoid null arrays. */ if ( x_hi - x_lo + 1 == 1 ) { r = r8vec_uniform_01_new ( 2, seed ); x[x_hi-1] = sqrt ( - 2.0 * log ( r[0] ) ) * cos ( 2.0 * r8_pi * r[1] ); 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, seed ); 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] ); } 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, seed ); 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] ); free ( r ); } return x; } /******************************************************************************/ void r8vec_normal_ab ( int n, double b, double c, int *seed, double x[] ) /******************************************************************************/ /* Purpose: R8VEC_NORMAL_AB returns a scaled pseudonormal R8VEC. Discussion: The scaled normal probability distribution function (PDF) has mean A and standard deviation B. Licensing: This code is distributed under the GNU LGPL license. Modified: 06 August 2013 Author: John Burkardt Parameters: Input, int N, the number of values desired. Input, double B, C, the mean and standard deviation. Input/output, int *SEED, a seed for the random number generator. Output, double X[N], a sample of the standard normal PDF. Local parameters: Local, 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. Local, int X_LO, X_HI, records the range of entries of X that we need to compute. */ { int i; int m; double *r; const double r8_pi = 3.141592653589793; int x_hi; int x_lo; /* Record the range of X we need to fill in. */ x_lo = 1; x_hi = n; /* If we need just one new value, do that here to avoid null arrays. */ if ( x_hi - x_lo + 1 == 1 ) { r = r8vec_uniform_01_new ( 2, seed ); x[x_hi-1] = sqrt ( - 2.0 * log ( r[0] ) ) * cos ( 2.0 * r8_pi * r[1] ); 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, seed ); 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] ); } 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, seed ); 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] ); free ( r ); } for ( i = 0; i < n; i++ ) { x[i] = b + c * x[i]; } return; } /******************************************************************************/ double *r8vec_normal_ab_new ( int n, double b, double c, int *seed ) /******************************************************************************/ /* Purpose: R8VEC_NORMAL_AB_NEW returns a scaled pseudonormal R8VEC. Discussion: The scaled normal probability distribution function (PDF) has mean A and standard deviation B. Licensing: This code is distributed under the GNU LGPL license. Modified: 06 August 2013 Author: John Burkardt Parameters: Input, int N, the number of values desired. Input, double B, C, the mean and standard deviation. Input/output, int *SEED, a seed for the random number generator. Output, double R8VEC_NORMAL_AB_NEW[N], a sample of the standard normal PDF. Local parameters: Local, 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. Local, int X_LO, X_HI, records the range of entries of X that we need to compute. */ { int i; int m; double *r; const double r8_pi = 3.141592653589793; double *x; int x_hi; int x_lo; x = ( double * ) malloc ( n * sizeof ( double ) ); /* Record the range of X we need to fill in. */ x_lo = 1; x_hi = n; /* If we need just one new value, do that here to avoid null arrays. */ if ( x_hi - x_lo + 1 == 1 ) { r = r8vec_uniform_01_new ( 2, seed ); x[x_hi-1] = sqrt ( - 2.0 * log ( r[0] ) ) * cos ( 2.0 * r8_pi * r[1] ); 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, seed ); 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] ); } 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). */ else { x_hi = x_hi - 1; m = ( x_hi - x_lo + 1 ) / 2 + 1; r = r8vec_uniform_01_new ( 2*m, seed ); 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] ); free ( r ); } for ( i = 0; i < n; i++ ) { x[i] = b + c * x[i]; } return x; } /******************************************************************************/ 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 GNU LGPL license. Modified: 08 April 2009 Author: John Burkardt Parameters: Input, int N, the number of components of the vector. Input, double A[N], the vector to be printed. Input, char *TITLE, a title. */ { int i; fprintf ( stdout, "\n" ); fprintf ( stdout, "%s\n", title ); fprintf ( stdout, "\n" ); for ( i = 0; i < n; i++ ) { fprintf ( stdout, " %8d: %14g\n", i, a[i] ); } return; } /******************************************************************************/ double *r8vec_uniform_01_new ( int n, int *seed ) /******************************************************************************/ /* Purpose: R8VEC_UNIFORM_01_NEW returns a unit pseudorandom R8VEC. Discussion: This routine implements the recursion seed = 16807 * seed mod ( 2^31 - 1 ) unif = seed / ( 2^31 - 1 ) The integer arithmetic never requires more than 32 bits, including a sign bit. Licensing: This code is distributed under the GNU LGPL license. Modified: 19 August 2004 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. Peter Lewis, Allen Goodman, James Miller, A Pseudo-Random Number Generator for the System/360, IBM Systems Journal, Volume 8, pages 136-143, 1969. Parameters: Input, int N, the number of entries in the vector. Input/output, int *SEED, a seed for the random number generator. Output, double R8VEC_UNIFORM_01_NEW[N], the vector of pseudorandom values. */ { int i; const int i4_huge = 2147483647; int k; double *r; r = ( double * ) malloc ( n * sizeof ( double ) ); for ( i = 0; i < n; i++ ) { k = *seed / 127773; *seed = 16807 * ( *seed - k * 127773 ) - k * 2836; if ( *seed < 0 ) { *seed = *seed + i4_huge; } r[i] = ( double ) ( *seed ) * 4.656612875E-10; } return r; } /******************************************************************************/ void timestamp ( ) /******************************************************************************/ /* Purpose: TIMESTAMP prints the current YMDHMS date as a time stamp. Example: May 31 2001 09:45:54 AM Licensing: This code is distributed under the GNU LGPL license. Modified: 02 October 2003 Author: John Burkardt Parameters: None */ { # define TIME_SIZE 40 static char time_buffer[TIME_SIZE]; const struct tm *tm; time_t now; now = time ( NULL ); tm = localtime ( &now ); strftime ( time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm ); printf ( "%s\n", time_buffer ); return; # undef TIME_SIZE }