#include <iostream>
#include <cstdlib>
#include <cmath>
#include <stdio.h>

#include "fastmath.h"

fastmath* fastmath::staticthis = NULL;

const double fastmath::atanxmax=200.;

const double fastmath::expxrang=100.;
const double fastmath::expxmin=-90.;
const double fastmath::expshiftrgt=std::exp(fastmath::expxrang);
const double fastmath::expshiftlft=1./fastmath::expshiftrgt;

const double fastmath::logxrang=100.;
const double fastmath::logxrangrecp=1./fastmath::logxrang;
const double fastmath::logxmin=0.1;
const double fastmath::logshift=std::log(fastmath::logxrang);

double fastmath::atanbinwrecp;
double fastmath::atanTbl[fastmath_nxbins];

double fastmath::expbinwrecp;
double fastmath::expTbl[fastmath_nxbins];

double fastmath::logbinwrecp;
double fastmath::logTbl[fastmath_nxbins];

fastmath::fastmath(){
  
  std::cout << "Initializing fastmath...";
  
  //Initializing atan
  atanbinwrecp=atanxmax/(double)(fastmath_nxbins-1);// atan bin width
  for (int i=0; i<fastmath_nxbins; i++) {
    atanTbl[i]=std::atan(i*atanbinwrecp);
  }
  atanbinwrecp=1./atanbinwrecp;// reciprocal of atan bin width
  
  //Initializing exp
  expbinwrecp=expxrang/(double)(fastmath_nxbins-1);// exp bin width
  for (int i=0; i<fastmath_nxbins; i++) {
    expTbl[i]=std::exp(i*expbinwrecp+expxmin);
  }
  expbinwrecp=1./expbinwrecp;// reciprocal of exp bin width
  
  //Initializing log
  logbinwrecp=(logxrang-1.)*logxmin/(double)(fastmath_nxbins-1);// log bin width
  for (int i=0; i<fastmath_nxbins; i++) {
    logTbl[i]=std::log(i*logbinwrecp+logxmin);
  }
  logbinwrecp=1./logbinwrecp;// reciprocal of log bin width
  
  std::cout << "complete!" << std::endl;
  
}

double fastmath::atan(double x) {
  
//  fastmath* dum=Get();
  
  double fsgn=1.;
  if (x<0.) {
    fsgn=-1.;
    x=-x;
  }
  
  if (!(x < atanxmax)) {// out of range!
    return fsgn*std::atan(x);
  }
  else {
    double dix=x*atanbinwrecp;
    int ix=(int)dix;
    double xd=dix-(double)ix;
    return fsgn*(atanTbl[ix]*(1.-xd)+atanTbl[ix+1]*xd);
  }
  
}

double fastmath::atan2(double y, double x) {
  const double pinum=3.1415926535897932;
  double phi=fastmath::atan(y/x);//[-pi/2,+pi/2]
  if (x<0.) {
    if (y<0.) {
      phi-=pinum;
    }
    else {
      phi+=pinum;
    }
  }
  
  return phi;//[-pi,+pi]
}

double fastmath::exp(double xtmp) {
  
//  fastmath* dum=Get();
  
  double x=xtmp-expxmin;
  
  int icount=5;
  double mulfact=1.;
  while (1) {
    icount--;
    if (icount) {
      if (x < 0.) {
        x+=expxrang;
        mulfact*=expshiftlft;
        continue;
      }
      else if (!(x < expxrang)) {
        x-=expxrang;
        mulfact*=expshiftrgt;
        continue;
      }
    }
    else {
//      std::cout << "fastmath::exp: too many iterations! x=" << xtmp << std::endl;
      return std::exp(xtmp);
    }
    break;
  }
  
  double dix=x*expbinwrecp;
  
  int ix=(int)dix;
  double xd=dix-(double)ix;
  return mulfact*(expTbl[ix]*(1.-xd)+expTbl[ix+1]*xd);
}

double fastmath::log(double xtmp) {
  
  //  fastmath* dum=Get();
  
  double x=xtmp;
  
  int icount=5;
  double shiftfact=0.;
  while (1) {
    icount--;
    if (icount) {
      if (x < logxmin) {
        x*=logxrang;
        shiftfact-=logshift;
        continue;
      }
      else if (!(x < logxmin*logxrang)) {
        x*=logxrangrecp;
        shiftfact+=logshift;
        continue;
      }
    }
    else {
//      std::cout << "fastmath::log x=" << xtmp << std::endl;
      return std::log(xtmp);
    }
    break;
  }
  
  double dix=(x-logxmin)*logbinwrecp;
  
  int ix=(int)dix;
  double xd=dix-(double)ix;
  return shiftfact+(logTbl[ix]*(1.-xd)+logTbl[ix+1]*xd);
}