pzreal.h

Go to the documentation of this file.

/******************************************************************************
 *
 * Data definition:     REAL
 *
 * Purpose:             Define what type of data the package is being compiled
 *                      for.  The choices are the following:
 *
 *                          double   real, double precision
 *                          float    real, single precision
 *
 *****************************************************************************/

#ifndef REALH
#define REALH

#include <pz_config.h>
#include <limits>

#ifndef _USE_MATH_DEFINES
#define _USE_MATH_DEFINES
#include <cmath>
#endif // _USE_MATH_DEFINES

#include <iostream>
#include <complex>
#include "fpo_exceptions.h"

#ifdef _AUTODIFF
template <int Num, class T> class TFad;
#endif

template <typename Enumeration>
typename std::underlying_type<Enumeration>::type as_integer(const Enumeration value) {
    return static_cast<typename std::underlying_type<Enumeration>::type>(value);
}


/*structs used for help identifying fundamental types in template parameters.
 For instance,
 
 template <class T,
 typename std::enable_if<(is_arithmetic_pz::value), int>::type* = nullptr>
 void Write(const TPZVec<T> &vec){
 //stuff here
 }
 
 This template would only match with T as char, int, long, float, double, 
 * std::complex<float> etc... (Not composite types).*/

template<class T>
struct is_complex_or_floating_point : std::is_floating_point<T> { };


template<class T>
struct is_complex_or_floating_point<std::complex<T>> : std::integral_constant<bool,
        std::is_integral<T>::value ||
        std::is_floating_point<T>::value> { };

template<class T>
struct is_arithmetic_pz : std::integral_constant<bool,
        std::is_integral<T>::value ||
        is_complex_or_floating_point<T>::value> { };

#ifndef MAX
#define MAX( a, b )   ( (a) > (b) ? (a) : (b) )
#endif

#ifndef MIN
#define MIN( a, b )   ( (a) < (b) ? (a) : (b) )
#endif

#ifdef WIN32
#define  __PRETTY_FUNCTION__ __FILE__
#endif


extern int gPrintLevel;

enum MOptype {ESum,EProd,EDiv,ESqrt,EPow,ECos,ESin,EAcos,EAsin,EAtan,EExp,ELog};

const int gNumOp = 12;

struct TPZCounter {
        uint64_t fCount[gNumOp];
                
        TPZCounter() 
        {
                  clear();
                }

                inline TPZCounter operator-(const TPZCounter &other)
                {
                  TPZCounter result;
                  for(int i=0; i<gNumOp; i++) result.fCount[i] = fCount[i]-other.fCount[i];
                  return result;
                }

                inline TPZCounter& operator+=(const TPZCounter &other)
                {
                  for(int i=0; i<gNumOp; i++) fCount[i] += other.fCount[i];
                  return *this;
                }

                inline TPZCounter& operator-=(const TPZCounter &other)
                {
                  for(int i=0; i<gNumOp; i++) fCount[i] -= other.fCount[i];
                  return *this;
                }

                inline TPZCounter& copy(const TPZCounter &other)
                {
                  for (int i=0; i<gNumOp; i++) fCount[i] = other.fCount[i];
                  return *this;
                }

        inline TPZCounter& clear()
                {
                  for(int i=0; i<gNumOp; i++) fCount[i] = 0 ;
                  return *this;
                }
                
                void Print(std::ostream &out = std::cout) const;
                
};

std::ostream &operator<<(std::ostream &out,const TPZCounter &count);

class TPZFlopCounter;

#ifdef REALfloat
typedef float REAL;
#endif // REALfloat
#ifdef REALdouble
typedef double REAL; //This is the default configuration
#endif // REALdouble
#ifdef REALlongdouble
typedef long double REAL;
#endif // REALlongdouble
#ifdef REALpzfpcounter
typedef TPZFlopCounter REAL;
#endif // REALpzfpcounter

#ifdef STATEfloat
typedef float STATE;
#endif // STATEfloat
#ifdef STATEdouble
typedef double STATE; //This is the default configuration
#endif // STATEdouble
#ifdef STATElongdouble
typedef long double STATE;
#endif // STATElongdouble
#ifdef STATEcomplexf
typedef std::complex<float> STATE;
#endif // STATEcomplexf
#ifdef STATEcomplexd
typedef std::complex<double> STATE;
#endif //STATEcomplexd
#ifdef STATEcomplexld
typedef std::complex<long double> STATE;
#endif //STATEcomplexld

