//#include "callfactory.h"
#include <zeta/module.h>
#include <cmath>
#include <cstdlib>
typedef double (*d_fun)();
typedef double (*d_fun_d)(double);
typedef double (*d_fun_dd)(double,double);
typedef double (*d_fun_ddd)(double,double,double);

typedef double (*d_fun_ud)(unsigned,double);
typedef double (*d_fun_id)(int,double);
typedef double (*d_fun_uud)(unsigned,unsigned,double);
typedef double (*d_fun_iid)(int,int,double);

#include <boost/math/special_functions/gamma.hpp>
#include <boost/math/special_functions/digamma.hpp>
#include <boost/math/special_functions/beta.hpp>
#include <boost/math/special_functions/erf.hpp>
#include <boost/math/special_functions/legendre.hpp>
#include <boost/math/special_functions/laguerre.hpp>
#include <boost/math/special_functions/hermite.hpp>

//Gamma functions
void gamma_functions(Module & mod){
  vector<string> X; X.push_back("x");
  vector<string> AB; AB.push_back("a"); AB.push_back("b");
  vector<string> AX; AX.push_back("a"); AX.push_back("x");
  vector<string> ABX; ABX.push_back("a"); ABX.push_back("b"); ABX.push_back("x");
  mod.def_function<d_fun_d>("tgamma",
			    "The tgamma() function returns the \\\"true gamma\\\" of x",
			    boost::math::tgamma,
			    X); 

  mod.def_function<d_fun_d>("tgamma1pm1",
			    "Return tgamma(x+1) - 1. Gives very accurate results for small x",
			    boost::math::tgamma1pm1,
			    X);

  mod.def_function<d_fun_d>("lgamma","Returns the logarithm of the absolute valaue of the gamma function.",
			    boost::math::lgamma,
			    X); 

  mod.def_function<d_fun_d>("digamma","Returns the digamma or psi function of x. Digamma is defined as the logarithmic derivative of the gamma function.",
			    boost::math::digamma,
			    X);

  mod.def_function<d_fun_dd>("tgamma_ratio","Returns the ratio of two gamma functions tgamma(a)/tgamma(b)",
			     boost::math::tgamma_ratio,AB);


  mod.def_function<d_fun_dd>("gamma_p",
			     "Returns the normalised lower incomplete gamma function of a and x",
			     boost::math::gamma_p,
			     AX);


  mod.def_function<d_fun_dd>("gamma_q","Returns the normalised upper incomplete gamma function of a and x",
			     boost::math::gamma_q,
			     AX);


  mod.def_function<d_fun_dd>("tgamma_lower",
			     "Returns the full (non-normalised) lower incomplete gamma function of a and x",
			     boost::math::tgamma_lower,
			     AX);


  mod.def_function<d_fun_dd>("tgamma_upper",
			     "Returns the full (non-normalised) upper incomplete gamma function of a and x",
			     boost::math::tgamma,
			     AX);


  mod.def_function<d_fun_dd>("gamma_q_inv",
			     "Returns a value x such that: q = gamma_q(a, x)",
			     boost::math::gamma_q_inv,
			     AX);

  mod.def_function<d_fun_dd>("gamma_p_inv",
			     "Returns a value x such that: p = gamma_p(a, x)",
			     boost::math::gamma_p_inv,
			     AX);

  mod.def_function<d_fun_dd>("gamma_q_inva",
			     "Returns a value a such that: q = gamma_q(a, x)",
			     boost::math::gamma_q_inva,
			     AX);

  mod.def_function<d_fun_dd>("gamma_p_inva",
			     "Returns a value a such that: p = gamma_p(a, x)",
			     boost::math::gamma_p_inva,
			     AX);

  mod.def_function<d_fun_dd>("gamma_p_derivative",
			     "Returns the partial derivative with respect to x of the incomplete gamma function gamma_p",
			     boost::math::gamma_p_derivative,
			     AX);

  //Beta
  mod.def_function<d_fun_dd>("beta",
			     "Returns the beta function of a nd b beta(a,b)=tgamma(a)tgamma(b)/tgamma(a+b)",
			     boost::math::beta,
			     AB);

  //Incomplete beta
  mod.def_function<d_fun_ddd>("ibeta",
			      "Returns the normalised incomplete beta function of a, b and x",
			      boost::math::ibeta,
			      ABX);

  mod.def_function<d_fun_ddd>("ibetac",
			      "Returns the normalised complement of the incomplete beta function of a, b and x",
			      boost::math::ibetac,
			      ABX);

  mod.def_function<d_fun_ddd>("nibeta",
			      "Returns the full (non-normalised) incomplete beta function of a, b and x",
			      boost::math::beta,
			      ABX);

  mod.def_function<d_fun_ddd>("nibetac",
			      "Returns the full (non-normalised) complement of the incomplete beta function of a, b and x",
			      boost::math::betac,
			      ABX);

  mod.def_function<d_fun_ddd>("ibeta_derivative",
			      "Return the partial derivative with respect to x of the incomplete beta function ibeta",
			      boost::math::ibeta_derivative,
			      ABX);
}



