numbers.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2006 - 2019 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.md at
// the top level directory of deal.II.
//
// ---------------------------------------------------------------------

#ifndef dealii_numbers_h
#define dealii_numbers_h


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

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

#ifdef DEAL_II_COMPILER_CUDA_AWARE
#  include <cuComplex.h>
#endif

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

#ifdef DEAL_II_COMPILER_CUDA_AWARE
#  define DEAL_II_CUDA_HOST_DEV __host__ __device__
#else
#  define DEAL_II_CUDA_HOST_DEV
#endif

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

// Declare / Import auto-differentiable math functions in(to) standard
// namespace before numbers::NumberTraits is defined
#ifdef DEAL_II_WITH_ADOLC
#  include <deal.II/differentiation/ad/adolc_math.h>

#  include <adolc/adouble.h> // Taped double
#endif
// Ideally we'd like to #include <deal.II/differentiation/ad/sacado_math.h>
// but header indirectly references numbers.h. We therefore simply
// import the whole Sacado header at this point to get the math
// functions imported into the standard namespace.
#ifdef DEAL_II_TRILINOS_WITH_SACADO
#  include <Sacado.hpp>
#endif

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> &);
  template <typename Number>
  ::VectorizedArray<Number>
  pow(const ::VectorizedArray<Number> &, const Number p);
  template <typename Number>
  ::VectorizedArray<Number>
  sin(const ::VectorizedArray<Number> &);
  template <typename Number>
  ::VectorizedArray<Number>
  cos(const ::VectorizedArray<Number> &);
  template <typename Number>
  ::VectorizedArray<Number>
  tan(const ::VectorizedArray<Number> &);
  template <typename Number>
  ::VectorizedArray<Number>
  exp(const ::VectorizedArray<Number> &);
  template <typename Number>
  ::VectorizedArray<Number>
  log(const ::VectorizedArray<Number> &);
} // namespace std

DEAL_II_NAMESPACE_OPEN

namespace numbers
{
  static constexpr double E = 2.7182818284590452354;

  static constexpr double LOG2E = 1.4426950408889634074;

  static constexpr double LOG10E = 0.43429448190325182765;

  static constexpr double LN2 = 0.69314718055994530942;

  static constexpr double LN10 = 2.30258509299404568402;

  static constexpr double PI = 3.14159265358979323846;

  static constexpr double PI_2 = 1.57079632679489661923;

  static constexpr double PI_4 = 0.78539816339744830962;

  static constexpr double SQRT2 = 1.41421356237309504880;

  static constexpr double SQRT1_2 = 0.70710678118654752440;

  template <typename Number, typename = void>
  struct is_cuda_compatible : std::true_type
  {};

  template <typename Number>
  struct is_cuda_compatible<std::complex<Number>, void> : std::false_type
  {};

  DEAL_II_DEPRECATED
  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 Number1, typename Number2>
  bool
  values_are_equal(const Number1 &value_1, const Number2 &value_2);

  template <typename Number1, typename Number2>
  bool
  values_are_not_equal(const Number1 &value_1, const Number2 &value_2);

  template <typename Number>
  bool
  value_is_zero(const Number &value);

  template <typename Number1, typename Number2>
  bool
  value_is_less_than(const Number1 &value_1, const Number2 &value_2);

  template <typename Number1, typename Number2>
  bool
  value_is_less_than_or_equal_to(const Number1 &value_1,
                                 const Number2 &value_2);



  template <typename Number1, typename Number2>
  bool
  value_is_greater_than(const Number1 &value_1, const Number2 &value_2);

  template <typename Number1, typename Number2>
  bool
  value_is_greater_than_or_equal_to(const Number1 &value_1,
                                    const Number2 &value_2);

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

    using real_type = number;

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

    template <typename Dummy = number>
    static constexpr DEAL_II_CUDA_HOST_DEV
      typename std::enable_if<std::is_same<Dummy, number>::value &&
                                is_cuda_compatible<Dummy>::value,
                              real_type>::type
      abs_square(const number &x);

