#ifndef __METAPROG_H__
#define __METAPROG_H__

namespace Metaprog
{

//Sont fournis / are provided
//Pi (calcule PI / computes PI)
//Pow, Fact, FactEven, FactOdd, Exp,
//Sin, Cos, Tan, Asin, Acos, Atan, Sinh, Cosh, Tanh, Atanh
//FOR

//Metaprog_sqr
//Utilitaire, calcule x*x / utility, computes x*x
template<typename T> static inline T Metaprog_sqr(T x) {return x*x;}

//Pow<P,T>(T x) = x^P.
//Pas de restriction particuliere/There is no special restriction
template <int P,typename T>
inline T Pow(T x)
{
  return (P == 0) ? 1
                  : (
                     (P<0) ? 1/Pow<(P<0)?-P:0,T>(x)
                           : ((P%2) ? x*Metaprog_sqr(Pow<(P%2) ? (P-1)/2 : 0,T>(x))
                                    : Metaprog_sqr(Pow<(P%2) ? 0 : P/2,T>(x))
                             )
                    );
}
//end Pow<>()

//Fact<N,T>() = N*(N-1)*(N-2)*(N-3)*...*1 (result of type T)
//Attention : Si N est grand, un type integral pour T peut subir le debordement
//Caution : for a great N, an integral type for T may overflow 
template<unsigned long N,typename T>
inline T Fact()
{
  return (N==0) ? 1 : N*Fact<(N > 0 ? N-1 : 0),T>();
}
//end Fact<>()

//FactOdd<2N+1>() = (2N+1)*(2N-1)*(2N-3)*...*1
template<unsigned long N,typename T>
inline T FactOdd()
{
  return (N<=1) ? 1 : (N%2) ? N*FactOdd<(N>2 ? N-2 : 0),T>() : FactOdd<(N > 1 ? N-1 : 0),T>();
}
//end FactOdd<>()

//FactEven<2N>() = (2N)*(2N-2)*(2N-4)*...*2
template<unsigned long N,typename T>
inline T FactEven()
{
  return (N<=1) ? 1 : (!(N%2)) ? N*FactEven<(N>2 ? N-2 : 0),T>() : FactEven<(N > 1 ? N-1 : 0),T>();
}
//end FactEven<>()

//Pi
//Si on s'en contente, une constante de bonne precision/
//If it is enough, a constant with good precision
static const double PI = 3.14159265358979323846;

//Pi_ est juste un intermediaire de calcul/
//Pi_ is just an intermediary of calculation
template<unsigned int N>
inline double Pi_(void)
{
  return Pi_<N-1>()+Fact<2*N,double>()/
                    (Pow<4*N+1,double>(2)*
                    Metaprog_sqr(Fact<N,double>())*(2*N+1));
}
template<> inline double Pi_<0>(void) {return 1/double(2);}

//Pi
//Pas de restriction particuliere/There is no special restriction
template<unsigned int N> inline double Pi(void) {return 6*Pi_<N>();}

//Exp
//Exp_ est juste un intermediaire de calcul/
//Exp_ is just an intermediary of calculation
template<unsigned int N>
inline double Exp_(double x)
{
  return (N==0) ? 1 : Exp_<(N>0)?(N-1):0>(x)+Pow<N,double>(x)/Fact<N,double>();
}

template<unsigned int N>
inline float Exp_(float x)
{
  return (N==0) ? 1 : Exp_<(N>0)?(N-1):0>(x)+Pow<N,float>(x)/Fact<N,float>();
}

//Exp
//Pas de restriction particuliere/No special restriction
template<unsigned int N>
inline double Exp(double x) {return (x>=0) ? Exp_<N>(x) : (1/Exp_<N>(-x));}

template<unsigned int N>
inline float Exp(float x) {return (x>=0) ? Exp_<N>(x) : (1/Exp_<N>(-x));}

//Sin
//Pas de restriction particuliere/No special restriction
//Pour une bonne precision, garder x dans [-2*PI,2*PI]/
//For a good precision, keep x in [-2*PI,2*PI]
template<unsigned int N>
inline double Sin(double x)
{
  return (!(N%2)) ? Sin<(!(N%2))?N-1:0>(x)
                  : Sin<(N%2)?N-2:0>(x) +
                    ((((N-1)/2)%2)?-1:1)*Pow<(N%2)?N:0,double>(x)
                                        /Fact<(N%2)?N:0,double>();
}
template<> inline double Sin<1>(double x) {return x;}
template<> inline double Sin<0>(double x) {return 0;}

template<unsigned int N>
inline float Sin(float x)
{
  return (!(N%2)) ? Sin<(!(N%2))?N-1:0>(x)
                  : Sin<(N%2)?N-2:0>(x) +
                    ((((N-1)/2)%2)?-1:1)*Pow<(N%2)?N:0,float>(x)
                                        /Fact<(N%2)?N:0,float>();
}
template<> inline float Sin<1>(float x) {return x;}
template<> inline float Sin<0>(float x) {return 0;}
//end Sin

//Asin
//Requis : |x|<1 / Requires : |x|<1
//INCORRECT avec |x|=1 / INCORRECT for |x|=1
template<unsigned int N>
inline double Asin(double x)
{
  return (!(N%2)) ? Asin<(!(N%2))?N-1:0>(x)
                  : Asin<(N%2)?N-2:0>(x)+(Fact<(N%2)?N-1:0,double>()*
                                         Pow<(N%2)?N:0,double>(x))/
                                        (Pow<(N%2)?N-1:0,double>(2)*
                                         Metaprog_sqr(
                                           Fact<(N%2)?(N-1)/2:0,double>())*N);
}
template<> inline double Asin<1>(double x) {return x;}
template<> inline double Asin<0>(double x) {return 0;}

template<unsigned int N>
inline float Asin(float x)
{
  return (!(N%2)) ? Asin<(!(N%2))?N-1:0>(x)
                  : Asin<(N%2)?N-2:0>(x)+(Fact<(N%2)?N-1:0,float>()*
                                         Pow<(N%2)?N:0,float>(x))/
                                        (Pow<(N%2)?N-1:0,float>(2)*
                                         Metaprog_sqr(
                                           Fact<(N%2)?(N-1)/2:0,float>())*N);
}
template<> inline float Asin<1>(float x) {return x;}
template<> inline float Asin<0>(float x) {return 0;}
//end Asin

//Cos
//Pas de restriction particuliere/No special restriction
//Pour une bonne precision, garder x dans [-2*PI,2*PI]/
//For a good precision, keep x in [-2*PI,2*PI]
template<unsigned int N>
inline double Cos(double x)
{
  return (N%2) ? Cos<(N%2)?N-1:0>(x)
               : Cos<(N%2)?0:N-2>(x) + 
                 (((N/2)%2)?-1:1)*Pow<(N%2)?0:N,double>(x)/Fact<(N%2)?0:N,double>();
}
template<> inline double Cos<1>(double x) {return x;}
template<> inline double Cos<0>(double x) {return 1;}

template<unsigned int N>
inline float Cos(float x)
{
  return (N%2) ? Cos<(N%2)?N-1:0>(x)
               : Cos<(N%2)?0:N-2>(x) + 
                 (((N/2)%2)?-1:1)*Pow<(N%2)?0:N,float>(x)/Fact<(N%2)?0:N,float>();
}
template<> inline float Cos<1>(float x) {return x;}
template<> inline float Cos<0>(float x) {return 1;}
//end Cos

//Acos
//Requis : |x|<1 / Requires : |x|<1
//INCORRECT avec |x|=1 / INCORRECT for |x|=1
template<unsigned int N>
inline double Acos(double x)
{
  return (PI/2)-Asin<N>(x);
}

template<unsigned int N>
inline float Acos(float x)
{
  return (PI/2)-Asin<N>(x);
}
//end Acos

//Tan
//Pas de restriction particuliere (a part x != n*Pi/2) /
//No special restriction (but x != n*Pi/2)
//Pour une bonne precision, garder x dans [-2*PI,2*PI]/
//For a good precision, keep x in [-2*PI,2*PI]
template<unsigned int N> inline double Tan(double x)
{
  return Sin<N>(x)/Cos<N>(x);
}

template<unsigned int N> inline float Tan(float x)
{
  return Sin<N>(x)/Cos<N>(x);
}
//end Tan

//Atan
//Eviter les valeurs poches de +-Pi/4 / Avoid values close to +-Pi/4
template<unsigned int N> inline double Atan(double x)
{
  return ((x >= -1-(1e-10)) && (x<=-1+(1e-10))) ? -PI/4
                        : ((x >= 1-(1e-10)) && (x<=1+(1e-10))) ? PI/4
                             : ((x<-1) || (x>1)) ? ((x<=-1)?-1:1)*PI/2-Atan<N>(1/x)
                                        : (!(N%2)) ? Atan<(!(N%2))?N-1:0>(x)
                                              : Atan<(!(N%2))?0:N-2>(x)+
                                                ((((N-1)/2)%2)?-1:1)*
                                                  Pow<(!(N%2))?0:N>(x)/N;
}
template<> inline double Atan<1>(double x) {return x;}
template<> inline double Atan<0>(double x) {return 0;}

template<unsigned int N> inline float Atan(float x)
{
  return ((x >= -1-(1e-10)) && (x<=-1+(1e-10))) ? -PI/4
                        : ((x >= 1-(1e-10)) && (x<=1+(1e-10))) ? PI/4
                             : ((x<-1) || (x>1)) ? ((x<=-1)?-1:1)*PI/2-Atan<N>(1/x)
                                        : (!(N%2)) ? Atan<(!(N%2))?N-1:0>(x)
                                              : Atan<(!(N%2))?0:N-2>(x)+
                                                ((((N-1)/2)%2)?-1:1)*
                                                  Pow<(!(N%2))?0:N>(x)/N;
}
template<> inline float Atan<1>(float x) {return x;}
template<> inline float Atan<0>(float x) {return 0;}
//end Atan

//Sinh
//Pas de restriction particuliere/No special restriction
template<unsigned int N>
inline double Sinh(double x)
{
  return (!(N%2)) ? Sinh<(!(N%2))?N-1:0>(x)
                  : Sinh<(!(N%2))?0:N-2>(x) +
                    Pow<(!(N%2))?0:N,double>(x)/Fact<(!(N%2))?0:N,double>();
}
template<> inline double Sinh<1>(double x) {return x;}
template<> inline double Sinh<0>(double x) {return 0;}

template<unsigned int N>
inline float Sinh(float x)
{
  return (!(N%2)) ? Sinh<(!(N%2))?N-1:0>(x)
                  : Sinh<(!(N%2))?0:N-2>(x) +
                    Pow<(!(N%2))?0:N,float>(x)/Fact<(!(N%2))?0:N,float>();
}
template<> inline float Sinh<1>(float x) {return x;}
template<> inline float Sinh<0>(float x) {return 0;}
//end Sinh

//Cosh
//Pas de restriction particuliere/No special restriction
template<unsigned int N>
inline double Cosh(double x)
{
  return (N%2) ? Cosh<(N%2)?N-1:0>(x)
               : Cosh<(N%2)?0:N-2>(x) +
                 Pow<(N%2)?0:N,double>(x)/Fact<(N%2)?0:N,double>();
}
template<> inline double Cosh<1>(double x) {return 0;}
template<> inline double Cosh<0>(double x) {return 1;}

template<unsigned int N>
inline float Cosh(float x)
{
  return (N%2) ? Cosh<(N%2)?N-1:0>(x)
               : Cosh<(N%2)?0:N-2>(x) +
                 Pow<(N%2)?0:N,float>(x)/Fact<(N%2)?0:N,float>();
}
template<> inline float Cosh<1>(float x) {return 0;}
template<> inline float Cosh<0>(float x) {return 1;}
//end Cosh

//Tanh
template<unsigned int N> inline double Tanh(double x)
{
  return Sinh<N>(x)/Cosh<N>(x);
}

template<unsigned int N> inline float Tanh(float x)
{
  return Sinh<N>(x)/Cosh<N>(x);
}
//end Tanh

//Atanh
//Pas de restriction particuliere/No special restriction
template<unsigned int N> inline double Atanh(double x)
{
  return (!(N%2)) ? Atanh<(!(N%2))?N-1:0>(x)
                  : Atanh<(!(N%2))?0:N-2>(x)+Pow<(!(N%2))?0:N>(x)/N;
}
template<> inline double Atanh<1>(double x) {return x;}
template<> inline double Atanh<0>(double x) {return 0;}

template<unsigned int N> inline float Atanh(float x)
{
  return (!(N%2)) ? Atanh<(!(N%2))?N-1:0>(x)
                  : Atanh<(!(N%2))?0:N-2>(x)+Pow<(!(N%2))?0:N>(x)/N;
}
template<> inline float Atanh<1>(float x) {return x;}
template<> inline float Atanh<0>(float x) {return 0;}
//end Atanh

} //end namespace Metaprog

#endif