eonc::exact Namespace Reference

Functions
double	dot (const double *x, const double *y, std::size_t n)
	Round-to-nearest-even value of the exact sum of x[i] * y[i].
double	orderedDot (const double *x, const double *y, std::size_t n)
	A fixed summation tree: four lanes (term i goes to lane i mod 4, each lane a left-to-right sum of rounded products), then (l0 + l1) + (l2 + l3).
double	sum (const double *x, std::size_t n)
	Round-to-nearest-even value of the exact sum of x[i].
double	sumStrided (const double *x, std::size_t n, std::size_t stride)
	Same sum for x[i * stride].
template<typename A, typename B>
double	dot (const Eigen::MatrixBase< A > &x, const Eigen::MatrixBase< B > &y)
template<typename A, typename B>
double	orderedDot (const Eigen::MatrixBase< A > &x, const Eigen::MatrixBase< B > &y)
template<typename A>
double	squaredNorm (const Eigen::MatrixBase< A > &x)
template<typename A>
double	norm (const Eigen::MatrixBase< A > &x)
	sqrt of the correctly rounded squared norm; IEEE sqrt is itself correctly rounded, so this is reproducible as well.
template<typename A>
double	orderedNorm (const Eigen::MatrixBase< A > &x)
	sqrt of orderedDot(x, x): reproducible, not correctly rounded.
template<typename A>
A::PlainObject	normalized (const Eigen::MatrixBase< A > &x, double min_norm=1e-300)
	x / norm(x), or zero below min_norm (the safemath convention).

Variables
constexpr std::size_t	kLanes = 4

Function Documentation

dot() [1/2]

double eonc::exact::dot ( const double * x, const double * y, std::size_t n )

nodiscard

Round-to-nearest-even value of the exact sum of x[i] * y[i].

An exact zero returns +0. A non-finite input returns the plain IEEE sum, so NaN and infinity propagate as they would without this function.

Definition at line 451 of file ExactReduction.cpp.

                                                          {
  double certified;
  if (certifiedDot(x, y, n, certified)) {
    return certified;
  }
  std::array<Accumulator, kLanes> acc;
  for (std::size_t i = 0; i < n; ++i) {
    Decoded a;
    Decoded b;
    if (!decode(x[i], a) || !decode(y[i], b)) [[unlikely]] {
      double plain = 0.0;
      for (std::size_t j = 0; j < n; ++j) {
        plain += x[j] * y[j];
      }
      return plain;
    }
    uint64_t hi;
    uint64_t lo;
    mul64(a.m, b.m, hi, lo);
    acc[i % kLanes].add(hi, lo, a.e + b.e, a.neg != b.neg);
    if ((i + 1) % kNormalizeEvery == 0) [[unlikely]] {
      for (auto &lane : acc) {
        lane.fold();
      }
    }
  }
  for (std::size_t k = 1; k < kLanes; ++k) {
    acc[0].absorb(acc[k]);
  }
  return acc[0].round();
}

dot() [2/2]

template<typename A, typename B>

double eonc::exact::dot ( const Eigen::MatrixBase< A > & x, const Eigen::MatrixBase< B > & y )

nodiscard

Definition at line 67 of file ExactReduction.h.

                                                      {
  eigen_assert(x.size() == y.size());
  const auto &xe = x.derived().eval();
  const auto &ye = y.derived().eval();
  return dot(xe.data(), ye.data(), static_cast<std::size_t>(xe.size()));
}

norm()

template<typename A>

double eonc::exact::norm ( const Eigen::MatrixBase< A > & x )

nodiscard

sqrt of the correctly rounded squared norm; IEEE sqrt is itself correctly rounded, so this is reproducible as well.

Definition at line 92 of file ExactReduction.h.

                                                                             {
  return std::sqrt(squaredNorm(x));
}

normalized()

template<typename A>

A::PlainObject eonc::exact::normalized ( const Eigen::MatrixBase< A > & x, double min_norm = 1e-300 )

nodiscard

x / norm(x), or zero below min_norm (the safemath convention).

