/* -*- mode:c++; coding: koi8-r -*- */

/* $Id: Random.cpp,v 1.1 2003/04/14 20:53:07 cher Exp $ */
/* Copyright (C) 2003 Alexander Chernov <cher@unicorn.cmc.msu.ru> */

/*
 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 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.

 See the `COPYING' file for the full terms and conditions.
*/

#include "Random.hpp"
#include "../double/cephes.h"

#include <math.h>

Random::Random(EnthropySource *_src)
  throw (BadArgsError)
{
  if (!_src) throw BadArgsError();
  src = _src;
}

Random::~Random()
  throw ()
{
  delete src;
}

double
Random::uniform()
  throw ()
{
  return src->getDouble();
}

double
Random::uniform(double a, double b)
  throw (BadArgsError)
{
  if (a >= b) throw BadArgsError();
  return this->uniform() * (b - a) + a;
}

int
Random::uniform(int a, int b)
  throw (BadArgsError)
{
  if (a >= b) throw BadArgsError();
  return int(this->uniform() * (b - a) + a);
}

double
Random::normal()
  throw ()
{
  return cephes_ndtri(src->getDouble());
}

double
Random::normal(double m, double sigma)
  throw (BadArgsError)
{
  if (sigma <= 0) throw BadArgsError();

  return this->normal() * sigma + m;
}

double
Random::exponent(double lambda)
  throw (BadArgsError)
{
  if (lambda <= 0) throw BadArgsError();
  return -lambda * log(src->getDouble());
}

double
Random::beta(double a, double b)
  throw (BadArgsError)
{
  double r1, r2, s1, s2;

  if (a <= 0 || b <= 0) throw BadArgsError();

  do {
    r1 = src->getDouble();
    r2 = src->getDouble();
    s1 = pow(r1, 1/a);
    s2 = pow(r2, 1/b);
  } while (s1 + s2 > 1);

  return s1 / (s1 + s2);
}

int
Random::bernoully(double p)
  throw (BadArgsError)
{
  if (p < 0.0 || p > 1.0) throw BadArgsError();
  if (p == 0.0) return 0;
  if (p == 1.0) return 1;
  if (src->getDouble() < p) return 1;
  return 0;
}

int
Random::binomial(int n, double p)
  throw (BadArgsError)
{
  int k = 0;

  if (p < 0.0 || p > 1.0 || n <= 0) throw BadArgsError();
  if (p == 0.0) return 0;
  if (p == 1.0) return 1;
  for (int i = 0; i < n; i++)
    k += this->bernoully(p);
  return k;
}

/*
 * Local variables:
 *  compile-command: "make -C .."
 *  c++-font-lock-extra-types: ("[A-Z]\\sw*[a-z]\\sw*")
 * End:
 */