#ifdef VC
#include <io.h>
#ifndef NOMINMAX
#define NOMINMAX // Preventing the redefinition of min and max as macros
#endif // NOMINMAX
#include <Windows.h>
// sqrt function adapted to int numbers. required for VC
inline double
sqrt(int __x)
{
  return sqrt ((double)__x);
}
// fabs function adapted to int numbers. required for VC
inline int
fabs(int __x)
{
  return (int)fabs((double) __x);
}
inline double
sin(int __x)
{
  return sin((double) __x);
}
inline double
cos(int __x)
{
  return cos((double) __x);
}
inline double
atan(int __x)
{
  return atan((double) __x);
}
#endif //VC

// fabs function adapted to complex numbers.
inline float
fabs(std::complex <float> __x)
{
  std::complex <float> ret = abs(__x);
  return std::abs(ret.real());
}
inline double
fabs(std::complex <double> __x)
{
  std::complex <double> ret = abs(__x);
  return std::abs(ret.real());
}
inline long double
fabs(std::complex <long double> __x)
{
  std::complex <long double> ret = abs(__x);
  return std::abs(ret.real());
}


class TPZFlopCounter {
public:
#ifdef REALpzfpcounter
                double fVal;
#else
                REAL fVal;
#endif //REALpzfpcounter

                static TPZCounter gCount;

                inline TPZFlopCounter()
                {
                }
                inline TPZFlopCounter(const double &val)
                {
                                fVal = (REAL)val;
                }
    
    inline REAL val() const
    {
        return  fVal;
    }
    
    operator REAL() const{
        return fVal;
    }
    
    bool operator<=(const REAL &val) const
    {
        return fVal <= val;
    }
    
    bool operator<(const REAL &val) const
    {
        return fVal < val;
    }
    bool operator>(const REAL &val) const
    {
        return fVal > val;
    }
    bool operator>=(const REAL &val) const
    {
        return fVal >= val;
    }
    bool operator==(const REAL &val) const
    {
        return fVal == val;
    }
                inline TPZFlopCounter operator*(const TPZFlopCounter &oth) const
                {
                                TPZFlopCounter result;
                                result.fVal = fVal*oth.fVal;
                                gCount.fCount[EProd]++;
                                return result;
                }
                inline TPZFlopCounter operator/(const TPZFlopCounter &oth) const
                {
                                TPZFlopCounter result;
                                result.fVal = fVal/oth.fVal;
                                gCount.fCount[EDiv]++;
                                return result;
                }
                
                inline TPZFlopCounter operator*(const REAL &oth) const
                {
                                TPZFlopCounter result;
                                result.fVal = fVal*oth;
                                gCount.fCount[EProd]++;
                                return result;
                }
                inline TPZFlopCounter operator/(const double &oth) const
                {
                                TPZFlopCounter result;
                                result.fVal = fVal/((REAL)oth);
                                gCount.fCount[EDiv]++;
                                return result;
                }
                
                inline TPZFlopCounter operator+(const TPZFlopCounter &oth) const
                {
                                TPZFlopCounter result;
                                result.fVal = fVal+oth.fVal;
                                gCount.fCount[ESum]++;
                                return result;
                }
                
                inline TPZFlopCounter operator-(const TPZFlopCounter &oth) const
                {
                                TPZFlopCounter result;
                                result.fVal = fVal-oth.fVal;
                                gCount.fCount[ESum]++;
                                return result;
                }
                inline TPZFlopCounter &operator+=(const TPZFlopCounter &oth) 
                {
                                fVal += oth.fVal;
                                gCount.fCount[ESum]++;
                                return *this;
                }
                inline TPZFlopCounter &operator-=(const TPZFlopCounter &oth) 
                {
                                fVal -= oth.fVal;
                                gCount.fCount[ESum]++;
                                return *this;
                }
                inline TPZFlopCounter &operator*=(const TPZFlopCounter &oth) 
                {
                                fVal *= oth.fVal;
                                gCount.fCount[EProd]++;
                                return *this;
                }
                inline TPZFlopCounter &operator/=(const TPZFlopCounter &oth) 
                {
                                fVal /= oth.fVal;
                                gCount.fCount[EDiv]++;
                                return *this;
                }
                inline TPZFlopCounter operator-() const
                {
                                TPZFlopCounter result;
                                result.fVal = -fVal;
                                gCount.fCount[ESum]++;
                                return result;
                }
                inline TPZFlopCounter operator+() const
                {
                                return *this;
                }
                inline bool operator<(const TPZFlopCounter &sec) const
                {
                                return (fVal < sec.fVal);
                }
                inline bool operator>(const TPZFlopCounter &sec) const
                {
                                return (fVal > sec.fVal);
                }
                inline bool operator<=(const TPZFlopCounter &sec) const
                {
                                return (fVal <= sec.fVal);
                }
                inline bool operator>=(const TPZFlopCounter &sec) const
                {
                                return (fVal >= sec.fVal);
                }
                inline bool operator==(const TPZFlopCounter &sec) const
                {
                                return (fVal == sec.fVal);
                }
                inline bool operator!=(const TPZFlopCounter &sec) const
                {
                                return (fVal != sec.fVal);
                }
                
