// C-code C-programs to compute (approximations) of
// two-sided Brownian motion (+ drift) and its successive k-1 integrals.
// Converted to C by Karim Filali from S-programs written by F. Balabdaoui 

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <string.h>

#define M_SQRT2 (1.414213562373095)

#define VERBOSE 0
#define NO_RAND 0
// if you set the below to 1 problems arise because apparently factorials become too big.
#define USE_OLD_FACT 0

double SchauderFunc(double);
double IntSchauderFunc(int l, int p, int i, double x);
void IntBrownFunc(double* IntBrown, int K, int m);
void IntBrown0C(double* Output, int K, double C,int m);
void IntBrownCCdrift(double* Output, int K, double C,int m);
void matMulAdd(double *Bi,double *Wi,double* BiMinus1,double* Matgrid, int k, int Lg,int stride);
long factorial(long n);
double factRatio(int k,int K);
double inverse_normal_func(double p);
double myPow(long x, long a);

int main(int argc, char*argv[]){

  char filename[50]="-";
  FILE*ofp;
  int y,k;
  int K,m;
  double C;
  int Lg;
  double* Output;

  srand(time(NULL));

  if(argc<4) {
    fprintf(stderr,"Need 3 arguments: K, C and m\n");
    exit(-1);
  }
  else{
    K=atoi(argv[1]);
    C=atof(argv[2]);
    m=atoi(argv[3]);
    if(argc >4)
      strcpy(filename,argv[4]);
  }

  if(strcmp(filename,"-")==0) {
    ofp=stdout;
  }
  else {
    if((ofp=fopen(filename,"w"))==NULL) {
      fprintf(stderr,"Could not open %s for writing\n",filename);
      exit(-1);
    }
  }

#if NO_RAND
  m=4;
#endif

#ifdef VERBOSE
  fprintf(stderr,"K=%d, C= %f, m=%d\n",K,C,m);
#endif

   Lg = (int)(pow(2.0,(double)(m))*2*C+1);
   Output = (double*)calloc((Lg+1)*(K+1),sizeof(double)); 
   if(Output == NULL) {
    fprintf(stderr,"Could not allocate memory for buffer\n");
    exit(-1);
  } 


  IntBrownCCdrift(Output, K, C,m);

  for(y=1;y<=Lg;++y) {
    for(k=1;k<=K;++k) {
      fprintf(ofp,"%e ",*(Output + k*(Lg+1) + y));
    }
    fprintf(ofp,"\n");
  }


  if(strcmp(filename,"-")!=0) {
    fclose(ofp);
  }

  free(Output);

  return 0;
}

void IntBrownCCdrift(double* Output, int K, double C,int m) {
  int k,y,i;
  int Lg = (int)(pow(2.0,(double)(m))*C+1);
  int stride = (int)(pow(2.0,(double)(m))*2*C+1+1);
  double* IntBrown = (double*)calloc((Lg+1)*(K+1),sizeof(double)); 
  double* grid = (double*)calloc(stride,sizeof(double)); 
  double val;
  double twoMinp = pow(2.0,(double)(-m));

  double fact_ratio;
  double yKk;

  i=1; ////////
  val=-C;
  while(val<=C) {
    grid[i]=val;
    val += twoMinp;  // equivalent to val = val+twoMinp;
    i++;  // i=i+1;
  }

  IntBrown0C(IntBrown,K,C,m);

  for(k=1;k<=K;++k)
    for(y=1;y<=Lg;++y)
      *(Output + k*(stride) + y) = *(IntBrown + k*(Lg+1) + (Lg-y+1));

  IntBrown0C(IntBrown,K,C,m);

  for(k=1;k<=K;++k)
    for(y=Lg;y<=stride-1;++y)
      *(Output + k*(stride) + y) = *(IntBrown + k*(Lg+1) + (y-Lg+1));

  for(k=1;k<=K;++k) {
#if USE_OLD_FACT
    fact_ratio=(double)factorial(K)/factorial(k+K);
#else
    fact_ratio=factRatio(k,K);
#endif
    for(y=1;y<=stride-1;++y) {
      yKk=pow(grid[y],(double)(K+k));
      *(Output + k*(stride) + y) += fact_ratio * yKk;
      }
  }

  free(IntBrown);
  free(grid);
}

