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