                friend TPZFlopCounter sqrt(const TPZFlopCounter &orig);
                friend TPZFlopCounter pow(const TPZFlopCounter &orig,const TPZFlopCounter &xp);
                friend TPZFlopCounter fabsFlop(const TPZFlopCounter &orig);
                friend REAL fabs(const TPZFlopCounter &orig);
                friend TPZFlopCounter acos(const TPZFlopCounter &orig);
                friend TPZFlopCounter cos(const TPZFlopCounter &orig);
                friend TPZFlopCounter asin(const TPZFlopCounter &orig);
                friend TPZFlopCounter sin(const TPZFlopCounter &orig);
                friend TPZFlopCounter atan(const TPZFlopCounter &orig);
                friend TPZFlopCounter atan2(const TPZFlopCounter &val1,const TPZFlopCounter &val2);
                friend TPZFlopCounter exp(const TPZFlopCounter &val);
                friend TPZFlopCounter log(const TPZFlopCounter &val);
                friend TPZFlopCounter log10(const TPZFlopCounter &val);
};

inline TPZFlopCounter sqrt(const TPZFlopCounter &orig)
{
                TPZFlopCounter result;
                result.fVal = ::sqrt(orig.fVal);
                TPZFlopCounter::gCount.fCount[ESqrt]++;
                return result;
}

inline TPZFlopCounter fabsFlop(const TPZFlopCounter &orig)
{
                TPZFlopCounter result;
                result.fVal = std::abs(orig.fVal);
                return result;
}
inline REAL fabs(const TPZFlopCounter &orig)
{
    return std::abs(orig.fVal);
}

inline TPZFlopCounter pow(const TPZFlopCounter &orig,const TPZFlopCounter &xp)
{
                TPZFlopCounter result;
                result.fVal = ::pow(orig.fVal,xp.fVal);
                TPZFlopCounter::gCount.fCount[EPow]++;
                return result;
}

inline TPZFlopCounter acos(const TPZFlopCounter &orig)
{
                TPZFlopCounter result;
                result.fVal = ::acos(orig.fVal);
                TPZFlopCounter::gCount.fCount[EAcos]++;
                return result;
}

inline TPZFlopCounter asin(const TPZFlopCounter &orig)
{
                TPZFlopCounter result;
                result.fVal = ::asin(orig.fVal);
                TPZFlopCounter::gCount.fCount[EAsin]++;
                return result;
}

inline TPZFlopCounter cos(const TPZFlopCounter &orig)
{
                TPZFlopCounter result;
                result.fVal = ::cos(orig.fVal);
                TPZFlopCounter::gCount.fCount[ECos]++;
                return result;
}

inline TPZFlopCounter sin(const TPZFlopCounter &orig)
{
                TPZFlopCounter result;
                result.fVal = ::sin(orig.fVal);
                TPZFlopCounter::gCount.fCount[ESin]++;
                return result;
}

inline TPZFlopCounter atan(const TPZFlopCounter &orig)
{
                TPZFlopCounter result;
                result.fVal = ::atan(orig.fVal);
                TPZFlopCounter::gCount.fCount[EAtan]++;
                return result;
}

inline TPZFlopCounter atan2(const TPZFlopCounter &val1,const TPZFlopCounter &val2)
{
                TPZFlopCounter result;
                result.fVal = ::atan2(val1.fVal,val2.fVal);
                TPZFlopCounter::gCount.fCount[EAtan]++;
                return result;
}

inline TPZFlopCounter exp(const TPZFlopCounter &orig)
{
                TPZFlopCounter result;
                result.fVal = ::exp(orig.fVal);
                TPZFlopCounter::gCount.fCount[EExp]++;
                return result;
}

inline TPZFlopCounter log(const TPZFlopCounter &orig)
{
                TPZFlopCounter result;
                result.fVal = ::log(orig.fVal);
                TPZFlopCounter::gCount.fCount[ELog]++;
                return result;
}

