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