/* Copyright (c) 1997-2007
   Ewgenij Gawrilow, Michael Joswig (Technische Universitaet Berlin, Germany)
   http://www.math.tu-berlin.de/polymake,  mailto:polymake@math.tu-berlin.de

   This program is free software; you can redistribute it and/or modify it
   under the terms of the GNU General Public License as published by the
   Free Software Foundation; either version 2, or (at your option) any
   later version: http://www.gnu.org/licenses/gpl.txt.

   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 General Public License for more details.
*/

#ifndef _POLYMAKE_RANDOM_GENERATORS_H
#define _POLYMAKE_RANDOM_GENERATORS_H "$Project: polymake $$Id: random_generators.h 7579 2007-01-22 11:34:46Z gawrilow $"

#include <cmath>
#include <cstdlib>

#ifndef _POLYMAKE_RATIONAL_H
#  include <Rational.h>
#endif
#ifndef _POLYMAKE_BITSET_H
#  include <Bitset.h>
#endif
#ifndef _POLYMAKE_VECTOR_H
#  include <Vector.h>
#endif

namespace pm {

unsigned long random_seed();

template <typename Generator>
class random_get_iterator {
protected:
   Generator* generator;
public:
   typedef input_iterator_tag iterator_category;
   typedef typename Generator::value_type value_type;
   typedef typename Generator::reference reference;
   typedef typename arrow_helper<Generator>::pointer pointer;
   typedef ptrdiff_t difference_type;

   random_get_iterator(const Generator* gen_arg=0) : generator(const_cast<Generator*>(gen_arg)) { }

   reference operator* () const { return generator->get(); }
   pointer operator-> () const { return arrow_helper<Generator>::arrow(this); }

   random_get_iterator& operator++() { return *this; }
   random_get_iterator& operator++(int) { return *this; }
   bool operator== (const random_get_iterator& it) const { return generator==it.generator; }
   bool operator!= (const random_get_iterator& it) const { return !operator==(it); }
};

/*
 *  The classes defined below comply with the STL container concept.
 *  To obtain the random numbers, you should either call the get() method,
 *  or create an iterator and dereference it.
 */

/// the general case is undefined
template <typename T> class UniformRNG;

/// Wrapper for lrand48
template <>
class UniformRNG<long> {
protected:
   unsigned short state[3];
public:
   typedef long value_type;
   typedef const long const_reference;
   typedef const_reference reference;

   UniformRNG() { reset(random_seed()); }
   explicit UniformRNG(unsigned long seed) { reset(seed); }

   void reset(unsigned long seed)
   {
      if (sizeof(short)*3 < sizeof(long)) {
	 state[0]=(reinterpret_cast<const unsigned short*>(&seed))[0];
	 state[1]=(reinterpret_cast<const unsigned short*>(&seed))[1];
	 state[2]=(reinterpret_cast<const unsigned short*>(&seed))[2];
      } else {
	 state[2]=12345;
	 reinterpret_cast<long&>(state[0])=seed;
      }
   }

   const_reference get() { return nrand48(state); }

   typedef random_get_iterator<UniformRNG> iterator;
   typedef iterator const_iterator;

   iterator begin() const { return iterator(this); }
   iterator end() const { return iterator(); }
};

/// Wrapper for drand48
template <>
class UniformRNG<double> : public UniformRNG<long> {
public:
   typedef double value_type;
   typedef const double const_reference;
   typedef const_reference reference;

   UniformRNG() {}
   explicit UniformRNG(unsigned long seed) : UniformRNG<long>(seed) {}

   const_reference get() { return erand48(state); }

   typedef random_get_iterator<UniformRNG> iterator;
   typedef iterator const_iterator;

   iterator begin() const { return iterator(this); }
   iterator end() const { return iterator(); }
};

/** Generator of floating-point numbers normally distributed in (-1,1).

    The algorithm is taken from Donald E. Knuth, The Art of Computer Programming, vol. II, p117.
*/
class NormalRNG : public UniformRNG<double> {
   typedef UniformRNG<double> super;
public:
   typedef const double& const_reference;
   typedef const_reference reference;

   NormalRNG() { fill(); }
   explicit NormalRNG(unsigned long seed) : super(seed) { fill(); }

   reference get()
   {
      if (++index==2) fill();
      return x[index];
   }

   typedef random_get_iterator<NormalRNG> iterator;
   typedef iterator const_iterator;

   iterator begin() const { return iterator(this); }
   iterator end() const { return iterator(); }
private:
   double x[2];
   int index;
   void fill();
};

/** Generator of random integer numbers uniformly distributed in [0, 2^k).

    This is a wrapper around GMP (GNU Multiple Precision Library) random functions.
    Developed and tested with GMP Version 3.1.1 and 4.0
    See the GMP Home Pages at `http://www.swox.com/gmp/'.
*/
template <>
class UniformRNG<Integer> {
protected:
   gmp_randstate_t state;
   const unsigned long bitlength;

public:
   typedef Integer value_type;
   typedef Integer reference;
   typedef reference const_reference;

   explicit UniformRNG(unsigned long bitlength_arg=48)
      : bitlength(bitlength_arg)
   {
      gmp_randinit(state, GMP_RAND_ALG_LC, bitlength);
      reset(random_seed());
   }

   UniformRNG(unsigned long bitlength_arg, unsigned long seed)
      : bitlength(bitlength_arg)
   {
      gmp_randinit(state, GMP_RAND_ALG_LC, bitlength);
      reset(seed);
   }

