distrib.cpp

#include <math.h>
#include "distrib.h" 
#include <cstdlib>
#include <iostream>

using namespace std;

int32_t Sample(const Distrib& Chi){
    if (Chi.max) {
        double r = static_cast <double> (rand()) / static_cast <double> (RAND_MAX);
        for (int i = 0; i < Chi.max; ++i)
            if (r<= Chi.table[i]){
                return i - Chi.offset;                
            }
        exit(1);
    }
    
    
    double r, s = Chi.std_dev;
    if (s < 500){
        int32_t x, maxx = ceil(s*8);
        while(true) {
            x = rand() % (2*maxx + 1)  - maxx;
            r = static_cast <float> (rand()) / static_cast <float> (RAND_MAX);
            if (r < exp(- x*x / (2*s*s))){
                return x;
            }
        }
    }    
    
    // For some reason unknown to us, the previous implementation provides a bad distribution for large s...
    // We switch from "discrete gaussian" to rounded gaussian when s gets larger

    double x;
    
    while(true){
        x = static_cast <float> (rand()) / static_cast <float> (RAND_MAX);
        x = 16 *x -8;
        r = static_cast <float> (rand()) / static_cast <float> (RAND_MAX);
        
        if (r < exp(- x*x / 2 )){
            return floor(.5 + x*s);
        }
    }
}