void IntBrown0C(double* Output, int K, double C,int m) {
  int k,y,i,j;
  double val;
  double twoMinp = pow(2.0,(double)(-m));
  int Lg = (int)pow(2.0,(double)(m))+1;
  double* grid = (double*)calloc((Lg+1),sizeof(double)); 
  double* Wi = (double*)calloc((Lg+1)*(K+1),sizeof(double)); 
  double* Matgrid = (double*)calloc((Lg+1)*(K+1),sizeof(double)); 
  double* BiMinus1 = (double*)calloc((K+1),sizeof(double)); 
  int stride = (int)(pow(2.0,(double)(m))*C+1+1);

  double* Bi=Output; 

  i=1; ////////
  val=0.0;
  while(val<=1) {
    grid[i]=val;
    val += twoMinp;  // equivalent to val = val+twoMinp;
    i++;  // i=i+1;
  }

  for(i=1;i<=C;++i) {
    IntBrownFunc(Wi,K,m);

    for(k=1;k<=K;++k) {
      for(j=1;j<=k;++j) {
	for(y=1;y<=Lg;++y) {
	  Matgrid[j*(Lg+1)+y]=pow(grid[y],(double)(k-j))/factorial(k-j);
	}
      }
      matMulAdd(Bi,Wi,BiMinus1,Matgrid,k,Lg,stride);
    }
    for(k=1;k<=K;++k) {
      BiMinus1[k]=Bi[k*(stride)+Lg];
    } 

    Bi+=Lg-1;
  }

  free(grid);
  free(Wi);
  free(Matgrid);
  free(BiMinus1);

}

void IntBrownFunc(double* IntBrown, int K, int m) {

  int i,y,k,n,j;
  double val;
  int twon;
  double twoMinp = pow(2.0,(double)(-m));
  int Lg = (int)pow(2.0,(double)(m))+1;
  double* grid = (double*)calloc(Lg+1,sizeof(double)); 
  double* IntUm = (double*)calloc((Lg+1)*(K+1),sizeof(double)); 
  double Z;

#if NO_RAND
double Zv[32]={9999999,0.315220356754406,-0.693934891037645,1.00410727769503,0.637379833073138,0.447755881608347,-0.959890939543513,0.0282058839228757,-0.850544262911404,-1.46369503097069,1.50253713491473,0.137841998852249,-0.401605497652938,0.965872946920573,-2.10800275305106,0.626552100563295,-1.02595529519296,-1.5290861653951,0.283879864780989,0.434617530200291,-0.755713880659146,-0.974956035422474,0.547855181900798,1.24478086115561,-1.00849395831567,-1.17673133386132,1.06082946495696,-0.535660164002144,-1.23955707460312,0.829774474538486,1.3322405161178,0.904988973718513};
Z=1.1463916094853;
#else
 int twomPlus1Min1 = (int) pow(2.0,(double)(m+1))-1; 
 double* Zv = (double*)calloc(twomPlus1Min1+1,sizeof(double)); 
 for (i = 0; i < twomPlus1Min1+1; i++ ) ///
   Zv[i] = inverse_normal_func((double)rand()/RAND_MAX);
 Z=inverse_normal_func((double)rand()/RAND_MAX); 
#endif

  i=1; ////////
  val=0;
  while(val<=1) {
    grid[i]=val;
    val += twoMinp;
    i++;
  }

  for(k=0;k<=K;++k)
    for(y=0;y<=Lg;++y)
    IntUm[k*(Lg+1)+y]=0;

  for(y=2;y<=Lg;++y)
    for(k=1;k<=K;++k)
      for(n=0;n<=m;++n) {
	twon=(int)pow(2.0,(double)n);
	for(j=0;j<=twon-1;++j) {
	  *(IntUm + k*(Lg+1) + y) += (Zv[twon + j] * IntSchauderFunc(k,n,j,grid[y]));
	}
      }

  for(k=1;k<=K;++k)
    for(y=1;y<=Lg;++y)
      IntBrown[k*(Lg+1) + y]=  IntUm[k*(Lg+1) + y] + Z*pow(grid[y],(double)k)/factorial(k);

  free(grid);
  free(Zv);
  free(IntUm);
}

