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