// Standard Precision GS routines adapted from C# from http://www.codeproject.com/KB/recipes/LaplaceTransforms.aspx   Included here because they are much more readable than their MP counterparts (Also, InverseTRansform() can be used for SP computation.)

#include <cmath>
#include <sstream>

#include "gmp.h"
#include "mpfr.h"

const int DefaultStehfest=14;

class Laplace {
  bool debug_;
  bool sp_;
  int prec_;
 public:

  typedef double (*fctptr_t)(double, long);
  typedef void (*mpfctptr_t)(mpfr_t&, const mpfr_t&, long);
  double   * V_; // Stehfest coefficients
  mpfr_t    * MPV_;
  double     ln2_; // log of 2
  mpfr_t      mpln2_;
  int        NV_;
  double     cumulV_;
  mpfr_t      mpCumulV_;

  double sumXsi() {
    if(sp_) {
      return cumulV_;
    } else {
      return mpfr_get_d(mpCumulV_,GMP_RNDN);
    }
  }

  std::string xsiStr() {
    std::stringstream ss;
    std::string res = ss.str();

    if(sp_) {
      for(int i = 0; i < NV_; ++i)  {
        if(i != 0) ss << " ";
        ss << V_[i];
      }
      res = ss.str();
    } else {
      for(int i = 0; i < NV_; ++i)  {
        if(i != 0) ss << " ";
        ss <<  mpfr_get_d(MPV_[i],GMP_RNDN);
      }
      res = ss.str();
    }
    return res;
  }

  void updV(mpfr_t& V, unsigned long M2, unsigned long j, unsigned i) {
    mpfr_t mpf1,mpf6,mpf7;
    mpz_t mpf2,mpf3,mpf4,mpf5;
    mpfr_init2(mpf1,prec_);
    mpz_init2(mpf2,prec_);
    mpz_init2(mpf3,prec_);
    mpz_init2(mpf4,prec_);
    mpz_init2(mpf5,prec_);
    mpfr_init2(mpf6,prec_);
    mpfr_init2(mpf7,prec_);

    mpfr_set_si(mpf1,j,GMP_RNDN);
    mpfr_pow_ui(mpf1,mpf1,M2+1,GMP_RNDN);
    mpz_fac_ui(mpf2,M2);
    mpz_bin_uiui(mpf3,M2,j);
    mpz_bin_uiui(mpf4,2*j,j);
    mpz_bin_uiui(mpf5,j,i+1-j);
    mpz_mul(mpf3,mpf3,mpf4);
    mpz_mul(mpf3,mpf3,mpf5);

    mpfr_set_z(mpf6,mpf2,GMP_RNDN);
    mpfr_div(mpf1,mpf1,mpf6,GMP_RNDN);
    mpfr_set_z(mpf7,mpf3,GMP_RNDN);
    mpfr_mul(mpf1,mpf1,mpf7,GMP_RNDN);

    mpfr_add(V,V,mpf1,GMP_RNDN);
    mpfr_clear(mpf1);
    mpz_clear(mpf2);
    mpz_clear(mpf3);
    mpz_clear(mpf4);
    mpz_clear(mpf5);
    mpfr_clear(mpf6);
    mpfr_clear(mpf7);
  }

  void updV(double& V, int M2, int j, int i) {
    V += (pow((long double)j, (long double)M2) / Factorial(j)) * (Factorial(2 * j) / Factorial(2 * j - i - 1)) / Factorial(M2 - j) / Factorial(j - 1) / Factorial(i + 1 - j);
  }

  void InitStehfest(int M)  {
    ln2_ = log(2.0);
    int M2 = M;
    NV_ = 2 * M2;
    V_ = new double[NV_];
    cumulV_=0.0;
    int sign = 1;
    if ((M2 % 2) != 0) sign = -1;
    for (int i = 0; i < NV_; i++) {
      int jmin = (i + 2) / 2;
      int jmax = i + 1;
      if (jmax > M2) jmax = M2;
      V_[i] = 0;
      sign = -sign;
      for (int j = jmin; j <= jmax; ++j) {
        updV(V_[i],M2,j,i);
      }
      V_[i] = sign * V_[i];
      cumulV_+=V_[i];
      if(debug_) {
        printf("V_[%d]=%lf; cumulV_=%lf\n",i,V_[i],cumulV_);
      }
    }
  }