Definition at line 106 of file ExactReduction.h.

                                                                         {
  const double n = norm(x);
  if (n < min_norm) [[unlikely]] {
    return A::PlainObject::Zero(x.rows(), x.cols());
  }
  return (x / n).eval();
}

orderedDot() [1/2]

double eonc::exact::orderedDot ( const double * x, const double * y, std::size_t n )

nodiscard

A fixed summation tree: four lanes (term i goes to lane i mod 4, each lane a left-to-right sum of rounded products), then (l0 + l1) + (l2 + l3).

Not exact, but every operation is a single IEEE multiply or add in an order this translation unit fixes (it is compiled without contraction), so the bits match on every compiler, Eigen release and CPU. For bulk products that set a vector, not a branch.

Definition at line 483 of file ExactReduction.cpp.

                                                                 {
  double l0 = 0.0;
  double l1 = 0.0;
  double l2 = 0.0;
  double l3 = 0.0;
  std::size_t i = 0;
  for (; i + 4 <= n; i += 4) {
    const double p0 = x[i] * y[i];
    const double p1 = x[i + 1] * y[i + 1];
    const double p2 = x[i + 2] * y[i + 2];
    const double p3 = x[i + 3] * y[i + 3];
    l0 += p0;
    l1 += p1;
    l2 += p2;
    l3 += p3;
  }
  if (i < n) {
    const double p = x[i] * y[i];
    l0 += p;
  }
  if (i + 1 < n) {
    const double p = x[i + 1] * y[i + 1];
    l1 += p;
  }
  if (i + 2 < n) {
    const double p = x[i + 2] * y[i + 2];
    l2 += p;
  }
  return (l0 + l1) + (l2 + l3);
}

orderedDot() [2/2]

template<typename A, typename B>

double eonc::exact::orderedDot ( const Eigen::MatrixBase< A > & x, const Eigen::MatrixBase< B > & y )

nodiscard

Definition at line 76 of file ExactReduction.h.

                                                             {
  eigen_assert(x.size() == y.size());
  const auto &xe = x.derived().eval();
  const auto &ye = y.derived().eval();
  return orderedDot(xe.data(), ye.data(), static_cast<std::size_t>(xe.size()));
}

orderedNorm()

template<typename A>

double eonc::exact::orderedNorm ( const Eigen::MatrixBase< A > & x )

nodiscard

sqrt of orderedDot(x, x): reproducible, not correctly rounded.

For the per-atom lengths of a step, where an exact sum per atom costs more than the step.

Definition at line 100 of file ExactReduction.h.

                                                              {
  return std::sqrt(orderedDot(x, x));
}

squaredNorm()

template<typename A>

double eonc::exact::squaredNorm ( const Eigen::MatrixBase< A > & x )

nodiscard

Definition at line 85 of file ExactReduction.h.

                                                              {
  const auto &xe = x.derived().eval();
  return dot(xe.data(), xe.data(), static_cast<std::size_t>(xe.size()));
}

sum()

double eonc::exact::sum ( const double * x, std::size_t n )

nodiscard

Round-to-nearest-even value of the exact sum of x[i].

Definition at line 514 of file ExactReduction.cpp.

{ return sumStrided(x, n, 1); }

sumStrided()

double eonc::exact::sumStrided ( const double * x, std::size_t n, std::size_t stride )

nodiscard

Same sum for x[i * stride].

A row-major atom matrix stores one column this way. stride 1 is sum().

Definition at line 516 of file ExactReduction.cpp.

                                                                  {
  Accumulator acc;
  for (std::size_t i = 0; i < n; ++i) {
    Decoded a;
    if (!decode(x[i * stride], a)) [[unlikely]] {
      double plain = 0.0;
      for (std::size_t j = 0; j < n; ++j) {
        plain += x[j * stride];
      }
      return plain;
    }
    acc.add(0, a.m, a.e, a.neg);
    if ((i + 1) % kNormalizeEvery == 0) [[unlikely]] {
      acc.fold();
    }
  }
  return acc.round();
}

Variable Documentation

kLanes

std::size_t eonc::exact::kLanes = 4

constexpr

Definition at line 278 of file ExactReduction.cpp.