 |
|
 |
To address a problem with seeding I have now re-coded FastRandom to use a completely different algorithm, namely an xor-shift method. Although there was a much simpler fix for the problem, I found this new algorithm a little faster so I decided to switch.
The algorithm itself is explained in a paper linked from the top of the article.
|
|
|
|
 |
|
 |
Thanks for the reference.
I have one remaining suggestion: inherit from System.Random and override methods as appropriate. This will truly make your class a direct replacement for System.Random.
Jeffrey
Everything should be as simple as possible, but not simpler. -- Albert Einstein
Numerical components for C# and VB.NET
|
|
|
|
 |
|
 |
Perhaps I'll leave that as an exercise for the reader to do if and when they need a direct replacement e.g. if some binary library method only accepts System.Random objects.
Non-virtual methods have ever-so-slightly more overhead than virtual ones, so the current arrangement gives the fastest possible speed.
Thanks for the ongoing feedback BTW. And I did spot the 'Marsaglia effect' in the link you posted
|
|
|
|
 |
|
 |
Hello
Why do you not use one of the standard algorithms that are proofed to behave randomly.
(I can't remember the name of THE standard algorithm that uses simple multiplication and modulo)
Nevertheless I think it's a good and efficient algorithm for random numbers
Aloha
Urs
|
|
|
|
 |
|
 |
Hi,
As I understand it the algorithm used is proven to be random and to have a long period. With regard to the short period of 1 problem, George Marsaglia stipulates some restrictions on the seeds values that should be used. However, I have now re-coded the FastRandom class to use another of Marsaglias's algorithms that is even faster than multiply-and-carry. This new algorithm is based on XOR and bit shifting only and does not exhibit the same problem with the Reinitialise methods.
Just waiting for the Codeproject folks to update the article
|
|
|
|
 |
|
 |
Please define "proven to be random".
I would like to note some references on rngs:
D. E. Knuth's "The Art of Computer Programming" gives an extensive discussion on (mostly linear) random number generators.
Theory and Practice of Random Number Generation by Peter Hellekalek at the Mathematics Department of the University of Salzburg is a good starting point for practitioneers.
Jens
--
Jens Scheidtmann
|
|
|
|
 |
|
 |
Well ok. The algorithm is proven to provide a good pseudo-random sequence of long period. Please see George Marsaglia's paper (linked in the article) and general body of work on RNGs for more info on tests of randomness, and mathemetical basis of this RNG, etc.
|
|
|
|
 |
|
 |
Urs Enzler wrote:
Why do you not use one of the standard algorithms that are proofed to behave randomly. (I can't remember the name of THE standard algorithm that uses simple multiplication and modulo)
I believe you're talking about the Linear Congruential method. The algorithm in the original article used a combination of two such RNG's. Whatever method you use, you must be careful to verify that the sequence you generate works well using the parameters you chose. A bad choice of parameters can destroy the appearance of randomness.
'Standard' is not always better. The RANDU random number generator was published by IBM in the 60's and something of a standard back then. It is of the linear congruential type. However, it was later discovered that it wasn't all that good. See this page[^] for an illustration (Java applet).
Constructing pseudo-random number generators is a bit of an art, it seems. There is a whole range of algorithms out there, with a large range of speed and quality.
Which algorithm you use depends on your application. For example, cryptography requires much better randomness than a simple game because the whole point of cryptography is to destroy any appearance of a pattern or relationship between the encrypted data and the original. In simulations, a pattern like that produced by the RANDU algorithm can introduce a bias that taints the results.
Jeffrey
Everything should be as simple as possible, but not simpler. -- Albert Einstein
Numerical components for C# and VB.NET
|
|
|
|
 |
|
 |
You are absolutely right. With 'standard' I meant an algorithm that is proven to be good.
The algorithm you mentioned is indeed very famous because a lot of scientific work had to be thrown away because of it On the other hand some simulation proofed to be very stable even against a very bad random number generator.
|
|
|
|
 |
|
 |
colgreen wrote:
That simulation requires that the RNG be reset with a given seed 1000s of times per second. FastRandom's Reinitialise() methods therefore provide a nice performance boost over System.Random in that case. I then discovered that a number of further performance improvements could be made to the Next*() methods.
How many random numbers do actually generate within the loop? Wouldn't it be a lot faster if you created the sequence outside the loop, put them in an array, and inside the loop just picked off numbers from that array?? You may be able to do some pre-processing of these random values outside the loop, leading to an even greater speed-up.
Jeffrey
Everything should be as simple as possible, but not simpler. -- Albert Einstein
Numerical components for C# and VB.NET
|
|
|
|
 |
|
 |
Hi, Yes that is an option but some simulations are open ended with regards to their end condition and so could require millions or billions of rnd numbers. But yes, for the more common requirement of say a few thousand numbers caching numbers would be faster, in fact I might try this out.
Thanks for the suggestion.
Colin
|
|
|
|
 |
|
 |
I can see why this routine is very fast, but this is an extremely basic mechanism for generating random numbers which can result in a very short period for particular values
|
|
|
|
 |
|
 |
Jabes wrote:
this is an extremely basic mechanism for generating random numbers which can result in a very short period for particular values
Indeed it can! The period can be as short as 1:
rng = new FastRandom(-2006912572, 1179647999);
if (rng.Next() == rng.Next())
Console.WriteLine("'Random' sample is always the same.");
Now, it is indeed very fast, and for most seeds it will produce a pretty good sequence. This example shows that caution is always a good thing.
The creator of this random number generator is one George Marsaglia, who wrote this message[^] to the sci.stat.math newsgroup. He also offers replacements for the parameters 36969 and 18000. Marsaglia also wrote DIEHARD[^], a series of programs for testing the quality of random number generators.
Jeffrey
Everything should be as simple as possible, but not simpler. -- Albert Einstein
Numerical components for C# and VB.NET
|
|
|
|
 |
|
 |
I think this problem can be addressed with some extra logic in the Reinitialise methods that ensure that bad seeds are avoided. Note that I already added conditions to avoid zero seeds, which will just result in a sequence of zero's. And as you rightly point out (Jabes) there are other seeds that also perform poorly - which is a problem I hadn't spotted, but the mechanism itself is capable of generating good sequences. That is, it isn't inherently flawed AFAIK.
|
|
|
|
 |
|
 |
If you are after an excellent random number generator, then the mersenne twister is excellent. I have used it in the past with very good, fast results.
I haven't seen a C# implementation yet, but I'm sure there must be one kicking around
|
|
|
|
 |
|
 |
There is one as part of the Netron project: http://netron.sourceforge.net/wp/[^]
The file would be C:\Program Files\The Netron Project\Cobalt\NetronGraphLibrary\Utils\Maths\MersenneTwister.cs
using System;
namespace Netron.GraphLib.Maths
{
/////////////////////////////////////////////////////////////////////////////
// C# Version Copyright (c) 2003 CenterSpace Software, LLC //
// //
// This code is free software under the Artistic license. //
// //
// CenterSpace Software //
// 2098 NW Myrtlewood Way //
// Corvallis, Oregon, 97330 //
// USA //
// http://www.centerspace.net //
/////////////////////////////////////////////////////////////////////////////
/*
A C-program for MT19937, with initialization improved 2002/1/26.
Coded by Takuji Nishimura and Makoto Matsumoto.
Before using, initialize the state by using init_genrand(seed)
or init_by_array(init_key, key_length).
Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura,
All rights reserved.
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions
are met:
1. Redistributions of source code must retain the above copyright
notice, this list of conditions and the following disclaimer.
2. Redistributions in binary form must reproduce the above copyright
notice, this list of conditions and the following disclaimer in the
documentation and/or other materials provided with the distribution.
3. The names of its contributors may not be used to endorse or promote
products derived from this software without specific prior written
permission.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
Any feedback is very welcome.
http://www.math.keio.ac.jp/matumoto/emt.html
email: matumoto@math.keio.ac.jp
*/
///
/// Class MersenneTwister generates random numbers from a uniform distribution using
/// the Mersenne Twister algorithm.
///
/// Caution: MT is for MonteCarlo, and is NOT SECURE for CRYPTOGRAPHY
/// as it is.
[Serializable()]public class MersenneTwister
{
#region Constants -------------------------------------------------------
// Period parameters.
private const int N = 624;
private const int M = 397;
private const uint MATRIX_A = 0x9908b0dfU; // constant NetronVector a
private const uint UPPER_MASK = 0x80000000U; // most significant w-r bits
private const uint LOWER_MASK = 0x7fffffffU; // least significant r bits
private const int MAX_RAND_INT = 0x7fffffff;
#endregion Constants
#region Instance Variables ----------------------------------------------
// mag01[x] = x * MATRIX_A for x=0,1
private uint[] mag01 = {0x0U, MATRIX_A};
// the array for the state NetronVector
private uint[] mt = new uint[N];
// mti==N+1 means mt[N] is not initialized
private int mti = N+1;
#endregion Instance Variables
#region Constructors ----------------------------------------------------
///
/// Creates a random number generator using the time of day in milliseconds as
/// the seed.
///
public MersenneTwister()
{
init_genrand( (uint)DateTime.Now.Millisecond );
}
///
/// Creates a random number generator initialized with the given seed.
///
/// The seed.
public MersenneTwister( int seed )
{
init_genrand( (uint)seed );
}
///
/// Creates a random number generator initialized with the given array.
///
/// The array for initializing keys.
public MersenneTwister( int[] init )
{
uint[] initArray = new uint[init.Length];
for ( int i = 0; i < init.Length; ++i )
initArray[i] = (uint)init[i];
init_by_array( initArray, (uint)initArray.Length );
}
#endregion Constructors
#region Properties ------------------------------------------------------
///
/// Gets the maximum random integer value. All random integers generated
/// by instances of this class are less than or equal to this value. This
/// value is 0x7fffffff (2,147,483,647).
///
public static int MaxRandomInt
{
get
{
return 0x7fffffff;
}
}
#endregion Properties
#region Member Functions ------------------------------------------------
///
/// Returns a random integer greater than or equal to zero and
/// less than or equal to MaxRandomInt.
///
/// The next random integer.
public int Next()
{
return genrand_int31();
}
///
/// Returns a positive random integer less than the specified maximum.
///
/// The maximum value. Must be greater than zero.
/// A positive random integer less than or equal to maxValue.
public int Next( int maxValue )
{
return Next( 0, maxValue );
}
///
/// Returns a random integer within the specified range.
///
/// The lower bound.
/// The upper bound.
/// A random integer greater than or equal to minValue, and less than
/// or equal to maxValue.
public int Next( int minValue, int maxValue )
{
if ( minValue > maxValue )
{
int tmp = maxValue;
maxValue = minValue;
minValue = tmp;
}
return (int)( Math.Floor((maxValue-minValue+1)*genrand_real1() + minValue) );
}
///
/// Returns a random number between 0.0 and 1.0.
///
/// A single-precision floating point number greater than or equal to 0.0,
/// and less than 1.0.
public float NextFloat()
{
return (float) genrand_real2();
}
///
/// Returns a random number greater than or equal to zero, and either strictly
/// less than one, or less than or equal to one, depending on the value of the
/// given boolean parameter.
///
///
/// If true, the random number returned will be
/// less than or equal to one; otherwise, the random number returned will
/// be strictly less than one.
///
///
/// If includeOne is true, this method returns a
/// single-precision random number greater than or equal to zero, and less
/// than or equal to one. If includeOne is false, this method
/// returns a single-precision random number greater than or equal to zero and
/// strictly less than one.
///
public float NextFloat( bool includeOne )
{
if ( includeOne )
{
return (float) genrand_real1();
}
return (float) genrand_real2();
}
///
/// Returns a random number greater than 0.0 and less than 1.0.
///
/// A random number greater than 0.0 and less than 1.0.
public float NextFloatPositive()
{
return (float) genrand_real3();
}
///
/// Returns a random number between 0.0 and 1.0.
///
/// A double-precision floating point number greater than or equal to 0.0,
/// and less than 1.0.
public double NextDouble()
{
return genrand_real2();
}
///
/// Returns a random number greater than or equal to zero, and either strictly
/// less than one, or less than or equal to one, depending on the value of the
/// given boolean parameter.
///
///
/// If true, the random number returned will be
/// less than or equal to one; otherwise, the random number returned will
/// be strictly less than one.
///
///
/// If includeOne is true, this method returns a
/// single-precision random number greater than or equal to zero, and less
/// than or equal to one. If includeOne is false, this method
/// returns a single-precision random number greater than or equal to zero and
/// strictly less than one.
///
public double NextDouble( bool includeOne )
{
if ( includeOne )
{
return genrand_real1();
}
return genrand_real2();
}
///
/// Returns a random number greater than 0.0 and less than 1.0.
///
/// A random number greater than 0.0 and less than 1.0.
public double NextDoublePositive()
{
return genrand_real3();
}
///
/// Generates a random number on [0,1) with 53-bit resolution.
///
/// A random number on [0,1) with 53-bit resolution
public double Next53BitRes()
{
return genrand_res53();
}
///
/// Reinitializes the random number generator using the time of day in
/// milliseconds as the seed.
///
public void Initialize()
{
init_genrand( (uint)DateTime.Now.Millisecond );
}
///
/// Reinitializes the random number generator with the given seed.
///
/// The seed.
public void Initialize( int seed )
{
init_genrand( (uint)seed );
}
///
/// Reinitializes the random number generator with the given array.
///
/// The array for initializing keys.
public void Initialize( int[] init )
{
uint[] initArray = new uint[init.Length];
for ( int i = 0; i < init.Length; ++i )
initArray[i] = (uint)init[i];
init_by_array( initArray, (uint)initArray.Length );
}
#region Methods ported from C -------------------------------------------
// initializes mt[N] with a seed
private void init_genrand( uint s)
{
mt[0]= s & 0xffffffffU;
for (mti=1; mti> 30)) + mti);
// See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier.
// In the previous versions, MSBs of the seed affect
// only MSBs of the array mt[].
// 2002/01/09 modified by Makoto Matsumoto
mt[mti] &= 0xffffffffU;
// for >32 bit machines
}
}
// initialize by an array with array-length
// init_key is the array for initializing keys
// key_length is its length
private void init_by_array(uint[] init_key, uint key_length)
{
int i, j, k;
init_genrand(19650218U);
i=1; j=0;
k = (int)(N>key_length ? N : key_length);
for (; k>0; k--)
{
mt[i] = (uint)((uint)(mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1664525U)) + init_key[j] + j); /* non linear */
mt[i] &= 0xffffffffU; // for WORDSIZE > 32 machines
i++; j++;
if (i>=N) { mt[0] = mt[N-1]; i=1; }
if (j>=key_length) j=0;
}
for (k=N-1; k>0; k--)
{
mt[i] = (uint)((uint)(mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1566083941U))- i); /* non linear */
mt[i] &= 0xffffffffU; // for WORDSIZE > 32 machines
i++;
if (i>=N) { mt[0] = mt[N-1]; i=1; }
}
mt[0] = 0x80000000U; // MSB is 1; assuring non-zero initial array
}
// generates a random number on [0,0xffffffff]-interval
uint genrand_int32()
{
uint y;
if (mti >= N)
{ /* generate N words at one time */
int kk;
if (mti == N+1) /* if init_genrand() has not been called, */
init_genrand(5489U); /* a default initial seed is used */
for (kk=0;kk> 1) ^ mag01[y & 0x1U];
}
for (;kk> 1) ^ mag01[y & 0x1U];
}
y = (mt[N-1]&UPPER_MASK)|(mt[0]&LOWER_MASK);
mt[N-1] = mt[M-1] ^ (y >> 1) ^ mag01[y & 0x1U];
mti = 0;
}
y = mt[mti++];
// Tempering
y ^= (y >> 11);
y ^= (y << 7) & 0x9d2c5680U;
y ^= (y << 15) & 0xefc60000U;
y ^= (y >> 18);
return y;
}
// generates a random number on [0,0x7fffffff]-interval
private int genrand_int31()
{
return (int)(genrand_int32()>>1);
}
// generates a random number on [0,1]-real-interval
double genrand_real1()
{
return genrand_int32()*(1.0/4294967295.0);
// divided by 2^32-1
}
// generates a random number on [0,1)-real-interval
double genrand_real2()
{
return genrand_int32()*(1.0/4294967296.0);
// divided by 2^32
}
// generates a random number on (0,1)-real-interval
double genrand_real3()
{
return (((double)genrand_int32()) + 0.5)*(1.0/4294967296.0);
// divided by 2^32
}
// generates a random number on [0,1) with 53-bit resolution
double genrand_res53()
{
uint a=genrand_int32()>>5, b=genrand_int32()>>6;
return(a*67108864.0+b)*(1.0/9007199254740992.0);
}
// These real versions are due to Isaku Wada, 2002/01/09 added
#endregion Methods ported from C
#endregion Member Functions
}
}
|
|
|
|
 |
|
 |
Correct.
Worst ever random generator!
slightly better than:
public void nextBytes(byte[] buffer)
{
// Why bother ?
}
Passes FIPS 140-1 with approximately 0.7 %.
|
|
|
|
 |
|
 |
I second your opinion Jabes.
|
|
|
|
 |