dd_real.h

#pragma once

#include <cmath>
#include <compare>
#include <limits>
#include <type_traits>

// Double-double arithmetic: represents x = hi + lo where |lo| <= eps * |hi|.
// This extends double precision to roughly 30 decimal digits.
class dd_real {
 public:
  using value_type = double;

  constexpr dd_real() noexcept = default;
  constexpr dd_real(double v) noexcept : hi_(v) {}
  constexpr dd_real(long double v) noexcept : hi_(static_cast<double>(v)) {}

  static constexpr auto from_parts(double hi, double lo) noexcept -> dd_real {
    return {raw_parts_t{}, hi, lo};
  }

  static auto normalized(double hi, double lo) noexcept -> dd_real {
    return quick_two_sum(hi, lo);
  }

  [[nodiscard]] constexpr auto hi() const noexcept -> double { return hi_; }
  [[nodiscard]] constexpr auto lo() const noexcept -> double { return lo_; }

  explicit constexpr operator float() const noexcept { return static_cast<float>(hi_); }
  explicit constexpr operator double() const noexcept { return hi_ + lo_; }

  auto operator-() const noexcept -> dd_real { return from_parts(-hi_, -lo_); }

  friend auto operator+(dd_real a, dd_real b) noexcept -> dd_real {
    auto s1 = two_sum(a.hi_, b.hi_);
    const double e1 = s1.lo_ + a.lo_ + b.lo_;
    return quick_two_sum(s1.hi_, e1);
  }
  friend auto operator-(dd_real a, dd_real b) noexcept -> dd_real { return a + (-b); }
  friend auto operator*(dd_real a, dd_real b) noexcept -> dd_real {
    auto p = two_prod(a.hi_, b.hi_);
    const double e = p.lo_ + a.hi_ * b.lo_ + a.lo_ * b.hi_;
    return quick_two_sum(p.hi_, e);
  }
  friend auto operator/(dd_real a, dd_real b) noexcept -> dd_real {
    const double q1 = a.hi_ / b.hi_;
    const dd_real r = a - b * q1;
    const double q2 = r.hi_ / b.hi_;
    return quick_two_sum(q1, q2);
  }

  friend auto operator+(dd_real a, double b) noexcept -> dd_real { return a + dd_real(b); }
  friend auto operator-(dd_real a, double b) noexcept -> dd_real { return a + (-dd_real(b)); }
  friend auto operator*(dd_real a, double b) noexcept -> dd_real {
    auto p = two_prod(a.hi_, b);
    const double e = p.lo_ + a.lo_ * b;
    return quick_two_sum(p.hi_, e);
  }
  friend auto operator/(dd_real a, double b) noexcept -> dd_real {
    const double q1 = a.hi_ / b;
    auto p = two_prod(q1, b);
    const double r = ((a.hi_ - p.hi_) + (a.lo_ - p.lo_)) / b;
    return quick_two_sum(q1, r);
  }

  friend auto operator+(double a, dd_real b) noexcept -> dd_real { return dd_real(a) + b; }
  friend auto operator-(double a, dd_real b) noexcept -> dd_real { return dd_real(a) - b; }
  friend auto operator*(double a, dd_real b) noexcept -> dd_real { return b * a; }
  friend auto operator/(double a, dd_real b) noexcept -> dd_real { return dd_real(a) / b; }

  auto operator+=(dd_real b) noexcept -> dd_real& { return *this = *this + b; }
  auto operator-=(dd_real b) noexcept -> dd_real& { return *this = *this - b; }
  auto operator*=(dd_real b) noexcept -> dd_real& { return *this = *this * b; }
  auto operator/=(dd_real b) noexcept -> dd_real& { return *this = *this / b; }
  auto operator+=(double b) noexcept -> dd_real& { return *this = *this + b; }
  auto operator-=(double b) noexcept -> dd_real& { return *this = *this - b; }
  auto operator*=(double b) noexcept -> dd_real& { return *this = *this * b; }
  auto operator/=(double b) noexcept -> dd_real& { return *this = *this / b; }

  friend auto operator==(dd_real a, dd_real b) noexcept -> bool {
    if (std::isnan(a.hi_) || std::isnan(a.lo_) || std::isnan(b.hi_) || std::isnan(b.lo_)) {
      return false;
    }
    return a.hi_ == b.hi_ && a.lo_ == b.lo_;
  }

  friend auto operator<=>(dd_real a, dd_real b) noexcept -> std::partial_ordering {
    if (std::isnan(a.hi_) || std::isnan(a.lo_) || std::isnan(b.hi_) || std::isnan(b.lo_)) {
      return std::partial_ordering::unordered;
    }
    if (a.hi_ < b.hi_) return std::partial_ordering::less;
    if (a.hi_ > b.hi_) return std::partial_ordering::greater;
    if (a.lo_ < b.lo_) return std::partial_ordering::less;
    if (a.lo_ > b.lo_) return std::partial_ordering::greater;
    return std::partial_ordering::equivalent;
  }

