---

math/random.h

Go to the documentation of this file.

#ifndef MATH_RANDOM_H
#define MATH_RANDOM_H

/*

  GLT OpenGL C++ Toolkit (LGPL)
  Copyright (C) 2000-2002 Nigel Stewart  

  Email: nigels.com@gmail.com   
  WWW:   http://www.nigels.com/glt/

  This library is free software; you can redistribute it and/or
  modify it under the terms of the GNU Lesser General Public
  License as published by the Free Software Foundation; either
  version 2.1 of the License, or (at your option) any later version.

  This library 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
  Lesser General Public License for more details.

  You should have received a copy of the GNU Lesser General Public
  License along with this library; if not, write to the Free Software
  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

*/

#include <iosfwd>
#include <cassert>
#include <cmath>

#include <glt/config.h>

#include <math/real.h>
#include <math/vector3.h>
#include <math/matrix4.h>

class GltRandomLCG 
{
public:
    GltRandomLCG(const uint32 seed = 199125);
    GltRandomLCG(const GltRandomLCG &rng);
    ~GltRandomLCG();

    uint32 rand() const;
    uint32 max() const;

    void seed(const uint32 seed = 199125);

    static GltRandomLCG rng;

          GltRandomLCG &base()       { return *this; }
    const GltRandomLCG &base() const { return *this; }

private:
    
    mutable uint32 _idum;

    const static int32 a;
    const static int32 m;
    const static int32 q;
    const static int32 r;
};

class GltRandomLFSRMix 
{
public:
    GltRandomLFSRMix(const uint32 seed1 = 199125,const uint32 seed2 = 90618,const uint32 seed3 = 189419);
    GltRandomLFSRMix(const GltRandomLFSRMix &rng);
    ~GltRandomLFSRMix();

    uint32 rand() const;
    uint32 max() const;

    static GltRandomLFSRMix rng;

          GltRandomLFSRMix &base()       { return *this; }
    const GltRandomLFSRMix &base() const { return *this; }

private:
    
    mutable uint32 _regA, _regB, _regC;
};

template<class R = GltRandomLFSRMix>
class GltRandomDouble
{
public:

    GltRandomDouble(R &random,const double min = 0.0,const double max = 1.0)
    : _random(random),
      _min(min),
      _max(max),
      _mult((max-min)/random.max())
    {
    }

    GltRandomDouble(const double min = 0.0,const double max = 1.0)
    : _random(R::rng),
      _min(min),
      _max(max),
      _mult((max-min)/R::rng.max())
    {
    }

    GltRandomDouble(const GltRandomDouble<R> &gen) 
    : _random(gen._random),
      _min(gen._min),
      _max(gen._max),
      _mult(gen._mult)
    {
    }

    ~GltRandomDouble() 
    {
    }

    double rand() const 
    {
        return _min + _mult*double(_random.rand());
    }

          R &base()       { return _random; }
    const R &base() const { return _random; }
    
private:

    R            &_random;
    const double  _min;
    const double  _max;
    const double  _mult;
};

template<class R = GltRandomLFSRMix>
class GltRandomInteger
{
public:

    GltRandomInteger(R &random,const uint32 min,const uint32 max)
    : _random(random),
      _min(min),
      _max(double(max)+1.0-1.0e-6),
      _mult((double(max)+1.0-1.0e-6-double(min))/random.max())
    {
        assert(min<max);
        assert(max<=random.max());
    }

    GltRandomInteger(const uint32 min,const uint32 max)
    : _random(R::rng),
      _min(min),
      _max(double(max)+1.0-1.0e-6),
      _mult((double(max)+1.0-1.0e-6-double(min))/R::rng.max())
    {
        assert(min<max);
        assert(max<=_random.max());
    }

    GltRandomInteger(const GltRandomInteger<R> &gen)
    : _random(gen._random),
      _min(gen._min),
      _max(gen._max),
      _mult(gen._mult)
    {
    }

    ~GltRandomInteger()
    {
    }

    uint32 rand() const
    {
        return _min + (uint32) floor(_mult*double(_random.rand()));
    }

          R &base()       { return _random; }
    const R &base() const { return _random; }

private:

    R            &_random;
    const uint32  _min;
    const double  _max;
    const double  _mult;
};

template<class R = GltRandomLFSRMix>
class GltRandomSphere
{
public:

    GltRandomSphere(R &random,const double radius = 1.0)
    : _z(random,-1.0,1.0),
      _t(random,0,2*M_PI),
      _radius(radius)
    {
    }

    GltRandomSphere(const double radius = 1.0)
    : _z(R::rng,-1.0,1.0),
      _t(R::rng, 0.0,2*M_PI),
      _radius(radius)
    {
    }

    GltRandomSphere(const GltRandomSphere<R> &gen) 
    : _z(gen._z._random),
          _t(gen._t._random),
      _radius(gen._radius)
    {
    }

    ~GltRandomSphere() 
    {
    }

    Vector rand() const
    {
        const double z = _z.rand();
        const double t = _t.rand();
        const double r = _radius*sqrt(1.0-z*z);
        return Vector(r*cos(t),r*sin(t),_radius*z);
    }

          R &base()       { return _z.base(); }
    const R &base() const { return _z.base(); }

private:

    GltRandomDouble<R> _z,_t;
    const double    _radius;
};

template<class R = GltRandomLFSRMix>
class GltRandomOrientation
{
public:

    GltRandomOrientation(R &random)
    : _s(random), _a(random,0.0,2*M_PI)
    {
    }

    GltRandomOrientation()
    : _s(R::rng), _a(R::rng,0.0,2*M_PI)
    {
    }

    GltRandomOrientation(const GltRandomOrientation<R> &gen) 
    : _s(gen._s), _a(gen._a)
    {
    }

    ~GltRandomOrientation() 
    {
    }

    Matrix rand() const
    {
        // A is a vector distributed 
        // randomly on the unit sphere
        
        const Vector a(_s.rand());

        // B is orthogonal to A
        
        Vector b;

        switch (a.dominant())
        {
        case 0: b = xProduct(a,VectorZ); break; // a in x direction
        case 1: b = xProduct(a,VectorZ); break; // a in y direction
        case 2: b = xProduct(a,VectorX); break; // a in z direction
        }

        // C is orthogonal to both A and B

        Vector c(xProduct(b,a));

        // TODO - Possible to avoid this step?
        b.normalize();
        c.normalize();

        // Choose an angle in the AB plane

        const real angle = _a.rand();
        const real sina  = sin(angle);
        const real cosa  = cos(angle);

        return matrixOrient
            (
                a,
                sina*b +  cosa*c,
                cosa*b -  sina*c
            );
    }

          R &base()       { return _s.base(); }
    const R &base() const { return _s.base(); }

private:
    GltRandomSphere<R> _s;
    GltRandomDouble<R> _a;
};

#endif 

---