#include "../rng.h"
#include "./incomplete_gamma.h"
#include <iostream>
#define PI 3.1415926535897932384626433832795028841971693993751

// using the power series proves to be very fast and very accurate.
double igam(double s, double z){


  if( (s<=0 || z<=0)){
      return 0.0;
  }

  int rounds = 50;

  double sum = 0;

  for(int k=0; k<rounds; k++){
    double toadd = (std::pow(z,s)*std::exp(-z)*std::pow(z,k));
    for(int i=0; i<=k; i++){
      toadd/=(s+i);
    }
    sum+=toadd;
  }


  return sum / std::tgamma(s);

}

// because we are dealing with regularized incomplete gamma function, we can minus one from the lower
double igamc(double s, double z){
  return 1-igam(s,z);
}

// alternative expansion found on wikipedia
double igam_chg(double a, double x){

  if( (x<=0 || a<=0)){
      return 0.0;
  }

  double sum = 0;

  int rounds = 10000000;

  for(int i=0; i<rounds; i++){
    sum+=(std::pow(-1,i)/std::tgamma(i+1))*(std::pow(x,a+i)/(a+i));
  }
  return sum;

}

// This uses simpsons composite rule and passed at a 0.01 accuracy test but I want better.
double igam_simpsons(double a, double x){
  if( (x<=0 || a<=0)){
      return 0.0;
  }
  int parts = 10000000;

  double sum = 0;

  auto func = [a](double t) {
    double ret = std::pow(t, a-1) * std::exp(-t);
    if(std::isinf(ret)) ret=0;
    return ret; 
  };

  double h = x/parts;
  for(int i=1; i<=parts/2; ++i){
    sum+= (func((2*i-2)*h) + 4*func(((2*i)-1)*h) + func((2*i)*h));
  }
  return (sum*h)/(std::tgamma(a)*3.0);
}


// This uses the trapezium approximation which is less accurate and slower
double igam_trapezoid(double a, double x){

    if( (x<=0 || a<=0)){
        return 0.0;
    }
  double parts = 100000000;

  // reimann sum the lower incomplete gamma function
  double sum = 0;
  auto func = [a](double t) {
    double ret = std::pow(t, a-1) * std::exp(-t);
    if(std::isinf(ret)) ret=0;
//    std::cout << std::format("func({}) -> {}", t,ret);
    return ret; 
  };


  sum+=func(0);
  sum+=func(x);
  for(int i=1; i<parts; i++){
    sum+=2*func((x/parts)*i);
  }
  double ans = ((x/parts)/2)*sum;

  return (1/std::tgamma(a))*ans;
}