#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <getopt.h>
#include <iostream>
#include <ctime>

#include "gmp.h"
#include "mpfr.h"
#include "gs.h"
#include "bessel.h"

using namespace std;

int    M              = 50;
int    Num            = 60;
int    Prec           = 400;
bool   SP             = false;  // single precision
bool   Verbose        = false;
bool   Debug          = false;
bool   Test           = false;
double X              = -1.0;  // If negative, use [ GridStart, GridEnd ];
double GridStart      = 0.33;
double GridEnd        = 4.0;
double GridStep       = 0.01;
bool   NoInstanceTime = false;
bool   NoOverallTime  = false;
bool   PrnFuncValOnly = false; //true;
bool   help           = false;

int toint (const char* const s) {
  size_t len = strlen(s);
  char *ptr;
  long val;
  val = strtol(s, &ptr, 0);
  if (ptr != (s+len)) {
    fprintf(stderr,"Not an integer argument.");
    exit(-1);
  }
  return val;
}
double todouble(const char*const s) {
  size_t len = strlen(s);
  char *ptr;
  double val;
  val = strtod(s, &ptr);
  if (ptr != (s+len)) {
    fprintf(stderr,"Not an floating point argument.");
    exit(-1);
  }
  return val;
}

void mp_g(mpfr_t& out, const mpfr_t& mpx, long N) {
  mpfr_set_default_prec(Prec);
  if(mpfr_cmp_d(mpx,1.0) < 0) {
    N=int(10/mpfr_get_d(mpx,GMP_RNDN))+N;
  }

  ////// 2.0 * sqrt(2.0) * x   ///
  mpfr_t mp_sq2; mpfr_init_set_str(mp_sq2,"2.0",10,GMP_RNDN);
  mpfr_sqrt(mp_sq2,mp_sq2,GMP_RNDN);
  mpfr_t mp_twoSq2X; mpfr_init_set_str(mp_twoSq2X,"2.0",10,GMP_RNDN);
  mpfr_mul(mp_twoSq2X,mp_twoSq2X,mp_sq2,GMP_RNDN);
  mpfr_mul(mp_twoSq2X,mp_twoSq2X,mpx,GMP_RNDN);

  mpfr_t mp_twoSq2XN; mpfr_init2(mp_twoSq2XN,Prec);

  mpfr_t b1;  mpfr_init2(b1,Prec);
  mpfr_t b0;  mpfr_init2(b0,Prec);
  mpfr_t ans; mpfr_init2(ans,Prec);

  mpfr_t mpH; mpfr_init_set_str(mpH,"0.0",10,GMP_RNDN);

  mpfr_prec_t p = mpfr_get_prec (mpx);

  for(int n=1;n < N;++n) {

    std::stringstream ss;

    mpfr_mul_ui(mp_twoSq2XN,mp_twoSq2X,n,GMP_RNDN);

    char* str_twoSq2XN = new char[Prec+1];
    mp_exp_t exponent = 0;
    mp_exp_t* expptr = &exponent;
    mpfr_get_str (str_twoSq2XN, expptr, 10, Prec, mp_twoSq2XN, GMP_RNDN);
    if(str_twoSq2XN == NULL ){
      printf("ERROR: Could not convert mp_twoSq2XN to string\n");
      exit(-1);
    }

    // Ugly workaround to initialize ALGLIB data structures from MPFR
    std::string expStr(str_twoSq2XN);
    ss << "0." << expStr<<"e"<< (long)(*expptr) ;
    std::string newstr = ss.str();

    amp::ampf<4000> amp_twoSq2XN(ss.str());
    amp::ampf<4000> amp_b1 = bessel::besselk1(amp_twoSq2XN);

    mpfr_set(b1,amp_b1.getReadPtr(),GMP_RNDN);

    mpfr_mul(ans,b1,mp_twoSq2XN,GMP_RNDN);

    amp::ampf<4000> amp_b0 = bessel::besselk0(amp_twoSq2XN);
    mpfr_set(b0,amp_b0.getReadPtr(),GMP_RNDN);

    delete [] str_twoSq2XN;

    mpfr_sub(ans,ans,b0,GMP_RNDN);
    mpfr_mul_ui(ans,ans,4,GMP_RNDN);
    mpfr_add(mpH,mpH,ans,GMP_RNDN);

  }

  mpfr_neg(mpH,mpH,GMP_RNDN);
  mpfr_exp(out,mpH,GMP_RNDN);

  mpfr_clear(b0);
  mpfr_clear(b1);
  mpfr_clear(ans);
  mpfr_clear(mpH);
  mpfr_clear(mp_twoSq2X);
  mpfr_clear(mp_twoSq2XN);

}




