1 //===-- Floating point divsion and remainder operations ---------*- 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_DIVISIONANDREMAINDEROPERATIONS_H
10 #define LLVM_LIBC_SRC___SUPPORT_FPUTIL_DIVISIONANDREMAINDEROPERATIONS_H
11
12 #include "FPBits.h"
13 #include "ManipulationFunctions.h"
14 #include "NormalFloat.h"
15
16 #include "src/__support/CPP/type_traits.h"
17 #include "src/__support/common.h"
18 #include "src/__support/macros/config.h"
19
20 namespace LIBC_NAMESPACE_DECL {
21 namespace fputil {
22
23 static constexpr int QUOTIENT_LSB_BITS = 3;
24
25 // The implementation is a bit-by-bit algorithm which uses integer division
26 // to evaluate the quotient and remainder.
27 template <typename T, cpp::enable_if_t<cpp::is_floating_point_v<T>, int> = 0>
remquo(T x,T y,int & q)28 LIBC_INLINE T remquo(T x, T y, int &q) {
29 FPBits<T> xbits(x), ybits(y);
30 if (xbits.is_nan())
31 return x;
32 if (ybits.is_nan())
33 return y;
34 if (xbits.is_inf() || ybits.is_zero())
35 return FPBits<T>::quiet_nan().get_val();
36
37 if (xbits.is_zero()) {
38 q = 0;
39 return LIBC_NAMESPACE::fputil::copysign(T(0.0), x);
40 }
41
42 if (ybits.is_inf()) {
43 q = 0;
44 return x;
45 }
46
47 const Sign result_sign =
48 (xbits.sign() == ybits.sign() ? Sign::POS : Sign::NEG);
49
50 // Once we know the sign of the result, we can just operate on the absolute
51 // values. The correct sign can be applied to the result after the result
52 // is evaluated.
53 xbits.set_sign(Sign::POS);
54 ybits.set_sign(Sign::POS);
55
56 NormalFloat<T> normalx(xbits), normaly(ybits);
57 int exp = normalx.exponent - normaly.exponent;
58 typename NormalFloat<T>::StorageType mx = normalx.mantissa,
59 my = normaly.mantissa;
60
61 q = 0;
62 while (exp >= 0) {
63 unsigned shift_count = 0;
64 typename NormalFloat<T>::StorageType n = mx;
65 for (shift_count = 0; n < my; n <<= 1, ++shift_count)
66 ;
67
68 if (static_cast<int>(shift_count) > exp)
69 break;
70
71 exp -= shift_count;
72 if (0 <= exp && exp < QUOTIENT_LSB_BITS)
73 q |= (1 << exp);
74
75 mx = n - my;
76 if (mx == 0) {
77 q = result_sign.is_neg() ? -q : q;
78 return LIBC_NAMESPACE::fputil::copysign(T(0.0), x);
79 }
80 }
81
82 NormalFloat<T> remainder(Sign::POS, exp + normaly.exponent, mx);
83
84 // Since NormalFloat to native type conversion is a truncation operation
85 // currently, the remainder value in the native type is correct as is.
86 // However, if NormalFloat to native type conversion is updated in future,
87 // then the conversion to native remainder value should be updated
88 // appropriately and some directed tests added.
89 T native_remainder(remainder);
90 T absy = ybits.get_val();
91 int cmp = remainder.mul2(1).cmp(normaly);
92 if (cmp > 0) {
93 q = q + 1;
94 if (x >= T(0.0))
95 native_remainder = native_remainder - absy;
96 else
97 native_remainder = absy - native_remainder;
98 } else if (cmp == 0) {
99 if (q & 1) {
100 q += 1;
101 if (x >= T(0.0))
102 native_remainder = -native_remainder;
103 } else {
104 if (x < T(0.0))
105 native_remainder = -native_remainder;
106 }
107 } else {
108 if (x < T(0.0))
109 native_remainder = -native_remainder;
110 }
111
112 q = result_sign.is_neg() ? -q : q;
113 if (native_remainder == T(0.0))
114 return LIBC_NAMESPACE::fputil::copysign(T(0.0), x);
115 return native_remainder;
116 }
117
118 } // namespace fputil
119 } // namespace LIBC_NAMESPACE_DECL
120
121 #endif // LLVM_LIBC_SRC___SUPPORT_FPUTIL_DIVISIONANDREMAINDEROPERATIONS_H
122