#include <cstdlib>
#include <ctime>
#include <cmath>
#include <iostream>
#include <fstream>
#include <sstream>
#include <string>
#include <assert.h>

using namespace std;
//
// Input for the program RMD part
//
int N = 5;			            // The size of the cube in power of 2
double Hurst_exponent = 0.1;	// The input of the exponent for the sponge
double start_value = 1.;        // The inital values inside the grid

int N_0_p_glob = int(pow(2.,N) + 1);
int N_1_p_glob = int(pow(2.,N) + 1);
int N_2_p_glob = int(pow(2.,N) + 1);					

int N_tot_p_glob = int(N_0_p_glob*N_1_p_glob*N_2_p_glob);
//
// Functions for the program
//
int main(int argc, char *argv[]); // The main program function
double gaussian_random(double m, double s); // Gaussian random number using Box Muller transformation
double delta_n(int n, double Hurst_exponent, double n_p); // Calculating the Gaussian random variable
bool AlmostEqual2sComplement(double A, double B, int maxUlps); // Compare two values 
//
// Function calculating the midpoint values - Line, Face, Volume
//
double midpoint(int n, double p_1, double p_2, double Hurst_exponent);
double midpoint(int n, double p_1, double p_2, double p_3, double p_4, double Hurst_exponent);
double midpoint(int n, double p_1, double p_2, double p_3, double p_4, double p_5, double p_6, double p_7, double p_8, double Hurst_exponent);