  void MPInitStehfest(int M)  {
    ln2_ = log(2.0);
    mpfr_init_set_str(mpln2_,".6931471805599453094172321214581765680755001343602552541206800094933936219696947156058633269964186875",10,GMP_RNDN);
    int M2 = M; //M / 2;
    NV_ = 2 * M2;
    MPV_ = new mpfr_t[NV_];
    mpfr_init_set_str(mpCumulV_,"0.0",10,GMP_RNDN);
    int sign = 1;
    if ((M2 % 2) != 0) sign = -1;
    for (int i = 0; i < NV_; i++) {
      int jmin = (i + 2) / 2;
      int jmax = i + 1;
      if (jmax > M2) jmax = M2;
      mpfr_init_set_str(MPV_[i],"0.0",10,GMP_RNDN);
      sign = -sign;
      for (int j = jmin; j <= jmax; ++j) {
        updV(MPV_[i],M2,j,i);
      }
      if(sign<0) mpfr_neg(MPV_[i],MPV_[i],GMP_RNDN);
      mpfr_add(mpCumulV_,mpCumulV_,MPV_[i],GMP_RNDN);
    }
  }


  double InverseTransform(fctptr_t f, double x, int Num) {
    double ln2x = ln2_ / x;
    double xx = 0;
    double y = 0;
    for (int i = 0; i < NV_; i++) {
      xx += ln2_;
      double s=sqrt(xx)*x;
      double fval=(*f)(s,Num);
      double prod =      V_[i] / (i+1) *fval;
      y += prod;
    }
    return y;
  }

  void MPInverseTransform(mpfr_t& out, mpfctptr_t f, double x, int Num) {
    mpfr_t mpx; mpfr_init2(mpx,prec_);
    mpfr_t s; mpfr_init2(s,prec_);
    mpfr_set_d(mpx,x,GMP_RNDN);
    mpfr_t nln2; mpfr_init_set_str(nln2,"0.0",10,GMP_RNDN);
    mpfr_t sqrtnln2; mpfr_init_set_str(sqrtnln2,"0.0",10,GMP_RNDN);
    mpfr_t y;
    mpfr_init_set_str(y,"0.0",10,GMP_RNDN);
    mpfr_t mpval;
    mpfr_init2(mpval,prec_);
    mpfr_t prod;
    mpfr_init2(prod,prec_);

    for (int k = 0; k < NV_; ++k) {
      if(debug_) printf("m\t%d\n",k);
      //  xx += ln2_;   //
      mpfr_add(nln2,nln2,mpln2_,GMP_RNDN);
      //  s  <-- sqrt(xx)*x; //
      mpfr_sqrt(sqrtnln2,nln2,GMP_RNDN);
      mpfr_mul(s,mpx,sqrtnln2,GMP_RNDN);

      // Evaluate function
      (*f)(mpval,s,Num);

      mpfr_div_ui(prod,MPV_[k],(unsigned long)(k+1),GMP_RNDN);
      mpfr_mul(prod,prod,mpval,GMP_RNDN);
      double sp_prod = mpfr_get_d(prod,GMP_RNDN);
      mpfr_add(y,y,prod,GMP_RNDN);
      double sp_y = mpfr_get_d(y,GMP_RNDN);
      sp_y += sp_prod;
    }

    // Cleanup
    mpfr_clear(mpval);
    mpfr_clear(prod);
    mpfr_set(out,y,GMP_RNDN);
  }

  static double Factorial(int N) {
    double x = 1;
    if (N > 1) {
      for (int i = 2; i <= N; i++) x = i * x;
    }
    return x;
  }

 Laplace(int M=DefaultStehfest, int prec=256, bool debug=false, bool sp=false):
  debug_(debug), sp_(sp) , prec_(prec) {
    if(sp_) {
      InitStehfest(M);
    } else {
      mpfr_set_default_prec(prec);
      MPInitStehfest(M);
    }
  }

  ~Laplace() {
    if(!sp_) {
      for (int i = 0; i < NV_; i++) {
        mpfr_clear(MPV_[i]);
      }
      delete [] MPV_;
      mpfr_clear(mpCumulV_);
      mpfr_clear(mpln2_);
    }
  }


};