libc/src/__support/FPUtil/DivisionAndRemainderOperations.h

   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