#include <stdlib.h>
#include <stdio.h>
#include <time.h>
#include <math.h>
#include "omp.h"

#define P 4


const int M=365; // nombre de pas de temps
const double two_pi =6.28318530718;


double gauss(){
double x,y;
x= (1.+rand())/(1.+RAND_MAX);
y= (1.+rand())/(1.+RAND_MAX);
return sqrt( -2 *log(x) )*cos(two_pi*y);
}

double EDOstoch(double S0, double dt, double sdt, double rdt){
double S= S0;
int i; for(i=1;i<M;i++)
S= S*(1.+gauss()*sdt+rdt);
return S;
}



int main(){
    double S[P]; srand(time(NULL));
    const double K=110, S0=100, T=1., r=0.03, sigma=0.2;
    int k,q, kmax=20000; // nb de realisations
    double Sa;
    for(q=0; q<P;q++) S[q]= S0; 
    double time0=omp_get_wtime();
    double dt=T/M, sdt=sigma*sqrt(dt), rdt=r*dt, P0=0, t0=clock();
    
    
#pragma omp parallel for num_threads(P) 
        for(q=0; q<P;q++) { 
                srand(q+time(NULL));
                for(k=0; k<kmax/P;k++){ 
                        Sa= EDOstoch(S0, dt, sdt, rdt);
                        if(K>Sa) S[q] += K-Sa;
                }
        }
        time0-=omp_get_wtime();

double PT = S[0];
for(q=1; q<P;q++) PT += S[q];
printf("P_T = %f CPUtime=%f\n",PT*exp(-r*T)/kmax, -time0);
while(!getchar());
return 0;}