un peu de nettoyage

5 files changed, 254 insertions, 28 deletions

diff --git a/Makefile b/Makefile
index 3e093eb..474ef18 100644
--- a/Makefile
+++ b/Makefile
@@ -17,7 +17,7 @@ test: test.o mt19937.o
 rqmc: rqmc.o mt19937.o
 	$(CXX) $^ -o $@ $(GSL_FLAGS) 
 
-option: option.o mt19937.o
+option: option.o mt19937.o p_adic.o
 	$(CXX) $^ -o $@ $(GSL_FLAGS) 
 
 clean:
diff --git a/low_discrepancy.hpp b/low_discrepancy.hpp
new file mode 100644
index 0000000..6b3c1bd
--- /dev/null
+++ b/low_discrepancy.hpp
@@ -0,0 +1,205 @@
+#ifndef LOW_DISCREPANCY_HPP
+#define LOW_DISCREPANCY_HPP
+#include <iostream>
+#include <cmath>
+#include <climits>
+#include <list>
+#include <vector>
+#include <algorithm>
+#include <numeric>
+#include <gsl/gsl_qrng.h>
+
+class p_adic {
+	public: 
+		typedef std::list<int> coeff;
+		p_adic(coeff ak, coeff pk) : ak(ak), pk(pk), p(*(++pk.begin())) {};
+		p_adic(int n, int p = 2);
+		p_adic(double x, int p = 2);
+		operator double() {
+			return std::inner_product(ak.begin(), ak.end(), ++pk.begin(), 
+					0.0, std::plus<double>(), std::divides<double>());
+		}
+		operator int() {
+			return std::inner_product(ak.begin(), ak.end(), pk.begin(), 0); 
+		}
+		friend p_adic operator+(const p_adic &m, const p_adic &o);
+		p_adic operator++(int) {
+			p_adic copie = *this;
+			increment();
+			return copie;
+		};
+		p_adic& operator++() {
+			increment();
+			return (*this);
+		};
+		friend struct halton;
+		friend struct kakutani;
+		friend struct faure;
+	protected:
+		void increment();
+		p_adic add(const p_adic &x) const;
+		coeff ak, pk;
+		int p;
+};
+
+static int primes[255] = {
+             2,    3,    5,    7,   11,   13,   17,   19,   23,
+      29,   31,   37,   41,   43,   47,   53,   59,   61,   67,
+      71,   73,   79,   83,   89,   97,  101,  103,  107,  109,
+     113,  127,  131,  137,  139,  149,  151,  157,  163,  167,
+     173,  179,  181,  191,  193,  197,  199,  211,  223,  227,
+     229,  233,  239,  241,  251,  257,  263,  269,  271,  277,
+     281,  283,  293,  307,  311,  313,  317,  331,  337,  347,
+     349,  353,  359,  367,  373,  379,  383,  389,  397,  401, 
+     409,  419,  421,  431,  433,  439,  443,  449,  457,  461,
+     463,  467,  479,  487,  491,  499,  503,  509,  521,  523,
+     541,  547,  557,  563,  569,  571,  577,  587,  593,  599,
+     601,  607,  613,  617,  619,  631,  641,  643,  647,  653,
+     659,  661,  673,  677,  683,  691,  701,  709,  719,  727,
+     733,  739,  743,  751,  757,  761,  769,  773,  787,  797,
+     809,  811,  821,  823,  827,  829,  839,  853,  857,  859,
+     863,  877,  881,  883,  887,  907,  911,  919,  929,  937,
+     941,  947,  953,  967,  971,  977,  983,  991,  997, 1009,
+    1013, 1019, 1021, 1031, 1033, 1039, 1049, 1051, 1061, 1063,
+    1069, 1087, 1091, 1093, 1097, 1103, 1109, 1117, 1123, 1129,
+    1151, 1153, 1163, 1171, 1181, 1187, 1193, 1201, 1213, 1217,
+    1223, 1229, 1231, 1237, 1249, 1259, 1277, 1279, 1283, 1289,
+    1291, 1297, 1301, 1303, 1307, 1319, 1321, 1327, 1361, 1367,
+    1373, 1381, 1399, 1409, 1423, 1427, 1429, 1433, 1439, 1447,
+    1451, 1453, 1459, 1471, 1481, 1483, 1487, 1489, 1493, 1499,
+    1511, 1523, 1531, 1543, 1549, 1553, 1559, 1567, 1571, 1579,
+    1583, 1597, 1601, 1607, 1609, 1613
+};
+
+struct halton {
+	typedef std::vector<double> result_type;
+	typedef std::vector<p_adic> list_p_adic;
+	halton(const list_p_adic &x) : nk(x), result(x.size()) {};
+	halton(int dimension) : result(dimension) {
+		for (int k = 0; k < dimension; k++)
+			nk.push_back(p_adic((int) 1, primes[k]));
+	};
+	halton(int dimension, double x_init[]) : result(dimension) {
+		for (int k = 0; k < dimension; k++)
+			nk.push_back(p_adic(x_init[k], primes[k]));
+	};
+	result_type operator()() {
+		result_type::iterator ir = result.begin();
+		list_p_adic::iterator ink = nk.begin();
+		while (ink != nk.end()) {
+			*ir++ = (double) (*ink++)++;
+		}
+		return result;
+	}
+	protected:
+		list_p_adic nk;
+		result_type result;
+};
+
+struct kakutani {
+	typedef std::vector<double> result_type;
+	typedef std::vector<p_adic> list_p_adic;
+	kakutani(const list_p_adic &x, const list_p_adic &y) : xk(x), yk(y), result(x.size()) {};
+	result_type operator()() {
+		std::transform(xk.begin(), xk.end(), yk.begin(), result.begin(), op);
+		return result;
+	}
+	protected:
+		list_p_adic xk;
+		list_p_adic yk;
+		result_type result;
+		struct t_op {
+			double operator()(p_adic &x, const p_adic &y) { 
+				x = x + y;
+				return (double) x;
+			}; 
+		} op;
+};
+
+struct faure {
+	typedef std::vector<double> result_type;
+	typedef std::vector< std::vector<int> > coeff_binom;
+	faure(int dimension, const p_adic &x) 
+		: x(x), result(dimension), Comb(32, std::vector<int>(32,0))  {
+		for (int n = 0; n < 32; n++) {
+			Comb[0][n] = 1;
+		}
+		for (int k = 1; k < 32; k++) {
+			Comb[k][0] = 1;
+			for (int n = 1; n <= 32-k; n++){
+				Comb[k][n] = Comb[k][n-1] + Comb[k-1][n];
+			}
+		}
+		std::cout << x.p << std::endl;
+	};
+	result_type operator()() {
+		p_adic xi = x++;
+		result_type::iterator it = result.begin();
+		while (it != result.end()) {
+			*it++ = (double) xi;
+			xi = T(xi);
+		}
+		return result;
+	}
+	p_adic T(const p_adic &y) {
+		p_adic::coeff bk;
+		coeff_binom::const_iterator it_Comb = Comb.begin();
+		p_adic::coeff::const_iterator it = y.ak.begin();
+		while (it != y.ak.end()) {
+			bk.push_back(std::inner_product(it, y.ak.end(), (*it_Comb).begin(), 0) % y.p);
+			it++; it_Comb++;
+		}
+		return p_adic(bk, y.pk);
+	}
+	protected:
+		p_adic x;
+		result_type result;
+		coeff_binom Comb; 
+};
+
+
+struct sobol {
+	typedef std::vector<double> result_type;
+	sobol(int dimension) : dimension(dimension), result(dimension) {
+		q = gsl_qrng_alloc(gsl_qrng_sobol, dimension);
+	}
+	sobol(sobol const & o) : dimension(o.dimension) {
+		q = gsl_qrng_clone(o.q);
+		result = o.result;
+	}
+    // constructeur move uniquement en C++11
+    sobol(sobol && o) : dimension(o.dimension), q(o.q), result(std::move(o.result)) {
+        o.dimension = 0;
+        o.q = nullptr;
+    }
+	sobol & operator=(sobol const & o) {
+		if (this != &o) {
+			dimension = o.dimension;
+			gsl_qrng_memcpy(q, o.q);
+			result = o.result;
+		}
+		return *this;
+	}
+    // operateur move uniquement en C++11
+    sobol & operator=(sobol && o) {
+        if (this != &o) {
+            dimension = o.dimension;
+            q = o.q;
+            result = std::move(o.result);
+            o.dimension = 0;
+            o.q = nullptr;
+        }
+        return *this;
+    }
+	~sobol() { gsl_qrng_free(q); }
+	result_type operator()() {
+		gsl_qrng_get(q, &(*result.begin()));
+		return result;
+	}
+	protected: 
+		int dimension;
+		gsl_qrng * q;
+		result_type result;
+};
+
+#endif
diff --git a/option.cpp b/option.cpp
index 7fb5ce3..0e5a56c 100644
--- a/option.cpp
+++ b/option.cpp
@@ -1,6 +1,7 @@
 #include <algorithm>
 #include <iostream>
 #include "rqmc.hpp"
+#include "p_adic.o"
 
 double frac_part(double x){
     return x - floor(x);
@@ -37,36 +38,24 @@ struct asian_option : public std::unary_function<std::vector<double>, double>
         
 
 
-template <typename Fct>
-struct quasi_option : public generator<double>
+template <typename Fct, typename LDS>
+struct quasi_option : public generator<typename Fct::result_type>
 {
-    quasi_option(int n, int d, Fct payoff) : n(n), d(d), payoff(payoff), U(0,1), s(d), seed(d) {};
+    quasi_option(int n, int d, Fct payoff) : n(n), d(d), payoff(payoff), G(d) {};
         
-    double operator()() {
-        std::vector<double> X(d);
+    typename Fct::result_type operator()() {
         double sum =0;
-        for(int i=0; i<d; i++){
-                seed[i]=U();
-            }
-        for (int i=0; i<n; i++){
-            X=s();
-            //std::cout<<X[i]<<std::endl;
-            for(int i=0; i<d; i++){
-                X[i] = gsl_cdf_gaussian_Pinv(frac_part(seed[i]+X[i]), 1);
-            }
-            sum += payoff(X);
-            //std::cout<<X[i]<<std::endl;
+        for(int i=0; i<n; i++){
+            sum += payoff(G());
+            
         }
-        
         return sum/n;
     };
     
     private:
         int n, d;
         Fct payoff;
-        uniform U;
-        sobol s;
-        std::vector<double> seed;
+        quasi_gaussian<LDS> G;
     };
 
 int main(){
@@ -76,10 +65,15 @@ int main(){
     
     int d =16; 
 
-    typedef asian_option payoff_type;
     
     std::vector<double> result(3);
-    result = monte_carlo(100, quasi_option<payoff_type> (N, d, A));
+    result = monte_carlo(100, quasi_option<asian_option, sobol> (N, d, A));
+    for(int i =0; i<3; i++){
+        std::cout<<result[i]<<std::endl;
+    }
+    
+    std::vector<double> result2(3);
+    result = monte_carlo(100, quasi_option<asian_option, halton> (N, d, A));
     for(int i =0; i<3; i++){
         std::cout<<result[i]<<std::endl;
     }
diff --git a/p_adic.cpp b/p_adic.cpp
new file mode 100644
index 0000000..30080ac
--- /dev/null
+++ b/p_adic.cpp
@@ -0,0 +1,23 @@
+#include "low_discrepancy.hpp"
+
+p_adic::p_adic(int n, int p) : p(p) {
+    int puiss = 1;
+    while (n > 0) {
+        ak.push_back(n % p);
+        pk.push_back(puiss);
+        puiss *= p;
+        n -= ak.back();
+        n /= p;
+    }
+    pk.push_back(puiss);
+};
+
+void p_adic::increment() {
+    coeff::iterator i = ak.begin();
+    while ((i != ak.end()) && ((*i)+1 == p)) { (*i) = 0; i++; }
+    if (i == ak.end()) {
+        ak.push_back(1);
+        pk.push_back(pk.back()*p);
+    }
+    else (*i) += 1;
+};
diff --git a/rqmc.hpp b/rqmc.hpp
index 415ab00..1b33713 100644
--- a/rqmc.hpp
+++ b/rqmc.hpp
@@ -79,22 +79,26 @@ std::vector<double> monte_carlo(int n, Fct f, L 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), sob(d) {};
+        : 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 = sob();
+        std::vector<double> result = s();
         for(int i=0;i<d; i++){
-            result[i] = mean + gsl_cdf_gaussian_Pinv(result[i], std);
+            result[i] = mean + gsl_cdf_gaussian_Pinv(frac_part(result[i]+seed[i]), std);
         }
         return value = result;
     }
     private: 
     int d;
     double mean, std;
-    sobol sob;
+    LDS s;
+    uniform U;
+    std::vector<double> seed;
 };