    template <typename Dummy = number>
    static constexpr
      typename std::enable_if<std::is_same<Dummy, number>::value &&
                                !is_cuda_compatible<Dummy>::value,
                              real_type>::type
      abs_square(const number &x);

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


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

    using real_type = number;

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

    static constexpr 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)
  {
    return std::isnan(x);
  }



  inline bool
  is_finite(const double x)
  {
    return std::isfinite(x);
  }



  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>
  constexpr DEAL_II_CUDA_HOST_DEV const number &
                                        NumberTraits<number>::conjugate(const number &x)
  {
    return x;
  }



  template <typename number>
  template <typename Dummy>
  constexpr DEAL_II_CUDA_HOST_DEV
    typename std::enable_if<std::is_same<Dummy, number>::value &&
                              is_cuda_compatible<Dummy>::value,
                            typename NumberTraits<number>::real_type>::type
    NumberTraits<number>::abs_square(const number &x)
  {
    return x * x;
  }



  template <typename number>
  template <typename Dummy>
  constexpr
    typename std::enable_if<std::is_same<Dummy, number>::value &&
                              !is_cuda_compatible<Dummy>::value,
                            typename NumberTraits<number>::real_type>::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>
  constexpr 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>
  constexpr typename NumberTraits<std::complex<number>>::real_type
  NumberTraits<std::complex<number>>::abs_square(const std::complex<number> &x)
  {
    return std::norm(x);
  }

} // namespace numbers


// Forward declarations
namespace Differentiation
{
  namespace AD
  {
    namespace internal
    {
      // Defined in differentiation/ad/ad_number_traits.h
      template <typename T>
      struct NumberType;
    } // namespace internal

    // Defined in differentiation/ad/ad_number_traits.h
    template <typename NumberType>
    struct is_ad_number;
  } // namespace AD
} // namespace Differentiation


namespace internal
{
  template <typename From, typename To>
  struct is_explicitly_convertible
  {
    // Source: https://stackoverflow.com/a/16944130
  private:
    template <typename T>
    static void f(T);

    template <typename F, typename T>
    static constexpr auto
    test(int) -> decltype(f(static_cast<T>(std::declval<F>())), true)
    {
      return true;
    }

    template <typename F, typename T>
    static constexpr auto
    test(...) -> bool
    {
      return false;
    }

  public:
    static bool const value = test<From, To>(0);
  };

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

    // Below are generic functions that allows an overload for any
    // type U that is transformable to type T. This is particularly
    // useful when needing to cast exotic number types
    // (e.g. auto-differentiable or symbolic numbers) to a floating
    // point one, such as might  happen when converting between tensor
    // types.

    // Type T is constructible from F.
    template <typename F>
    static constexpr DEAL_II_CUDA_HOST_DEV T
                                           value(const F &f,
                                                 typename std::enable_if<
            !std::is_same<typename std::decay<T>::type,
                          typename std::decay<F>::type>::value &&
            std::is_constructible<T, F>::value>::type * = nullptr)
    {
      return T(f);
    }

    // Type T is explicitly convertible (but not constructible) from F.
    template <typename F>
    static constexpr T
    value(const F &f,
          typename std::enable_if<
            !std::is_same<typename std::decay<T>::type,
                          typename std::decay<F>::type>::value &&
            !std::is_constructible<T, F>::value &&
            is_explicitly_convertible<const F, T>::value>::type * = nullptr)
    {
      return static_cast<T>(f);
    }

