random.h

#pragma once
#include <algorithm>
#include <ctime>
#include <cmath>
#include <set>
 
#include "constants.h"
#include "functions.h"
#include "vec.h"
#include "fvec.h"
#include "mat.h"
#include "eig.h"
#include "diag_mat.h"
#include "perm_mat.h"
#include "up_tri_mat.h"
#include "low_tri_mat.h"
#include "quaternion.h"
 
namespace cgv {
namespace math {
 
struct random
{
private:
    unsigned long long u,v,w;
    double storedval;
public:
    void set_seed(unsigned long long seed)
    {
        w = 1;
        v = 4101842887655102017LL;
        u = seed^v;
        unsigned long long ull;
        uniform(ull);
        v = u; uniform(ull);
        w = v; uniform(ull);
    }
 
    random() : storedval(0.0)
    {
        unsigned long long seed = clock();
        set_seed(seed);
    }
 
    random(unsigned long long seed) : storedval(0.0)
    {
        set_seed(seed);
    }
    
    void uniform(unsigned long long& rv)
    {
        u = u*2862933555777941757LL + 7046029254386353087LL;
        v ^= v >> 17; v ^= v <<31; v ^= v >> 8;
        w = 4294957665U*(w & 0xffffffff) + (w >> 32);
        unsigned long long x = u ^(u << 21); x ^= x >>35; x^=x << 4;
        rv= (x+v)^w;
    }
 
    void uniform(unsigned int& rv)
    {
        unsigned long long v;
        uniform(v);
        rv= (unsigned int) v;
    }
 
    void uniform(int& rv)
    {
        unsigned long long v;
        uniform(v);
        rv = (int) v;
    }
 
    void uniform(long& rv)
    {
        unsigned long long v;
        uniform(v);
        rv = (long) v;
    }
 
    void uniform(long long& rv)
    {
        unsigned long long v;
        uniform(v);
        rv = (long long) v;
    }
 
    void uniform(bool& rv)
    {
        unsigned long long v;
        uniform(v);
        rv= v%2 == 0;
    }
 
    void uniform(const unsigned min,const unsigned max, unsigned& rv)
    {
        uniform(rv);
        rv =rv%(max-min+1)+min;
    }
 
    void uniform(const int min,const int max, int& rv)
    {
        uniform(rv);
        rv =(int)(rv%(unsigned)(max-min+1))+min;
    }
 
 
    void uniform(double& rv)
    {
        unsigned long long v;
        uniform(v);
        rv= 5.42101086242752217E-20 * v;
    }
 
    void uniform(const double min,const double max, double &rv)
    {
        double v;
        uniform(v);
        rv= (max-min)*v + min;
    }
    
    void uniform(float & rv)
    {
        double v;
        uniform(v);
        rv =(float) v;
    }
 
    void uniform(const float min,const float max, float& rv)
    {
        float f;
        uniform(f);
        rv= (max-min)*f + min;
    }
 
    void uniform( vec<float>& rv)
    {
        for(unsigned i = 0; i < rv.size();i++)
            uniform(rv(i)); 
    }
    
    void uniform(const float min,const float max, vec<float>& rv)
    {
        for(unsigned i = 0; i < rv.size();i++)
            uniform(min,max,rv(i));
        
    }
 
    
    void uniform(vec<double>& rv)
    {
        for(unsigned i = 0; i < rv.size();i++)
            uniform(rv(i)); 
    }
    
    void uniform(const double min,const double max, vec<double>& rv)
    {
        for(unsigned i = 0; i < rv.size();i++)
            uniform(min,max,rv(i)); 
    }
 
    template<typename T, cgv::type::uint32_type N>
    void uniform(cgv::math::fvec<T, N>& rv)
    {
        for(unsigned i = 0; i < N; i++)
            uniform(rv(i));
    }
 
    template<typename T, cgv::type::uint32_type N>
    void uniform(const float min, const float max, cgv::math::fvec<T, N>& rv)
    {
        for(unsigned i = 0; i < N; i++)
            uniform(min, max, rv(i));
    }
 
