1 files changed, 96 insertions, 0 deletions
diff --git a/rqmc.hpp b/rqmc.hpp
new file mode 100644
index 0000000..62efbaa
--- /dev/null
+++ b/rqmc.hpp
@@ -0,0 +1,96 @@
+#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 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)); 
+};
+
+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), sob(d) {};
+    std::vector<double> operator ()(){
+        std::vector<double> result = sob();
+        for(int i=0;i<d; i++){
+            result[i] = mean + gsl_cdf_gaussian_Pinv(result[i], std);
+        }
+        return value = result;
+    }
+    private: 
+    int d;
+    double mean, std;
+    sobol sob;
+};
+
+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;
+};
+
+
+
+
+
+