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
119
|
//===-- 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_UTILS_FPUTIL_DIVISION_AND_REMAINDER_OPERATIONS_H
#define LLVM_LIBC_UTILS_FPUTIL_DIVISION_AND_REMAINDER_OPERATIONS_H
#include "FPBits.h"
#include "ManipulationFunctions.h"
#include "NormalFloat.h"
#include "utils/CPP/TypeTraits.h"
namespace __llvm_libc {
namespace fputil {
static constexpr int quotientLSBBits = 3;
// The implementation is a bit-by-bit algorithm which uses integer division
// to evaluate the quotient and remainder.
template <typename T,
cpp::EnableIfType<cpp::IsFloatingPointType<T>::Value, int> = 0>
static inline T remquo(T x, T y, int &q) {
FPBits<T> xbits(x), ybits(y);
if (xbits.isNaN())
return x;
if (ybits.isNaN())
return y;
if (xbits.isInf() || ybits.isZero())
return FPBits<T>::buildNaN(1);
if (xbits.isZero()) {
q = 0;
return __llvm_libc::fputil::copysign(T(0.0), x);
}
if (ybits.isInf()) {
q = 0;
return x;
}
bool resultSign = (xbits.getSign() == ybits.getSign() ? 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.setSign(0);
ybits.setSign(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 shiftCount = 0;
typename NormalFloat<T>::UIntType n = mx;
for (shiftCount = 0; n < my; n <<= 1, ++shiftCount)
;
if (static_cast<int>(shiftCount) > exp)
break;
exp -= shiftCount;
if (0 <= exp && exp < quotientLSBBits)
q |= (1 << exp);
mx = n - my;
if (mx == 0) {
q = resultSign ? -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 nativeRemainder(remainder);
T absy = T(ybits);
int cmp = remainder.mul2(1).cmp(normaly);
if (cmp > 0) {
q = q + 1;
if (x >= T(0.0))
nativeRemainder = nativeRemainder - absy;
else
nativeRemainder = absy - nativeRemainder;
} else if (cmp == 0) {
if (q & 1) {
q += 1;
if (x >= T(0.0))
nativeRemainder = -nativeRemainder;
} else {
if (x < T(0.0))
nativeRemainder = -nativeRemainder;
}
} else {
if (x < T(0.0))
nativeRemainder = -nativeRemainder;
}
q = resultSign ? -q : q;
if (nativeRemainder == T(0.0))
return __llvm_libc::fputil::copysign(T(0.0), x);
return nativeRemainder;
}
} // namespace fputil
} // namespace __llvm_libc
#endif // LLVM_LIBC_UTILS_FPUTIL_DIVISION_AND_REMAINDER_OPERATIONS_H
|
