1//===-- A class to store high precision floating point numbers --*- 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_DYADIC_FLOAT_H
10#define LLVM_LIBC_SRC___SUPPORT_FPUTIL_DYADIC_FLOAT_H
11
12#include "FEnvImpl.h"
13#include "FPBits.h"
14#include "hdr/errno_macros.h"
15#include "hdr/fenv_macros.h"
16#include "multiply_add.h"
17#include "rounding_mode.h"
18#include "src/__support/CPP/type_traits.h"
19#include "src/__support/big_int.h"
20#include "src/__support/macros/attributes.h"
21#include "src/__support/macros/config.h"
22#include "src/__support/macros/optimization.h" // LIBC_UNLIKELY
23#include "src/__support/macros/properties/types.h"
24
25#include <stddef.h>
26
27namespace LIBC_NAMESPACE_DECL {
28namespace fputil {
29
30// Decide whether to round a UInt up, down or not at all at a given bit
31// position, based on the current rounding mode. The assumption is that the
32// caller is going to make the integer `value >> rshift`, and then might need
33// to round it up by 1 depending on the value of the bits shifted off the
34// bottom.
35//
36// `logical_sign` causes the behavior of FE_DOWNWARD and FE_UPWARD to
37// be reversed, which is what you'd want if this is the mantissa of a
38// negative floating-point number.
39//
40// Return value is +1 if the value should be rounded up; -1 if it should be
41// rounded down; 0 if it's exact and needs no rounding.
42template <size_t Bits>
43LIBC_INLINE LIBC_CONSTEXPR_DEFAULT int
44rounding_direction(const LIBC_NAMESPACE::UInt<Bits> &value, size_t rshift,
45 [[maybe_unused]] Sign logical_sign) {
46 // logical_sign only affects FE_DOWNWARD and FE_UPWARD rounding modes. In the
47 // case of LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY being enabled, this option
48 // is no-op
49
50 if (rshift == 0 || (rshift < Bits && (value << (Bits - rshift)) == 0) ||
51 (rshift >= Bits && value == 0))
52 return 0; // exact
53
54#ifdef LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
55 if (rshift > 0 && rshift <= Bits && value.get_bit(rshift - 1)) {
56 // We round up, unless the value is an exact halfway case and
57 // the bit that will end up in the units place is 0, in which
58 // case tie-break-to-even says round down.
59 bool round_bit = rshift < Bits ? value.get_bit(rshift) : 0;
60 return round_bit != 0 || (value << (Bits - rshift + 1)) != 0 ? +1 : -1;
61 } else {
62 return -1;
63 }
64#else // !LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
65 switch (quick_get_round()) {
66 case FE_TONEAREST:
67 if (rshift > 0 && rshift <= Bits && value.get_bit(rshift - 1)) {
68 // We round up, unless the value is an exact halfway case and
69 // the bit that will end up in the units place is 0, in which
70 // case tie-break-to-even says round down.
71 bool round_bit = rshift < Bits ? value.get_bit(rshift) : 0;
72 return round_bit != 0 || (value << (Bits - rshift + 1)) != 0 ? +1 : -1;
73 } else {
74 return -1;
75 }
76 case FE_TOWARDZERO:
77 return -1;
78 case FE_DOWNWARD:
79 return logical_sign.is_neg() &&
80 (rshift < Bits && (value << (Bits - rshift)) != 0)
81 ? +1
82 : -1;
83 case FE_UPWARD:
84 return logical_sign.is_pos() &&
85 (rshift < Bits && (value << (Bits - rshift)) != 0)
86 ? +1
87 : -1;
88 default:
89 __builtin_unreachable();
90 }
91#endif // LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
92}
93
94// A generic class to perform computations of high precision floating points.
95// We store the value in dyadic format, including 3 fields:
96// sign : boolean value - false means positive, true means negative
97// exponent: the exponent value of the least significant bit of the mantissa.
98// mantissa: unsigned integer of length `Bits`.
99// So the real value that is stored is:
100// real value = (-1)^sign * 2^exponent * (mantissa as unsigned integer)
101// The stored data is normal if for non-zero mantissa, the leading bit is 1.
102// The outputs of the constructors and most functions will be normalized.
103// To simplify and improve the efficiency, many functions will assume that the
104// inputs are normal.
105template <size_t Bits> struct DyadicFloat {
106 using MantissaType = LIBC_NAMESPACE::UInt<Bits>;
107
108 Sign sign = Sign::POS;
109 int exponent = 0;
110 MantissaType mantissa = MantissaType(0);
111
112 LIBC_INLINE constexpr DyadicFloat() = default;
113
114 template <typename T, cpp::enable_if_t<cpp::is_floating_point_v<T>, int> = 0>
115 LIBC_INLINE LIBC_BIT_CAST_CONSTEXPR DyadicFloat(T x) {
116 static_assert(FPBits<T>::FRACTION_LEN < Bits);
117 FPBits<T> x_bits(x);
118 sign = x_bits.sign();
119 exponent = x_bits.get_explicit_exponent() - FPBits<T>::FRACTION_LEN;
120 mantissa = MantissaType(x_bits.get_explicit_mantissa());
121 normalize();
122 }
123
124 LIBC_INLINE constexpr DyadicFloat(Sign s, int e, const MantissaType &m)
125 : sign(s), exponent(e), mantissa(m) {
126 normalize();
127 }
128
129 // Normalizing the mantissa, bringing the leading 1 bit to the most
130 // significant bit.
131 LIBC_INLINE constexpr DyadicFloat &normalize() {
132 if (!mantissa.is_zero()) {
133 int shift_length = cpp::countl_zero(mantissa);
134 exponent -= shift_length;
135 mantissa <<= static_cast<size_t>(shift_length);
136 }
137 return *this;
138 }
139
140 // Used for aligning exponents. Output might not be normalized.
141 LIBC_INLINE constexpr DyadicFloat &shift_left(unsigned shift_length) {
142 exponent -= static_cast<int>(shift_length);
143 if (shift_length < Bits)
144 mantissa <<= shift_length;
145 else
146 mantissa = MantissaType(0);
147 return *this;
148 }
149
150 // Used for aligning exponents. Output might not be normalized.
151 LIBC_INLINE constexpr DyadicFloat &shift_right(unsigned shift_length) {
152 exponent += static_cast<int>(shift_length);
153 if (shift_length < Bits)
154 mantissa >>= shift_length;
155 else
156 mantissa = MantissaType(0);
157 return *this;
158 }
159
160 // Assume that it is already normalized. Output the unbiased exponent.
161 LIBC_INLINE constexpr int get_unbiased_exponent() const {
162 return exponent + (Bits - 1);
163 }
164
165 // Produce a correctly rounded DyadicFloat from a too-large mantissa,
166 // by shifting it down and rounding if necessary.
167 template <size_t MantissaBits>
168 LIBC_INLINE LIBC_CONSTEXPR_DEFAULT static DyadicFloat<Bits>
169 round(Sign result_sign, int result_exponent,
170 const LIBC_NAMESPACE::UInt<MantissaBits> &input_mantissa,
171 size_t rshift) {
172 MantissaType result_mantissa(input_mantissa >> rshift);
173 if (rounding_direction(input_mantissa, rshift, result_sign) > 0) {
174 ++result_mantissa;
175 if (result_mantissa == 0) {
176 // Rounding up made the mantissa integer wrap round to 0,
177 // carrying a bit off the top. So we've rounded up to the next
178 // exponent.
179 result_mantissa.set_bit(Bits - 1);
180 ++result_exponent;
181 }
182 }
183 return DyadicFloat(result_sign, result_exponent, result_mantissa);
184 }
185
186 template <typename T, bool ShouldSignalExceptions = true>
187 LIBC_INLINE LIBC_CONSTEXPR_DEFAULT cpp::enable_if_t<
188 cpp::is_floating_point_v<T> && (FPBits<T>::FRACTION_LEN < Bits), T>
189 generic_as() const {
190 using FPBits = FPBits<T>;
191 using StorageType = typename FPBits::StorageType;
192
193 constexpr int EXTRA_FRACTION_LEN = Bits - 1 - FPBits::FRACTION_LEN;
194
195 if (mantissa == 0)
196 return FPBits::zero(sign).get_val();
197
198 int unbiased_exp = get_unbiased_exponent();
199
200 if (unbiased_exp + FPBits::EXP_BIAS >= FPBits::MAX_BIASED_EXPONENT) {
201 if constexpr (ShouldSignalExceptions) {
202 set_errno_if_required(ERANGE);
203 raise_overflow_except_if_required<T>();
204 }
205
206#ifdef LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
207 return FPBits::inf(sign).get_val();
208#else // !LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
209 switch (quick_get_round()) {
210 case FE_TONEAREST:
211 return FPBits::inf(sign).get_val();
212 case FE_TOWARDZERO:
213 return FPBits::max_normal(sign).get_val();
214 case FE_DOWNWARD:
215 if (sign.is_pos())
216 return FPBits::max_normal(Sign::POS).get_val();
217 return FPBits::inf(Sign::NEG).get_val();
218 case FE_UPWARD:
219 if (sign.is_neg())
220 return FPBits::max_normal(Sign::NEG).get_val();
221 return FPBits::inf(Sign::POS).get_val();
222 default:
223 __builtin_unreachable();
224 }
225#endif // LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
226 }
227
228 StorageType out_biased_exp = 0;
229 StorageType out_mantissa = 0;
230 bool round = false;
231 bool sticky = false;
232 bool underflow = false;
233
234 if (unbiased_exp < -FPBits::EXP_BIAS - FPBits::FRACTION_LEN) {
235 sticky = true;
236 underflow = true;
237 } else if (unbiased_exp == -FPBits::EXP_BIAS - FPBits::FRACTION_LEN) {
238 round = true;
239 // underflow is detected pre-rounding FE_UNDERFLOW may be raised
240 // even if rounding produces a non-underflow result
241 underflow = true;
242 MantissaType sticky_mask = (MantissaType(1) << (Bits - 1)) - 1;
243 sticky = (mantissa & sticky_mask) != 0;
244 } else {
245 int extra_fraction_len = EXTRA_FRACTION_LEN;
246
247 if (unbiased_exp < 1 - FPBits::EXP_BIAS) {
248 underflow = true;
249 extra_fraction_len += 1 - FPBits::EXP_BIAS - unbiased_exp;
250 } else {
251 out_biased_exp =
252 static_cast<StorageType>(unbiased_exp + FPBits::EXP_BIAS);
253 }
254
255 if (extra_fraction_len > 0) {
256 MantissaType round_mask = MantissaType(1) << (extra_fraction_len - 1);
257 round = (mantissa & round_mask) != 0;
258 MantissaType sticky_mask = round_mask - 1;
259 sticky = (mantissa & sticky_mask) != 0;
260 }
261
262 out_mantissa = static_cast<StorageType>(mantissa >> extra_fraction_len);
263 }
264
265 bool lsb = (out_mantissa & 1) != 0;
266
267 StorageType result =
268 FPBits::create_value(sign, out_biased_exp, out_mantissa).uintval();
269
270#ifdef LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
271 if (round && (lsb || sticky))
272 ++result;
273#else
274 switch (quick_get_round()) {
275 case FE_TONEAREST:
276 if (round && (lsb || sticky))
277 ++result;
278 break;
279 case FE_DOWNWARD:
280 if (sign.is_neg() && (round || sticky))
281 ++result;
282 break;
283 case FE_UPWARD:
284 if (sign.is_pos() && (round || sticky))
285 ++result;
286 break;
287 default:
288 break;
289 }
290#endif // LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
291
292 if (ShouldSignalExceptions && (round || sticky)) {
293 if (FPBits(result).is_inf()) {
294 set_errno_if_required(ERANGE);
295 raise_overflow_except_if_required<T>();
296 } else if (underflow) {
297 set_errno_if_required(ERANGE);
298 raise_underflow_except_if_required<T>();
299 } else {
300 raise_except_if_required(FE_INEXACT);
301 }
302 }
303
304 return FPBits(result).get_val();
305 }
306
307 template <typename T, bool ShouldSignalExceptions = true,
308 typename = cpp::enable_if_t<cpp::is_floating_point_v<T> &&
309 (FPBits<T>::FRACTION_LEN < Bits),
310 void>>
311 LIBC_INLINE LIBC_CONSTEXPR_DEFAULT T fast_as() const {
312 if (LIBC_UNLIKELY(mantissa.is_zero()))
313 return FPBits<T>::zero(sign).get_val();
314
315 // Assume that it is normalized, and output is also normal.
316 constexpr uint32_t PRECISION = FPBits<T>::FRACTION_LEN + 1;
317 using output_bits_t = typename FPBits<T>::StorageType;
318 constexpr output_bits_t IMPLICIT_MASK =
319 FPBits<T>::SIG_MASK - FPBits<T>::FRACTION_MASK;
320
321 int exp_hi = exponent + static_cast<int>((Bits - 1) + FPBits<T>::EXP_BIAS);
322
323 if (LIBC_UNLIKELY(exp_hi > 2 * FPBits<T>::EXP_BIAS)) {
324 // Results overflow.
325 T d_hi =
326 FPBits<T>::create_value(sign, 2 * FPBits<T>::EXP_BIAS, IMPLICIT_MASK)
327 .get_val();
328 // volatile prevents constant propagation that would result in infinity
329 // always being returned no matter the current rounding mode.
330 volatile T two = static_cast<T>(2.0);
331 T r = two * d_hi;
332
333 // TODO: Whether rounding down the absolute value to max_normal should
334 // also raise FE_OVERFLOW and set ERANGE is debatable.
335 if (ShouldSignalExceptions && FPBits<T>(r).is_inf())
336 set_errno_if_required(ERANGE);
337
338 return r;
339 }
340
341 bool denorm = false;
342 uint32_t shift = Bits - PRECISION;
343 if (LIBC_UNLIKELY(exp_hi <= 0)) {
344 // Output is denormal.
345 denorm = true;
346 shift = (Bits - PRECISION) + static_cast<uint32_t>(1 - exp_hi);
347
348 exp_hi = FPBits<T>::EXP_BIAS;
349 }
350
351 int exp_lo = exp_hi - static_cast<int>(PRECISION) - 1;
352
353 MantissaType m_hi =
354 shift >= MantissaType::BITS ? MantissaType(0) : mantissa >> shift;
355
356 T d_hi = FPBits<T>::create_value(
357 sign, static_cast<output_bits_t>(exp_hi),
358 (static_cast<output_bits_t>(m_hi) & FPBits<T>::SIG_MASK) |
359 IMPLICIT_MASK)
360 .get_val();
361
362 MantissaType round_mask =
363 shift - 1 >= MantissaType::BITS ? 0 : MantissaType(1) << (shift - 1);
364 MantissaType sticky_mask = round_mask - MantissaType(1);
365
366 bool round_bit = !(mantissa & round_mask).is_zero();
367 bool sticky_bit = !(mantissa & sticky_mask).is_zero();
368 int round_and_sticky = int(round_bit) * 2 + int(sticky_bit);
369
370 T d_lo;
371
372 if (LIBC_UNLIKELY(exp_lo <= 0)) {
373 // d_lo is denormal, but the output is normal.
374 int scale_up_exponent = 1 - exp_lo;
375 T scale_up_factor =
376 FPBits<T>::create_value(Sign::POS,
377 static_cast<output_bits_t>(
378 FPBits<T>::EXP_BIAS + scale_up_exponent),
379 IMPLICIT_MASK)
380 .get_val();
381 T scale_down_factor =
382 FPBits<T>::create_value(Sign::POS,
383 static_cast<output_bits_t>(
384 FPBits<T>::EXP_BIAS - scale_up_exponent),
385 IMPLICIT_MASK)
386 .get_val();
387
388 d_lo = FPBits<T>::create_value(
389 sign, static_cast<output_bits_t>(exp_lo + scale_up_exponent),
390 IMPLICIT_MASK)
391 .get_val();
392
393 return multiply_add(d_lo, T(round_and_sticky), d_hi * scale_up_factor) *
394 scale_down_factor;
395 }
396
397 d_lo = FPBits<T>::create_value(sign, static_cast<output_bits_t>(exp_lo),
398 IMPLICIT_MASK)
399 .get_val();
400
401 // Still correct without FMA instructions if `d_lo` is not underflow.
402 T r = multiply_add(d_lo, T(round_and_sticky), d_hi);
403
404 if (LIBC_UNLIKELY(denorm)) {
405 // Exponent before rounding is in denormal range, simply clear the
406 // exponent field.
407 output_bits_t clear_exp = static_cast<output_bits_t>(
408 output_bits_t(exp_hi) << FPBits<T>::SIG_LEN);
409 output_bits_t r_bits = FPBits<T>(r).uintval() - clear_exp;
410
411 if (!(r_bits & FPBits<T>::EXP_MASK)) {
412 // Output is denormal after rounding, clear the implicit bit for
413 // 80-bit long double.
414 r_bits -= IMPLICIT_MASK;
415 }
416
417 // Underflow exception and ERANGE are signaled when an unrounded result
418 // in the denormal range is inexact, even if destination rounding rounds
419 // it up to min_normal.
420 if (ShouldSignalExceptions && round_and_sticky) {
421 set_errno_if_required(ERANGE);
422 raise_underflow_except_if_required<T>();
423 }
424
425 return FPBits<T>(r_bits).get_val();
426 }
427
428 return r;
429 }
430
431 // Assume that it is already normalized.
432 // Output is rounded correctly with respect to the current rounding mode.
433 template <typename T, bool ShouldSignalExceptions = true,
434 typename = cpp::enable_if_t<cpp::is_floating_point_v<T> &&
435 (FPBits<T>::FRACTION_LEN < Bits),
436 void>>
437 LIBC_INLINE LIBC_CONSTEXPR_DEFAULT T as() const {
438 if constexpr (cpp::is_same_v<T, bfloat16> || cpp::is_same_v<T, Float128> ||
439 cpp::is_same_v<T, Float80>
440#if defined(LIBC_TYPES_HAS_FLOAT16) && !defined(__LIBC_USE_FLOAT16_CONVERSION)
441 || cpp::is_same_v<T, float16>
442#endif
443#if defined(LIBC_TYPES_HAS_NATIVE_FLOAT128)
444 || cpp::is_same_v<T, float128>
445#endif
446 )
447 return generic_as<T, ShouldSignalExceptions>();
448 else
449 return fast_as<T, ShouldSignalExceptions>();
450 }
451
452 template <typename T,
453 typename = cpp::enable_if_t<cpp::is_floating_point_v<T> &&
454 (FPBits<T>::FRACTION_LEN < Bits),
455 void>>
456 LIBC_INLINE explicit constexpr operator T() const {
457 return as<T, /*ShouldSignalExceptions=*/true>();
458 }
459
460 LIBC_INLINE constexpr MantissaType as_mantissa_type() const {
461 if (mantissa.is_zero())
462 return 0;
463
464 MantissaType new_mant = mantissa;
465 if (exponent > 0) {
466 new_mant <<= exponent;
467 } else {
468 // Cast the exponent to size_t before negating it, rather than after,
469 // to avoid undefined behavior negating INT_MIN as an integer (although
470 // exponents coming in to this function _shouldn't_ be that large). The
471 // result should always end up as a positive size_t.
472 size_t shift = -static_cast<size_t>(exponent);
473 size_t limit = cpp::numeric_limits<MantissaType>::digits;
474 if (shift >= limit)
475 new_mant = 0;
476 else
477 new_mant >>= shift;
478 }
479
480 if (sign.is_neg()) {
481 new_mant = (~new_mant) + 1;
482 }
483
484 return new_mant;
485 }
486
487 LIBC_INLINE LIBC_CONSTEXPR_DEFAULT MantissaType
488 as_mantissa_type_rounded(int *round_dir_out = nullptr) const {
489 int round_dir = 0;
490 MantissaType new_mant;
491 if (mantissa.is_zero()) {
492 new_mant = 0;
493 } else {
494 new_mant = mantissa;
495 if (exponent > 0) {
496 new_mant <<= exponent;
497 } else if (exponent < 0) {
498 // Cast the exponent to size_t before negating it, rather than after,
499 // to avoid undefined behavior negating INT_MIN as an integer
500 // (although exponents coming in to this function _shouldn't_ be that
501 // large). The result should always end up as a positive size_t.
502 size_t shift = -static_cast<size_t>(exponent);
503 if (shift >= Bits)
504 new_mant = 0;
505 else
506 new_mant >>= shift;
507 round_dir = rounding_direction(mantissa, shift, sign);
508 if (round_dir > 0)
509 ++new_mant;
510 }
511
512 if (sign.is_neg()) {
513 new_mant = (~new_mant) + 1;
514 }
515 }
516
517 if (round_dir_out)
518 *round_dir_out = round_dir;
519
520 return new_mant;
521 }
522
523 LIBC_INLINE constexpr DyadicFloat operator-() const {
524 return DyadicFloat(sign.negate(), exponent, mantissa);
525 }
526};
527
528// Quick add - Add 2 dyadic floats with rounding toward 0 and then normalize
529// the output:
530// - Align the exponents so that:
531// new a.exponent = new b.exponent = max(a.exponent, b.exponent)
532// - Add or subtract the mantissas depending on the signs.
533// - Normalize the result.
534// The absolute errors compared to the mathematical sum is bounded by:
535// | quick_add(a, b) - (a + b) | < MSB(a + b) * 2^(-Bits + 2),
536// i.e., errors are up to 2 ULPs.
537// Assume inputs are normalized (by constructors or other functions) so that
538// we don't need to normalize the inputs again in this function. If the
539// inputs are not normalized, the results might lose precision significantly.
540template <size_t Bits>
541LIBC_INLINE constexpr DyadicFloat<Bits> quick_add(DyadicFloat<Bits> a,
542 DyadicFloat<Bits> b) {
543 if (LIBC_UNLIKELY(a.mantissa.is_zero()))
544 return b;
545 if (LIBC_UNLIKELY(b.mantissa.is_zero()))
546 return a;
547
548 // Align exponents
549 if (a.exponent > b.exponent)
550 b.shift_right(static_cast<unsigned>(a.exponent - b.exponent));
551 else if (b.exponent > a.exponent)
552 a.shift_right(static_cast<unsigned>(b.exponent - a.exponent));
553
554 DyadicFloat<Bits> result;
555
556 if (a.sign == b.sign) {
557 // Addition
558 result.sign = a.sign;
559 result.exponent = a.exponent;
560 result.mantissa = a.mantissa;
561 if (result.mantissa.add_overflow(b.mantissa)) {
562 // Mantissa addition overflow.
563 result.shift_right(1);
564 result.mantissa.val[DyadicFloat<Bits>::MantissaType::WORD_COUNT - 1] |=
565 (uint64_t(1) << 63);
566 }
567 // Result is already normalized.
568 return result;
569 }
570
571 // Subtraction
572 if (a.mantissa >= b.mantissa) {
573 result.sign = a.sign;
574 result.exponent = a.exponent;
575 result.mantissa = a.mantissa - b.mantissa;
576 } else {
577 result.sign = b.sign;
578 result.exponent = b.exponent;
579 result.mantissa = b.mantissa - a.mantissa;
580 }
581
582 return result.normalize();
583}
584
585template <size_t Bits>
586LIBC_INLINE constexpr DyadicFloat<Bits> quick_sub(DyadicFloat<Bits> a,
587 DyadicFloat<Bits> b) {
588 return quick_add(a, -b);
589}
590
591// Quick Mul - Slightly less accurate but efficient multiplication of 2 dyadic
592// floats with rounding toward 0 and then normalize the output:
593// result.exponent = a.exponent + b.exponent + Bits,
594// result.mantissa = quick_mul_hi(a.mantissa + b.mantissa)
595// ~ (full product a.mantissa * b.mantissa) >> Bits.
596// The errors compared to the mathematical product is bounded by:
597// 2 * errors of quick_mul_hi = 2 * (UInt<Bits>::WORD_COUNT - 1) in ULPs.
598// Assume inputs are normalized (by constructors or other functions) so that
599// we don't need to normalize the inputs again in this function. If the
600// inputs are not normalized, the results might lose precision significantly.
601template <size_t Bits>
602LIBC_INLINE constexpr DyadicFloat<Bits> quick_mul(const DyadicFloat<Bits> &a,
603 const DyadicFloat<Bits> &b) {
604 DyadicFloat<Bits> result;
605 result.sign = (a.sign != b.sign) ? Sign::NEG : Sign::POS;
606 result.exponent = a.exponent + b.exponent + static_cast<int>(Bits);
607
608 if (!(a.mantissa.is_zero() || b.mantissa.is_zero())) {
609 result.mantissa = a.mantissa.quick_mul_hi(b.mantissa);
610 // Check the leading bit directly, should be faster than using clz in
611 // normalize().
612 if (result.mantissa.val[DyadicFloat<Bits>::MantissaType::WORD_COUNT - 1] >>
613 (DyadicFloat<Bits>::MantissaType::WORD_SIZE - 1) ==
614 0)
615 result.shift_left(1);
616 } else {
617 result.mantissa = (typename DyadicFloat<Bits>::MantissaType)(0);
618 }
619 return result;
620}
621
622// Correctly rounded multiplication of 2 dyadic floats, assuming the
623// exponent remains within range.
624template <size_t Bits>
625LIBC_INLINE constexpr DyadicFloat<Bits>
626rounded_mul(const DyadicFloat<Bits> &a, const DyadicFloat<Bits> &b) {
627 using DblMant = LIBC_NAMESPACE::UInt<(2 * Bits)>;
628 Sign result_sign = (a.sign != b.sign) ? Sign::NEG : Sign::POS;
629 int result_exponent = a.exponent + b.exponent + static_cast<int>(Bits);
630 auto product = DblMant(a.mantissa) * DblMant(b.mantissa);
631 // As in quick_mul(), renormalize by 1 bit manually rather than countl_zero
632 if (product.get_bit(2 * Bits - 1) == 0) {
633 product <<= 1;
634 result_exponent -= 1;
635 }
636
637 return DyadicFloat<Bits>::round(result_sign, result_exponent, product, Bits);
638}
639
640// Approximate reciprocal - given a nonzero a, make a good approximation to
641// 1/a. The method is Newton-Raphson iteration, based on quick_mul.
642template <size_t Bits, typename = cpp::enable_if_t<(Bits >= 32)>>
643LIBC_INLINE constexpr DyadicFloat<Bits>
644approx_reciprocal(const DyadicFloat<Bits> &a) {
645 // Given an approximation x to 1/a, a better one is x' = x(2-ax).
646 //
647 // You can derive this by using the Newton-Raphson formula with the function
648 // f(x) = 1/x - a. But another way to see that it works is to say: suppose
649 // that ax = 1-e for some small error e. Then ax' = ax(2-ax) = (1-e)(1+e) =
650 // 1-e^2. So the error in x' is the square of the error in x, i.e. the
651 // number of correct bits in x' is double the number in x.
652
653 // An initial approximation to the reciprocal
654 DyadicFloat<Bits> x(Sign::POS, -32 - a.exponent - int(Bits),
655 uint64_t(0xFFFFFFFFFFFFFFFF) /
656 static_cast<uint64_t>(a.mantissa >> (Bits - 32)));
657
658 // The constant 2, which we'll need in every iteration
659 DyadicFloat<Bits> two(Sign::POS, 1, 1);
660
661 // We expect at least 31 correct bits from our 32-bit starting approximation
662 size_t ok_bits = 31;
663
664 // The number of good bits doubles in each iteration, except that rounding
665 // errors introduce a little extra each time. Subtract a bit from our
666 // accuracy assessment to account for that.
667 while (ok_bits < Bits) {
668 x = quick_mul(x, quick_sub(two, quick_mul(a, x)));
669 ok_bits = 2 * ok_bits - 1;
670 }
671
672 return x;
673}
674
675// Correctly rounded division of 2 dyadic floats, assuming the
676// exponent remains within range.
677template <size_t Bits>
678LIBC_INLINE constexpr DyadicFloat<Bits>
679rounded_div(const DyadicFloat<Bits> &af, const DyadicFloat<Bits> &bf) {
680 using DblMant = LIBC_NAMESPACE::UInt<(Bits * 2 + 64)>;
681
682 // Make an approximation to the quotient as a * (1/b). Both the
683 // multiplication and the reciprocal are a bit sloppy, which doesn't
684 // matter, because we're going to correct for that below.
685 auto qf = fputil::quick_mul(af, fputil::approx_reciprocal(bf));
686
687 // Switch to BigInt and stop using quick_add and quick_mul: now
688 // we're working in exact integers so as to get the true remainder.
689 DblMant a = af.mantissa, b = bf.mantissa, q = qf.mantissa;
690 q <<= 2; // leave room for a round bit, even if exponent decreases
691 a <<= af.exponent - bf.exponent - qf.exponent + 2;
692 DblMant qb = q * b;
693 if (qb < a) {
694 DblMant too_small = a - b;
695 while (qb <= too_small) {
696 qb += b;
697 ++q;
698 }
699 } else {
700 while (qb > a) {
701 qb -= b;
702 --q;
703 }
704 }
705
706 DyadicFloat<(Bits * 2)> qbig(qf.sign, qf.exponent - 2, q);
707 return DyadicFloat<Bits>::round(qbig.sign, qbig.exponent + Bits,
708 qbig.mantissa, Bits);
709}
710
711// Simple polynomial approximation.
712template <size_t Bits>
713LIBC_INLINE constexpr DyadicFloat<Bits>
714multiply_add(const DyadicFloat<Bits> &a, const DyadicFloat<Bits> &b,
715 const DyadicFloat<Bits> &c) {
716 return quick_add(c, quick_mul(a, b));
717}
718
719// Simple exponentiation implementation for printf. Only handles positive
720// exponents, since division isn't implemented.
721template <size_t Bits>
722LIBC_INLINE constexpr DyadicFloat<Bits> pow_n(const DyadicFloat<Bits> &a,
723 uint32_t power) {
724 DyadicFloat<Bits> result = 1.0;
725 DyadicFloat<Bits> cur_power = a;
726
727 while (power > 0) {
728 if ((power % 2) > 0) {
729 result = quick_mul(result, cur_power);
730 }
731 power = power >> 1;
732 cur_power = quick_mul(cur_power, cur_power);
733 }
734 return result;
735}
736
737template <size_t Bits>
738LIBC_INLINE constexpr DyadicFloat<Bits> mul_pow_2(const DyadicFloat<Bits> &a,
739 int32_t pow_2) {
740 DyadicFloat<Bits> result = a;
741 result.exponent += pow_2;
742 return result;
743}
744
745} // namespace fputil
746} // namespace LIBC_NAMESPACE_DECL
747
748#endif // LLVM_LIBC_SRC___SUPPORT_FPUTIL_DYADIC_FLOAT_H
749