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_DIVISION_AND_REMAINDER_OPERATIONS_H
10 #define LLVM_LIBC_SRC_SUPPORT_FPUTIL_DIVISION_AND_REMAINDER_OPERATIONS_H
11 
12 #include "FPBits.h"
13 #include "ManipulationFunctions.h"
14 #include "NormalFloat.h"
15 
16 #include "src/__support/CPP/TypeTraits.h"
17 
18 namespace __llvm_libc {
19 namespace fputil {
20 
21 static constexpr int QUOTIENT_LSB_BITS = 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>
27 static inline T remquo(T x, T y, int &q) {
28   FPBits<T> xbits(x), ybits(y);
29   if (xbits.is_nan())
30     return x;
31   if (ybits.is_nan())
32     return y;
33   if (xbits.is_inf() || ybits.is_zero())
34     return FPBits<T>::build_nan(1);
35 
36   if (xbits.is_zero()) {
37     q = 0;
38     return __llvm_libc::fputil::copysign(T(0.0), x);
39   }
40 
41   if (ybits.is_inf()) {
42     q = 0;
43     return x;
44   }
45 
46   bool result_sign = (xbits.get_sign() == ybits.get_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.set_sign(0);
52   ybits.set_sign(0);
53 
54   NormalFloat<T> normalx(xbits), normaly(ybits);
55   int exp = normalx.exponent - normaly.exponent;
56   typename NormalFloat<T>::UIntType mx = normalx.mantissa,
57                                     my = normaly.mantissa;
58 
59   q = 0;
60   while (exp >= 0) {
61     unsigned shift_count = 0;
62     typename NormalFloat<T>::UIntType n = mx;
63     for (shift_count = 0; n < my; n <<= 1, ++shift_count)
64       ;
65 
66     if (static_cast<int>(shift_count) > exp)
67       break;
68 
69     exp -= shift_count;
70     if (0 <= exp && exp < QUOTIENT_LSB_BITS)
71       q |= (1 << exp);
72 
73     mx = n - my;
74     if (mx == 0) {
75       q = result_sign ? -q : q;
76       return __llvm_libc::fputil::copysign(T(0.0), x);
77     }
78   }
79 
80   NormalFloat<T> remainder(exp + normaly.exponent, mx, 0);
81 
82   // Since NormalFloat to native type conversion is a truncation operation
83   // currently, the remainder value in the native type is correct as is.
84   // However, if NormalFloat to native type conversion is updated in future,
85   // then the conversion to native remainder value should be updated
86   // appropriately and some directed tests added.
87   T native_remainder(remainder);
88   T absy = T(ybits);
89   int cmp = remainder.mul2(1).cmp(normaly);
90   if (cmp > 0) {
91     q = q + 1;
92     if (x >= T(0.0))
93       native_remainder = native_remainder - absy;
94     else
95       native_remainder = absy - native_remainder;
96   } else if (cmp == 0) {
97     if (q & 1) {
98       q += 1;
99       if (x >= T(0.0))
100         native_remainder = -native_remainder;
101     } else {
102       if (x < T(0.0))
103         native_remainder = -native_remainder;
104     }
105   } else {
106     if (x < T(0.0))
107       native_remainder = -native_remainder;
108   }
109 
110   q = result_sign ? -q : q;
111   if (native_remainder == T(0.0))
112     return __llvm_libc::fputil::copysign(T(0.0), x);
113   return native_remainder;
114 }
115 
116 } // namespace fputil
117 } // namespace __llvm_libc
118 
119 #endif // LLVM_LIBC_SRC_SUPPORT_FPUTIL_DIVISION_AND_REMAINDER_OPERATIONS_H
120