8514ecb02d
The implementation is not fully standards compliant in the sense that errno is not set on error, and floating point exceptions are not raised. Subnormal range and normal range are tested separately in the tests. Reviewed By: lntue Differential Revision: https://reviews.llvm.org/D86666
112 lines
3.3 KiB
C++
112 lines
3.3 KiB
C++
//===-- 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() || ybits.isInf()) {
|
|
q = 0;
|
|
return __llvm_libc::fputil::copysign(T(0.0), x);
|
|
}
|
|
|
|
bool resultSign = (xbits.sign == ybits.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.sign = ybits.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 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)
|
|
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
|