mirror of
https://port.numenaute.org/aleajactaest/khanat-opennel-code.git
synced 2024-12-17 06:38:43 +00:00
153 lines
4.2 KiB
C++
153 lines
4.2 KiB
C++
// Ryzom - MMORPG Framework <http://dev.ryzom.com/projects/ryzom/>
|
|
// Copyright (C) 2010 Winch Gate Property Limited
|
|
//
|
|
// This program is free software: you can redistribute it and/or modify
|
|
// it under the terms of the GNU Affero General Public License as
|
|
// published by the Free Software Foundation, either version 3 of the
|
|
// License, or (at your option) any later version.
|
|
//
|
|
// This program is distributed in the hope that it will be useful,
|
|
// but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
|
// GNU Affero General Public License for more details.
|
|
//
|
|
// You should have received a copy of the GNU Affero General Public License
|
|
// along with this program. If not, see <http://www.gnu.org/licenses/>.
|
|
|
|
#include "simulation_random.h"
|
|
#include "nel/misc/debug.h"
|
|
|
|
//
|
|
// Implementation of a random number generator that draws from a
|
|
// gaussian distribution. The underlying (uniform deviate)
|
|
// pseudorandom number generator is the Park and Miller "Minimum
|
|
// standard" generator with Bayes-Durham shuffle
|
|
//
|
|
// Copyright (c) 2001 Wesley H. Huang. All rights reserved.
|
|
//
|
|
#include <math.h>
|
|
//#include <rand.h>
|
|
|
|
|
|
// default random generator seed (nothing special about this number
|
|
// except that it's not 0)
|
|
const long defaultSeed = 18846224;
|
|
long pmRandGenSeed = defaultSeed;
|
|
|
|
// set the seed (optional)
|
|
void pmRandSeed(long s)
|
|
{
|
|
if (s == 0) // can't use 0 as a seed
|
|
pmRandGenSeed = defaultSeed; // use the default
|
|
else
|
|
pmRandGenSeed = s;
|
|
}
|
|
|
|
// parameters for the random number generator
|
|
const long a = 16807;
|
|
const long m = 2147483647;
|
|
const long q = m/a; // 127773
|
|
const long r = m % a; // 2836
|
|
const long pmRandMax = m;
|
|
|
|
// Park and Miller's "Minimal Standard" generator using Schrage's
|
|
// method to do the computation in 32 bit arithmetic without overflow.
|
|
//
|
|
// This is a multiplicative congruential generator: I_{j+1} = a I_j mod m
|
|
// using the above constants (which are very carefully picked).
|
|
//
|
|
// The result is an integer between 1 and m-1 inclusive.
|
|
//
|
|
// Reference: Numerical Recipes in C, 2nd ed, pg 278.
|
|
long pmRand()
|
|
{
|
|
long k = pmRandGenSeed/q;
|
|
pmRandGenSeed = a*(pmRandGenSeed - k*q) - r*k;
|
|
if (pmRandGenSeed < 0)
|
|
pmRandGenSeed += m;
|
|
return pmRandGenSeed;
|
|
}
|
|
|
|
// the Bayes-Durham shuffle
|
|
long pmRandShuffle()
|
|
{
|
|
const int size = 32;
|
|
const long div = 1 + (pmRandMax-1)/size;
|
|
static bool initialized = false;
|
|
static long rn[size];
|
|
static long next;
|
|
int index;
|
|
|
|
if (!initialized)
|
|
{
|
|
int k;
|
|
for (k= 0; k < size; k++) // warm up
|
|
pmRand();
|
|
for (k= 0; k < size; k++) // fill table
|
|
rn[k] =pmRand();
|
|
next = pmRand();
|
|
initialized = true;
|
|
}
|
|
index = next/div; // find out which table entry to use next
|
|
next = rn[index]; // that random nember determines the next index
|
|
rn[index] = pmRand(); // refill the table slot
|
|
return next; // return the random number
|
|
}
|
|
|
|
//
|
|
// This is a simple but pretty good random number generator.
|
|
//
|
|
// Returns a floating point number in (0.0, 1.0)
|
|
float myRand()
|
|
{
|
|
return pmRandShuffle()/(pmRandMax + 1.0f);
|
|
}
|
|
|
|
void myRandSeed(long s)
|
|
{
|
|
pmRandSeed(s);
|
|
}
|
|
|
|
// returns a number drawn from the exponential distribution with the given mean
|
|
float exponential( float mean )
|
|
{
|
|
float uniform = myRand();
|
|
|
|
nlassert( uniform < 1.0f );
|
|
float exponential = -logf( 1.0f - uniform ) * mean;
|
|
|
|
// nldebug( "uniform rand: %.2f, exponential rand: %.2f", uniform, exponential );
|
|
return exponential;
|
|
}
|
|
|
|
// returns a number drawn from a gaussian distribution with the
|
|
// specified mean and standard deviation
|
|
float gaussian(float mean, float stdev)
|
|
{
|
|
static bool cached = false;
|
|
static float extra; // the extra number from a 0 mean unit stdev gaussian
|
|
|
|
if (cached)
|
|
{
|
|
cached = false;
|
|
return extra*stdev + mean;
|
|
}
|
|
else
|
|
{
|
|
// pick a random point in the unit circle (excluding the origin)
|
|
float a,b,c;
|
|
do
|
|
{
|
|
a = 2.0f*myRand()-1.0f;
|
|
b = 2.0f*myRand()-1.0f;
|
|
c = a*a + b*b;
|
|
}
|
|
while (c >= 1.0f || c == 0.0);
|
|
|
|
// transform it into two values drawn from a gaussian
|
|
float t = sqrtf(-2.0f*logf(c)/c);
|
|
extra = t*a;
|
|
cached = true;
|
|
return t*b*stdev + mean;
|
|
}
|
|
}
|