RhodeCode

    /// \todo Consider using the "ziggurat" method instead.

    double gauss() 

    {

      double V1,V2,S;

      do {

	V1=2*real<double>()-1;

	V2=2*real<double>()-1;

	S=V1*V1+V2*V2;

      } while(S>=1);

      return std::sqrt(-2*std::log(S)/S)*V1;

    }

    /// Gauss distribution with given mean and standard deviation

    /// Gauss distribution with given mean and standard deviation.

    /// \sa gauss()

    double gauss(double mean,double std_dev)

    {

      return gauss()*std_dev+mean;

    }

    /// Exponential distribution with given mean

    /// This function generates an exponential distribution random number

    /// with mean <tt>1/lambda</tt>.

    ///

    double exponential(double lambda=1.0)

    {

      return -std::log(1.0-real<double>())/lambda;

    }

    /// Gamma distribution with given integer shape

    /// This function generates a gamma distribution random number.

    /// 

    ///\param k shape parameter (<tt>k>0</tt> integer)

    double gamma(int k) 

    {

      double s = 0;

      for(int i=0;i<k;i++) s-=std::log(1.0-real<double>());

      return s;

    }

    /// Gamma distribution with given shape and scale parameter

    /// This function generates a gamma distribution random number.

    /// 

    ///\param k shape parameter (<tt>k>0</tt>)

    ///\param theta scale parameter

    ///

    double gamma(double k,double theta=1.0)

    {

      double xi,nu;

      const double delta = k-std::floor(k);

      const double v0=E/(E-delta);

      do {

	double V0=1.0-real<double>();

	double V1=1.0-real<double>();

	double V2=1.0-real<double>();

	if(V2<=v0) 

	  {

	    xi=std::pow(V1,1.0/delta);

	    nu=V0*std::pow(xi,delta-1.0);

	  }

	else 

	  {

	    xi=1.0-std::log(V1);

	    nu=V0*std::exp(-xi);

	  }

      } while(nu>std::pow(xi,delta-1.0)*std::exp(-xi));

      return theta*(xi-gamma(int(std::floor(k))));

    }

    /// Weibull distribution

    /// This function generates a Weibull distribution random number.

    /// 

    ///\param k shape parameter (<tt>k>0</tt>)

    ///\param lambda scale parameter (<tt>lambda>0</tt>)

    ///

    double weibull(double k,double lambda)

    {

      return lambda*pow(-std::log(1.0-real<double>()),1.0/k);

    }  

    /// Pareto distribution

    /// This function generates a Pareto distribution random number.

    /// 

    ///\param k shape parameter (<tt>k>0</tt>)

    ///\param x_min location parameter (<tt>x_min>0</tt>)

    ///

    double pareto(double k,double x_min)

    {

      return exponential(gamma(k,1.0/x_min));

    }  

    ///@}

    ///\name Two dimensional distributions

    ///

    ///@{

    /// Uniform distribution on the full unit circle

    /// Uniform distribution on the full unit circle.

    ///

    dim2::Point<double> disc() 

    {

      double V1,V2;

      do {

	V1=2*real<double>()-1;

	V2=2*real<double>()-1;

      } while(V1*V1+V2*V2>=1);

      return dim2::Point<double>(V1,V2);

    }

    /// A kind of two dimensional Gauss distribution

    /// This function provides a turning symmetric two-dimensional distribution.

    /// Both coordinates are of standard normal distribution, but they are not

    /// independent.

    ///

    /// \note The coordinates are the two random variables provided by

    /// the Box-Muller method.

    dim2::Point<double> gauss2()

    {

      double V1,V2,S;

      do {

	V1=2*real<double>()-1;

	V2=2*real<double>()-1;

	S=V1*V1+V2*V2;

      } while(S>=1);

      double W=std::sqrt(-2*std::log(S)/S);

      return dim2::Point<double>(W*V1,W*V2);

    }

    /// A kind of two dimensional exponential distribution

    /// This function provides a turning symmetric two-dimensional distribution.

    /// The x-coordinate is of conditionally exponential distribution

    /// with the condition that x is positive and y=0. If x is negative and 

    /// y=0 then, -x is of exponential distribution. The same is true for the

    /// y-coordinate.

    dim2::Point<double> exponential2() 

    {

      double V1,V2,S;

      do {

	V1=2*real<double>()-1;

	V2=2*real<double>()-1;

	S=V1*V1+V2*V2;

      } while(S>=1);

      double W=-std::log(S)/S;

      return dim2::Point<double>(W*V1,W*V2);

    }

    ///@}    

  };

  extern Random rnd;

}

#endif

/* -*- C++ -*-

 *

 * This file is a part of LEMON, a generic C++ optimization library

 *

 * Copyright (C) 2003-2008

 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport

 * (Egervary Research Group on Combinatorial Optimization, EGRES).

 *

 * Permission to use, modify and distribute this software is granted

 * provided that this copyright notice appears in all copies. For

 * precise terms see the accompanying LICENSE file.

 *

 * This software is provided "AS IS" with no warranty of any kind,

 * express or implied, and with no claim as to its suitability for any

 * purpose.

 *

 */

#include <lemon/random.h>

#include "test_tools.h"

///\file \brief Test cases for random.h

///

///\todo To be extended

///

int main()

{

  double a=lemon::rnd();

  check(a<1.0&&a>0.0,"This should be in [0,1)");

  a=lemon::rnd.gauss();

  a=lemon::rnd.gamma(3.45,0);

  a=lemon::rnd.gamma(4);

  //Does gamma work with integer k?

  a=lemon::rnd.gamma(4.0,0);

}