int main(int argc,char*argv[]) {
  while (1) {
    static struct option long_options[] = {
      {"help",    no_argument,        0,  'h' },
      {"M",       required_argument,  0,  'm' },
      {"Num",     required_argument,  0,  'n' },
      {"Prec",    required_argument,  0,  'p' },
      {"X",       required_argument,  0,  'x' },
      {"verbose", optional_argument,  0,  'v' },
      {"debug", optional_argument,    0,  'd' },
      {"test", optional_argument,     0,  't' },
      {"noInstanceTime", optional_argument,     0,  'i' },
      {"noOverallTime", optional_argument,     0,  'o' },
      {"prnFuncValOnly",   optional_argument,  0,  'y' },
      {"!prnFuncValOnly",   optional_argument,  0,  'z' },
      {"SP",       optional_argument, 0,  's' },
      {0,0,0,0} /* This is a filler for -1 */
    };
    int option_index = 0;
    int c = getopt_long (argc, argv, "0:f:h:m:n:p:x:vdtsioyz", long_options, &option_index);
    if (c == -1) break;
    switch (c) {
    case '0':
      GridStart   = todouble(optarg);
      break;
    case 'f':
      GridEnd   = todouble(optarg);
      break;
    case 'h':
      GridStep   = todouble(optarg);
      break;
    case 'm':
      M = toint(optarg);
      break;
    case 'n':
      Num = toint(optarg);
      break;
    case 'p':
      Prec = toint(optarg);
      break;
    case 'v' :
      Verbose = true;
      break;
    case 'd' :
      Debug   = true;
      break;
    case 't' :
      Test    = true;
      break;
    case 's' :
      SP = true;
      break;
    case 'x' :
      X = todouble(optarg);
      break;
    case 'i' :
      NoInstanceTime   = true;
      break;
    case 'o' :
      NoOverallTime   = true;
      break;
    case 'y':
      PrnFuncValOnly = true;
    case 'z':
      PrnFuncValOnly = false;
    default:
      fprintf(stderr,"Unknown option '%c'\n",c);
      return 1;
      break;
    }
  }

  Laplace laInst(M,Prec,Debug,SP);
  if(Verbose) {
    printf("Precision = %d\n",Prec);
    printf("M = %d\n",M);
    printf("Num = %d\n",Num);
    printf("GridStart (-0) = %lf\n",GridStart);
    printf("GridEnd   (-f) = %lf\n",GridEnd);
    printf("GridStep  (-h) = %lf\n",GridStep);
    printf("SumXsi = %lf\n", laInst.sumXsi());
    printf("Xsi = %s\n", laInst.xsiStr().c_str());
  }

  mpfr_t laVal; mpfr_init2(laVal,Prec);
  if(X!=-1.0) {
    if(!PrnFuncValOnly) {
      printf("%.2f\t",X);
    }
    laInst.MPInverseTransform(laVal,mp_g, X, Num);
    mpfr_out_str(NULL,10,0,laVal,GMP_RNDN);
    printf("\n");
  } else { // use grid
    time_t startT, endT;
    time(&startT);
    double epsilon = 1e-10;
    for(double x = GridStart; x <= GridEnd + epsilon; x +=GridStep) {
      time_t i_startT, i_endT;
      time(&i_startT);
      laInst.MPInverseTransform(laVal,mp_g, x, Num);
      printf("%lf\t",x);
      mpfr_out_str(NULL,10,100,laVal,GMP_RNDN);
      time(&i_endT);
      int i_sec = i_endT - i_startT;
        if(!NoInstanceTime) printf("\t%d", i_sec);
        printf("\n");
    }
    time(&endT);
    int sec = endT - startT;
    if(!NoOverallTime) printf("Number of seconds elapsed = %d\n", sec);
  }
  mpfr_clear(laVal);

  return 0;
}