/*@@@@@@@@@@@@@@@@@@@@@@@@@@@ utilities.cpp @@@@@@@@@@@@@@@@@@@@@@@@@@@
 Random Option Pricing
 Summer 2002
 Authors: A. Korobchevsky
          S. Ferrando
 Ryerson University ALL RIGHTS RESERVED
 
 Purpose: Implements all utility functions required for the application 
 @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*/

#include <iostream>
#include "utilities.h"




/*************************** cumulative *****************************
 * Function: cumulative
 * Purpose: Computes normal cumulative function by using numerical approximation
 * Param: x - any real number
 * Returns: cumulative in range [0..1]
 ********************************************************************/

// BEGIN Cum
double cumulative(double x)
{  double width = x / 10000;
   if (width < 0) 
      width *= (-1);
   double sum = 0;

   for (int i = 0; i < 10000; i++)
   {  double value = 1 / sqrt(2 * M_PI) * exp((-1) * (width*i) * (width*i) / 2);
      sum += value * width;
   }
   if (x > 0)
      return 0.5 + sum;
   return 0.5 - sum;
}


// END Cum


/************************ NDApproximationProbability ****************
 * Function: NDApproximationProbability
 * Purpose: Computes normal approximation of binomial probability
 * Param: probability_up - probability of success
                       N - number of trials
                       n - number of successes
 * Returns: probability in range [0..1]
 ********************************************************************/
double NDApproximationProbability(double probability_up, int N, int n)
{ double mean, sd;
  
  mean = N * probability_up;
  sd = sqrt(N * probability_up * (1 - probability_up));

  return cumulative(( n + 0.5 - mean) / sd) - cumulative(( n - 0.5 - mean) / sd);
}




/************************ factorial ********************************
 * Function: factorial
 * Purpose: Computes the factorial (iterative non-recursive)
 * Param: i - argument
 * Returns: value of the factorial
 * Note: inpredictable for i >> 20
 ********************************************************************/
long factorial(int i)
{ long f = 1;
  for (int j = 1; j <= i; j++)
    f *= j;
  return f;
}




/************************ combinations ******************************
 * Function: factorial
 * Purpose: Computes combinations (recursive)
 * Param: N, k - combination standard parameters
 * Returns: number of combinations
 ********************************************************************/
long combinations(int N, int k) 
{ double comb = 1;
  for (int i = 0; i < k; i++)
    comb *= (double)(N - i) / (i + 1);
  return (long)comb;
}





/************************ directProbability *************************
 * Function: directProbability
 * Purpose: Computes binomial probability directly
 * Param: pU - probability of success
           N - number of trials
           n - number of successes
 * Returns: probability in range [0..1]
 * Note: Unreliable when N is large
 ********************************************************************/
double directProbability(double probability_up, int N, int n)
{ return combinations(N, n) * pow(probability_up, n) * pow(1 - probability_up, N - n);
}




/***************************** permutePath *************************
 * Function: directProbability
 * Author: Jihong Cai
 * Purpose: Computes the next path with same number of zeros in lexigraphical order 
 * Param: *path - pointer to path 
              N - length of the path
 * Returns: N/A
 ********************************************************************/
void permutePath(int *path,int N)
{
  int i,k,j=0;

  for(i=N-1; i>=0;i--)
    if(path[i]==1)
    {	j=i;
	   break;
    }
  
  for(i=j-1;i>=0;i--)
    if(path[i]==0)
    {	//swap
	   path[i] = 1;
	   path[i+1] = 0;
	
      //move back
  	   if( (i!=j-1) && (j!=N-1) )
	   { //move 1's to right
	     for(k=1; k < j-i; k++)
	       path[N-k] = 1;
	     //put 0's on left
	     for(k=1;k<N-j;k++)
	       path[i+1+k] = 0;
	   }
	   break;
    }
}




/************************ logProbability ****************************
 * Function: logProbability
 * Purpose: Computes logarithmic approximation of binomial probability
 * Param: pU - probability of success
           N - number of trials
           n - number of successes
 * Returns: probability in range [0..1]
 ********************************************************************/

// BEGIN Log
double logProbability(double pU, int N, int n)
{ double S = 0;
  for (int i = 0; i < n; i++)
    S += log ((double)(N - i) / (i + 1));
  return exp(N * log (1 - pU) + n * log (pU / (1 - pU)) + S); 
}
// END Log