double IntSchauderFunc(int l, int p, int i, double x) {
  double IntSchauder=0.0;
  double twop = pow((double)2,(double)p);
  double twopMin1 = twop -1;
  double twoMinp = pow((double)2,-(double)p);
  double twoHalfp = pow((double)2,(double)p/2.0);
  double twoMinHalfp = pow((double)2,-(double)p/2.0);
  double twoMinpPlus1l = pow((double)2,-(double)(p+1)*l);
  double twoMinpPlus1lMin1 = pow((double)2,-(double)(p+1)*(l-1));
  double twoMinHalfpPlus1 = twoMinHalfp/2.0;
  double factlMin1 = factorial(l-1);
  double factl = factlMin1*l;

  if(i < 0 || i > twopMin1) {
    fprintf(stderr,"i (%d) has to be between 0 and 2^p-1 (%d)\n",i, (int)twopMin1);	
    exit(-1); 
  }
  if(l==1) {
    IntSchauder = twoMinHalfp*SchauderFunc((double) twop*x-i);
  }
  else {
    if(x >= twoMinp * i) {  
      if( x <= twoMinp*(i + 1/2.0))	
	IntSchauder = twoHalfp /factl * pow((double)(x- twoMinp*i),(double)l); 
      else {
	if( x <=  twoMinp*(i + 1))
	  IntSchauder = twoHalfp/factl*(twoMinpPlus1l - pow((double) x - twoMinp*(i+1/2.0),(double) l)) + (twoMinHalfpPlus1/factlMin1)*pow((double) x-twoMinp*(i+1/2.0),(double) l-1);
	else
	  IntSchauder = twoMinHalfpPlus1*twoMinpPlus1lMin1/factlMin1;
      } 
    }
  }

  return IntSchauder;
}

double SchauderFunc(double x) {
  double Schauder= 0.0;
  if(x >= 0 && x <= 0.5) 
    Schauder = x;
  else if(x > 0.5 && x <= 1) 
    Schauder = 1-x;
  return Schauder ;
}

void matMulAdd(double *Bi,double *Wi,double* BiMinus1,double* Matgrid, int k, int Lg,int stride) {
  int y,j;
  double tmp=0.0;
  for(y=1;y<=Lg;++y) {
    tmp=0.0;
    for(j=1;j<=k;++j) {
      tmp+=(BiMinus1[j]*Matgrid[j*(Lg+1)+y]);
    }
    Bi[k*(stride)+y]=Wi[k*(Lg+1)+y]+tmp;
  }

}

double myPow(long x, long a) {
  long i;
  int invert=0;
  double result=(double)x;
  if(a==0) return 1;
  else if(a < 0) {
    invert=1;
    a=-a;
  }
  for(i=1;i<a;++i)
    result*=x;
  if(invert)
    return 1.0/result;
  else
    return (double)result;

}

long factorial(long n) {
  long fact=1;
  long i;

  for(i=2; i <= n; ++i) 
    fact *= i;

  return fact;

}


double factRatio(int k,int K) {
  double fact=1;
  long i;

  for(i=K+1; i <= k+K; ++i)
    fact *= (1.0/(double)i);

  return fact;
}



