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_DOUBLE_DOUBLE_H
10 #define LLVM_LIBC_SRC___SUPPORT_FPUTIL_DOUBLE_DOUBLE_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 LIBC_NAMESPACE::fputil {
18
19 using DoubleDouble = LIBC_NAMESPACE::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.
exact_add(double a,double b)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|
add(const DoubleDouble & a,const DoubleDouble & b)34 LIBC_INLINE constexpr DoubleDouble add(const DoubleDouble &a,
35 const DoubleDouble &b) {
36 DoubleDouble r = exact_add(a.hi, b.hi);
37 double lo = a.lo + b.lo;
38 return exact_add(r.hi, r.lo + lo);
39 }
40
41 // Assumption: |a.hi| >= |b|
add(const DoubleDouble & a,double b)42 LIBC_INLINE constexpr DoubleDouble add(const DoubleDouble &a, double b) {
43 DoubleDouble r = exact_add(a.hi, b);
44 return exact_add(r.hi, r.lo + a.lo);
45 }
46
47 // Velkamp's Splitting for double precision.
split(double a)48 LIBC_INLINE constexpr DoubleDouble split(double a) {
49 DoubleDouble r{0.0, 0.0};
50 // Splitting constant = 2^ceil(prec(double)/2) + 1 = 2^27 + 1.
51 constexpr double C = 0x1.0p27 + 1.0;
52 double t1 = C * a;
53 double t2 = a - t1;
54 r.hi = t1 + t2;
55 r.lo = a - r.hi;
56 return r;
57 }
58
exact_mult(double a,double b)59 LIBC_INLINE DoubleDouble exact_mult(double a, double b) {
60 DoubleDouble r{0.0, 0.0};
61
62 #ifdef LIBC_TARGET_CPU_HAS_FMA
63 r.hi = a * b;
64 r.lo = fputil::multiply_add(a, b, -r.hi);
65 #else
66 // Dekker's Product.
67 DoubleDouble as = split(a);
68 DoubleDouble bs = split(b);
69 r.hi = a * b;
70 double t1 = as.hi * bs.hi - r.hi;
71 double t2 = as.hi * bs.lo + t1;
72 double t3 = as.lo * bs.hi + t2;
73 r.lo = as.lo * bs.lo + t3;
74 #endif // LIBC_TARGET_CPU_HAS_FMA
75
76 return r;
77 }
78
quick_mult(double a,const DoubleDouble & b)79 LIBC_INLINE DoubleDouble quick_mult(double a, const DoubleDouble &b) {
80 DoubleDouble r = exact_mult(a, b.hi);
81 r.lo = multiply_add(a, b.lo, r.lo);
82 return r;
83 }
84
quick_mult(const DoubleDouble & a,const DoubleDouble & b)85 LIBC_INLINE DoubleDouble quick_mult(const DoubleDouble &a,
86 const DoubleDouble &b) {
87 DoubleDouble r = exact_mult(a.hi, b.hi);
88 double t1 = multiply_add(a.hi, b.lo, r.lo);
89 double t2 = multiply_add(a.lo, b.hi, t1);
90 r.lo = t2;
91 return r;
92 }
93
94 // Assuming |c| >= |a * b|.
95 template <>
96 LIBC_INLINE DoubleDouble multiply_add<DoubleDouble>(const DoubleDouble &a,
97 const DoubleDouble &b,
98 const DoubleDouble &c) {
99 return add(c, quick_mult(a, b));
100 }
101
102 } // namespace LIBC_NAMESPACE::fputil
103
104 #endif // LLVM_LIBC_SRC___SUPPORT_FPUTIL_DOUBLE_DOUBLE_H
105