/* Copyright (c) 1997-2004
   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,v 1.14 2004/06/22 16:22:42 gawrilow Exp $"

#if defined(__GNUC__) && __GNUC_MINOR__<3 && !defined(__APPLE__) && !defined(__COMO__)
#pragma interface
#endif

#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 {

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(Generator* gen_arg=0) : 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 <class 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;

   explicit UniformRNG(unsigned long seed=0L) { reset(seed); }

   void reset(unsigned long seed=0L) {
      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() { return iterator(this); }
   iterator end() { 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;

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

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

   typedef random_get_iterator<UniformRNG> iterator;
   typedef iterator const_iterator;

   iterator begin() { return iterator(this); }
   iterator end() { 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;

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

   typedef random_get_iterator<NormalRNG> iterator;
   typedef iterator const_iterator;

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

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

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;

   /** @param seed start value for the generator
       @param bitlength maximal size in bits of random numbers to be produced
   */
   explicit UniformRNG(unsigned long seed=0L, unsigned long bitlength_arg=48)
      : bitlength(bitlength_arg) {
      gmp_randinit(state, GMP_RAND_ALG_LC, bitlength);
      reset(seed);
   }

   /// The same with an \\Integer\\ seed value.
   explicit UniformRNG(const Integer& seed, unsigned long bitlength_arg=48)
      : 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() { return iterator(this); }
   iterator end() { 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;

   /** @param seed start value for the generator
       @param bitlength maximal size in bits of numerators and denominators
   */
   explicit UniformRNG(unsigned long seed=0L, unsigned long bitlength_arg=48)
      : super(seed,bitlength_arg) { }

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

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

   typedef random_get_iterator<UniformRNG> iterator;
   typedef iterator const_iterator;
   iterator begin() { return iterator(this); }
   iterator end() { 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, unsigned long seed=0L)
      : super(seed,max_elem) { }
   UniformRNG(int max_elem, const Integer& seed)
      : super(seed,max_elem) { }

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

   typedef random_get_iterator<UniformRNG> iterator;
   typedef iterator const_iterator;
   iterator begin() { return iterator(this); }
   iterator end() { 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, unsigned long seed=0L)
      : NormalRNG(seed), point(dim) { }

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

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

protected:
   value_type point;

   void fill_vector();
};

template <class 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;
}

#endif // _POLYMAKE_RANDOM_GENERATORS_H

// Local Variables:
// mode:C++
// End:
