//#####################################################################
// Copyright 2002-2009, Robert Bridson, Ronald Fedkiw, Eran Guendelman, Geoffrey Irving, Neil Molino, Igor Neverov, Duc Nguyen, Craig Schroeder, Andrew Selle.
// This file is part of PhysBAM whose distribution is governed by the license contained in the accompanying file PHYSBAM_COPYRIGHT.txt.
//#####################################################################
// Class RANDOM_NUMBERS
//#####################################################################
#ifndef __RANDOM_NUMBERS__
#define __RANDOM_NUMBERS__

#include <Core/Arrays/ARRAYS_FORWARD.h>
#include <complex>
#include <random>
namespace PhysBAM{

template<class TV> class RANGE;
template<class TV> class ROTATION;
template<class TV> class FRAME;
template<class T,int d> class VECTOR;
template<class T,int d> class DIAGONAL_MATRIX;
template<class T,int d> class SYMMETRIC_MATRIX;
template<class T,class T_MATRIX> class MATRIX_BASE;
class TYPED_ISTREAM;
class TYPED_OSTREAM;

template<class T>
class RANDOM_NUMBERS
{
public:
    int gaussian_iset; // Used to force Get_Gaussian to reset
    T gset; // used internally by Get_Gaussian
    std::mt19937 random_number_generator;

    T Get_Number()
    {return (T)random_number_generator()*(1/(T)4294967296.0);} // in [0,1)

    int Get_Uniform_Integer(const int a,const int b) // in [a,b]
    {return random_number_generator()%(b+1-a)+a;}

    T Get_Uniform_Number(const T a,const T b)
    {return a+(b-a)*Get_Number();} // in [a,b)

    RANDOM_NUMBERS();
    explicit RANDOM_NUMBERS(const unsigned int seed);
    RANDOM_NUMBERS(const RANDOM_NUMBERS&) = delete;
    void operator=(const RANDOM_NUMBERS&) = delete;
    ~RANDOM_NUMBERS() = default;

    template<class T2,class T_ARRAY,class ID> void Random_Shuffle(ARRAY_BASE<T2,T_ARRAY,ID>& array)
    {T_ARRAY& derived=array.Derived();for(ID i(derived.Size()-1);i>=ID(0);i--) exchange(array(i),array(Get_Uniform_Integer(0,i)));}

    void Fill_Uniform(T& x,T a,T b)
    {x=Get_Uniform_Number(a,b);}

    void Fill_Uniform(std::complex<T>& x,T a,T b)
    {T r=Get_Uniform_Number(a,b),c=Get_Uniform_Number(a,b);x={r,c};}

    void Fill_Uniform(T& x,const INTERVAL<T>& a)
    {x=Get_Uniform_Number(a.min_corner,a.max_corner);}

    template<class T_OBJECT,class ...Args> auto
    Fill_Uniform(T_OBJECT& obj,Args&&...args)
        -> enable_if_t<sizeof((obj.Fill_Uniform(*this,std::forward<Args>(args)...),1))!=0>
    {obj.Fill_Uniform(*this,std::forward<Args>(args)...);}

    template<class T_OBJECT,class TV=typename T_OBJECT::VECTOR_T> auto
    Fill_Uniform(T_OBJECT& obj,T a,T b)
        -> enable_if_t<sizeof(Fill_Uniform(obj,TV(),TV()),1)!=0>
    {Fill_Uniform(obj,TV()+a,TV()+b);}

    template<class T_OBJECT,int d> auto
    Fill_Uniform(T_OBJECT& obj,const VECTOR<T,d>& a,const VECTOR<T,d>& b)
        -> enable_if_t<sizeof(Fill_Uniform(obj,RANGE<VECTOR<T,d>>()),1)!=0>
    {Fill_Uniform(obj,RANGE<VECTOR<T,d>>(a,b));}

//#####################################################################
    void Set_Seed(const unsigned int seed_input=time(0));
    T Get_Gaussian();
    template<class TV> TV Get_Vector_In_Unit_Sphere();
    template<class TV> TV Get_Direction();
    template<class RW> void Read(UNTYPED_ISTREAM<RW>& input);
    template<class RW> void Write(UNTYPED_OSTREAM<RW>& output) const;
//#####################################################################
};
}
#endif
