122 lines · c
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-exception6//7//===----------------------------------------------------------------------===//8 9#ifndef LLVM_LIBC_SRC___SUPPORT_FPUTIL_DIVISIONANDREMAINDEROPERATIONS_H10#define LLVM_LIBC_SRC___SUPPORT_FPUTIL_DIVISIONANDREMAINDEROPERATIONS_H11 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 20namespace LIBC_NAMESPACE_DECL {21namespace fputil {22 23static constexpr int QUOTIENT_LSB_BITS = 3;24 25// The implementation is a bit-by-bit algorithm which uses integer division26// to evaluate the quotient and remainder.27template <typename T, cpp::enable_if_t<cpp::is_floating_point_v<T>, int> = 0>28LIBC_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 absolute51 // values. The correct sign can be applied to the result after the result52 // 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 operation85 // 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 updated88 // 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 else97 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 fputil119} // namespace LIBC_NAMESPACE_DECL120 121#endif // LLVM_LIBC_SRC___SUPPORT_FPUTIL_DIVISIONANDREMAINDEROPERATIONS_H122