   UniformRNG(unsigned long bitlength_arg, const Integer& seed)
      : bitlength(bitlength_arg)
   {
      gmp_randinit(state, GMP_RAND_ALG_LC, bitlength);
      reset(seed);
   }

   UniformRNG(const UniformRNG& rg)
      : bitlength(rg.bitlength)
   {
      gmp_randinit(state, GMP_RAND_ALG_LC, bitlength);
#if __GNU_MP_VERSION >= 4
      gmp_randseed(state, rg.state[0]._mp_seed);
#else
      gmp_randseed(state, const_cast<__mpz_struct*>(rg.state[0].seed));	// cast due to an inaccurate declaration in GMP
#endif
   }

   void reset(unsigned long seed) { gmp_randseed_ui(state, seed); }

   void reset(const Integer& seed)
   {
#if __GNU_MP_VERSION >= 4
      gmp_randseed(state, seed.get_rep());
#else
      gmp_randseed(state, const_cast<__mpz_struct*>(seed.get_rep()));	// cast due to an inaccurate declaration in GMP
#endif
   }

   ~UniformRNG() { gmp_randclear(state); }

   reference get() { return Integer(mpz_urandomb, state, bitlength); }

   typedef random_get_iterator<UniformRNG> iterator;
   typedef iterator const_iterator;
   iterator begin() const { return iterator(this); }
   iterator end() const { return iterator(); }
};

/** Generator of random rational numbers uniformly distributed in [0, 1).

    Only the numerators of the numbers are really random.
    The denominators are always some powers of 2.
*/
template <>
class UniformRNG<Rational> : public UniformRNG<Integer> {
   typedef UniformRNG<Integer> super;
public:
   typedef Rational value_type;
   typedef Rational reference;
   typedef reference const_reference;

   explicit UniformRNG(unsigned long bitlength_arg=48)
      : super(bitlength_arg) {}
   UniformRNG(unsigned long bitlength_arg, unsigned long seed)
      : super(bitlength_arg,seed) {}
   UniformRNG(unsigned long bitlength_arg, const Integer& seed)
      : super(bitlength_arg,seed) {}

   reference get()
   {
      Rational r=super::get();
      r>>=bitlength;
      return r;
   }

   typedef random_get_iterator<UniformRNG> iterator;
   typedef iterator const_iterator;
   iterator begin() const { return iterator(this); }
   iterator end() const { return iterator(); }
};

/// Generator of random Bitset of a given maximal cardinality
template <>
class UniformRNG<Bitset> : public UniformRNG<Integer> {
   typedef UniformRNG<Integer> super;
public:
   typedef Bitset value_type;
   typedef Bitset reference;
   typedef reference const_reference;

   explicit UniformRNG(int max_elem)
      : super(max_elem,random_seed()) {}
   UniformRNG(int max_elem, unsigned long seed)
      : super(max_elem,seed) {}
   UniformRNG(int max_elem, const Integer& seed)
      : super(max_elem,seed) {}

   reference get() { return Bitset(mpz_urandomb, state, bitlength); }

   typedef random_get_iterator<UniformRNG> iterator;
   typedef iterator const_iterator;
   iterator begin() const { return iterator(this); }
   iterator end() const { return iterator(); }
};

/// Generator of uniformly distributed random points on the unit sphere in R^d
template <>
class UniformRNG< Vector<double> >
   : protected NormalRNG {
public:
   typedef Vector<double> value_type;
   typedef const value_type& reference;
   typedef reference const_reference;

   UniformRNG(int dim)
      : NormalRNG(random_seed()), point(dim) {}
   UniformRNG(int dim, unsigned long seed)
      : NormalRNG(seed), point(dim) {}

   reference get()
   {
      fill_vector();
      return point;
   }

   typedef random_get_iterator<UniformRNG> iterator;
   typedef iterator const_iterator;
   iterator begin() const { return iterator(this); }
   iterator end() const { return iterator(); }

protected:
   value_type point;

   void fill_vector();
};

class DiscreteRNG {
public:
   typedef int value_type;
   typedef int reference;
   typedef int const_reference;

   template <typename Container>
   DiscreteRNG(const Container& distrib_src,
	       typename enable_if<void*, identical<typename object_traits<Container>::model,
	                                          is_container>::value>::type=0)
      : distribution(distrib_src)
   {
      normalize();
   }

   template <typename Container>
   DiscreteRNG(const Container& distrib_src, unsigned long seed,
	       typename enable_if<void*, identical<typename object_traits<Container>::model,
	                                          is_container>::value>::type=0)
      : rg(seed), distribution(distrib_src) 
   {
      normalize();
   }

   int get() const;
protected:
   mutable UniformRNG<double> rg;
   Vector<double> distribution;
private:
   void normalize();
};


template <typename Generator>
struct check_iterator_feature<random_get_iterator<Generator>, unlimited> : True { };

template <typename T>
struct spec_object_traits< UniformRNG<T> > : spec_object_traits<is_container> { };

template <>
struct spec_object_traits< NormalRNG > : spec_object_traits<is_container> { };

} // end namespace pm

namespace polymake {
   using pm::UniformRNG;
   using pm::NormalRNG;
   using pm::DiscreteRNG;
   using pm::random_seed;
}

#endif // _POLYMAKE_RANDOM_GENERATORS_H

// Local Variables:
// mode:C++
// c-basic-offset:3
// End:
