#include <algorithm>
#include <iostream>
#include "rqmc.hpp"

double frac_part(double x){
    return x - floor(x);
}

double pos (double x){
    return x>0?x:0;
}

struct asian_option : public std::unary_function<std::vector<double>, double>
{
    asian_option(double r, double T, double S0, double V, int d, double K)
        : r(r), T(T), S0(S0), V(V), d(d), K(K) {};
    
    double operator()(std::vector<double> X) const {
        std::vector<double> S(d);
        S[0]= S0*exp((r-V*V/2)*(T/d)+V*sqrt(T/d)*X[0]);
        for(int i=1;i<d;i++){
            S[i]=S[i-1]*exp((r-V*V/2)*(T/d)+V*sqrt(T/d)*X[i]);
        }
        double temp = std::accumulate(S.begin(), S.end(), 0.)/d;
        return exp(-r*T)*pos(temp-K);
        };
        
    private:
        double r;
        double T;
        double S0;
        double V;
        int d;
        double K;
    };
        
        


template <typename Fct>
struct quasi_option : public generator<double>
{
    quasi_option(int n, int d, Fct payoff) : n(n), d(d), payoff(payoff), U(0,1), s(d), seed(d) {};
        
    double operator()() {
        std::vector<double> X(d);
        double sum =0;
        for (int i=0; i<n; i++){
            X=s();
            //std::cout<<X[i]<<std::endl;
            for(int i=0; i<d; i++){
                seed[i]=U();
                X[i] = gsl_cdf_gaussian_Pinv(frac_part(seed[i]+X[i]), 1);
            }
            sum += payoff(X);
            //std::cout<<X[i]<<std::endl;
        }
        
        return sum/n;
    };
    
    private:
        int n, d;
        Fct payoff;
        uniform U;
        sobol s;
        std::vector<double> seed;
    };

int main(){
    init_alea(1);
    asian_option A(0.05, 1.0, 50.0, 0.1, 16, 45);
    int N= 10000;
    
    int d =16; 

    typedef asian_option payoff_type;
    for(int i =0; i<3; i++){
        std::cout<<monte_carlo(100, quasi_option<payoff_type> (N, d, A))[i]<<std::endl;
    }


    
    return 0;
}