1 files changed, 127 insertions, 0 deletions
diff --git a/src/rqmc.hpp b/src/rqmc.hpp
new file mode 100644
index 0000000..1b33713
--- /dev/null
+++ b/src/rqmc.hpp
@@ -0,0 +1,127 @@
+#include "low_discrepancy.hpp"
+#include <vector>
+#include <gsl/gsl_cdf.h>
+#include "var_alea.hpp"
+
+
+double frac_part(double x);
+
+//fonctions de compositions de monte_carlo.hpp
+
+template <class _Result>
+struct generator 
+{
+	typedef _Result result_type;
+};
+
+template <typename Fct, typename VA>
+struct compose_t : public generator< typename Fct::result_type >
+{
+	compose_t(Fct f, VA X) : f(f), X(X) {};
+	typename Fct::result_type operator()() { 
+		return f(X()); 
+	};
+	private:
+		Fct f; VA X;
+};
+
+template <typename Fct, typename VA>
+inline compose_t<Fct, VA>
+compose(Fct f, VA X) {
+    return compose_t<Fct, VA>(f, X);
+};
+
+template <typename VA>
+struct mean_t: public generator<double>
+{
+    mean_t(int n, VA X): n(n), X(X) {};
+    double operator()(){
+        double sum = 0;
+        for(int i=0; i<n; i++){
+            sum+=X();
+        }
+        return sum/n;
+    }
+    private :
+        int n; VA X;
+    };
+    
+template <typename VA>
+inline mean_t<VA>
+mean(int n, VA X){
+    return mean_t<VA>(n, X);
+};
+
+
+//Les classes de monte-Carlo
+
+template <typename L>
+std::vector<double> monte_carlo(int n, L X)
+{
+	std::vector<double> result(3,0);
+	double x;
+	for (int j = 0; j < n; j++) {
+		x = X();
+		result[0] += x;
+		result[1] += x*x;
+	}
+	result[0] /= (double) n;
+	result[1] = (result[1] - n*result[0]*result[0])/(double)(n-1);
+	result[2] = 1.96*sqrt(result[1]/(double) n);
+	return result;
+};
+
+template <typename L, typename Fct>
+std::vector<double> monte_carlo(int n, Fct f, L X)
+{
+	return monte_carlo(n, compose(f, X)); 
+};
+
+//Les classes de quasi-mean
+
+template <typename LDS>
+struct quasi_gaussian : public var_alea<std::vector<double> >
+{
+    quasi_gaussian (int d, double mean = 0, double std =1)
+        : d(d), mean(mean), std(std), s(d), U(0,1), seed(d) {
+            for(int i=0; i<d; i++) {seed[i]=U();}
+                };
+    std::vector<double> operator ()(){
+        std::vector<double> result = s();
+        for(int i=0;i<d; i++){
+            result[i] = mean + gsl_cdf_gaussian_Pinv(frac_part(result[i]+seed[i]), std);
+        }
+        return value = result;
+    }
+    private: 
+    int d;
+    double mean, std;
+    LDS s;
+    uniform U;
+    std::vector<double> seed;
+};
+
+
+
+struct gaussian_d : public var_alea<std::vector<double> >
+{
+    gaussian_d(int d, double mean = 0, double std =1)
+        : d(d), mean(mean), std(std), G(mean, std) {};
+    std::vector<double> operator ()(){
+        std::vector<double> result(d);
+        for(int i=0;i<d; i++){
+            result[i] = G();
+        }
+        return value = result;
+    }
+    private: 
+    int d;
+    double mean, std;
+    gaussian G;
+};
+
+
+
+
+
+