    void uniform(mat<float>& rv)
    {
        for(unsigned i = 0; i < rv.nrows(); i++)
            for(unsigned j = 0; j < rv.ncols(); j++)
            uniform(rv(i,j));   
    }
    
    void uniform(const float min,const float max, mat<float>& rv)
    {
        for(unsigned i = 0; i < rv.nrows(); i++)
            for(unsigned j = 0; j < rv.ncols(); j++)
                uniform(min,max,rv(i,j));       
    }
 
    void uniform(perm_mat& rv)
    {
        unsigned j;
        for(unsigned i = 0; i < rv.size();i++)
        {
            uniform(j);
            j=(i+j)%rv.size();
            rv.swap(i,j);
        }
    }
 
    void uniform(up_tri_mat<float>& rv)
    {
        for(unsigned i = 0; i < rv.nrows(); i++)
            for(unsigned j = i; j < rv.ncols(); j++)
            uniform(rv(i,j));   
    }
 
    
 
    void uniform(const float min, const float max,up_tri_mat<float>& rv)
    {
        for(unsigned i = 0; i < rv.nrows(); i++)
            for(unsigned j = i; j < rv.ncols(); j++)
            uniform(min,max,rv(i,j));   
    }
 
    void uniform(up_tri_mat<double>& rv)
    {
        for(unsigned i = 0; i < rv.nrows(); i++)
            for(unsigned j = i; j < rv.ncols(); j++)
            uniform(rv(i,j));   
    }
 
    void uniform(const double min, const double max,up_tri_mat<double>& rv)
    {
        for(unsigned i = 0; i < rv.nrows(); i++)
            for(unsigned j = i; j < rv.ncols(); j++)
            uniform(min,max,rv(i,j));   
    }
 
    void uniform( low_tri_mat<float>& rv)
    {
        for(unsigned i = 0; i < rv.nrows(); i++)
            for(unsigned j = 0; j <= i; j++)
                uniform(rv(i,j));       
    }
    
    void uniform(const float min,const float max, low_tri_mat<float>& rv)
    {
        for(unsigned i = 0; i < rv.nrows(); i++)
            for(unsigned j = 0; j <= i; j++)
                uniform(min,max,rv(i,j));       
    }
 
    void uniform( low_tri_mat<double>& rv)
    {
        for(unsigned j = 0; j < rv.ncols(); j++)
            for(unsigned i = 0; i < j; i++)
                uniform(rv(i,j));       
    }
    
    void uniform(const double min,const double max, low_tri_mat<double>& rv)
    {
        for(unsigned j = 0; j < rv.ncols(); j++)
            for(unsigned i = 0; i < j; i++)
                uniform(min,max,rv(i,j));       
    }
 
    void uniform(diag_mat<float>& rv)
    {
        for(unsigned i = 0; i < rv.nrows(); i++)
            uniform(rv(i)); 
    }
    
    void uniform(const float min,const float max, diag_mat<float>& rv)
    {
            for(unsigned i = 0; i < rv.nrows(); i++)
                uniform(min,max,rv(i));     
    }
 
    void uniform(diag_mat<double>& rv)
    {
        for(unsigned i = 0; i < rv.nrows(); i++)
            uniform(rv(i)); 
    }
    
    void uniform(const double min,const double max, diag_mat<double>& rv)
    {
            for(unsigned i = 0; i < rv.nrows(); i++)
                uniform(min,max,rv(i));     
    }
 
 
    void uniform(mat<double>& rv)
    {
        for(unsigned i = 0; i < rv.nrows(); i++)
            for(unsigned j = 0; j < rv.ncols(); j++)
                uniform(rv(i,j));   
    }
    
    void uniform(const double min,const double max, mat<double>& rv)
    {
        for(unsigned  i = 0; i < rv.nrows(); i++)
            for(unsigned j = 0; j < rv.ncols(); j++)
                uniform(min,max,rv(i,j));       
    }
 
 
 
