#ifndef DUAL_HXX
#define DUAL_HXX
#include "vect.hxx"

  }
};
template <typename T, typename U> struct zero<dual<T, U> > {
  constexpr dual<T, U> operator()() const {
    return dual<T, U>(zero<T>()(), zero<U>()());

constexpr Arith::dual<T, U> abs(const Arith::dual<T, U> &x);
template <typename T, typename U>

constexpr bool isnan(const Arith::dual<T, U> &x);
template <typename T, typename U>

template <typename T, typename U>
constexpr Arith::dual<T, U> abs(const Arith::dual<T, U> &x) {
  // return sqrt(pow2(x));
  return {abs(x.val), copysign(T{1}, x.val) * x.eps};
}

constexpr bool isnan(const Arith::dual<T, U> &x) {
  return isnan(x.val) || x.eps.isnan().any();
}
template <typename T, typename U>

////////////////////////////////////////////////////////////////////////////////
namespace Arith {
template <typename T> struct one;
template <> struct one<bool> {
  constexpr bool operator()() const { return 1; }
};
template <> struct one<short> {
  constexpr short operator()() const { return 1; }
};
template <> struct one<int> {
  constexpr int operator()() const { return 1; }
};
template <> struct one<long> {
  constexpr long operator()() const { return 1; }
};
template <> struct one<float> {
  constexpr float operator()() const { return 1; }
};
template <> struct one<double> {
  constexpr double operator()() const { return 1; }
};
template <typename T> struct one<complex<T> > {
  constexpr complex<T> operator()() const {
    return complex<T>(one<T>()(), zero<T>()());
  }
};
template <typename T, typename U> struct one<dual<T, U> > {
  constexpr dual<T, U> operator()() const {
    return dual<T, U>(one<T>()(), zero<U>()());
  }
};
} // namespace Arith
#endif // #ifndef DUAL_HXX

#include <complex>

template <> struct zero<bool> {
  constexpr bool operator()() const { return 0; }
};

};
template <typename T> struct zero<complex<T> > {
  constexpr complex<T> operator()() const {
    return complex<T>(zero<T>()(), zero<T>()());
  }

  return r;
}
template <typename T, int D>
constexpr CCTK_ATTRIBUTE_ALWAYS_INLINE Arith::vect<bool, D>
isnan(const Arith::vect<T, D> &x) {
  Arith::vect<bool, D> r;
  for (int d = 0; d < D; ++d)
    r.elts[d] = isnan(x.elts[d]);