numbers.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2006 - 2017 by the deal.II authors
//
// This file is part of the deal.II library.
//
// The deal.II library is free software; you can use it, redistribute
// it, and/or modify it under the terms of the GNU Lesser General
// Public License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
// The full text of the license can be found in the file LICENSE at
// the top level of the deal.II distribution.
//
// ---------------------------------------------------------------------

#ifndef dealii__numbers_h
#define dealii__numbers_h


#include <deal.II/base/config.h>
#include <deal.II/base/types.h>

#include <cmath>
#include <cstdlib>
#include <complex>

DEAL_II_NAMESPACE_OPEN

// forward declarations to support abs or sqrt operations on VectorizedArray
template <typename Number> class VectorizedArray;
template <typename T> struct EnableIfScalar;

DEAL_II_NAMESPACE_CLOSE

namespace std
{
  template <typename Number> DEAL_II_ALWAYS_INLINE ::VectorizedArray<Number>
  sqrt(const ::VectorizedArray<Number> &);
  template <typename Number> DEAL_II_ALWAYS_INLINE ::VectorizedArray<Number>
  abs(const ::VectorizedArray<Number> &);
  template <typename Number> DEAL_II_ALWAYS_INLINE ::VectorizedArray<Number>
  max(const ::VectorizedArray<Number> &, const ::VectorizedArray<Number> &);
  template <typename Number> DEAL_II_ALWAYS_INLINE ::VectorizedArray<Number>
  min (const ::VectorizedArray<Number> &, const ::VectorizedArray<Number> &);
}

DEAL_II_NAMESPACE_OPEN

namespace numbers
{
  static const double  E       = 2.7182818284590452354;

  static const double  LOG2E   = 1.4426950408889634074;

  static const double  LOG10E  = 0.43429448190325182765;

  static const double  LN2     = 0.69314718055994530942;

  static const double  LN10    = 2.30258509299404568402;

  static const double  PI      = 3.14159265358979323846;

  static const double  PI_2    = 1.57079632679489661923;

  static const double  PI_4    = 0.78539816339744830962;

  static const double  SQRT2   = 1.41421356237309504880;

  static const double  SQRT1_2 = 0.70710678118654752440;

  bool is_nan (const double x);

  bool is_finite (const double x);

  bool is_finite (const std::complex<double> &x);

  bool is_finite (const std::complex<float> &x);

  bool is_finite (const std::complex<long double> &x);

  template <typename number>
  struct NumberTraits
  {
    static const bool is_complex = false;

    typedef number real_type;

    static
    const number &conjugate (const number &x);

    static
    real_type abs_square (const number &x);

    static
    real_type abs (const number &x);
  };


  template <typename number>
  struct NumberTraits<std::complex<number> >
  {
    static const bool is_complex = true;

    typedef number real_type;

    static
    std::complex<number> conjugate (const std::complex<number> &x);

    static
    real_type abs_square (const std::complex<number> &x);


    static
    real_type abs (const std::complex<number> &x);
  };

  // --------------- inline and template functions ---------------- //

  inline bool is_nan (const double x)
  {
#ifdef DEAL_II_HAVE_STD_ISNAN
    return std::isnan(x);
#elif defined(DEAL_II_HAVE_ISNAN)
    return isnan(x);
#elif defined(DEAL_II_HAVE_UNDERSCORE_ISNAN)
    return _isnan(x);
#else
    return false;
#endif
  }

  inline bool is_finite (const double x)
  {
#ifdef DEAL_II_HAVE_ISFINITE
    return !is_nan(x) && std::isfinite (x);
#else
    // Check against infinities. Note
    // that if x is a NaN, then both
    // comparisons will be false
    return ((x >= -std::numeric_limits<double>::max())
            &&
            (x <= std::numeric_limits<double>::max()));
#endif
  }



  inline bool is_finite (const std::complex<double> &x)
  {
    // Check complex numbers for infinity
    // by testing real and imaginary part
    return ( is_finite (x.real())
             &&
             is_finite (x.imag()) );
  }



  inline bool is_finite (const std::complex<float> &x)
  {
    // Check complex numbers for infinity
    // by testing real and imaginary part
    return ( is_finite (x.real())
             &&
             is_finite (x.imag()) );
  }



  inline bool is_finite (const std::complex<long double> &x)
  {
    // Same for std::complex<long double>
    return ( is_finite (x.real())
             &&
             is_finite (x.imag()) );
  }


  template <typename number>
  const number &
  NumberTraits<number>::conjugate (const number &x)
  {
    return x;
  }



  template <typename number>
  typename NumberTraits<number>::real_type
  NumberTraits<number>::abs_square (const number &x)
  {
    return x * x;
  }



  template <typename number>
  typename NumberTraits<number>::real_type
  NumberTraits<number>::abs (const number &x)
  {
    return std::abs(x);
  }



  template <typename number>
  std::complex<number>
  NumberTraits<std::complex<number> >::conjugate (const std::complex<number> &x)
  {
    return std::conj(x);
  }



  template <typename number>
  typename NumberTraits<std::complex<number> >::real_type
  NumberTraits<std::complex<number> >::abs (const std::complex<number> &x)
  {
    return std::abs(x);
  }



  template <typename number>
  typename NumberTraits<std::complex<number> >::real_type
  NumberTraits<std::complex<number> >::abs_square (const std::complex<number> &x)
  {
    return std::norm (x);
  }

}

namespace internal
{
  template <typename T>
  struct NumberType
  {
    static T value (const T &t)
    {
      return t;
    }
  };

  template <typename T>
  struct NumberType<std::complex<T> >
  {
    static std::complex<T> value (const T &t)
    {
      return std::complex<T>(t);
    }
  };
}

DEAL_II_NAMESPACE_CLOSE

#endif