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
17namespace LIBC_NAMESPACE::fputil {
18
19using 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.
25LIBC_INLINE constexpr DoubleDouble exact_add(double a, double b) {
26 DoubleDouble r{.lo: 0.0, .hi: 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|
34LIBC_INLINE constexpr DoubleDouble add(const DoubleDouble &a,
35 const DoubleDouble &b) {
36 DoubleDouble r = exact_add(a: a.hi, b: b.hi);
37 double lo = a.lo + b.lo;
38 return exact_add(a: r.hi, b: r.lo + lo);
39}
40
41// Assumption: |a.hi| >= |b|
42LIBC_INLINE constexpr DoubleDouble add(const DoubleDouble &a, double b) {
43 DoubleDouble r = exact_add(a: a.hi, b);
44 return exact_add(a: r.hi, b: r.lo + a.lo);
45}
46
47// Velkamp's Splitting for double precision.
48LIBC_INLINE constexpr DoubleDouble split(double a) {
49 DoubleDouble r{.lo: 0.0, .hi: 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
59LIBC_INLINE DoubleDouble exact_mult(double a, double b) {
60 DoubleDouble r{.lo: 0.0, .hi: 0.0};
61
62#ifdef LIBC_TARGET_CPU_HAS_FMA
63 r.hi = a * b;
64 r.lo = fputil::multiply_add(x: a, y: b, z: -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
79LIBC_INLINE DoubleDouble quick_mult(double a, const DoubleDouble &b) {
80 DoubleDouble r = exact_mult(a, b: b.hi);
81 r.lo = multiply_add(x: a, y: b.lo, z: r.lo);
82 return r;
83}
84
85LIBC_INLINE DoubleDouble quick_mult(const DoubleDouble &a,
86 const DoubleDouble &b) {
87 DoubleDouble r = exact_mult(a: a.hi, b: b.hi);
88 double t1 = multiply_add(x: a.hi, y: b.lo, z: r.lo);
89 double t2 = multiply_add(x: a.lo, y: b.hi, z: t1);
90 r.lo = t2;
91 return r;
92}
93
94// Assuming |c| >= |a * b|.
95template <>
96LIBC_INLINE DoubleDouble multiply_add<DoubleDouble>(const DoubleDouble &a,
97 const DoubleDouble &b,
98 const DoubleDouble &c) {
99 return add(a: c, b: quick_mult(a, b));
100}
101
102} // namespace LIBC_NAMESPACE::fputil
103
104#endif // LLVM_LIBC_SRC___SUPPORT_FPUTIL_DOUBLE_DOUBLE_H
105

source code of libc/src/__support/FPUtil/double_double.h