#include <stdio.h>
#include <math.h>
#include "StatisticalFunctions.h"
#include "HypothesisTesting.h"

/* Complement of the error function approximation

    erfc(z) = t * exp(-z*z)    * exp( a0 + a1*t + a2*t^2 + a3*t^3 ...
 where t = 1/(1 + z/2)

  Error:  1.2e-7 maximal relative error for z >= 0
  Source:  Numerical Recipies, Chapter 6, Special functions  */

double erfc (double z) {

    const double a0 = -1.26551223;
    const double a1 =  1.00002368;
    const double a2 =  0.37409196;
    const double a3 =  0.09678418;
    const double a4 = -0.18628806;
    const double a5 =  0.27886807;
    const double a6 = -1.13520398;
    const double a7 =  1.48851587;
    const double a8 = -0.82215223;
    const double a9 =  0.17087277;
    const double t = 2./(2. + z);
    return t*exp(-z*z+a0+t*(a1+t*(a2+t*(a3+t*(a4+t*(a5+t*(a6+t*(a7+t*(a8+t*a9)))))))));

};




/* Gamma function approximation

    gamma(z) = ( (sqrt(2*pi)/z) * ( p0 + (p1/z+1) + ... + (p6/z+6) )) *
       ( (z+5.5)^(z+0.5) ) * exp( -1*(z+5.5) )
    
  Error:  Unknown
  Source:  Unknown  */

double gamma (double z) {

    const double p0 = 1.000000000190015;
    const double p1 = 76.18009172947146;
    const double p2 = -86.50532032941677;
    const double p3 = 24.01409824083091;
    const double p4 =  -1.231739572450155;
    const double p5 = 1.208650973866179*0.001;
    const double p6 = -5.395239384953*0.000001;

    double P = p0 + (p1/(z+1)) + (p2/(z+2)) + (p3/(z+3)) + (p4/(z+4)) + (p5/(z+5)) + (p6/(z+6));
    return (sqrt(2*M_PI)/z)*P*pow(z+5.5,z+0.5)*exp( -1*(z+5.5) );

};




/* Incomplete gamma function approximation

    incgamma(x,a) = (x^a)*exp(-x)*sum(n=0:infinity){(x^n)/(a*(a+1)*...*(a+n))}

  Error:  Arbitrary.  As accurate as double allows
  Source:  Unknown  */

double incgamma (double x, double a) {

    printf("Calculating incgamma\n");
    double sum = 0;
    double term = 1.0/a;
    int n = 1;
    while (term != 0) {
 sum = sum + term;
 term = term*(x/(a+n));
 n++;
    }
    return pow(x,a)*exp(-1*x)*sum;

}



/* Chi-square distribution approximation

    chidist(x,v) = ( exp((-1*x)/2) * x^( (v/2) - 1) ) / ( 2^(v/2) * gamma(v/2) )
 where x >= 0

 Error:  Dependant on the gamma and incomplete gamma functions
 Source:  Standard Statistical Function  */

double chidist (double x, double v) {
    return incgamma(x/2.0,v/2.0)/gamma(x/2.0);

}

import java.math.*;
//import java.lang.Math.*;

public class IncGamma {
    
    public static double calculate(double a, double x) {
 double sum = 0;
 long n = 1;
 double term = 1/a;
 while (term != 0) {
     sum = sum + term;
     term = term*(x/(a+n));
     n++;
 }
 return java.lang.Math.pow(x,a)*java.lang.Math.exp(-1*x)*sum;
    }
    
    public static void main(String args[]) {
 System.out.println("IncGamma(1,2) = " + calculate(1,2));
 System.out.println("IncGamma(2,2) = " + calculate(2,2));
 System.out.println("IncGamma(3,2) = " + calculate(3,2));
 System.out.println("IncGamma(4,2) = " + calculate(4,2));
 System.out.println("IncGamma(5,2) = " + calculate(5,2));
    }
}