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

source code of libc/src/__support/FPUtil/DivisionAndRemainderOperations.h