/usr/include/random/normal.h is in libblitz0-dev 1:0.10-3.3.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 | // -*- C++ -*-
// $Id$
/*
* This is a modification of the Kinderman + Monahan algorithm for
* generating normal random numbers, due to Leva:
*
* J.L. Leva, Algorithm 712. A normal random number generator, ACM Trans.
* Math. Softw. 18 (1992) 454--455.
*
* http://www.acm.org/pubs/citations/journals/toms/1992-18-4/p449-leva/
*
* Note: Some of the constants used below look like they have dubious
* precision. These constants are used for an approximate bounding
* region test (see the paper). If the approximate test fails,
* then an exact region test is performed.
*
* Only 0.012 logarithm evaluations are required per random number
* generated, making this method comparatively fast.
*
* Adapted to C++ by T. Veldhuizen.
*/
#ifndef BZ_RANDOM_NORMAL
#define BZ_RANDOM_NORMAL
#ifndef BZ_RANDOM_UNIFORM
#include <random/uniform.h>
#endif
BZ_NAMESPACE(ranlib)
template<typename T = double, typename IRNG = defaultIRNG,
typename stateTag = defaultState>
class NormalUnit : public UniformOpen<T,IRNG,stateTag>
{
public:
typedef T T_numtype;
NormalUnit() {}
explicit NormalUnit(unsigned int i) :
UniformOpen<T,IRNG,stateTag>(i) {};
T random()
{
const T s = 0.449871, t = -0.386595, a = 0.19600, b = 0.25472;
const T r1 = 0.27597, r2 = 0.27846;
T u, v;
for (;;) {
// Generate P = (u,v) uniform in rectangle enclosing
// acceptance region:
// 0 < u < 1
// - sqrt(2/e) < v < sqrt(2/e)
// The constant below is 2*sqrt(2/e).
u = this->getUniform();
v = 1.715527769921413592960379282557544956242L
* (this->getUniform() - 0.5);
// Evaluate the quadratic form
T x = u - s;
T y = fabs(v) - t;
T q = x*x + y*(a*y - b*x);
// Accept P if inside inner ellipse
if (q < r1)
break;
// Reject P if outside outer ellipse
if (q > r2)
continue;
// Between ellipses: perform exact test
if (v*v <= -4.0 * log(u)*u*u)
break;
}
return v/u;
}
};
template<typename T = double, typename IRNG = defaultIRNG,
typename stateTag = defaultState>
class Normal : public NormalUnit<T,IRNG,stateTag> {
public:
typedef T T_numtype;
Normal(T mean, T standardDeviation)
{
mean_ = mean;
standardDeviation_ = standardDeviation;
}
Normal(T mean, T standardDeviation, unsigned int i) :
NormalUnit<T,IRNG,stateTag>(i)
{
mean_ = mean;
standardDeviation_ = standardDeviation;
};
T random()
{
return mean_ + standardDeviation_
* NormalUnit<T,IRNG,stateTag>::random();
}
private:
T mean_;
T standardDeviation_;
};
BZ_NAMESPACE_END
#endif // BZ_RANDOM_NORMAL
|