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