Eigen/src/Core/arch/Default/GenericPacketMathDoubleWord.h

// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2018-2025 Rasmus Munk Larsen <rmlarsen@gmail.com>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
// SPDX-License-Identifier: MPL-2.0

#ifndef EIGEN_ARCH_GENERIC_PACKET_MATH_DOUBLE_WORD_H
#define EIGEN_ARCH_GENERIC_PACKET_MATH_DOUBLE_WORD_H

// IWYU pragma: private
#include "../../InternalHeaderCheck.h"

namespace Eigen {
namespace internal {

// This function splits x into the nearest integer n and fractional part r,
// such that x = n + r holds exactly.
template <typename Packet>
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void absolute_split(const Packet& x, Packet& n, Packet& r) {
  n = pround(x);
  r = psub(x, n);
}

// This function computes the sum {s_hi, s_lo}, such that x + y = s_hi + s_lo
// holds exactly, and s_hi = fl(x+y), if |x| >= |y|.
template <typename Packet>
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void fast_twosum(const Packet& x, const Packet& y, Packet& s_hi, Packet& s_lo) {
  s_hi = padd(x, y);
  const Packet t = psub(s_hi, x);
  s_lo = psub(y, t);
}

#ifdef EIGEN_VECTORIZE_FMA
// Given x, y, and xy = fl(x*y), return the residual such that x*y = xy + residual exactly.
template <typename Packet, std::enable_if_t<!is_scalar<Packet>::value, int> = 0>
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet twoprod_low(const Packet& x, const Packet& y, const Packet& xy) {
  return pmsub(x, y, xy);
}

template <typename Scalar, std::enable_if_t<is_scalar<Scalar>::value, int> = 0>
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Scalar twoprod_low(const Scalar& x, const Scalar& y, const Scalar& xy) {
  // Error-free products require FMA even when EIGEN_SCALAR_MADD_USE_FMA disables fusion in scalar madd.
  return numext::fma(x, y, Scalar(-xy));
}

// This function implements the extended precision product of
// a pair of floating point numbers. Given {x, y}, it computes the pair
// {p_hi, p_lo} such that x * y = p_hi + p_lo holds exactly and
// p_hi = fl(x * y).
template <typename Packet>
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void twoprod(const Packet& x, const Packet& y, Packet& p_hi, Packet& p_lo) {
  p_hi = pmul(x, y);
  p_lo = twoprod_low(x, y, p_hi);
}

#else

// Dekker's product requires separately rounded operations. GCC's C++ default -ffp-contract=fast may fuse
// fl(x*y) into its consumers, so use fma for scalar types that <cmath> advertises as fast.
// This does not prevent contraction in the remaining scalar or packet fallback paths.
template <typename Scalar>
struct has_fast_fma : std::false_type {};
#ifdef FP_FAST_FMAF
template <>
struct has_fast_fma<float> : std::true_type {};
#endif
#ifdef FP_FAST_FMA
template <>
struct has_fast_fma<double> : std::true_type {};
#endif
#ifdef FP_FAST_FMAL
template <>
struct has_fast_fma<long double> : std::true_type {};
#endif

template <typename Scalar, std::enable_if_t<has_fast_fma<Scalar>::value, int> = 0>
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Scalar twoprod_low(const Scalar& x, const Scalar& y, const Scalar& xy) {
  return numext::fma(x, y, Scalar(-xy));
}

template <typename Scalar, std::enable_if_t<has_fast_fma<Scalar>::value, int> = 0>
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void twoprod(const Scalar& x, const Scalar& y, Scalar& p_hi, Scalar& p_lo) {
  p_hi = x * y;
  p_lo = twoprod_low(x, y, p_hi);
}

// This function implements the Veltkamp splitting. Given a floating point
// number x it returns the pair {x_hi, x_lo} such that x_hi + x_lo = x holds
// exactly and that half of the significand of x fits in x_hi.
// This is Algorithm 3 from Jean-Michel Muller, "Elementary Functions",
// 3rd edition, Birkh\"auser, 2016.
template <typename Packet>
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void veltkamp_splitting(const Packet& x, Packet& x_hi, Packet& x_lo) {
  using Scalar = typename unpacket_traits<Packet>::type;
  constexpr int shift = (NumTraits<Scalar>::digits() + 1) / 2;
  const Scalar shift_scale = Scalar(uint64_t(1) << shift);  // Scalar constructor not necessarily constexpr.
  const Packet gamma = pmul(pset1<Packet>(shift_scale + Scalar(1)), x);
  Packet rho = psub(x, gamma);
  x_hi = padd(rho, gamma);
  x_lo = psub(x, x_hi);
}

// This function implements Dekker's algorithm for products x * y.
// Given floating point numbers {x, y} computes the pair
// {p_hi, p_lo} such that x * y = p_hi + p_lo holds exactly and
// p_hi = fl(x * y).
template <typename Packet, std::enable_if_t<!has_fast_fma<Packet>::value, int> = 0>
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void twoprod(const Packet& x, const Packet& y, Packet& p_hi, Packet& p_lo) {
  Packet x_hi, x_lo, y_hi, y_lo;
  veltkamp_splitting(x, x_hi, x_lo);
  veltkamp_splitting(y, y_hi, y_lo);

  p_hi = pmul(x, y);
  p_lo = pmadd(x_hi, y_hi, pnegate(p_hi));
  p_lo = pmadd(x_hi, y_lo, p_lo);
  p_lo = pmadd(x_lo, y_hi, p_lo);
  p_lo = pmadd(x_lo, y_lo, p_lo);
}

// A version of twoprod that takes x, y, and fl(x*y) as input and returns the p_lo such that
// x * y = xy + p_lo holds exactly.
template <typename Packet, std::enable_if_t<!has_fast_fma<Packet>::value, int> = 0>
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet twoprod_low(const Packet& x, const Packet& y, const Packet& xy) {
  Packet x_hi, x_lo, y_hi, y_lo;
  veltkamp_splitting(x, x_hi, x_lo);
  veltkamp_splitting(y, y_hi, y_lo);

  Packet p_lo = pmadd(x_hi, y_hi, pnegate(xy));
  p_lo = pmadd(x_hi, y_lo, p_lo);
  p_lo = pmadd(x_lo, y_hi, p_lo);
  p_lo = pmadd(x_lo, y_lo, p_lo);
  return p_lo;
}

#endif  // EIGEN_VECTORIZE_FMA

// This function implements the addition of two double word numbers
// represented by {x_hi, x_lo} and {y_hi, y_lo}, with operands of either sign
// and magnitude: Knuth's two-sum of the high parts is error-free, the low
// parts are added to its error term, and the result is renormalized.
// This is Algorithm 5 (SloppyDWPlusDW) from Joldes, Muller, & Popescu (2017),
// "Tight and rigorous error bounds for basic building blocks of double-word
// arithmetic", https://hal.science/hal-01351529: the error is a few u^2
// relative to |x| + |y|, not to the sum, which may cancel.
template <typename Packet>
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void twosum(const Packet& x_hi, const Packet& x_lo, const Packet& y_hi,
                                                  const Packet& y_lo, Packet& s_hi, Packet& s_lo) {
  Packet s = padd(x_hi, y_hi);
  Packet y_part = psub(s, x_hi);
  Packet err = padd(psub(x_hi, psub(s, y_part)), psub(y_hi, y_part));
  fast_twosum(s, padd(err, padd(x_lo, y_lo)), s_hi, s_lo);
}

// The difference {x_hi, x_lo} - {y_hi, y_lo} of two double word numbers,
// the same algorithm as twosum with the sign of y folded into the operations.
template <typename Packet>
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void twodiff(const Packet& x_hi, const Packet& x_lo, const Packet& y_hi,
                                                   const Packet& y_lo, Packet& s_hi, Packet& s_lo) {
  Packet s = psub(x_hi, y_hi);
  Packet y_part = psub(x_hi, s);
  Packet err = psub(psub(x_hi, padd(s, y_part)), psub(y_hi, y_part));
  fast_twosum(s, padd(err, psub(x_lo, y_lo)), s_hi, s_lo);
}

// This is a version of twosum for double word numbers,
// which assumes that |x_hi| >= |y_hi|.
template <typename Packet>
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void fast_twosum(const Packet& x_hi, const Packet& x_lo, const Packet& y_hi,
                                                       const Packet& y_lo, Packet& s_hi, Packet& s_lo) {
  Packet r_hi, r_lo;
  fast_twosum(x_hi, y_hi, r_hi, r_lo);
  const Packet s = padd(padd(y_lo, r_lo), x_lo);
  fast_twosum(r_hi, s, s_hi, s_lo);
}

// This is a version of twosum for adding a floating point number x to
// double word number {y_hi, y_lo} number, with the assumption
// that |x| >= |y_hi|.
template <typename Packet>
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void fast_twosum(const Packet& x, const Packet& y_hi, const Packet& y_lo,
                                                       Packet& s_hi, Packet& s_lo) {
  Packet r_hi, r_lo;
  fast_twosum(x, y_hi, r_hi, r_lo);
  const Packet s = padd(y_lo, r_lo);
  fast_twosum(r_hi, s, s_hi, s_lo);
}

// This function implements the multiplication of a double word
// number represented by {x_hi, x_lo} by a floating point number y.
// It returns the result as a pair {p_hi, p_lo} such that
// (x_hi + x_lo) * y = p_hi + p_lo holds with a relative error
// of less than 2*2^{-2p}, where p is the number of significand bits
// in the floating point type.
// This is Algorithm 7 from Jean-Michel Muller, "Elementary Functions",
// 3rd edition, Birkh\"auser, 2016.
template <typename Packet>
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void twoprod(const Packet& x_hi, const Packet& x_lo, const Packet& y,
                                                   Packet& p_hi, Packet& p_lo) {
  Packet c_hi, c_lo1;
  twoprod(x_hi, y, c_hi, c_lo1);
  const Packet c_lo2 = pmul(x_lo, y);
  Packet t_hi, t_lo1;
  fast_twosum(c_hi, c_lo2, t_hi, t_lo1);
  const Packet t_lo2 = padd(t_lo1, c_lo1);
  fast_twosum(t_hi, t_lo2, p_hi, p_lo);
}

// This function implements the multiplication of two double word
// numbers represented by {x_hi, x_lo} and {y_hi, y_lo}.
// For normalized inputs (x_hi = RN(x_hi + x_lo), likewise for y), the result
// {p_hi, p_lo} has relative error < 5*u^2 for p-bit significands, u = 2^{-p}, p >= 8.
// This assumes round-to-nearest, no overflow/underflow, an exact two-product
// residual, and separately rounded operations outside that residual.
// Each DW*FP product has relative error <= (3/2 + 4*u)*u^2: Joldes, Muller,
// and Popescu (2017), "Tight and rigorous error bounds for basic building
// blocks of double-word arithmetic", Theorem 4.1, https://hal.science/hal-01351529.
// The final sum adds absolute error <= u^2*((3+3*u+u^2)*abs(p_hi_hi)
// + (4+3*u+u^2)*abs(p_lo_hi)); using abs(y_lo) <= u*abs(y_hi) gives the bound above.
template <typename Packet>
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void twoprod(const Packet& x_hi, const Packet& x_lo, const Packet& y_hi,
                                                   const Packet& y_lo, Packet& p_hi, Packet& p_lo) {
  Packet p_hi_hi, p_hi_lo;
  twoprod(x_hi, x_lo, y_hi, p_hi_hi, p_hi_lo);
  Packet p_lo_hi, p_lo_lo;
  twoprod(x_hi, x_lo, y_lo, p_lo_hi, p_lo_lo);
  fast_twosum(p_hi_hi, p_hi_lo, p_lo_hi, p_lo_lo, p_hi, p_lo);
}

// This function implements the multiplication of two double word numbers
// {x_hi, x_lo} and {y_hi, y_lo} at a lower cost than twoprod above: the
// cross terms are accumulated onto the two-product residual and the x_lo*y_lo
// term is dropped. For normalized inputs the relative error is < 7*u^2.
// This is Algorithm 10 (DWTimesDW1) from Joldes, Muller, & Popescu (2017),
// "Tight and rigorous error bounds for basic building blocks of double-word
// arithmetic", https://hal.science/hal-01351529.
template <typename Packet>
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void fast_twoprod(const Packet& x_hi, const Packet& x_lo, const Packet& y_hi,
                                                        const Packet& y_lo, Packet& p_hi, Packet& p_lo) {
  Packet c_hi, c_lo;
  twoprod(x_hi, y_hi, c_hi, c_lo);
  c_lo = pmadd(x_hi, y_lo, pmadd(x_lo, y_hi, c_lo));
  fast_twosum(c_hi, c_lo, p_hi, p_lo);
}

// This function implements the division of double word {x_hi, x_lo}
// by float y. This is Algorithm 15 from "Tight and rigorous error bounds
// for basic building blocks of double-word arithmetic", Joldes, Muller, & Popescu,
// 2017. https://hal.archives-ouvertes.fr/hal-01351529
template <typename Packet>
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void doubleword_div_fp(const Packet& x_hi, const Packet& x_lo, const Packet& y,
                                                             Packet& z_hi, Packet& z_lo) {
  const Packet t_hi = pdiv(x_hi, y);
  Packet pi_hi, pi_lo;
  twoprod(t_hi, y, pi_hi, pi_lo);
  const Packet delta_hi = psub(x_hi, pi_hi);
  const Packet delta_t = psub(delta_hi, pi_lo);
  const Packet delta = padd(delta_t, x_lo);
  const Packet t_lo = pdiv(delta, y);
  fast_twosum(t_hi, t_lo, z_hi, z_lo);
}

}  // end namespace internal
}  // end namespace Eigen

#endif  // EIGEN_ARCH_GENERIC_PACKET_MATH_DOUBLE_WORD_H