math.h

#ifndef gridripper_math_h
#define gridripper_math_h

#include <gridripper/config.h>
#include <cmath>
#include <complex>
#include <cfloat>

namespace gridripper {

const GReal_t GNaN =       NAN;                                     /* Quiet NaN */
const GReal_t GE =         2.7182818284590452353602874713526625L;   /* e */
const GReal_t GLog2E =     1.4426950408889634073599246810018921L;   /* log_2 e */
const GReal_t GLlog10E =   0.4342944819032518276511289189166051L;   /* log_10 e */
const GReal_t GLn2 =       0.6931471805599453094172321214581766L;   /* log_e 2 */
const GReal_t GLn10 =      2.3025850929940456840179914546843642L;   /* log_e 10 */
const GReal_t GPi =        3.1415926535897932384626433832795029L;   /* pi */
const GReal_t GSqrtPi =    1.7724538509055160272981674833411452L;   /* sqrt(pi) */
const GReal_t GPi_2 =      1.5707963267948966192313216916397514L;   /* pi/2 */
const GReal_t GPi_4 =      0.7853981633974483096156608458198757L;   /* pi/4 */
const GReal_t G1_Pi =      0.3183098861837906715377675267450287L;   /* 1/pi */
const GReal_t G2_Pi =      0.6366197723675813430755350534900574L;   /* 2/pi */
const GReal_t G2_SqrtPi =  1.1283791670955125738961589031215452L;   /* 2/sqrt(pi) */
const GReal_t GSqrt2 =     1.4142135623730950488016887242096981L;   /* sqrt(2) */
const GReal_t GSqrt1_2 =   0.7071067811865475244008443621048490L;   /* 1/sqrt(2) */
const GReal_t G2Pi =       6.2831853071795864769252867665590058L;   /* 2pi */
const GReal_t G1_2Pi =     0.1591549430918953357688837633725144L;   /* 1/(2pi) */
const GReal_t GSqrt2Pi =   2.5066282746310005024157652848110452L;   /* sqrt(2pi) */
const GReal_t G1_Sqrt2Pi = 0.3989422804014326779399460599343819L;   /* 1.0/sqrt(2pi) */
const GReal_t G4Pi =       12.5663706143591729538505735331180116L;  /* 4pi */
const GReal_t G1_4Pi =     0.0795774715459476678844418816862572L;   /* 1/(4pi) */
const GReal_t GSqrt4Pi =   3.5449077018110320545963349666822903L;   /* sqrt(4pi) */
const GReal_t G1_Sqrt4Pi = 0.2820947917738781434740397257803863L;   /* 1.0/sqrt(4pi) */

const GReal_t GReal_EPSILON = sizeof(GReal_t) == 4? FLT_EPSILON :
                              sizeof(GReal_t) == 8? DBL_EPSILON :
                                                    LDBL_EPSILON;

typedef std::complex<GReal_t> GComplex_t;
const GComplex_t GI(0.0, 1.0);

inline GReal_t square(GReal_t x) {
    return x*x;
}

inline GComplex_t square(GComplex_t x) {
    return x*x;
}

inline GReal_t normsquare(GReal_t x) {
    return x*x;
}


inline GReal_t normsquare(GComplex_t x) {
    // This is the original normsquare function implemented with complex<>, 
    // it has wrong name 'norm'.
    return norm(x);
}

/*template <class T>
inline T diff_O2(T (*f)(GReal_t, void*),
                 GReal_t t, GReal_t dt, void* p) {
    T f0 = (*f)(t - dt, p);
    T f2 = (*f)(t + dt, p);
    return (f2 - f0)/(2*dt);
}

template <class T>
inline T diff_O4(T (*f)(GReal_t, void*),
                 GReal_t t, GReal_t dt, void* p) {
    T f0 = (*f)(t - 2*dt, p);
    T f1 = (*f)(t - dt, p);
    T f3 = (*f)(t + dt, p);
    T f4 = (*f)(t + 2*dt, p);
    return (f0 - 8*(f1 - f3) - f4) / (12*dt);
}*/

}

#endif /* gridripper_math_h */