HPReal.h

// -*- tab-width: 2; indent-tabs-mode: nil; coding: utf-8-with-signature -*-
//-----------------------------------------------------------------------------
// Copyright 2000-2026 CEA (www.cea.fr) IFPEN (www.ifpenergiesnouvelles.com)
// See the top-level COPYRIGHT file for details.
// SPDX-License-Identifier: Apache-2.0
//-----------------------------------------------------------------------------
/*---------------------------------------------------------------------------*/
/* HPReal.h                                                    (C) 2000-2019 */
/*                                                                           */
/* High-precision real number.                                               */
/*---------------------------------------------------------------------------*/
#ifndef ARCANE_UTILS_HPREAL_H
#define ARCANE_UTILS_HPREAL_H
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
#include "arcane/utils/Numeric.h"
 
#include <iosfwd>
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
namespace Arcane
{
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/

class ARCANE_UTILS_EXPORT HPReal
{
 public:

  HPReal() {}

  explicit HPReal(double avalue)
  : m_value(avalue)
  , m_correction(0.0)
  {}

  HPReal(double avalue, double acorrection)
  : m_value(avalue)
  , m_correction(acorrection)
  {}

  Real value() const { return m_value; }

  Real correction() const { return m_correction; }

  void operator+=(Real v)
  {
    *this = accumulate(v, *this);
  }

  inline void operator+=(HPReal v)
  {
    *this = reduce(*this, v);
  }

  inline Real toReal() const
  {
    return value() + correction();
  }

  HPReal reduce(HPReal b) const
  {
    return reduce(*this, b);
  }

  void operator*=(Real v)
  {
    *this = product(v, *this);
  }

  inline void operator*=(HPReal v)
  {
    *this = product(*this, v);
  }

  void operator/=(Real v)
  {
    *this = div2(v, *this);
  }

  inline void operator/=(HPReal v)
  {
    *this = div2(*this, v);
  }

  std::istream& assign(std::istream& i);
  std::ostream& print(std::ostream& o) const;
  std::ostream& printPretty(std::ostream& o) const;
 
 public:

  static HPReal zero() { return HPReal(0.0); }
 
 public:
 
  static HPReal accumulate(Real a, HPReal b)
  {
    HPReal x(_doTwoSum(a, b.value()));
    Real c = x.correction() + b.correction();
    return _doQuickTwoSum(x.value(), c);
  }
 
  // Mpi passes through
  static HPReal reduce(HPReal a, HPReal b)
  {
    HPReal x(_doTwoSum(a.value(), b.value()));
    Real c = x.correction() + a.correction() + b.correction();
    return _doQuickTwoSum(x.value(), c);
  }
 
  // algo 11 of AC_TWOPRODUCTS Ref 2
  static HPReal product(HPReal a, HPReal b)
  {
    HPReal x(_doTwoProducts(a.value(), b.value()));
    Real w = x.correction() + (a.value() * b.correction() + b.value() * a.correction());
    return HPReal(x.value(), w);
  }
  // algo 11 of AC_TWOPRODUCTS Ref 2
  static HPReal product(Real a, HPReal b)
  {
    HPReal x(_doTwoProducts(a, b.value()));
    Real w = x.correction() + (a * b.correction());
    return HPReal(x.value(), w);
  }
 
  // algo div2 Ref 6    div2 =  a/b
  // I put parentheses in the evaluation of w to force an order sometimes the compiler options
  // ATTENTION the last 2 lines of the div2 algorithm are not included (we do not renormalize)
  // The "mul2" algorithm also contains them and did not include them for "product" because algorithm 11 of Ref 2 does not have them (in my humble opinion)
  static HPReal div2(HPReal a, HPReal b)
  {
    Real c = a.value() / b.value();
    HPReal u(_doTwoProducts(c, b.value()));
    Real w = ((((a.value() - u.value()) + -u.correction()) + a.correction()) - c * b.correction()) / b.value();
    return HPReal(c, w);
  }
  // algo div2 Ref 6     div2 =  a/b
  static HPReal div2(Real b, HPReal a)
  {
    Real c = a.value() / b;
    HPReal u(_doTwoProducts(c, b));
    Real w = (((a.value() - u.value()) + -u.correction()) + a.correction()) / b;
    return HPReal(c, w);
  }
 
 private:
 
  Real m_value;
  Real m_correction;
 
 private:
 
  // Correction of ()
  static HPReal _doTwoSum(Real a, Real b)
  {
    Real value = a + b;
    Real approx_b = value - a;
    Real sum_error = (a - (value - approx_b)) + (b - approx_b);
    return HPReal(value, sum_error);
  }
  // correction tests of the values of a and b and we put absolute values
  static HPReal _doQuickTwoSum(Real a1, Real b1)
  {
    Real a = a1;
    Real b = b1;
    if (std::abs(b1) > std::abs(a1)) {
      a = b1;
      b = a1;
    }
    Real value = a + b;
    Real error_value = (b - (value - a));
    return HPReal(value, error_value);
  }
  // algorithm 4 ref 2 or algorithm 2.4 ref 3 or algorithm 6 p.4
  static HPReal _doTwoProducts(Real a, Real b)
  {
    Real x = a * b;
    HPReal aw = SPLIT(a);
    HPReal bw = SPLIT(b);
    Real y = aw.correction() * bw.correction() - (((x - aw.value() * bw.value()) - aw.correction() * bw.value()) - aw.value() * bw.correction());
    return HPReal(x, y);
  }
  static HPReal SPLIT(Real a)
  {
    const Real f = 134217729; // 1+ 2^ceil(52/2);
    Real c = f * a;
    Real x = c - (c - a);
    Real y = a - x;
    return HPReal(x, y);
  }
};
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
inline bool
operator<(const HPReal& a, const HPReal& b)
{
  return a.toReal() < b.toReal();
}
 
inline bool
operator>(const HPReal& a, const HPReal& b)
{
  return a.toReal() > b.toReal();
}
 
inline bool
operator==(const HPReal& a, const HPReal& b)
{
  return a.value() == b.value() && a.correction() == b.correction();
}
 
inline bool
operator!=(const HPReal& a, const HPReal& b)
{
  return !operator==(a, b);
}
 
inline HPReal
operator+(const HPReal& a, const HPReal& b)
{
  return HPReal::reduce(a, b);
}
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
inline std::ostream&
operator<<(std::ostream& o, HPReal t)
{
  return t.printPretty(o);
}
 
inline std::istream&
operator>>(std::istream& i, HPReal& t)
{
  return t.assign(i);
}
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
} // End namespace Arcane
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
#endif