#include <stdio.h>
#include <mpl.h>

/*#define INLINE_MATH */
#include <math.h>

#include <time.h>  /* timer */

extern float f_pow();
extern float f_sqrt();
extern plural float fp_sqrt();
extern plural float fp_exp();
extern float reduceAddf();
extern void dpuTimerStart();
extern double dpuTimerElapsed();
extern double dpuTimerConst();
extern unsigned long dpuTimerTicks();

main()
{
  float sz,x,r,t,v,psi,xmu,cor,div,call,tdiv;
  int reps, reps1, i;
  float american(),european();
  double time;

  sz = 110;
  reps = 17;
  psi = 1.0;
  xmu = 0.0 ;
  cor = 0.234;
  x = 100.00 ;
  t = .25;                                                  
  r = 0.0;                                                    
  v = 0.15;
  reps = 17; tdiv = 0.125; div = 1.0;

/*  printf("****** American ******\n");

  for (reps1=1; reps1<reps; reps1++) {
    time = 0.0;
    dpuTimerStart();
    for (i=1; i<=10; i++)
      call = american(sz,x,r,t,v,psi,xmu,cor,reps1,reps,div);

    time = dpuTimerElapsed();
    printf("Elapsed time = %f millisec.\n",time*100.0); 

    printf("reps1 = %i call price = %f \n",reps1,call);
  };
*/
  printf("****** European ******\n");

  for (reps=13; reps<18; reps++) {
    time = 0.0;
    dpuTimerStart();
    for (i=1; i<=10; i++)
      call = european(sz,x,r,t,v,psi,xmu,cor,reps);

    time = dpuTimerElapsed();
    printf("Elapsed time = %f millisec.\n",time*100.0); 

    printf("reps = %i call price = %f \n",reps,call);
  };

  return 0;
}

float american(s,x,r,t,v,psi,xmu,cor,reps1,reps,div)
  float s,x,r,t,v,psi,xmu,cor,div;
  int reps, reps1;
{
  plural float ss1,vv1,aatmp1,aatmp2,uux,ddx;
  plural float cearly, cback;

  plural float sss1[33],vvv1[33],cearly_big[33];

  float u,d,dt,r1,tmp1,tmp2,rx,xreps,call;
  int i,n,nreps1,start,end,repsdiv,reps2,nreps2,j;
  int m,eproc,ncol,nrow,early;

  xreps = (float)reps;
  dt = t/xreps;
  rx = f_pow((r + 1.0),dt);
  
  tmp1 = (xmu - 0.5*psi*psi)*dt;
  tmp2 = cor*psi*f_sqrt(dt);
  u = f_exp(tmp1 + tmp2);
  d = f_exp(tmp1 - tmp2);
  
  vv1 = ss1 = cearly = cback = 0.0;
  proc[0].ss1 = s; proc[0].vv1 = v*v;

  n = m = 1;
  eproc = 0;
  ncol = nxproc/2;
  nrow = nyproc/2;
  early = 0;
  for (i=1; i<=13; i++) {
     aatmp1 = (r - 0.5*vv1)*dt;
     aatmp2 = fp_sqrt(vv1*dt);
     uux = fp_exp(aatmp1 + aatmp2);
     ddx = fp_exp(aatmp1 - aatmp2);
     if (n <= ncol) {
	if (ixproc<n) {
	  xnetE[n].vv1 = d*vv1;
          vv1 = u*vv1;
	  xnetE[n].ss1 = ddx*ss1;
          ss1 = uux*ss1;
        }; 
        if (i == reps1) {
            eproc = 2*n;
            if (ixproc < eproc) {
               cearly = ss1 - x;
               ss1 = ss1 - div;
            }
        };         
        n = 2*n;
     }
     else {
	if (iyproc<m) {
	  xnetS[m].vv1 = d*vv1;
          vv1 = u*vv1;
	  xnetS[m].ss1 = ddx*ss1;
          ss1 = uux*ss1;
        }; 
	if (i == reps1) {
          eproc = 2*m;
          early = 1;
          if (iyproc < eproc) {
            cearly = ss1 - x;
            ss1 = ss1 - div;
          }
	};
        m = 2*m;
     }        
  };

  sss1[1] = ss1; vvv1[1] = vv1;

  reps2 = reps - 13; 
  nreps1 = (int)f_pow(2.0,(float)(reps2+1)) - 1;
  repsdiv = reps1 - 14;
  start = (int)f_pow(2.0,(float)repsdiv);
  end = 2*start - 1;
  
  for (i=1; i<=nreps1/2; i++) {
     aatmp1 = (r - 0.5*vvv1[i])*dt;
     aatmp2 = fp_sqrt(vvv1[i]*dt);
     uux = fp_exp(aatmp1 + aatmp2);
     ddx = fp_exp(aatmp1 - aatmp2);
     vvv1[2*i] =  u*vvv1[i];  
     vvv1[2*i+1] = d*vvv1[i];
     if (i>=start && i<= end) {
       aatmp1 = uux*sss1[i];
       aatmp2 = ddx*sss1[i];
       cearly_big[2*(i-start)+1] = aatmp1 - x ;
       cearly_big[2*(i-start+1)] = aatmp2 - x ;
       sss1[2*i] =  aatmp1 - div ;
       sss1[2*i+1] = aatmp2 - div ;
     }
     else {
       sss1[2*i]   = uux*sss1[i];
       sss1[2*i+1] = ddx*sss1[i];
     };
     sss1[i] = 0.0;
   };
 
  for (i=nreps1/2; i<=nreps1; i++) {  
     sss1[i] -= x ;
     if (sss1[i] < 0.0) sss1[i] = 0.0;
   };

  if (reps1 <= 13) {
    for (i=nreps1/2+1; i<=nreps1; i++) 
      sss1[1] += sss1[i];
    sss1[1] = sss1[1]/f_pow(2.0*rx,(float)(reps-13));
    if (early == 0) {
      for (i = nyproc/2; i>=1; i = i/2) 
	if (iyproc <i) sss1[1] += xnetS[i].(sss1[1]) ;
      for (i = nxproc/2; i>=eproc; i = i/2)
	if (ixproc <i) sss1[1] += xnetE[i].(sss1[1]) ;
      if (iproc >= eproc) sss1[1] = 0.0;
    }
    else {
      for (i = nyproc/2; i>=eproc; i = i/2) 
	if (iyproc <i) sss1[1] += xnetS[i].(sss1[1]) ;
      if (iyproc >= eproc) sss1[1] = 0.0;
    };
    
    cback = sss1[1] / f_pow(2.0*rx,(float)(13-reps1));
    if (cearly > cback) cback = cearly;
    if (cback < 0.0) cback = 0.0;
    call = reduceAddf(cback)/f_pow(2.0*rx,(float)reps1);
  }
  else {
    nreps2 = (int) f_pow(2.0,(float)(reps-reps1));
    call = 0.0;
    for (i=1; i<=(2*start); i++) {
      aatmp1 = 0.0;
      n = nreps1/2 + (i-1)*nreps2;
      for (j=1; j<=nreps2; j++)
	aatmp1 += sss1[n+j];
      aatmp1 = aatmp1/f_pow(2.0*rx,(float)(reps-reps1));
      if (cearly_big[i] > aatmp1) aatmp1 = cearly_big[i] ;
      if (aatmp1 < 0.0) aatmp1 = 0.0;
      call = call + reduceAddf(aatmp1);
    };
    call = call/f_pow(2.0*rx,(float)reps1);
  };

  return call;
}