int main(int argc, char *argv[]){
//
// Initialize a random seed
//
  srand(time(NULL));
//
// Program variables
//
  // RMD data file - Start
  ofstream rmd_data_file;
  stringstream N_number_rmd, H_number_rmd;
  
  string rmd_data = "rmd_data_";

  N_number_rmd << N_0_p_glob;
  H_number_rmd << Hurst_exponent;
  rmd_data = rmd_data + N_number_rmd.str() + "_";
  rmd_data = rmd_data + H_number_rmd.str() + ".dat";
  // RMD data file - End
  
  rmd_data_file.open(rmd_data.c_str());
  
  cout << "----                 RMD ALGORITHM                    ----" << endl;
  cout << "                                                          " << endl;
  cout << "                 THREE DIMENSIONAL CASE                   " << endl;
  cout << "----                                                  ----" << endl;
  cout << "                                                          " << endl;
  cout << "Hurst Exponent: " << Hurst_exponent << endl;
  cout << "Size of cube: (" << N_0_p_glob << "," << N_1_p_glob << "," << N_2_p_glob << ")" << endl;
  cout << endl;
  
  rmd_data_file << "# Hurst Exponent: " << Hurst_exponent << endl;
  rmd_data_file << "# Size of cube: (" << N_0_p_glob << "," << N_1_p_glob << "," << N_2_p_glob << ")" << endl;
  
//
// RMD algorithm for generating the structure of the grid
//
  double *point = new double[N_tot_p_glob]; // The global structure of the grid

  double corner_points[8] = {1.}; // The inital corner points of the grid

  for (int i_N = 0; i_N <= N; i_N++){
    int n_H_step = i_N; // Value for the Hurst calculation
	
    if (i_N == 0){ // Init: All the points get assigned a start value according to delta function
      
	  for (int i = 0; i < N_tot_p_glob; i++){ *(point+i) = start_value; } // Assign all the points a start value

	  // The corner points - Start
      for (int i_1 = 0; i_1 < 8; i_1++){ corner_points[i_1] = delta_n(n_H_step, Hurst_exponent,3.); } 
      
      *(point+(0*N_2_p_glob*N_1_p_glob + 0*N_1_p_glob + 0)) = corner_points[0];
      *(point+(0*N_2_p_glob*N_1_p_glob + 0*N_1_p_glob + (N_0_p_glob-1))) = corner_points[1];
      *(point+(0*N_2_p_glob*N_1_p_glob + (N_1_p_glob-1)*N_1_p_glob + 0)) = corner_points[2];
      *(point+(0*N_2_p_glob*N_1_p_glob + (N_1_p_glob-1)*N_1_p_glob + (N_0_p_glob-1))) = corner_points[3];
      *(point+((N_2_p_glob-1)*N_2_p_glob*N_1_p_glob + 0*N_1_p_glob + 0)) = corner_points[4];
      *(point+((N_2_p_glob-1)*N_2_p_glob*N_1_p_glob + 0*N_1_p_glob + (N_0_p_glob-1))) = corner_points[5];
      *(point+((N_2_p_glob-1)*N_2_p_glob*N_1_p_glob + (N_1_p_glob-1)*N_1_p_glob + 0)) = corner_points[6];
      *(point+((N_2_p_glob-1)*N_2_p_glob*N_1_p_glob + (N_1_p_glob-1)*N_1_p_glob + (N_0_p_glob-1))) = corner_points[7];
	  // The corner points - End
	  
    }
	else if (i_N > 0){ // Continue: After the init we proceed with the RMD algorithm
	  int p_1_nn, p_2_nn, p_3_nn, p_4_nn, p_5_nn, p_6_nn, p_7_nn, p_8_nn, p_9_nn;
      double p_1, p_2, p_3, p_4, p_5, p_6, p_7, p_8, p_9;
	  
      int p_0_n[3], p_1_n[3], p_2_n[3];
      int p_step_0 = 0, p_step_1 = 0, p_step_2 = 0;
	  
	  // Determine the step size - Start
      p_step_0 = int((N_0_p_glob - 1)/(pow(2.,i_N)));
      p_step_1 = int((N_1_p_glob - 1)/(pow(2.,i_N)));
      p_step_2 = int((N_2_p_glob - 1)/(pow(2.,i_N)));
      // Determine the step size - End

	  // Step through the grid - Start
      for (int cube_0 = 0; (cube_0 < int(pow(2,(i_N - 1.)))); cube_0++){
        for (int cube_1 = 0; (cube_1 < int(pow(2,(i_N - 1.)))); cube_1++){
          for (int cube_2 = 0; (cube_2 < int(pow(2,(i_N - 1.)))); cube_2++){
			
			p_0_n[0] = int(0 + cube_0*2*p_step_0);
			p_0_n[1] = int(p_0_n[0] + p_step_0);
			p_0_n[2] = int(p_0_n[1] + p_step_0);
			
			p_1_n[0] = int(0 + cube_1*2*p_step_1);
			p_1_n[1] = int(p_1_n[0] + p_step_1);
			p_1_n[2] = int(p_1_n[1] + p_step_1);
			
			p_2_n[0] = int(0 + cube_2*2*p_step_2);
			p_2_n[1] = int(p_2_n[0] + p_step_2);
			p_2_n[2] = int(p_2_n[1] + p_step_2);

			for (int ll = 0; ll < 3; ll = ll + 2){
			  
              // Faces 1 and 3
			  p_1_nn = int(p_2_n[ll]*N_2_p_glob*N_1_p_glob + p_1_n[0]*N_1_p_glob + p_0_n[0]); p_1 = *(point+p_1_nn);
			  p_2_nn = int(p_2_n[ll]*N_2_p_glob*N_1_p_glob + p_1_n[2]*N_1_p_glob + p_0_n[0]); p_2 = *(point+p_2_nn);
			  p_3_nn = int(p_2_n[ll]*N_2_p_glob*N_1_p_glob + p_1_n[2]*N_1_p_glob + p_0_n[2]); p_3 = *(point+p_3_nn);
			  p_4_nn = int(p_2_n[ll]*N_2_p_glob*N_1_p_glob + p_1_n[0]*N_1_p_glob + p_0_n[2]); p_4 = *(point+p_4_nn);

			  p_5_nn = int((p_2_n[ll])*N_2_p_glob*N_1_p_glob + (p_1_n[1])*N_1_p_glob + (p_0_n[1])); p_5 = *(point+p_5_nn);

			  if (AlmostEqual2sComplement(p_5, start_value, 1)){
				p_5 = midpoint(n_H_step,p_1,p_2,p_3,p_4,Hurst_exponent);
				*(point+p_5_nn) = p_5;
			  }
			  else{
                p_5 = 0.5*(p_5 + midpoint(n_H_step,p_1,p_2,p_3,p_4,Hurst_exponent));
				*(point+p_5_nn) = p_5;
			  }

			  // Faces 2 and 4
			  p_1_nn = int(p_2_n[0]*N_2_p_glob*N_1_p_glob + p_1_n[0]*N_1_p_glob + p_0_n[ll]); p_1 = *(point+p_1_nn);
			  p_2_nn = int(p_2_n[0]*N_2_p_glob*N_1_p_glob + p_1_n[2]*N_1_p_glob + p_0_n[ll]); p_2 = *(point+p_2_nn);
			  p_3_nn = int(p_2_n[2]*N_2_p_glob*N_1_p_glob + p_1_n[2]*N_1_p_glob + p_0_n[ll]); p_3 = *(point+p_3_nn);
			  p_4_nn = int(p_2_n[2]*N_2_p_glob*N_1_p_glob + p_1_n[0]*N_1_p_glob + p_0_n[ll]); p_4 = *(point+p_4_nn);

			  p_5_nn = int((p_2_n[1])*N_2_p_glob*N_1_p_glob + (p_1_n[1])*N_1_p_glob + (p_0_n[ll])); p_5 = *(point+p_5_nn);

			  if (AlmostEqual2sComplement(p_5, start_value, 1)){
				p_5 = midpoint(n_H_step,p_1,p_2,p_3,p_4,Hurst_exponent);
				*(point+p_5_nn) = p_5;
			  }
			  else{
                p_5 = 0.5*(p_5 + midpoint(n_H_step,p_1,p_2,p_3,p_4,Hurst_exponent));
				*(point+p_5_nn) = p_5;
			  }

			}
			
			for (int ll_1 = 0; ll_1 < 3; ll_1 = ll_1 + 2){

			  // Faces 1 and 3 - Bottom
			  p_1_nn = int(p_2_n[ll_1]*N_2_p_glob*N_1_p_glob + p_1_n[0]*N_1_p_glob + p_0_n[0]); p_1 = *(point+p_1_nn);
			  p_2_nn = int(p_2_n[ll_1]*N_2_p_glob*N_1_p_glob + p_1_n[0]*N_1_p_glob + p_0_n[2]); p_2 = *(point+p_2_nn);

			  p_3_nn = int((p_2_n[ll_1])*N_2_p_glob*N_1_p_glob + (p_1_n[0])*N_1_p_glob + (p_0_n[1])); p_3 = *(point+p_3_nn);

			  if (AlmostEqual2sComplement(p_3, start_value, 1)){
                p_3 = midpoint(n_H_step,p_1,p_2,Hurst_exponent);
				*(point+p_3_nn) = p_3;
			  }
			  else{
                p_3 = 0.5*(p_3 + midpoint(n_H_step,p_1,p_2,Hurst_exponent));
				*(point+p_3_nn) = p_3;
			  }
			  
			  // Faces 1 and 3 - Right
			  p_1_nn = int(p_2_n[ll_1]*N_2_p_glob*N_1_p_glob + p_1_n[0]*N_1_p_glob + p_0_n[2]); p_1 = *(point+p_1_nn);
			  p_2_nn = int(p_2_n[ll_1]*N_2_p_glob*N_1_p_glob + p_1_n[2]*N_1_p_glob + p_0_n[2]); p_2 = *(point+p_2_nn);

			  p_3_nn = int((p_2_n[ll_1])*N_2_p_glob*N_1_p_glob + (p_1_n[1])*N_1_p_glob + (p_0_n[2])); p_3 = *(point+p_3_nn);

			  if (AlmostEqual2sComplement(p_3, start_value, 1)){
                p_3 = midpoint(n_H_step,p_1,p_2,Hurst_exponent);
				*(point+p_3_nn) = p_3;
			  }
			  else{
                p_3 = 0.5*(p_3 + midpoint(n_H_step,p_1,p_2,Hurst_exponent));
				*(point+p_3_nn) = p_3;
			  }

			  // Faces 1 and 3 - Top
			  p_1_nn = int(p_2_n[ll_1]*N_2_p_glob*N_1_p_glob + p_1_n[2]*N_1_p_glob + p_0_n[2]); p_1 = *(point+p_1_nn);
			  p_2_nn = int(p_2_n[ll_1]*N_2_p_glob*N_1_p_glob + p_1_n[2]*N_1_p_glob + p_0_n[0]); p_2 = *(point+p_2_nn);

			  p_3_nn = int((p_2_n[ll_1])*N_2_p_glob*N_1_p_glob + (p_1_n[2])*N_1_p_glob + (p_0_n[1])); p_3 = *(point+p_3_nn);
			  
			  if (AlmostEqual2sComplement(p_3, start_value, 1)){
                p_3 = midpoint(n_H_step,p_1,p_2,Hurst_exponent);
				*(point+p_3_nn) = p_3;
			  }
			  else{
                p_3 = 0.5*(p_3 + midpoint(n_H_step,p_1,p_2,Hurst_exponent));
				*(point+p_3_nn) = p_3;
			  }

			  // Faces 1 and 3 - Left
			  p_1_nn = int(p_2_n[ll_1]*N_2_p_glob*N_1_p_glob + p_1_n[2]*N_1_p_glob + p_0_n[0]); p_1 = *(point+p_1_nn);
			  p_2_nn = int(p_2_n[ll_1]*N_2_p_glob*N_1_p_glob + p_1_n[0]*N_1_p_glob + p_0_n[0]); p_2 = *(point+p_2_nn);
			  
			  p_3_nn = int((p_2_n[ll_1])*N_2_p_glob*N_1_p_glob + (p_1_n[1])*N_1_p_glob + (p_0_n[0])); p_3 = *(point+p_3_nn);

			  if (AlmostEqual2sComplement(p_3, start_value, 1)){
                p_3 = midpoint(n_H_step,p_1,p_2,Hurst_exponent);
				*(point+p_3_nn) = p_3;
			  }
			  else{
                p_3 = 0.5*(p_3 + midpoint(n_H_step,p_1,p_2,Hurst_exponent));
				*(point+p_3_nn) = p_3;
			  }
			  
			  // Faces 2 and 4 - Bottom
			  p_1_nn = int(p_2_n[0]*N_2_p_glob*N_1_p_glob + p_1_n[0]*N_1_p_glob + p_0_n[ll_1]); p_1 = *(point+p_1_nn);
			  p_2_nn = int(p_2_n[2]*N_2_p_glob*N_1_p_glob + p_1_n[0]*N_1_p_glob + p_0_n[ll_1]); p_2 = *(point+p_2_nn);
			  
			  p_3_nn = int((p_2_n[1])*N_2_p_glob*N_1_p_glob + (p_1_n[0])*N_1_p_glob + (p_0_n[ll_1])); p_3 = *(point+p_3_nn);
			  
			  if (AlmostEqual2sComplement(p_3, start_value, 1)){
                p_3 = midpoint(n_H_step,p_1,p_2,Hurst_exponent);
				*(point+p_3_nn) = p_3;
			  }
			  else{
                p_3 = 0.5*(p_3 + midpoint(n_H_step,p_1,p_2,Hurst_exponent));
				*(point+p_3_nn) = p_3;
			  }
			  
			  // Faces 2 and 4 - Top
			  p_1_nn = int(p_2_n[2]*N_2_p_glob*N_1_p_glob + p_1_n[2]*N_1_p_glob + p_0_n[ll_1]); p_1 = *(point+p_1_nn);
			  p_2_nn = int(p_2_n[0]*N_2_p_glob*N_1_p_glob + p_1_n[2]*N_1_p_glob + p_0_n[ll_1]); p_2 = *(point+p_2_nn);
			  
			  p_3_nn = int((p_2_n[1])*N_2_p_glob*N_1_p_glob + (p_1_n[2])*N_1_p_glob + (p_0_n[ll_1])); p_3 = *(point+p_3_nn);
			  
			  if (AlmostEqual2sComplement(p_3, start_value, 1)){
                p_3 = midpoint(n_H_step,p_1,p_2,Hurst_exponent);
				*(point+p_3_nn) = p_3;
			  }
			  else{
                p_3 = 0.5*(p_3 + midpoint(n_H_step,p_1,p_2,Hurst_exponent));
				*(point+p_3_nn) = p_3;
			  }
			  
			}
			
			for (int ll_2 = 0; ll_2 < 3; ll_2 = ll_2 + 2){
			  
			  // Faces 5 and 6			
			  p_1_nn = int(p_2_n[0]*N_2_p_glob*N_1_p_glob + p_1_n[ll_2]*N_1_p_glob + p_0_n[0]); p_1 = *(point+p_1_nn);
			  p_2_nn = int(p_2_n[0]*N_2_p_glob*N_1_p_glob + p_1_n[ll_2]*N_1_p_glob + p_0_n[2]); p_2 = *(point+p_2_nn);
			  p_3_nn = int(p_2_n[2]*N_2_p_glob*N_1_p_glob + p_1_n[ll_2]*N_1_p_glob + p_0_n[2]); p_3 = *(point+p_3_nn);
			  p_4_nn = int(p_2_n[2]*N_2_p_glob*N_1_p_glob + p_1_n[ll_2]*N_1_p_glob + p_0_n[0]); p_4 = *(point+p_4_nn);
			  
			  p_5_nn = int((p_2_n[1])* N_2_p_glob*N_1_p_glob+ (p_1_n[ll_2])*N_1_p_glob + (p_0_n[1])); p_5 = *(point+p_5_nn);
			  
			  if (AlmostEqual2sComplement(p_5, start_value, 1)){
				p_5 = midpoint(n_H_step,p_1,p_2,p_3,p_4,Hurst_exponent);
				*(point+p_5_nn) = p_5;
			  }
			  else{
                p_5 = 0.5*(p_5 + midpoint(n_H_step,p_1,p_2,p_3,p_4,Hurst_exponent));
				*(point+p_5_nn) = p_5;
			  }
			  
			}
			
			// Volume
			p_1_nn = int(p_2_n[0]*N_2_p_glob*N_1_p_glob + p_1_n[0]*N_1_p_glob + p_0_n[0]); p_1 = *(point+p_1_nn);
			p_2_nn = int(p_2_n[0]*N_2_p_glob*N_1_p_glob + p_1_n[2]*N_1_p_glob + p_0_n[0]); p_2 = *(point+p_2_nn);
			p_3_nn = int(p_2_n[0]*N_2_p_glob*N_1_p_glob + p_1_n[2]*N_1_p_glob + p_0_n[2]); p_3 = *(point+p_3_nn);
			p_4_nn = int(p_2_n[0]*N_2_p_glob*N_1_p_glob + p_1_n[0]*N_1_p_glob + p_0_n[2]); p_4 = *(point+p_4_nn);
			p_5_nn = int(p_2_n[2]*N_2_p_glob*N_1_p_glob + p_1_n[0]*N_1_p_glob + p_0_n[0]); p_5 = *(point+p_5_nn);
			p_6_nn = int(p_2_n[2]*N_2_p_glob*N_1_p_glob + p_1_n[2]*N_1_p_glob + p_0_n[0]); p_6 = *(point+p_6_nn);
			p_7_nn = int(p_2_n[2]*N_2_p_glob*N_1_p_glob + p_1_n[0]*N_1_p_glob + p_0_n[2]); p_7 = *(point+p_7_nn);
			p_8_nn = int(p_2_n[2]*N_2_p_glob*N_1_p_glob + p_1_n[2]*N_1_p_glob + p_0_n[2]); p_8 = *(point+p_8_nn);
			
			p_9_nn = int((p_2_n[1])*N_2_p_glob*N_1_p_glob + (p_1_n[1])*N_1_p_glob + (p_0_n[1])); p_9 = *(point+p_9_nn);
			
			if (AlmostEqual2sComplement(p_9, start_value, 1)){
			  p_9 = midpoint(n_H_step,p_1,p_2,p_3,p_4,p_5,p_6,p_7,p_8,Hurst_exponent);
			  *(point+p_9_nn) = p_9;
			}
			else{
              p_9 = 0.5*(p_9 + midpoint(n_H_step,p_1,p_2,p_3,p_4,p_5,p_6,p_7,p_8,Hurst_exponent));
			  *(point+p_9_nn) = p_9;
			}
			
		  }
		}
      }
	  // Step through the grid - End
    }
  }
  
  for (int i = 0; i < N_tot_p_glob; i++){
	rmd_data_file << *(point+i) << endl; 
  }

  rmd_data_file.close();
  
  return 0;
}