inline TPZFlopCounter log10(const TPZFlopCounter &orig)
{
                TPZFlopCounter result;
                result.fVal = ::log10(orig.fVal);
                TPZFlopCounter::gCount.fCount[ELog]++;
                return result;
}
inline TPZFlopCounter operator*(double val1, const TPZFlopCounter &val2)
{
                TPZFlopCounter result;
                result = TPZFlopCounter(val1)*val2;
                return result;
                
}
inline TPZFlopCounter operator/(double val1, const TPZFlopCounter &val2)
{
                TPZFlopCounter result;
                result = TPZFlopCounter(val1)/val2;
                return result;
                
}
inline TPZFlopCounter operator+(double val1, const TPZFlopCounter &val2)
{
                TPZFlopCounter result;
                result = TPZFlopCounter(val1)+val2;
                return result;
                
}
inline TPZFlopCounter operator+(const TPZFlopCounter &val2, double val1 )
{
                TPZFlopCounter result;
                result = TPZFlopCounter(val1)+val2;
                return result;
                
}
inline TPZFlopCounter operator-(double val1, const TPZFlopCounter &val2)
{
                TPZFlopCounter result;
                result = TPZFlopCounter(val1)-val2;
                return result;
                
}

inline std::ostream &operator<<(std::ostream &out, const TPZFlopCounter &val)
{
                return out << val.fVal;
}
inline std::istream &operator>>(std::istream &out, /*const*/ TPZFlopCounter &val)
{
                return out >> val.fVal;
}



inline REAL ZeroTolerance() {
    typedef std::numeric_limits< REAL > dbl;
                return (REAL)pow(10,(-1 * (dbl::max_digits10- 5)));
}

template<typename T,
        typename std::enable_if< is_arithmetic_pz<T>::value, int>::type* = nullptr>
inline void ZeroTolerance(T &tol) {
    typedef std::numeric_limits< T > dbl;
    tol = (T)pow(10,(-1 * (dbl::max_digits10- 5)));
}

template<typename T,
        typename std::enable_if<std::is_same<T,TPZFlopCounter>::value>::type* = nullptr>
inline void ZeroTolerance(T &tol) {
    typedef std::numeric_limits< REAL > dbl;
    tol = (T)pow(10,(-1 * (dbl::max_digits10- 5)));
}

#ifdef _AUTODIFF
template<int Num, typename T> 
inline void ZeroTolerance(TFad<Num,T> &Tol) {
    ZeroTolerance(Tol.val());
}
#endif

#ifdef _AUTODIFF

template<class T>
inline bool IsZero( T a ) {
                return ( std::abs( a.val() ) < ZeroTolerance() );
}
#endif

//template<>
inline bool IsZero( long double a ) {
#ifdef WIN32
                return ( fabs( a ) < 1.e-12 );
#else
                return ( std::fabs( a ) < 1.e-16 );
#endif
}
//template<>
inline bool IsZero( double a ) {
#ifdef WIN32
                return ( fabs( a ) < 1.e-10 );
#else
                return ( fabs( a ) < 1.e-12 );
#endif
}
//template<>
inline bool IsZero( float a ) {
#ifdef WIN32
                return ( fabs( a ) < 1.e-6 );
#else
                return ( fabs( a ) < 1.e-6 );
#endif
}
//template<>
inline bool IsZero( std::complex<long double> a ) {
#ifdef WIN32
                return ( fabs( a ) < 1.e-12 );
#else
                return ( fabs( a ) < 1.e-16 );
#endif
}
//template<>
inline bool IsZero( std::complex<double> a ) {
#ifdef WIN32
                return ( fabs( a ) < 1.e-9 );
#else
                return ( fabs( a ) < 1.e-12 );
#endif
}
//template<>
inline bool IsZero( std::complex<float> a ) {
#ifdef WIN32
                return ( fabs( a ) < 1.e-5 );
#else
                return ( fabs( a ) < 1.e-7 );
#endif
}
//template<>
inline bool IsZero( int a ) {
                return ( a==0 );
}
inline bool IsZero( int64_t a ) {
                return ( a==0L );
}
template <class T>
inline const T& Max( const T &a, const T &b ) {
                return a > b ? a : b;
}
template <class T>
inline const T& Min( const T & a, const T &b ) {
                return a < b ? a : b;
}

// In the math library (cmath.h) don't exist some overloading for some functions

// SPECIAL FUNCTIONS NON STANDARD IN WINDOWS SYSTEM
#if (!defined(__cplusplus) || __cplusplus < 201103L) && (!defined(_MSC_VER) || _MSC_VER < 1900)// If we aren't using C++11.

REAL erf(REAL arg);

#endif // not C++11

#if defined(_MSC_VER) && _MSC_VER < 1900 // Microsoft Visual Studio < 2015

#include <cfloat>
#define isnan(x) _isnan(x)

#endif

#endif