  friend auto operator==(dd_real a, double b) noexcept -> bool {
    return (a <=> dd_real(b)) == std::partial_ordering::equivalent;
  }
  friend auto operator==(double a, dd_real b) noexcept -> bool { return b == a; }
  friend auto operator<=>(dd_real a, double b) noexcept -> std::partial_ordering { return a <=> dd_real(b); }
  friend auto operator<=>(double a, dd_real b) noexcept -> std::partial_ordering { return dd_real(a) <=> b; }

 private:
  struct raw_parts_t {
    explicit constexpr raw_parts_t() noexcept = default;
  };

  constexpr dd_real(raw_parts_t, double hi, double lo) noexcept : hi_(hi), lo_(lo) {}

  static auto two_sum(double a, double b) noexcept -> dd_real {
    const double s = a + b;
    const double v = s - a;
    return from_parts(s, (a - (s - v)) + (b - v));
  }
  static auto quick_two_sum(double a, double b) noexcept -> dd_real {
    const double s = a + b;
    return from_parts(s, b - (s - a));
  }
  static auto two_prod(double a, double b) noexcept -> dd_real {
    const double p = a * b;
    return from_parts(p, std::fma(a, b, -p));
  }

  double hi_ = 0.0;
  double lo_ = 0.0;
};

// Complex double-double for high-precision reference orbit computation.
struct cmp_dd {
  dd_real re, im;

  cmp_dd() = default;
  cmp_dd(dd_real r, dd_real i) : re(r), im(i) {}

  friend auto operator+(cmp_dd a, cmp_dd b) -> cmp_dd { return {a.re + b.re, a.im + b.im}; }
  friend auto operator*(cmp_dd a, cmp_dd b) -> cmp_dd {
    return {a.re * b.re - a.im * b.im, a.re * b.im + a.im * b.re};
  }
  [[nodiscard]] auto norm() const -> dd_real { return re * re + im * im; }
};

inline auto isfinite(dd_real v) noexcept -> bool {
  return std::isfinite(v.hi()) && std::isfinite(v.lo());
}
inline auto isnan(dd_real v) noexcept -> bool {
  return std::isnan(v.hi()) || std::isnan(v.lo());
}
inline auto isinf(dd_real v) noexcept -> bool {
  return std::isinf(v.hi()) || std::isinf(v.lo());
}

inline auto to_double(dd_real v) noexcept -> double {
  return static_cast<double>(v);
}

namespace std {

template <>
class numeric_limits<dd_real> {
 public:
  static constexpr bool is_specialized = true;

  static constexpr int digits = 106;
  static constexpr int digits10 = 31;
  static constexpr int max_digits10 = 33;

  static constexpr bool is_signed = true;
  static constexpr bool is_integer = false;
  static constexpr bool is_exact = false;
  static constexpr int radix = 2;

  static constexpr int min_exponent = numeric_limits<double>::min_exponent;
  static constexpr int min_exponent10 = numeric_limits<double>::min_exponent10;
  static constexpr int max_exponent = numeric_limits<double>::max_exponent;
  static constexpr int max_exponent10 = numeric_limits<double>::max_exponent10;

  static constexpr bool has_infinity = numeric_limits<double>::has_infinity;
  static constexpr bool has_quiet_NaN = numeric_limits<double>::has_quiet_NaN;
  static constexpr bool has_signaling_NaN = numeric_limits<double>::has_signaling_NaN;
  static constexpr bool is_iec559 = numeric_limits<double>::is_iec559;
  static constexpr bool is_bounded = numeric_limits<double>::is_bounded;
  static constexpr bool is_modulo = numeric_limits<double>::is_modulo;

  static constexpr bool traps = numeric_limits<double>::traps;
  static constexpr bool tinyness_before = numeric_limits<double>::tinyness_before;
  static constexpr float_round_style round_style = numeric_limits<double>::round_style;

  static constexpr auto min() noexcept -> dd_real { return {numeric_limits<double>::min()}; }
  static constexpr auto lowest() noexcept -> dd_real { return {numeric_limits<double>::lowest()}; }
  static constexpr auto max() noexcept -> dd_real { return {numeric_limits<double>::max()}; }
  static constexpr auto epsilon() noexcept -> dd_real {
    return dd_real::from_parts(numeric_limits<double>::epsilon(), 0.0);
  }
  static constexpr auto round_error() noexcept -> dd_real { return dd_real::from_parts(0.5, 0.0); }
  static constexpr auto infinity() noexcept -> dd_real { return {numeric_limits<double>::infinity()}; }
  static constexpr auto quiet_NaN() noexcept -> dd_real { return {numeric_limits<double>::quiet_NaN()}; }
  static constexpr auto signaling_NaN() noexcept -> dd_real {
    return {numeric_limits<double>::signaling_NaN()};
  }
  static constexpr auto denorm_min() noexcept -> dd_real { return {numeric_limits<double>::denorm_min()}; }
};

template <>
struct common_type<dd_real, double> {
  using type = dd_real;
};

template <>
struct common_type<double, dd_real> {
  using type = dd_real;
};

template <>
struct common_type<dd_real, float> {
  using type = dd_real;
};

template <>
struct common_type<float, dd_real> {
  using type = dd_real;
};

}  // namespace std