#include <iostream>
#include <cstdlib>
#include <ctime>
#include <cmath>
#include "PoissonRandom.h"
using namespace std;
// simulation de v.a.r. uniforme sur ]0,1[ :
inline double unif_rand() {
return (rand() + 0.5)/(RAND_MAX + 1.0);
}
class Alea {
public:
virtual double rand() const = 0;
};
class AleaDiscret : public Alea {
public:
virtual double loi(int n) const = 0;
virtual double rand() const;
};
double AleaDiscret::rand() const {
double y = unif_rand();
int x = 0;
double Fx = loi(0);
while ( y > Fx ) {
Fx += loi(++x);
}
return double(x);
}
class Poisson : public AleaDiscret {
double lambda;
double e_lbda;
public :
Poisson(double l = 1.0) : lambda(l) , e_lbda(exp(-l)) {}
virtual double loi(int n) const {
return e_lbda*pow(lambda,n)/tgamma(n+1);
}
// plus rapide que AleaDiscret::rand()
virtual double rand() const;
};
double Poisson::rand() const {
double x = 1;
int n = -1;
do {
x *= unif_rand();
++n;
} while ( x > e_lbda );
return double(n);
}
double moy(const Alea & a) {
const int N = 1000;
double s = 0.0;
for (int n = 0 ; n < N ; ++n)
s += a.rand();
return s/N;
}
void PoissonRandom(int durationNum, int poiNum, int totalNum) {
double total=0;
// double duration[100];
double *duration = new double[durationNum];
int tmp=0, num;
// Poisson poi(12);
Poisson poi(poiNum);
srand(int(time(0)));
do{
duration[tmp]=poi.rand();
total+=duration[tmp];
tmp++;
// }while(total<=50);
}while(total<=totalNum);
for(int i=1; i<tmp; i++){
duration[i]+=duration[i-1];
}
cout<<tmp<<endl;
for(int i=0; i<tmp; i++){
cout<<duration[i]<<endl;
}
delete[] duration;
}