/*************************************************************************** * Rand.C (Random Number Generators) * * * * Source file for pseudo-random number utilities. Rand is the * * base class for several different algorithms for generating pseudo-random * * numbers. Any method can generate individual samples or arrays of * * samples using "Eval". The random seed can be reset at any time by * * calling "Seed" with any integer. Random permutations of the integers * * 0,1,...(n-1) are generated by "Perm(n,P)". * * * * HISTORY * * Name Date Description * * * * arvo 08/04/97 Changed to virtual functions. * * arvo 06/06/93 Optimization, especially for array evaluators. * * arvo 10/06/91 Converted to C++ * * arvo 11/20/89 Added "gen_seed" function to handle. * * arvo 10/30/89 "state" allocation now done in rand_alloc. * * arvo 07/08/89 Initial coding. * * * *--------------------------------------------------------------------------* * Copyright (C) 1989, James Arvo * * * * This program is free software; you can redistribute it and/or modify it * * under the terms of the GNU General Public License as published by the * * Free Software Foundation. See http://www.fsf.org/copyleft/gpl.html * * * * This program is distributed in the hope that it will be useful, but * * WITHOUT EXPRESS OR IMPLIED WARRANTY of merchantability or fitness for * * any particular purpose. See the GNU General Public License for more * * details. * * * ***************************************************************************/ #include #include #include #ifndef ABS #define ABS( x ) ((x) > 0 ? (x) : -(x)) #endif /*-------------------------------------------------------------------------* * M E T H O D 1 * * * * From "Numerical Recipes," by William H. Press, Brian P. Flannery, * * Saul A. Teukolsky, and William T. Vetterling, p. 197. * * * *-------------------------------------------------------------------------*/ static const long M1 = 714025; static const long IA = 1366; static const long IC = 150889; static const double RM = 1.400512E-6; float RandGen_1::Eval() { register long *elem; register long offset; register float rand; offset = 1 + ( 97 * index ) / M1; if( offset > 97 ) offset = 97; if( offset < 1 ) offset = 1; elem = shuffle + offset; rand = ( index = *elem ) * RM; *elem = ( seed = ( IA * seed + IC ) % M1 ); return rand; } void RandGen_1::Eval( int n, float *array ) { register long *shfl = shuffle; register long *elem; register long offset; for( int i = 0; i < n; i++ ) { offset = 1 + ( 97 * index ) / M1; if( offset > 97 ) offset = 97; if( offset < 1 ) offset = 1; elem = shfl + offset; *array++ = ( index = *elem ) * RM; *elem = ( seed = ( IA * seed + IC ) % M1 ); } } void RandGen_1::Seed( long seed ) { long t = ( IC + ABS( seed ) + 1 ) % M1; for( register int k = 1; k <= 97; k++ ) { t = ( IA * t + IC ) % M1; shuffle[k] = ABS( t ); } t = ( IA * t + IC ) % M1; seed = ABS( t ); index = ABS( t ); } /*-------------------------------------------------------------------------* * M E T H O D 2 * * * * From "The Multiple Prime Random Number Generator," by Alexander Haas, * * ACM Transactions on Mathematical Software, Vol. 13, No. 4, December * * 1987, pp. 368-381. * * * *-------------------------------------------------------------------------*/ float RandGen_2::Eval() { if( (m += 7 ) >= 9973 ) m -= 9871; if( (i += 1907 ) >= 99991 ) i -= 89989; if( (j += 73939) >= 224729 ) j -= 96233; r = ((r * m + i + j) % 100000) / 10; return r * 1.00010001E-4; } void RandGen_2::Eval( int n, float *array ) { for( register int k = 0; k < n; k++ ) { if( (m += 7 ) >= 9973 ) m -= 9871; if( (i += 1907 ) >= 99991 ) i -= 89989; if( (j += 73939) >= 224729 ) j -= 96233; r = ((r * m + i + j) % 100000) / 10; *array++ = r * 1.00010001E-4; } } void RandGen_2::Seed( long seed ) { r = ABS( seed ); m = ABS( seed * 7 ); i = ABS( seed * 11 ); j = ABS( seed * 13 ); if( m < 100 ) m += 100; if( i < 10000 ) i += 10000; if( j < 128000 ) j += 128000; } /*-------------------------------------------------------------------------* * M E T H O D 3 * * * * From "A More Portable Fortran Random Number Generator," by Linus * * Schrage, ACM Transactions on Mathematical Software, Vol. 5, No, 2, * * June 1979, pp. 132-138. * * * *-------------------------------------------------------------------------*/ static const long A3 = 16807; static const long P3 = 2147483647; float RandGen_3::Eval() { long xhi = ix >> 16; long xalo = ( ix & 0xFFFF ) * A3; long leftlo = xalo >> 16; long fhi = xhi * A3 + leftlo; long k = fhi >> 15; ix = ( ((xalo - (leftlo << 16)) - P3) + ((fhi - (k << 15)) << 16) ) + k; if( ix < 0 ) ix += P3; return ix * 4.656612875E-10; } void RandGen_3::Eval( int n, float *array ) { register long xhi, xalo, leftlo; register long fhi, k; for( register int i = 0; i < n; i++ ) { xhi = ix >> 16; xalo = ( ix & 0xFFFF ) * A3; leftlo = xalo >> 16; fhi = xhi * A3 + leftlo; k = fhi >> 15; ix = ( ((xalo - (leftlo << 16)) - P3) + ((fhi - (k << 15)) << 16) ) + k; if( ix < 0 ) ix += P3; *array++ = ix * 4.656612875E-10; } } void RandGen_3::Seed( long seed ) { ix = ABS( seed ); } /*-------------------------------------------------------------------------* * R A N D : : P E R M (Permutation) * * * * This routine fills an integer array of length "len" with a random * * permutation of the integers 0, 1, 2, ... (len-1). * * * * For efficiency, the random numbers are generated in batches of up to * * "Nmax" at a time. The constant Nmax can be set to any value >= 1. * * * *-------------------------------------------------------------------------*/ static const int Nmax = 20; void RandGen::Perm( int len, int perm[] ) { float R[ Nmax ]; // A buffer for getting random numbers. int L = len - 1; // Total number of random numbers needed. int N = 0; // How many to generate when we call Eval. int n = 0; // The array index into R. // First initialize the array "perm" to the identity permutation. for( int j = 0; j < len; j++ ) perm[j] = j; // Now swap a random element in the front with the i'th element. // When i gets down to 0, we're done. for( int i = len - 1; i > 0; i-- ) // Element i is a swap candidate. { if( n == N ) // Generate more random numbers. { N = ( L < Nmax ) ? L : Nmax; // Can't get more than "Nmax". Eval( N, R ); // Generate N random numbers. L -= N; // Decrement total counter. n = 0; // Start index at beginning of R. } float r = ( i + 1 ) * R[ n++ ]; // Pick a float in [0,i+1]. int k = (int)r; // Truncate r to an integer. if( k < i ) // Disregard k == i and k == i+1. { int tmp = perm[i]; // Swap elements i and k. perm[i] = perm[k]; perm[k] = tmp; } } }