    void normal(double& rv)
    {
        double v1,v2,rsq,fac;
        if(storedval == 0.)
        {
            do{
                uniform(v1);
                uniform(v2);
 
                v1 = 2.0 * v1-1.0;
                v2 = 2.0 * v2-1.0;
                rsq = v1*v1+v2*v2;
            }while(rsq >= 1.0 || rsq == 0.0);
            fac = sqrt(-2.0*std::log(rsq)/rsq);
            storedval = v1*fac;
            rv= v2*fac;
 
        }else
        {
            fac = storedval;
            storedval = 0.;
            rv= fac;
        }
    }
 
    void normal(float& rv)
    {
        double d;
        normal(d);
        rv = (float)d;
    }
 
    void normal(const double mu, const double sigma, double &rv)
    {   
        double d;
        normal(d);
        rv =mu + sigma * d;
    }
 
    void normal(const float mu, const float sigma, float& rv)
    {   float f;
        normal(f);
        rv= mu + sigma * f;
    }
 
    void uniform_direction(vec<float>& v)
    {
        
        for(unsigned i = 0;i < v.size(); i++)
             normal(v(i));
        v.normalize();
    }
 
    
    void uniform_direction(vec<double>& v)
    {
        for(unsigned i = 0;i < v.size();i++)
             normal(v(i));
        v.normalize();
    }
    
    template<typename T>
    void uniform_quat_orientation(quaternion<T>& q) {
        T x[3];
        uniform(x[0]);
        uniform(x[1]);
        uniform(x[2]);
        T a, b, c, d, s;
        T z, r, theta, omega;
 
        z = x[0];
        r = sqrt(1 - z * z);
        theta = static_cast<T>(2.0 * cgv::math::constants::pi) * x[1];
        omega = static_cast<T>(cgv::math::constants::pi) * x[2];
 
        s = sin(omega);
        a = cos(omega);
        b = s * cos(theta) * r;
        c = s * sin(theta) * r;
        d = s * z;
 
        q.set(b, c, d, a);
    }
    
    template<typename T>
    void uniform_orientation(fmat<T, 3, 3>& m)
    {
        quaternion<T> q;
        uniform_quat_orientation(q);
        q.put_matrix(m);
        return;
    }
 
    template<typename T>
    void uniform_orientation(fmat<T, 4, 4>& m)
    {
        quaternion<T> q;
        uniform_quat_orientation(q);
        q.put_homogeneous_matrix(m);
        return;
    }
 
    void uniform_point_in_unit_box(vec<double>& p)
    {
        if (p.size() == 0)
            p.resize(3);
        for (unsigned i=0; i<p.size(); ++i)
            uniform(p(i));
    }
 
    void uniform_point_in_box(const vec<double>& minp,const vec<double>& maxp, vec<double>& p)
    {
        assert(minp.size() == maxp.size());
        if (p.size() != minp.size())
            p.resize(minp.size());
        for (unsigned i=0; i<p.size(); ++i)
            uniform(minp(i),maxp(i),p(i));
    }
 
    void uniform_point_in_unit_ball(vec<double>& p)
    {
        uniform_direction(p);
        double r;
        uniform(r);
        r = std::pow(r, 1.0/p.size());
        p *= r;
    }
    void uniform_point_in_triangle(const vec<double>& p1, const vec<double>& p2, const vec<double>& p3,vec<double> p)
    {
        double u,v;
        uniform(u);
        uniform(v);
        if(u+v > 1.0)
        {
            u=1.0-u;
            v=1.0-v;
        }
        p= u*p1 + v*p2 +(1.0-u-v)*p3;
    }
 
    void uniform_point_on_sphere(const vec<double>& center,const double& radius,
        vec<double>& p)
    {
        for(unsigned i = 0;i < p.size();i++)
             normal(p(i));
        p.normalize();
        p=center+radius*p;
    }
 
 
    void uniform_nchoosek(unsigned n, unsigned k, vec<unsigned>& indices)
    {
        std::set<unsigned> s;   
        unsigned v;
        while(s.size() < k)
        {
            uniform(0,n-1,v);
            s.insert(v);
        }
 
        indices.resize(k);
        unsigned i=0;
        for(std::set<unsigned>::iterator it =s.begin(); it != s.end();it++,i++)
            indices(i)=*it;
            
    }
 
};
 
 
} // namespace math
} // namespace cgv