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/config.h"
15#include "src/__support/macros/properties/cpu_features.h" // LIBC_TARGET_CPU_HAS_FMA
16#include "src/__support/number_pair.h"
17
18namespace LIBC_NAMESPACE_DECL {
19namespace fputil {
20
21template <typename T> struct DefaultSplit;
22template <> struct DefaultSplit<float> {
23 static constexpr size_t VALUE = 12;
24};
25template <> struct DefaultSplit<double> {
26 static constexpr size_t VALUE = 27;
27};
28
29using DoubleDouble = NumberPair<double>;
30using FloatFloat = NumberPair<float>;
31
32// Dekker's FastTwoSum:
33// r_hi + r_lo ~ a + b
34// In precision `p`, it's exact if `lsb(a) >= ulp(b)` and one of the following
35// conditions satisfies:
36// (i) rounding mode = RN
37// (ii) e_a - e_b <= p
38// (iii) rounding mode = RD and b >= 0
39// (iv) rounding mode = RU and b <= 0
40// (v) rounding mode = RZ and ab >= 0
41// Otherwise, the errors err = (a + b) - (r_hi + r_lo) is bounded by:
42// |err| <= 2^(-2p + 1) * ufp(a + b) <= 2^(-2p + 1) * ufp(r_hi)
43// when `lsb(a) >= ulp(b)`.
44// In particular,
45// |err| <= 2^-47 * ufp(r_hi) for single precision,
46// <= 2^-105 * ufp(r_hi) for double precision.
47// If the condition `lsb(a) >= ulp(b)` does NOT satisfy, then:
48// |err| < 3 * 2^(-p) * |r_hi|.
49// Reference:
50// Jeannerod, C.-P. and Zimmermann, P., "FastTwoSum revisited", ARITH 2025,
51// <https://www.arith2025.org/proceedings/215900a141.pdf>.
52template <bool FAST2SUM = true, typename T = double>
53LIBC_INLINE constexpr NumberPair<T> exact_add(T a, T b) {
54 NumberPair<T> r{0.0, 0.0};
55 if constexpr (FAST2SUM) {
56 r.hi = a + b;
57 T t = r.hi - a;
58 r.lo = b - t;
59 } else {
60 r.hi = a + b;
61 T t1 = r.hi - a;
62 T t2 = r.hi - t1;
63 T t3 = b - t1;
64 T t4 = a - t2;
65 r.lo = t3 + t4;
66 }
67 return r;
68}
69
70// Following the above analysis of FastTwoSum,
71// If `lsb(a.hi) >= ulp(b.hi)`, then the errors:
72// err = (a.hi + a.lo + b.hi + b.lo) - (r.hi - r.lo) is bounded by:
73// |err| < 2^(-2p) * ufp(r.hi) + 2^(-p + 2) * ufp(a.lo + b.lo).
74template <bool FAST2SUM = true, typename T>
75LIBC_INLINE constexpr NumberPair<T> add(const NumberPair<T> &a,
76 const NumberPair<T> &b) {
77 NumberPair<T> r = exact_add<FAST2SUM>(a.hi, b.hi);
78 T lo = a.lo + b.lo;
79 T r_lo = r.lo + lo;
80 return exact_add<FAST2SUM>(r.hi, r_lo);
81}
82
83// Assumption: when FAST2SUM = true, |a.hi| >= |b|
84template <bool FAST2SUM = true, typename T>
85LIBC_INLINE constexpr NumberPair<T> add(const NumberPair<T> &a, T b) {
86 NumberPair<T> r = exact_add<FAST2SUM>(a.hi, b);
87 T r_lo = r.lo + a.lo;
88 return exact_add<FAST2SUM>(r.hi, r_lo);
89}
90
91// Veltkamp's Splitting for double precision.
92// Note: This is proved to be correct for all rounding modes:
93// Zimmermann, P., "Note on the Veltkamp/Dekker Algorithms with Directed
94// Roundings," https://inria.hal.science/hal-04480440.
95// Default splitting constant = 2^ceil(prec(double)/2) + 1 = 2^27 + 1.
96template <typename T = double, size_t N = DefaultSplit<T>::VALUE>
97LIBC_INLINE constexpr NumberPair<T> split(T a) {
98 NumberPair<T> r{0.0, 0.0};
99 // CN = 2^N.
100 constexpr T CN = static_cast<T>(1 << N);
101 constexpr T C = CN + T(1);
102 T t1 = C * a;
103 T t2 = a - t1;
104 r.hi = t1 + t2;
105 r.lo = a - r.hi;
106 return r;
107}
108
109// Helper for non-fma exact mult where the first number is already split.
110template <typename T = double, size_t SPLIT_B = DefaultSplit<T>::VALUE>
111LIBC_INLINE constexpr NumberPair<T> exact_mult(const NumberPair<T> &as, T a,
112 T b) {
113 NumberPair<T> bs = split<T, SPLIT_B>(b);
114 NumberPair<T> r{0.0, 0.0};
115
116 r.hi = a * b;
117 T t1 = as.hi * bs.hi - r.hi;
118 T t2 = as.hi * bs.lo + t1;
119 T t3 = as.lo * bs.hi + t2;
120 r.lo = as.lo * bs.lo + t3;
121
122 return r;
123}
124
125// The templated exact multiplication needs template version of
126// LIBC_TARGET_CPU_HAS_FMA_* macro to correctly select the implementation.
127// These can be moved to "src/__support/macros/properties/cpu_features.h" if
128// other part of libc needed.
129template <typename T> struct TargetHasFmaInstruction {
130 static constexpr bool VALUE = false;
131};
132
133#ifdef LIBC_TARGET_CPU_HAS_FMA_FLOAT
134template <> struct TargetHasFmaInstruction<float> {
135 static constexpr bool VALUE = true;
136};
137#endif // LIBC_TARGET_CPU_HAS_FMA_FLOAT
138
139#ifdef LIBC_TARGET_CPU_HAS_FMA_DOUBLE
140template <> struct TargetHasFmaInstruction<double> {
141 static constexpr bool VALUE = true;
142};
143#endif // LIBC_TARGET_CPU_HAS_FMA_DOUBLE
144
145// Note: When FMA instruction is not available, the `exact_mult` function is
146// only correct for round-to-nearest mode. See:
147// Zimmermann, P., "Note on the Veltkamp/Dekker Algorithms with Directed
148// Roundings," https://inria.hal.science/hal-04480440.
149// Using Theorem 1 in the paper above, without FMA instruction, if we restrict
150// the generated constants to precision <= 51, and splitting it by 2^28 + 1,
151// then a * b = r.hi + r.lo is exact for all rounding modes.
152template <typename T = double, size_t SPLIT_B = DefaultSplit<T>::VALUE>
153LIBC_INLINE LIBC_CONSTEXPR NumberPair<T> exact_mult(T a, T b) {
154 NumberPair<T> r{0.0, 0.0};
155
156 if constexpr (TargetHasFmaInstruction<T>::VALUE) {
157 r.hi = a * b;
158 r.lo = fputil::multiply_add(a, b, -r.hi);
159 } else {
160 // Dekker's Product.
161 NumberPair<T> as = split(a);
162
163 r = exact_mult<T, SPLIT_B>(as, a, b);
164 }
165
166 return r;
167}
168
169template <typename T = double>
170LIBC_INLINE NumberPair<T> quick_mult(T a, const NumberPair<T> &b) {
171 NumberPair<T> r = exact_mult(a, b.hi);
172 r.lo = multiply_add(a, b.lo, r.lo);
173 return r;
174}
175
176template <size_t SPLIT_B = 27>
177LIBC_INLINE constexpr DoubleDouble quick_mult(const DoubleDouble &a,
178 const DoubleDouble &b) {
179 DoubleDouble r = exact_mult<double, SPLIT_B>(a.hi, b.hi);
180 double t1 = multiply_add(x: a.hi, y: b.lo, z: r.lo);
181 double t2 = multiply_add(x: a.lo, y: b.hi, z: t1);
182 r.lo = t2;
183 return r;
184}
185
186// Assuming |c| >= |a * b|.
187template <>
188LIBC_INLINE DoubleDouble multiply_add<DoubleDouble>(const DoubleDouble &a,
189 const DoubleDouble &b,
190 const DoubleDouble &c) {
191 return add(a: c, b: quick_mult(a, b));
192}
193
194// Accurate double-double division, following Karp-Markstein's trick for
195// division, implemented in the CORE-MATH project at:
196// https://gitlab.inria.fr/core-math/core-math/-/blob/master/src/binary64/tan/tan.c#L1855
197//
198// Error bounds:
199// Let a = ah + al, b = bh + bl.
200// Let r = rh + rl be the approximation of (ah + al) / (bh + bl).
201// Then:
202// (ah + al) / (bh + bl) - rh =
203// = ((ah - bh * rh) + (al - bl * rh)) / (bh + bl)
204// = (1 + O(bl/bh)) * ((ah - bh * rh) + (al - bl * rh)) / bh
205// Let q = round(1/bh), then the above expressions are approximately:
206// = (1 + O(bl / bh)) * (1 + O(2^-52)) * q * ((ah - bh * rh) + (al - bl * rh))
207// So we can compute:
208// rl = q * (ah - bh * rh) + q * (al - bl * rh)
209// as accurate as possible, then the error is bounded by:
210// |(ah + al) / (bh + bl) - (rh + rl)| < O(bl/bh) * (2^-52 + al/ah + bl/bh)
211template <typename T>
212LIBC_INLINE NumberPair<T> div(const NumberPair<T> &a, const NumberPair<T> &b) {
213 NumberPair<T> r;
214 T q = T(1) / b.hi;
215 r.hi = a.hi * q;
216
217#ifdef LIBC_TARGET_CPU_HAS_FMA
218 T e_hi = fputil::multiply_add(b.hi, -r.hi, a.hi);
219 T e_lo = fputil::multiply_add(b.lo, -r.hi, a.lo);
220#else
221 NumberPair<T> b_hi_r_hi = fputil::exact_mult(b.hi, -r.hi);
222 NumberPair<T> b_lo_r_hi = fputil::exact_mult(b.lo, -r.hi);
223 T e_hi = (a.hi + b_hi_r_hi.hi) + b_hi_r_hi.lo;
224 T e_lo = (a.lo + b_lo_r_hi.hi) + b_lo_r_hi.lo;
225#endif // LIBC_TARGET_CPU_HAS_FMA
226
227 r.lo = q * (e_hi + e_lo);
228 return r;
229}
230
231} // namespace fputil
232} // namespace LIBC_NAMESPACE_DECL
233
234#endif // LLVM_LIBC_SRC___SUPPORT_FPUTIL_DOUBLE_DOUBLE_H
235