double inverse_error_func(double p) {
	/* 
           Source: This routine was derived (using f2c) from the 
                   FORTRAN subroutine MERFI found in 
                   ACM Algorithm 602 obtained from netlib.

                   MDNRIS code contains the 1978 Copyright 
                   by IMSL, INC. .  Since MERFI has been 
                   submitted to netlib, it may be used with 
                   the restriction that it may only be 
                   used for noncommercial purposes and that
                   IMSL be acknowledged as the copyright-holder
                   of the code.
        */



	/* Initialized data */
	static double a1 = -.5751703;
	static double a2 = -1.896513;
	static double a3 = -.05496261;
	static double b0 = -.113773;
	static double b1 = -3.293474;
	static double b2 = -2.374996;
	static double b3 = -1.187515;
	static double c0 = -.1146666;
	static double c1 = -.1314774;
	static double c2 = -.2368201;
	static double c3 = .05073975;
	static double d0 = -44.27977;
	static double d1 = 21.98546;
	static double d2 = -7.586103;
	static double e0 = -.05668422;
	static double e1 = .3937021;
	static double e2 = -.3166501;
	static double e3 = .06208963;
	static double f0 = -6.266786;
	static double f1 = 4.666263;
	static double f2 = -2.962883;
	static double g0 = 1.851159e-4;
	static double g1 = -.002028152;
	static double g2 = -.1498384;
	static double g3 = .01078639;
	static double h0 = .09952975;
	static double h1 = .5211733;
	static double h2 = -.06888301;

	/* Local variables */
	static double a, b, f, w, x, y, z, sigma, z2, sd, wi, sn;

	x = p;

	/* determine sign of x */
	if (x > 0)
		sigma = 1.0;
	else
		sigma = -1.0;

	/* Note: -1.0 < x < 1.0 */

	z = fabs(x);

	/* z between 0.0 and 0.85, approx. f by a 
	   rational function in z  */

	if (z <= 0.85) {
		z2 = z * z;
		f = z + z * (b0 + a1 * z2 / (b1 + z2 + a2 
		    / (b2 + z2 + a3 / (b3 + z2))));

	/* z greater than 0.85 */
	} else {
		a = 1.0 - z;
		b = z;

		/* reduced argument is in (0.85,1.0), 
		   obtain the transformed variable */

		w = sqrt(-(double)log(a + a * b));

		/* w greater than 4.0, approx. f by a 
		   rational function in 1.0 / w */

		if (w >= 4.0) {
			wi = 1.0 / w;
			sn = ((g3 * wi + g2) * wi + g1) * wi;
			sd = ((wi + h2) * wi + h1) * wi + h0;
			f = w + w * (g0 + sn / sd);

		/* w between 2.5 and 4.0, approx. 
		   f by a rational function in w */

		} else if (w < 4.0 && w > 2.5) {
			sn = ((e3 * w + e2) * w + e1) * w;
			sd = ((w + f2) * w + f1) * w + f0;
			f = w + w * (e0 + sn / sd);

		/* w between 1.13222 and 2.5, approx. f by 
		   a rational function in w */
		} else if (w <= 2.5 && w > 1.13222) {
			sn = ((c3 * w + c2) * w + c1) * w;
			sd = ((w + d2) * w + d1) * w + d0;
			f = w + w * (c0 + sn / sd);
		}
	}
	y = sigma * f;
	return(y);
}


double inverse_normal_func(double p) {
	/* 
           Source: This routine was derived (using f2c) from the 
                   FORTRAN subroutine MDNRIS found in 
                   ACM Algorithm 602 obtained from netlib.

                   MDNRIS code contains the 1978 Copyright 
                   by IMSL, INC. .  Since MDNRIS has been 
                   submitted to netlib it may be used with 
                   the restriction that it may only be 
                   used for noncommercial purposes and that
                   IMSL be acknowledged as the copyright-holder
                   of the code.
        */

	/* Initialized data */
	static double eps = 1e-10;
	static double g0 = 1.851159e-4;
	static double g1 = -.002028152;
	static double g2 = -.1498384;
	static double g3 = .01078639;
	static double h0 = .09952975;
	static double h1 = .5211733;
	static double h2 = -.06888301;
	static double sqrt2 = M_SQRT2; /* 1.414213562373095; */

	/* Local variables */
	static double a, w, x;
	static double sd, wi, sn, y;

	double inverse_error_func(double p);

	/* Note: 0.0 < p < 1.0 */
	/* assert ( 0.0 < p && p < 1.0 ); */

	/* p too small, compute y directly */
	if (p <= eps) {
		a = p + p;
		w = sqrt(-(double)log(a + (a - a * a)));

		/* use a rational function in 1.0 / w */
		wi = 1.0 / w;
		sn = ((g3 * wi + g2) * wi + g1) * wi;
		sd = ((wi + h2) * wi + h1) * wi + h0;
		y = w + w * (g0 + sn / sd);
		y = -y * sqrt2;
	} else {
		x = 1.0 - (p + p);
		y = inverse_error_func(x);
		y = -sqrt2 * y;
	}
	return(y);
}