297 lines · c
1//===-- Common header for fmod implementations ------------------*- 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_GENERIC_FMOD_H10#define LLVM_LIBC_SRC___SUPPORT_FPUTIL_GENERIC_FMOD_H11 12#include "src/__support/CPP/bit.h"13#include "src/__support/CPP/limits.h"14#include "src/__support/CPP/type_traits.h"15#include "src/__support/FPUtil/FEnvImpl.h"16#include "src/__support/FPUtil/FPBits.h"17#include "src/__support/macros/config.h"18#include "src/__support/macros/optimization.h" // LIBC_UNLIKELY19 20namespace LIBC_NAMESPACE_DECL {21namespace fputil {22namespace generic {23 24// Objective:25// The algorithm uses integer arithmetic (max uint64_t) for general case.26// Some common cases, like abs(x) < abs(y) or abs(x) < 1000 * abs(y) are27// treated specially to increase performance. The part of checking special28// cases, numbers NaN, INF etc. treated separately.29//30// Objective:31// 1) FMod definition (https://cplusplus.com/reference/cmath/fmod/):32// fmod = numer - tquot * denom, where tquot is the truncated33// (i.e., rounded towards zero) result of: numer/denom.34// 2) FMod with negative x and/or y can be trivially converted to fmod for35// positive x and y. Therefore the algorithm below works only with36// positive numbers.37// 3) All positive floating point numbers can be represented as m * 2^e,38// where "m" is positive integer and "e" is signed.39// 4) FMod function can be calculated in integer numbers (x > y):40// fmod = m_x * 2^e_x - tquot * m_y * 2^e_y41// = 2^e_y * (m_x * 2^(e_x - e^y) - tquot * m_y).42// All variables in parentheses are unsigned integers.43//44// Mathematical background:45// Input x,y in the algorithm is represented (mathematically) like m_x*2^e_x46// and m_y*2^e_y. This is an ambiguous number representation. For example:47// m * 2^e = (2 * m) * 2^(e-1)48// The algorithm uses the facts that49// r = a % b = (a % (N * b)) % b,50// (a * c) % (b * c) = (a % b) * c51// where N is positive integer number. a, b and c - positive. Let's adopt52// the formula for representation above.53// a = m_x * 2^e_x, b = m_y * 2^e_y, N = 2^k54// r(k) = a % b = (m_x * 2^e_x) % (2^k * m_y * 2^e_y)55// = 2^(e_y + k) * (m_x * 2^(e_x - e_y - k) % m_y)56// r(k) = m_r * 2^e_r = (m_x % m_y) * 2^(m_y + k)57// = (2^p * (m_x % m_y) * 2^(e_y + k - p))58// m_r = 2^p * (m_x % m_y), e_r = m_y + k - p59//60// Algorithm description:61// First, let write x = m_x * 2^e_x and y = m_y * 2^e_y with m_x, m_y, e_x, e_y62// are integers (m_x amd m_y positive).63// Then the naive implementation of the fmod function with a simple64// for/while loop:65// while (e_x > e_y) {66// m_x *= 2; --e_x; // m_x * 2^e_x == 2 * m_x * 2^(e_x - 1)67// m_x %= m_y;68// }69// On the other hand, the algorithm exploits the fact that m_x, m_y are the70// mantissas of floating point numbers, which use less bits than the storage71// integers: 24 / 32 for floats and 53 / 64 for doubles, so if in each step of72// the iteration, we can left shift m_x as many bits as the storage integer73// type can hold, the exponent reduction per step will be at least 32 - 24 = 874// for floats and 64 - 53 = 11 for doubles (double example below):75// while (e_x > e_y) {76// m_x <<= 11; e_x -= 11; // m_x * 2^e_x == 2^11 * m_x * 2^(e_x - 11)77// m_x %= m_y;78// }79// Some extra improvements are done:80// 1) Shift m_y maximum to the right, which can significantly improve81// performance for small integer numbers (y = 3 for example).82// The m_x shift in the loop can be 62 instead of 11 for double.83// 2) For some architectures with very slow division, it can be better to84// calculate inverse value ones, and after do multiplication in the loop.85// 3) "likely" special cases are treated specially to improve performance.86//87// Simple example:88// The examples below use byte for simplicity.89// 1) Shift hy maximum to right without losing bits and increase iy value90// m_y = 0b00101100 e_y = 20 after shift m_y = 0b00001011 e_y = 22.91// 2) m_x = m_x % m_y.92// 3) Move m_x maximum to left. Note that after (m_x = m_x % m_y) CLZ in m_x93// is not lower than CLZ in m_y. m_x=0b00001001 e_x = 100, m_x=0b10010000,94// e_x = 100-4 = 96.95// 4) Repeat (2) until e_x == e_y.96//97// Complexity analysis (double):98// Converting x,y to (m_x,e_x),(m_y, e_y): CTZ/shift/AND/OR/if. Loop count:99// (m_x - m_y) / (64 - "length of m_y").100// max("length of m_y") = 53,101// max(e_x - e_y) = 2048102// Maximum operation is 186. For rare "unrealistic" cases.103//104// Special cases (double):105// Supposing that case where |y| > 1e-292 and |x/y|<2000 is very common106// special processing is implemented. No m_y alignment, no loop:107// result = (m_x * 2^(e_x - e_y)) % m_y.108// When x and y are both subnormal (rare case but...) the109// result = m_x % m_y.110// Simplified conversion back to double.111 112// Exceptional cases handler according to cppreference.com113// https://en.cppreference.com/w/cpp/numeric/math/fmod114// and POSIX standard described in Linux man115// https://man7.org/linux/man-pages/man3/fmod.3p.html116// C standard for the function is not full, so not by default (although it can117// be implemented in another handler.118// Signaling NaN converted to quiet NaN with FE_INVALID exception.119// https://www.open-std.org/JTC1/SC22/WG14/www/docs/n1011.htm120template <typename T> struct FModDivisionSimpleHelper {121 LIBC_INLINE constexpr static T execute(int exp_diff, int sides_zeroes_count,122 T m_x, T m_y) {123 while (exp_diff > sides_zeroes_count) {124 exp_diff -= sides_zeroes_count;125 m_x <<= sides_zeroes_count;126 m_x %= m_y;127 }128 m_x <<= exp_diff;129 m_x %= m_y;130 return m_x;131 }132};133 134template <typename T> struct FModDivisionInvMultHelper {135 LIBC_INLINE constexpr static T execute(int exp_diff, int sides_zeroes_count,136 T m_x, T m_y) {137 constexpr int LENGTH = sizeof(T) * CHAR_BIT;138 if (exp_diff > sides_zeroes_count) {139 T inv_hy = (cpp::numeric_limits<T>::max() / m_y);140 while (exp_diff > sides_zeroes_count) {141 exp_diff -= sides_zeroes_count;142 T hd = (m_x * inv_hy) >> (LENGTH - sides_zeroes_count);143 m_x <<= sides_zeroes_count;144 m_x -= hd * m_y;145 while (LIBC_UNLIKELY(m_x > m_y))146 m_x -= m_y;147 }148 T hd = (m_x * inv_hy) >> (LENGTH - exp_diff);149 m_x <<= exp_diff;150 m_x -= hd * m_y;151 while (LIBC_UNLIKELY(m_x > m_y))152 m_x -= m_y;153 } else {154 m_x <<= exp_diff;155 m_x %= m_y;156 }157 return m_x;158 }159};160 161template <typename T, typename U = typename FPBits<T>::StorageType,162 typename DivisionHelper = FModDivisionSimpleHelper<U>>163class FMod {164 static_assert(cpp::is_floating_point_v<T> &&165 is_unsigned_integral_or_big_int_v<U> &&166 (sizeof(U) * CHAR_BIT > FPBits<T>::FRACTION_LEN),167 "FMod instantiated with invalid type.");168 169private:170 using FPB = FPBits<T>;171 using StorageType = typename FPB::StorageType;172 173 LIBC_INLINE static bool pre_check(T x, T y, T &out) {174 using FPB = fputil::FPBits<T>;175 const T quiet_nan = FPB::quiet_nan().get_val();176 FPB sx(x), sy(y);177 if (LIBC_LIKELY(!sy.is_zero() && !sy.is_inf_or_nan() &&178 !sx.is_inf_or_nan()))179 return false;180 181 if (sx.is_nan() || sy.is_nan()) {182 if (sx.is_signaling_nan() || sy.is_signaling_nan())183 fputil::raise_except_if_required(FE_INVALID);184 out = quiet_nan;185 return true;186 }187 188 if (sx.is_inf() || sy.is_zero()) {189 fputil::raise_except_if_required(FE_INVALID);190 fputil::set_errno_if_required(EDOM);191 out = quiet_nan;192 return true;193 }194 195 out = x;196 return true;197 }198 199 LIBC_INLINE static constexpr FPB eval_internal(FPB sx, FPB sy) {200 201 if (LIBC_LIKELY(sx.uintval() <= sy.uintval())) {202 if (sx.uintval() < sy.uintval())203 return sx; // |x|<|y| return x204 return FPB::zero(); // |x|=|y| return 0.0205 }206 207 int e_x = sx.get_biased_exponent();208 int e_y = sy.get_biased_exponent();209 210 // Most common case where |y| is "very normal" and |x/y| < 2^EXP_LEN211 if (LIBC_LIKELY(e_y > int(FPB::FRACTION_LEN) &&212 e_x - e_y <= int(FPB::EXP_LEN))) {213 StorageType m_x = sx.get_explicit_mantissa();214 StorageType m_y = sy.get_explicit_mantissa();215 StorageType d = (e_x == e_y)216 ? (m_x - m_y)217 : static_cast<StorageType>(m_x << (e_x - e_y)) % m_y;218 if (d == 0)219 return FPB::zero();220 // iy - 1 because of "zero power" for number with power 1221 return FPB::make_value(d, e_y - 1);222 }223 // Both subnormal special case.224 if (LIBC_UNLIKELY(e_x == 0 && e_y == 0)) {225 FPB d;226 d.set_mantissa(sx.uintval() % sy.uintval());227 return d;228 }229 230 // Note that hx is not subnormal by conditions above.231 U m_x = static_cast<U>(sx.get_explicit_mantissa());232 e_x--;233 234 U m_y = static_cast<U>(sy.get_explicit_mantissa());235 constexpr int DEFAULT_LEAD_ZEROS =236 sizeof(U) * CHAR_BIT - FPB::FRACTION_LEN - 1;237 int lead_zeros_m_y = DEFAULT_LEAD_ZEROS;238 if (LIBC_LIKELY(e_y > 0)) {239 e_y--;240 } else {241 m_y = static_cast<U>(sy.get_mantissa());242 lead_zeros_m_y = cpp::countl_zero(m_y);243 }244 245 // Assume hy != 0246 int tail_zeros_m_y = cpp::countr_zero(m_y);247 int sides_zeroes_count = lead_zeros_m_y + tail_zeros_m_y;248 // n > 0 by conditions above249 int exp_diff = e_x - e_y;250 {251 // Shift hy right until the end or n = 0252 int right_shift = exp_diff < tail_zeros_m_y ? exp_diff : tail_zeros_m_y;253 m_y >>= right_shift;254 exp_diff -= right_shift;255 e_y += right_shift;256 }257 258 {259 // Shift hx left until the end or n = 0260 int left_shift =261 exp_diff < DEFAULT_LEAD_ZEROS ? exp_diff : DEFAULT_LEAD_ZEROS;262 m_x <<= left_shift;263 exp_diff -= left_shift;264 }265 266 m_x %= m_y;267 if (LIBC_UNLIKELY(m_x == 0))268 return FPB::zero();269 270 if (exp_diff == 0)271 return FPB::make_value(static_cast<StorageType>(m_x), e_y);272 273 // hx next can't be 0, because hx < hy, hy % 2 == 1 hx * 2^i % hy != 0274 m_x = DivisionHelper::execute(exp_diff, sides_zeroes_count, m_x, m_y);275 return FPB::make_value(static_cast<StorageType>(m_x), e_y);276 }277 278public:279 LIBC_INLINE static T eval(T x, T y) {280 if (T out; LIBC_UNLIKELY(pre_check(x, y, out)))281 return out;282 FPB sx(x), sy(y);283 Sign sign = sx.sign();284 sx.set_sign(Sign::POS);285 sy.set_sign(Sign::POS);286 FPB result = eval_internal(sx, sy);287 result.set_sign(sign);288 return result.get_val();289 }290};291 292} // namespace generic293} // namespace fputil294} // namespace LIBC_NAMESPACE_DECL295 296#endif // LLVM_LIBC_SRC___SUPPORT_FPUTIL_GENERIC_FMOD_H297