344 lines · c
1//===-- Utility class to manipulate fixed point numbers. --*- 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_FIXED_POINT_FX_BITS_H10#define LLVM_LIBC_SRC___SUPPORT_FIXED_POINT_FX_BITS_H11 12#include "include/llvm-libc-macros/stdfix-macros.h"13#include "src/__support/CPP/algorithm.h"14#include "src/__support/CPP/bit.h"15#include "src/__support/CPP/limits.h" // numeric_limits16#include "src/__support/CPP/type_traits.h"17#include "src/__support/libc_assert.h"18#include "src/__support/macros/attributes.h" // LIBC_INLINE19#include "src/__support/macros/config.h" // LIBC_NAMESPACE_DECL20#include "src/__support/macros/null_check.h" // LIBC_CRASH_ON_VALUE21#include "src/__support/macros/optimization.h" // LIBC_UNLIKELY22#include "src/__support/math_extras.h"23 24#include "fx_rep.h"25 26#include <stdio.h>27 28#ifdef LIBC_COMPILER_HAS_FIXED_POINT29 30namespace LIBC_NAMESPACE_DECL {31namespace fixed_point {32 33template <typename T> struct FXBits {34private:35 using fx_rep = FXRep<T>;36 using StorageType = typename fx_rep::StorageType;37 38 StorageType value;39 40 static_assert(fx_rep::FRACTION_LEN > 0);41 42 static constexpr size_t FRACTION_OFFSET = 0; // Just for completeness43 static constexpr size_t INTEGRAL_OFFSET =44 fx_rep::INTEGRAL_LEN == 0 ? 0 : fx_rep::FRACTION_LEN;45 static constexpr size_t SIGN_OFFSET =46 fx_rep::SIGN_LEN == 047 ? 048 : ((sizeof(StorageType) * CHAR_BIT) - fx_rep::SIGN_LEN);49 50 static constexpr StorageType FRACTION_MASK =51 mask_trailing_ones<StorageType, fx_rep::FRACTION_LEN>()52 << FRACTION_OFFSET;53 static constexpr StorageType INTEGRAL_MASK =54 mask_trailing_ones<StorageType, fx_rep::INTEGRAL_LEN>()55 << INTEGRAL_OFFSET;56 static constexpr StorageType SIGN_MASK =57 (fx_rep::SIGN_LEN == 0 ? 0 : StorageType(1) << SIGN_OFFSET);58 59 // mask for <integral | fraction>60 static constexpr StorageType VALUE_MASK = INTEGRAL_MASK | FRACTION_MASK;61 62 // mask for <sign | integral | fraction>63 static constexpr StorageType TOTAL_MASK = SIGN_MASK | VALUE_MASK;64 65public:66 LIBC_INLINE constexpr FXBits() = default;67 68 template <typename XType> LIBC_INLINE constexpr explicit FXBits(XType x) {69 using Unqual = typename cpp::remove_cv_t<XType>;70 if constexpr (cpp::is_same_v<Unqual, T>) {71 value = cpp::bit_cast<StorageType>(x);72 } else if constexpr (cpp::is_same_v<Unqual, StorageType>) {73 value = x;74 } else {75 // We don't want accidental type promotions/conversions, so we require76 // exact type match.77 static_assert(cpp::always_false<XType>);78 }79 }80 81 LIBC_INLINE constexpr StorageType get_fraction() {82 return (value & FRACTION_MASK) >> FRACTION_OFFSET;83 }84 85 LIBC_INLINE constexpr StorageType get_integral() {86 return (value & INTEGRAL_MASK) >> INTEGRAL_OFFSET;87 }88 89 // returns complete bitstring representation the fixed point number90 // the bitstring is of the form: padding | sign | integral | fraction91 LIBC_INLINE constexpr StorageType get_bits() {92 return (value & TOTAL_MASK) >> FRACTION_OFFSET;93 }94 95 // TODO: replace bool with Sign96 LIBC_INLINE constexpr bool get_sign() {97 return static_cast<bool>((value & SIGN_MASK) >> SIGN_OFFSET);98 }99 100 // This represents the effective negative exponent applied to this number101 LIBC_INLINE constexpr int get_exponent() { return fx_rep::FRACTION_LEN; }102 103 LIBC_INLINE constexpr void set_fraction(StorageType fraction) {104 value = (value & (~FRACTION_MASK)) |105 ((fraction << FRACTION_OFFSET) & FRACTION_MASK);106 }107 108 LIBC_INLINE constexpr void set_integral(StorageType integral) {109 value = (value & (~INTEGRAL_MASK)) |110 ((integral << INTEGRAL_OFFSET) & INTEGRAL_MASK);111 }112 113 // TODO: replace bool with Sign114 LIBC_INLINE constexpr void set_sign(bool sign) {115 value = (value & (~SIGN_MASK)) |116 ((static_cast<StorageType>(sign) << SIGN_OFFSET) & SIGN_MASK);117 }118 119 LIBC_INLINE constexpr T get_val() const { return cpp::bit_cast<T>(value); }120};121 122// Bit-wise operations are not available for fixed point types yet.123template <typename T>124LIBC_INLINE constexpr cpp::enable_if_t<cpp::is_fixed_point_v<T>, T>125bit_and(T x, T y) {126 using BitType = typename FXRep<T>::StorageType;127 BitType x_bit = cpp::bit_cast<BitType>(x);128 BitType y_bit = cpp::bit_cast<BitType>(y);129 // For some reason, bit_cast cannot deduce BitType from the input.130 return cpp::bit_cast<T, BitType>(x_bit & y_bit);131}132 133template <typename T>134LIBC_INLINE constexpr cpp::enable_if_t<cpp::is_fixed_point_v<T>, T>135bit_or(T x, T y) {136 using BitType = typename FXRep<T>::StorageType;137 BitType x_bit = cpp::bit_cast<BitType>(x);138 BitType y_bit = cpp::bit_cast<BitType>(y);139 // For some reason, bit_cast cannot deduce BitType from the input.140 return cpp::bit_cast<T, BitType>(x_bit | y_bit);141}142 143template <typename T>144LIBC_INLINE constexpr cpp::enable_if_t<cpp::is_fixed_point_v<T>, T>145bit_not(T x) {146 using BitType = typename FXRep<T>::StorageType;147 BitType x_bit = cpp::bit_cast<BitType>(x);148 // For some reason, bit_cast cannot deduce BitType from the input.149 return cpp::bit_cast<T, BitType>(static_cast<BitType>(~x_bit));150}151 152template <typename T> LIBC_INLINE constexpr T abs(T x) {153 using FXRep = FXRep<T>;154 if constexpr (FXRep::SIGN_LEN == 0)155 return x;156 else {157 if (LIBC_UNLIKELY(x == FXRep::MIN()))158 return FXRep::MAX();159 return (x < FXRep::ZERO() ? -x : x);160 }161}162 163// Round-to-nearest, tie-to-(+Inf)164template <typename T> LIBC_INLINE constexpr T round(T x, int n) {165 using FXRep = FXRep<T>;166 if (LIBC_UNLIKELY(n < 0))167 n = 0;168 if (LIBC_UNLIKELY(n >= FXRep::FRACTION_LEN))169 return x;170 171 T round_bit = FXRep::EPS() << (FXRep::FRACTION_LEN - n - 1);172 // Check for overflow.173 if (LIBC_UNLIKELY(FXRep::MAX() - round_bit < x))174 return FXRep::MAX();175 176 T all_ones = bit_not(FXRep::ZERO());177 178 int shift = FXRep::FRACTION_LEN - n;179 T rounding_mask =180 (shift == FXRep::TOTAL_LEN) ? FXRep::ZERO() : (all_ones << shift);181 return bit_and((x + round_bit), rounding_mask);182}183 184// count leading sign bits185// TODO: support fixed_point_padding186template <typename T>187LIBC_INLINE constexpr cpp::enable_if_t<cpp::is_fixed_point_v<T>, int>188countls(T f) {189 using FXRep = FXRep<T>;190 using BitType = typename FXRep::StorageType;191 using FXBits = FXBits<T>;192 193 if constexpr (FXRep::SIGN_LEN > 0) {194 if (f < 0)195 f = bit_not(f);196 }197 198 BitType value_bits = FXBits(f).get_bits();199 return cpp::countl_zero(value_bits) - FXRep::SIGN_LEN;200}201 202// fixed-point to integer conversion203template <typename T, typename XType>204LIBC_INLINE constexpr cpp::enable_if_t<cpp::is_fixed_point_v<T>, XType>205bitsfx(T f) {206 return cpp::bit_cast<XType, T>(f);207}208 209// divide the two fixed-point types and return an integer result210template <typename T, typename XType>211LIBC_INLINE constexpr cpp::enable_if_t<cpp::is_fixed_point_v<T>, XType>212idiv(T x, T y) {213 using FXBits = FXBits<T>;214 using FXRep = FXRep<T>;215 using CompType = typename FXRep::CompType;216 217 // If the value of the second operand of the / operator is zero, the218 // behavior is undefined. Ref: ISO/IEC TR 18037:2008(E) p.g. 16219 LIBC_CRASH_ON_VALUE(y, FXRep::ZERO());220 221 CompType x_comp = static_cast<CompType>(FXBits(x).get_bits());222 CompType y_comp = static_cast<CompType>(FXBits(y).get_bits());223 224 // If an integer result of one of these functions overflows, the behavior is225 // undefined. Ref: ISO/IEC TR 18037:2008(E) p.g. 16226 CompType result = x_comp / y_comp;227 228 return static_cast<XType>(result);229}230 231LIBC_INLINE long accum nrstep(long accum d, long accum x0) {232 auto v = x0 * (2.lk - (d * x0));233 return v;234}235 236// Divide the two integers and return a fixed_point value237//238// For reference, see:239// https://en.wikipedia.org/wiki/Division_algorithm#Newton%E2%80%93Raphson_division240// https://stackoverflow.com/a/9231996241 242template <typename XType> LIBC_INLINE constexpr XType divi(int n, int d) {243 // If the value of the second operand of the / operator is zero, the244 // behavior is undefined. Ref: ISO/IEC TR 18037:2008(E) p.g. 16245 LIBC_CRASH_ON_VALUE(d, 0);246 247 if (LIBC_UNLIKELY(n == 0)) {248 return FXRep<XType>::ZERO();249 }250 auto is_power_of_two = [](int n) { return (n > 0) && ((n & (n - 1)) == 0); };251 long accum max_val = static_cast<long accum>(FXRep<XType>::MAX());252 long accum min_val = static_cast<long accum>(FXRep<XType>::MIN());253 254 if (is_power_of_two(cpp::abs(d))) {255 int k = cpp::countr_zero<uint32_t>(static_cast<uint32_t>(cpp::abs(d)));256 constexpr int F = FXRep<XType>::FRACTION_LEN;257 int64_t scaled_n = static_cast<int64_t>(n) << F;258 int64_t res64 = scaled_n >> k;259 constexpr int TOTAL_BITS = sizeof(XType) * 8;260 const int64_t max_limit = (1LL << (TOTAL_BITS - 1)) - 1;261 const int64_t min_limit = -(1LL << (TOTAL_BITS - 1));262 if (res64 > max_limit) {263 return FXRep<XType>::MAX();264 } else if (res64 < min_limit) {265 return FXRep<XType>::MIN();266 }267 long accum res_accum =268 static_cast<long accum>(res64) / static_cast<long accum>(1 << F);269 res_accum = (d < 0) ? static_cast<long accum>(-1) * res_accum : res_accum;270 if (res_accum > max_val) {271 return FXRep<XType>::MAX();272 } else if (res_accum < min_val) {273 return FXRep<XType>::MIN();274 }275 return static_cast<XType>(res_accum);276 }277 278 bool result_is_negative = ((n < 0) != (d < 0));279 int64_t n64 = static_cast<int64_t>(n);280 int64_t d64 = static_cast<int64_t>(d);281 282 uint64_t nv = static_cast<uint64_t>(n64 < 0 ? -n64 : n64);283 uint64_t dv = static_cast<uint64_t>(d64 < 0 ? -d64 : d64);284 285 if (d == INT_MIN) {286 nv <<= 1;287 dv >>= 1;288 }289 290 uint32_t clz = cpp::countl_zero<uint32_t>(static_cast<uint32_t>(dv)) - 1;291 uint64_t scaled_val = dv << clz;292 // Scale denominator to be in the range of [0.5,1]293 FXBits<long accum> d_scaled{scaled_val};294 uint64_t scaled_val_n = nv << clz;295 // Scale the numerator as much as the denominator to maintain correctness of296 // the original equation297 FXBits<long accum> n_scaled{scaled_val_n};298 long accum n_scaled_val = n_scaled.get_val();299 long accum d_scaled_val = d_scaled.get_val();300 // x0 = (48/17) - (32/17) * d_n301 long accum a = 0x2.d89d89d8p0lk; // 48/17 = 2.8235294...302 long accum b = 0x1.e1e1e1e1p0lk; // 32/17 = 1.8823529...303 // Error of the initial approximation, as derived304 // from the wikipedia article is305 // E0 = 1/17 = 0.059 (5.9%)306 long accum initial_approx = a - (b * d_scaled_val);307 // Since, 0.5 <= d_scaled_val <= 1.0, 0.9412 <= initial_approx <= 1.88235308 LIBC_ASSERT((initial_approx >= 0x0.78793dd9p0lk) &&309 (initial_approx <= 0x1.f0f0d845p0lk));310 // Each newton-raphson iteration will square the error, due311 // to quadratic convergence. So,312 // E1 = (0.059)^2 = 0.0034313 long accum val = nrstep(d_scaled_val, initial_approx);314 if constexpr (FXRep<XType>::FRACTION_LEN > 8) {315 // E2 = 0.0000121316 val = nrstep(d_scaled_val, val);317 if constexpr (FXRep<XType>::FRACTION_LEN > 16) {318 // E3 = 1.468e−10319 val = nrstep(d_scaled_val, val);320 }321 }322 long accum res = n_scaled_val * val;323 324 if (result_is_negative) {325 res *= static_cast<long accum>(-1);326 }327 328 // Per clause 7.18a.6.1, saturate values on overflow329 if (res > max_val) {330 return FXRep<XType>::MAX();331 } else if (res < min_val) {332 return FXRep<XType>::MIN();333 } else {334 return static_cast<XType>(res);335 }336}337 338} // namespace fixed_point339} // namespace LIBC_NAMESPACE_DECL340 341#endif // LIBC_COMPILER_HAS_FIXED_POINT342 343#endif // LLVM_LIBC_SRC___SUPPORT_FIXED_POINT_FX_BITS_H344