void error_functions(Module & mod){
  vector<string> X; X.push_back("x");


  mod.def_function<d_fun_d>("erf",
			    "Returns the error function erf of x",
			    boost::math::erf,
			    X);

  mod.def_function<d_fun_d>("erfc",
			      "Returns the complement of the error function of x",
			      boost::math::erfc,
			      X);

  mod.def_function<d_fun_d>("erf_inv",
			      "Returns the inverse error function of x",
			      boost::math::erf_inv,
			      X);

  mod.def_function<d_fun_d>("erfc_inv",
			      "Returns the inverse of the complement of the error function of x",
			      boost::math::erfc_inv,
			      X);


}

void polynomial_functions(Module &mod){
  vector<string> LX; LX.push_back("l"); LX.push_back("x");
  vector<string> LMX; LMX.push_back("l"); LMX.push_back("m"); LMX.push_back("x");
  vector<string> NX; NX.push_back("N"); NX.push_back("X");
  vector<string> NMX; NMX.push_back("N"); NMX.push_back("M");  NMX.push_back("X");

  mod.def_function<d_fun_id>("legendre_pl",
			     "Returns the Legendre Polynomial of the first kind",
			     boost::math::legendre_p,
			     LX);

  mod.def_function<d_fun_iid>("legendre_plm",
			      "Returns the associated Legendre polynomial of the first kind",
			      boost::math::legendre_p,
			      LMX);

  mod.def_function<d_fun_ud>("legendre_ql",
			      "Returns the value of the Legendre polynomial that is the second solution to the Legendre differential equation",
			      boost::math::legendre_q,
			      LX);



  mod.def_function<d_fun_ud>("laguerre_n",
			     "Returns the value of the Laguerre Polynomial of order n at point x",
			     boost::math::laguerre,
			     NX);

  mod.def_function<d_fun_uud>("laguerre_nm",
			      "Returns the Associated Laguerre polynomial of degree n and order m at point x",
			      boost::math::laguerre,
			      NMX);

  mod.def_function<d_fun_ud>("hermite",
			     "Returns the value of the Hermite Polynomial of order n at point x",
			     boost::math::hermite,
			      NX);

}

extern "C" {
  void init_math(){
    
    Module math = Module("math","Mathematical functions");
    vector<string> args;
    vector<string> args2;
    args.push_back("x");
    args2.push_back("y");
    args2.push_back("x");

    
    math
      .def_function<d_fun_d>("sin",
			     "The sin() function returns the sine of x, where x is given in radians.",sin,args)
      .def_function<d_fun_d>("cos",
			     "The cos() function returns the cosine of x, where x is given in radians.",
			     cos,args)
      .def_function<d_fun_d>("asin",
"The asin() function calculates the arc sine of x; that is the value whose cosine is x.",
			     asin,args)
      .def_function<d_fun_d>("acos",
			     "The acos() function calculates the arc cosine of x; that is the value whose cosine is x.",
			     acos,args)
      .def_function<d_fun_d>("atan",
			     "The  atan()  function  calculates  the  arc tangent of x; that is the value whose tangent is x.",
			     atan,args)
      .def_function<d_fun_dd>("atan2",
			      "The  atan2() function calculates the arc tangent of the two variables x and y.",
			      atan2,args2)
      .def_function<d_fun_dd>("pow","",pow)
      .def_function<d_fun_d>("exp","The  exp() function returns the value of e (the base of natural logarithms) raised to the power of x.",exp,args)
      .def_function<d_fun_d>("log","The log() function returns the natural logarithm of x.",log,args)
      .def_function<d_fun_d>("log10","he log10() function returns the base-10 logarithm of x.",log10,args)
      .def_function<d_fun_d>("sqrt","he  sqrt()  function  returns the non-negative square root of x.",sqrt,args)
      .def_function("drand48","The  drand48() function return non-negative double-precision floating-point values uniformly distributed between [0.0, 1.0).",drand48)
      .def_function("srand48","The srand48() functions  is  the initialization  functions which should be called before using drand48().",srand48,args)
      .def_function("rand","The rand() function returns a pseudo-random integer between 0 and RAND_MAX.",rand)
      .def_function("srand","The srand() function sets its argument as the seed for a  new  sequence  of pseudo-random  integers  to  be  returned  by  rand().",srand,args)
      .def_const("pi",M_PI)
      .def_const("e",M_E);



    //Gamma functions
    gamma_functions(math);

    //Factorials??

    //Error functions
    error_functions(math);

    //Polynomials
    polynomial_functions(math);
  }
}