libc/src/__support/FPUtil/double_double.h

   1 //===-- Utilities for double-double data type. ------------------*- C++ -*-===//
   2 //
   3 // Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
   4 // See https://llvm.org/LICENSE.txt for license information.
   5 // SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
   6 //
   7 //===----------------------------------------------------------------------===//
   8
   9 #ifndef LLVM_LIBC_SRC_SUPPORT_FPUTIL_DOUBLEDOUBLE_H
  10 #define LLVM_LIBC_SRC_SUPPORT_FPUTIL_DOUBLEDOUBLE_H
  11
  12 #include "multiply_add.h"
  13 #include "src/__support/common.h"
  14 #include "src/__support/macros/properties/cpu_features.h" // LIBC_TARGET_CPU_HAS_FMA
  15 #include "src/__support/number_pair.h"
  16
  17 namespace __llvm_libc::fputil {
  18
  19 using DoubleDouble = __llvm_libc::NumberPair<double>;
  20
  21 // The output of Dekker's FastTwoSum algorithm is correct, i.e.:
  22 //   r.hi + r.lo = a + b exactly
  23 //   and |r.lo| < eps(r.lo)
  24 // if ssumption: |a| >= |b|, or a = 0.
  25 LIBC_INLINE constexpr DoubleDouble exact_add(double a, double b) {
  26   DoubleDouble r{0.0, 0.0};
  27   r.hi = a + b;
  28   double t = r.hi - a;
  29   r.lo = b - t;
  30   return r;
  31 }
  32
  33 // Assumption: |a.hi| >= |b.hi|
  34 LIBC_INLINE constexpr DoubleDouble add(DoubleDouble a, DoubleDouble b) {
  35   DoubleDouble r = exact_add(a.hi, b.hi);
  36   double lo = a.lo + b.lo;
  37   return exact_add(r.hi, r.lo + lo);
  38 }
  39
  40 // Assumption: |a.hi| >= |b|
  41 LIBC_INLINE constexpr DoubleDouble add(DoubleDouble a, double b) {
  42   DoubleDouble r = exact_add(a.hi, b);
  43   return exact_add(r.hi, r.lo + a.lo);
  44 }
  45
  46 // Velkamp's Splitting for double precision.
  47 LIBC_INLINE constexpr DoubleDouble split(double a) {
  48   DoubleDouble r{0.0, 0.0};
  49   // Splitting constant = 2^ceil(prec(double)/2) + 1 = 2^27 + 1.
  50   constexpr double C = 0x1.0p27 + 1.0;
  51   double t1 = C * a;
  52   double t2 = a - t1;
  53   r.hi = t1 + t2;
  54   r.lo = a - r.hi;
  55   return r;
  56 }
  57
  58 LIBC_INLINE DoubleDouble exact_mult(double a, double b) {
  59   DoubleDouble r{0.0, 0.0};
  60
  61 #ifdef LIBC_TARGET_CPU_HAS_FMA
  62   r.hi = a * b;
  63   r.lo = fputil::multiply_add(a, b, -r.hi);
  64 #else
  65   // Dekker's Product.
  66   DoubleDouble as = split(a);
  67   DoubleDouble bs = split(b);
  68   r.hi = a * b;
  69   double t1 = as.hi * bs.hi - r.hi;
  70   double t2 = as.hi * bs.lo + t1;
  71   double t3 = as.lo * bs.hi + t2;
  72   r.lo = as.lo * bs.lo + t3;
  73 #endif // LIBC_TARGET_CPU_HAS_FMA
  74
  75   return r;
  76 }
  77
  78 LIBC_INLINE DoubleDouble quick_mult(DoubleDouble a, DoubleDouble b) {
  79   DoubleDouble r = exact_mult(a.hi, b.hi);
  80   double t1 = fputil::multiply_add(a.hi, b.lo, r.lo);
  81   double t2 = fputil::multiply_add(a.lo, b.hi, t1);
  82   r.lo = t2;
  83   return r;
  84 }
  85
  86 } // namespace __llvm_libc::fputil
  87
  88 #endif // LLVM_LIBC_SRC_SUPPORT_FPUTIL_DOUBLEDOUBLE_H