    // Sacado doesn't provide any conversion operators, so we have
    // to extract the value and perform further conversions from there.
    // To be safe, we extend this to other possible AD numbers that
    // might fall into the same category.
    template <typename F>
    static T
    value(const F &f,
          typename std::enable_if<
            !std::is_same<typename std::decay<T>::type,
                          typename std::decay<F>::type>::value &&
            !std::is_constructible<T, F>::value &&
            !is_explicitly_convertible<const F, T>::value &&
            Differentiation::AD::is_ad_number<F>::value>::type * = nullptr)
    {
      return Differentiation::AD::internal::NumberType<T>::value(f);
    }
  };

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

    static constexpr std::complex<T>
    value(const T &t)
    {
      return std::complex<T>(t);
    }

    // Facilitate cast from complex<double> to complex<float>
    template <typename U>
    static constexpr std::complex<T>
    value(const std::complex<U> &t)
    {
      return std::complex<T>(NumberType<T>::value(t.real()),
                             NumberType<T>::value(t.imag()));
    }
  };

#ifdef DEAL_II_COMPILER_CUDA_AWARE
  template <>
  struct NumberType<cuComplex>
  {
    static cuComplex
    value(const float t)
    {
      return make_cuComplex(t, 0.f);
    }
  };

  template <>
  struct NumberType<cuDoubleComplex>
  {
    static cuDoubleComplex
    value(const double t)
    {
      return make_cuDoubleComplex(t, 0.);
    }
  };
#endif
} // namespace internal

namespace numbers
{
#ifdef DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING

  // Defined in differentiation/ad/adolc_number_types.cc
  bool
  values_are_equal(const adouble &value_1, const adouble &value_2);


  template <typename Number>
  bool
  values_are_equal(const adouble &value_1, const Number &value_2)
  {
    // Use the specialized definition for two ADOL-C taped types
    return values_are_equal(value_1,
                            internal::NumberType<adouble>::value(value_2));
  }


  template <typename Number>
  bool
  values_are_equal(const Number &value_1, const adouble &value_2)
  {
    // Use the above definition
    return values_are_equal(value_2, value_1);
  }

  // Defined in differentiation/ad/adolc_number_types.cc
  bool
  value_is_less_than(const adouble &value_1, const adouble &value_2);


  template <typename Number>
  bool
  value_is_less_than(const adouble &value_1, const Number &value_2)
  {
    // Use the specialized definition for two ADOL-C taped types
    return value_is_less_than(value_1,
                              internal::NumberType<adouble>::value(value_2));
  }


  template <typename Number>
  bool
  value_is_less_than(const Number &value_1, const adouble &value_2)
  {
    // Use the specialized definition for two ADOL-C taped types
    return value_is_less_than(internal::NumberType<adouble>::value(value_1),
                              value_2);
  }

#endif


  template <typename Number1, typename Number2>
  inline bool
  values_are_equal(const Number1 &value_1, const Number2 &value_2)
  {
    return (value_1 == internal::NumberType<Number1>::value(value_2));
  }


  template <typename Number1, typename Number2>
  inline bool
  values_are_not_equal(const Number1 &value_1, const Number2 &value_2)
  {
    return !(values_are_equal(value_1, value_2));
  }


  template <typename Number>
  inline bool
  value_is_zero(const Number &value)
  {
    return values_are_equal(value, 0.0);
  }


  template <typename Number1, typename Number2>
  inline bool
  value_is_less_than(const Number1 &value_1, const Number2 &value_2)
  {
    return (value_1 < internal::NumberType<Number1>::value(value_2));
  }


  template <typename Number1, typename Number2>
  inline bool
  value_is_less_than_or_equal_to(const Number1 &value_1, const Number2 &value_2)
  {
    return (value_is_less_than(value_1, value_2) ||
            values_are_equal(value_1, value_2));
  }


  template <typename Number1, typename Number2>
  bool
  value_is_greater_than(const Number1 &value_1, const Number2 &value_2)
  {
    return !(value_is_less_than_or_equal_to(value_1, value_2));
  }


  template <typename Number1, typename Number2>
  inline bool
  value_is_greater_than_or_equal_to(const Number1 &value_1,
                                    const Number2 &value_2)
  {
    return !(value_is_less_than(value_1, value_2));
  }
} // namespace numbers

DEAL_II_NAMESPACE_CLOSE

#endif