//
// Function : Calculating the Gaussian random variable
//
double delta_n(int n, double Hurst_exponent, double n_p){
  double delta_1 = 0.;
  double sigma_1 = 0.;
  
  sigma_1 = gaussian_random(0., 1.);
    
  delta_1 = sigma_1*(1./pow(n_p,(0.5*Hurst_exponent)))*(1./pow(2.,(1.*Hurst_exponent*n)))*sqrt(1. - pow(2.,(2.*Hurst_exponent - 2.*n_p)));
  
  return (delta_1);
}
//
// Function : Generating a Gaussian random number using Box Muller transformation
//
double gaussian_random(double m, double s){
  double x_1, x_2, w, y_1;
  static double y_2;
  static int use_last = 0;
  
  if (use_last){
    y_1 = y_2;
	use_last = 0;
  }
  else{
    do{
      x_1 = 2.*((double)rand()/(double)RAND_MAX) - 1.;
      x_2 = 2.*((double)rand()/(double)RAND_MAX) - 1.;
      w = x_1*x_1 + x_2*x_2;
    } while (w >= 1.);
	
    w = sqrt((-2.*log(w))/w);
    y_1 = x_1*w;
    y_2 = x_2*w;
    use_last = 1;
  }
  
  return (m + y_1*s);
}
//
// --- Function
//
double midpoint(int n, double p_1, double p_2, double Hurst_exponent){
  double p_new = 0;
  
  p_new = 0.5*(p_1 + p_2) + delta_n(n,Hurst_exponent,1.);
  
  return p_new;
}
//
// --- Function
//
double midpoint(int n, double p_1, double p_2, double p_3, double p_4, double Hurst_exponent){
  double p_new = 0;
  
  p_new = 0.25*(p_1 + p_2 + p_3 + p_4) + delta_n(n,Hurst_exponent,2.);
  
  return p_new;
}
//
// --- Function
//
double midpoint(int n, double p_1, double p_2, double p_3, double p_4, double p_5, double p_6, double p_7, double p_8, double Hurst_exponent){
  double p_new = 0;
  
  p_new = 0.125*(p_1 + p_2 + p_3 + p_4 + p_5 + p_6 + p_7 + p_8) + delta_n(n,Hurst_exponent,3.);
  
  return p_new;
}
//
// AlmostEqual function
//
bool AlmostEqual2sComplement(double A, double B, int maxUlps){
  // Make sure maxUlps is non-negative and small enough that the
  // default NAN won't compare as equal to anything.
  assert(maxUlps > 0 && maxUlps < 4 * 1024 * 1024);
  int aInt = *(int*)&A;
  // Make aInt lexicographically ordered as a twos-complement int
  if (aInt < 0)
    aInt = 0x80000000 - aInt;
  // Make bInt lexicographically ordered as a twos-complement int
  int bInt = *(int*)&B;
  if (bInt < 0)
    bInt = 0x80000000 - bInt;
  int intDiff = abs(aInt - bInt);
  if (intDiff <= maxUlps)
    return true;
  return false;
}