/* European Model */
float european(s,x,r,t,v,psi,xmu,cor,reps)
  float s,x,r,t,v,psi,xmu,cor;
  int reps;
{
  plural float ss1,vv1,aatmp1,aatmp2,uux,ddx;

  plural float sss1[33],vvv1[33];

  float u,d,dt,r1,tmp1,tmp2,rx,xreps,call;
  int i,n,nreps1,start,end,repsdiv,reps2,nreps2,j;
  int m,eproc,ncol,nrow;

  xreps = (float)reps;
  dt = t/xreps;
  rx = f_pow((r + 1.0),dt);
  
  tmp1 = (xmu - 0.5*psi*psi)*dt;
  tmp2 = cor*psi*f_sqrt(dt);
  u = f_exp(tmp1 + tmp2);
  d = f_exp(tmp1 - tmp2);
  
  vv1 = ss1 = 0.0;
  proc[0].ss1 = s; proc[0].vv1 = v*v;

  n = m = 1;
  ncol = nxproc/2;
  nrow = nyproc/2;
  for (i=1; i<=13; i++) {
     aatmp1 = (r - 0.5*vv1)*dt;
     aatmp2 = fp_sqrt(vv1*dt);
     uux = fp_exp(aatmp1 + aatmp2);
     ddx = fp_exp(aatmp1 - aatmp2);
     if (n <= ncol) {
	if (ixproc<n) {
	  xnetE[n].vv1 = d*vv1;
          vv1 = u*vv1;
	  xnetE[n].ss1 = ddx*ss1;
          ss1 = uux*ss1;
	};
        n = 2*n;
     }
     else {
	if (iyproc<m) {
	  xnetS[m].vv1 = d*vv1;
          vv1 = u*vv1;
	  xnetS[m].ss1 = ddx*ss1;
          ss1 = uux*ss1; 
	};
        m = 2*m;
     };
  };

  reps2 = reps - 13; 
  nreps1 = (int)f_pow(2.0,(float)(reps2+1)) - 1;

  sss1[1] = ss1; vvv1[1] = vv1;  

  for (i=1; i<=nreps1/2; i++) {
     aatmp1 = (r - 0.5*vvv1[i])*dt;
     aatmp2 = fp_sqrt(vvv1[i]*dt);
     uux = fp_exp(aatmp1 + aatmp2);
     ddx = fp_exp(aatmp1 - aatmp2);
     vvv1[2*i] =  u*vvv1[i];  
     vvv1[2*i+1] = d*vvv1[i];
     sss1[2*i]   = uux*sss1[i];
     sss1[2*i+1] = ddx*sss1[i];
     sss1[i] = 0.0; 
   };
 
  call = 0.0;
  for (i=nreps1/2; i<=nreps1; i++) {  
     sss1[i] -= x ;
     if (sss1[i] < 0.0) sss1[i] = 0.0;
     call += reduceAddf(sss1[i]);
   };

  call = call/f_pow(2.0*rx,(float)xreps);

  return call;
}