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