From fc7eb739e73d565bb2edf083ada78510cb6fb37d Mon Sep 17 00:00:00 2001 From: Mateusz Pusz Date: Sat, 23 May 2026 23:19:38 +0200 Subject: [PATCH] fix: complete `double_width_int` operator suite so it can stand in for `__int128` on MSVC (#795) `double_width_int` was previously missing a number of operators that generic numerical code expects from a 128-bit integer. On platforms with native `__int128` this was harmless, but on MSVC (where `int128_t` / `uint128_t` alias `double_width_int<(u)int64_t>`) the recent `UsesIntegerScaling` change wired the synthetic dwint into concept checks, and `compare_quantities` + magnitude folding into runtime paths, exposing every gap. This change rounds out the operator set so `double_width_int` truly behaves as an integer type: * binary `+`, `-`, `*`, `/` between two `double_width_int`s (alongside the existing narrow-rhs overloads). * unary `~` and binary `&`, `|`, `^`. * compound assignment for arithmetic, bitwise, and shift operations. * pre/post `++` and `--`. * `static_cast` (and `double`/`float`) with sign-preserving conversion that avoids catastrophic cancellation on platforms where `long double == double`. * `std::numeric_limits>` specialization so generic code probing `::max()`, `::min()`, `::digits`, `::is_signed`, etc. gets correct answers (needed by `checked_int_pow`, `compute_base_power`, `safe_int::operator-`). * fields `hi_` / `lo_` and the `(hi, lo)` ctor are public so cross-instance inline friend operators can access each other without further friend declarations. While here, also: * extract `double_width_int` (and its `std::numeric_limits` specialization) into a dedicated header `bits/double_width_int.h`; `fixed_point.h` now just includes it and keeps the `int128_t` aliases, `min_width_uint_t` / `double_width_int_for_t` / `wide_product_of` helpers, and the `fixed_point` class itself. * rename a local variable `m` in `wide_product_of` to `mid` to avoid shadowing `si::unit_symbols::m` (MSVC C4459). * rewrite the `lo_ > 0 ? -1 : 0` unary-minus body to use `Tl{0} - lo_` instead of `-lo_`, silencing MSVC C4146 about unary minus on an unsigned operand. Tests: * new `test/static/double_width_int_test.cpp` pins down each new operator at compile time (carry/borrow edges, schoolbook multiplication, narrow/wide division paths, bitwise identities, numeric_limits values, long-double round-trips). Co-Authored-By: Claude Opus 4.7 (1M context) --- src/core/CMakeLists.txt | 1 + .../include/mp-units/bits/double_width_int.h | 495 ++++++++++++++++++ src/core/include/mp-units/bits/fixed_point.h | 266 +--------- test/static/CMakeLists.txt | 1 + test/static/double_width_int_test.cpp | 140 +++++ test/static/fixed_point_test.cpp | 6 +- 6 files changed, 639 insertions(+), 270 deletions(-) create mode 100644 src/core/include/mp-units/bits/double_width_int.h create mode 100644 test/static/double_width_int_test.cpp diff --git a/src/core/CMakeLists.txt b/src/core/CMakeLists.txt index 0540dbffa..6a317600e 100644 --- a/src/core/CMakeLists.txt +++ b/src/core/CMakeLists.txt @@ -31,6 +31,7 @@ add_mp_units_module( core mp-units-core HEADERS include/mp-units/bits/constexpr_math.h include/mp-units/bits/core_gmf.h + include/mp-units/bits/double_width_int.h include/mp-units/bits/fixed_point.h include/mp-units/bits/get_associated_quantity.h include/mp-units/bits/hacks.h diff --git a/src/core/include/mp-units/bits/double_width_int.h b/src/core/include/mp-units/bits/double_width_int.h new file mode 100644 index 000000000..1b822441c --- /dev/null +++ b/src/core/include/mp-units/bits/double_width_int.h @@ -0,0 +1,495 @@ +// The MIT License (MIT) +// +// Copyright (c) 2018 Mateusz Pusz +// +// Permission is hereby granted, free of charge, to any person obtaining a copy +// of this software and associated documentation files (the "Software"), to deal +// in the Software without restriction, including without limitation the rights +// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +// copies of the Software, and to permit persons to whom the Software is +// furnished to do so, subject to the following conditions: +// +// The above copyright notice and this permission notice shall be included in all +// copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +// SOFTWARE. + +#pragma once + +#include // IWYU pragma: keep +#include // IWYU pragma: keep + +#ifndef MP_UNITS_IN_MODULE_INTERFACE +#ifdef MP_UNITS_IMPORT_STD +import std; +#else +#include +#include +#include +#include +#include +#include +#include +#include +#include +// becomes partially freestanding in C++26 +// before that, there is no guarantee about this header even existing +// (GCC 14 has it, but actively #error's out) +#if __STDC_HOSTED__ +#include +#endif +#endif +#endif + +namespace mp_units::detail { + +template +constexpr std::size_t integer_rep_width_v = std::numeric_limits>::digits; + +template +[[nodiscard]] consteval T int_power(T base, int exponent) +{ +#if __STDC_HOSTED__ && defined(__cpp_lib_constexpr_cmath) && __cpp_lib_constexpr_cmath >= 202202L + return std::ldexp(base, exponent); +#else + if (exponent < 0) { + base = T{1} / base; + exponent = -exponent; + } + T ret = 1; + while (exponent) { + if (exponent & 1) ret *= base; + exponent >>= 1; + base *= base; + } + return ret; + +#endif +} + +// A double-width integer synthesized from two base-width integers. On platforms without a +// native `__int128` (notably MSVC), `int128_t` / `uint128_t` are aliases for this template +// instantiated on (u)int64_t; the operator set is therefore comprehensive enough for it to +// stand in for a real 128-bit integer. See issue #795 for the original report from the +// MSVC team about missing operators. +template +struct double_width_int { + static constexpr bool is_signed = std::is_signed_v; + static constexpr std::size_t base_width = integer_rep_width_v; + static constexpr std::size_t width = 2 * base_width; + + using Th = T; + using Tl = std::make_unsigned_t; + + constexpr double_width_int() = default; + + constexpr double_width_int(Th hi, Tl lo) : hi_(hi), lo_(lo) {} + + static constexpr double_width_int from_hi_lo(Th hi, Tl lo) { return {hi, lo}; } + + explicit constexpr double_width_int(long double v) + { + constexpr auto scale = int_power(2, base_width); + constexpr auto iscale = 1.l / scale; + auto scaled = v * iscale; + hi_ = static_cast(scaled); + auto resid = (scaled - static_cast(hi_)); + if (resid < 0) { + --hi_; + resid += 1; + } + lo_ = static_cast(resid * scale); + } + template + requires(is_signed || !std::is_signed_v) + explicit(false) constexpr double_width_int(U v) + { + if constexpr (is_signed) { + hi_ = v < 0 ? Th{-1} : Th{0}; + } else { + hi_ = 0; + } + lo_ = static_cast(v); + } + + template + explicit constexpr operator U() const + { + if constexpr (integer_rep_width_v > base_width) { + return (static_cast(hi_) << base_width) + static_cast(lo_); + } else { + return static_cast(lo_); + } + } + + // Conversion to floating-point. Reverses the long-double constructor above; needed by + // generic numerical algorithms (e.g. constexpr_math.h root() does the initial bisection in + // long double) so dwint can stand in for a built-in __int128 there. Negative values are + // converted via their absolute value to avoid catastrophic cancellation on platforms + // where long double has fewer than `width` bits of mantissa (e.g. MSVC, where long double + // == double has only 53 bits, so `hi_ * 2^base_width + lo_` would round both terms before + // they nearly cancel). + explicit constexpr operator long double() const + { + if constexpr (is_signed) { + if (hi_ < Th{0}) return -static_cast(-*this); + } + constexpr auto scale = int_power(2, base_width); + return static_cast(hi_) * scale + static_cast(lo_); + } + explicit constexpr operator double() const { return static_cast(static_cast(*this)); } + explicit constexpr operator float() const { return static_cast(static_cast(*this)); } + + [[nodiscard]] constexpr auto operator<=>(const double_width_int&) const = default; + + // calculates the double-width product of two base-size integers; this implementation requires at least one of them to + // be unsigned + static constexpr double_width_int wide_product_of(Th lhs, Tl rhs) + { + constexpr std::size_t half_width = base_width / 2; + constexpr Tl msk = (Tl(1) << half_width) - 1u; + Th l1 = lhs >> half_width; + Tl l0 = static_cast(lhs) & msk; + Tl r1 = rhs >> half_width; + Tl r0 = rhs & msk; + Tl t00 = l0 * r0; + Tl t01 = l0 * r1; + Th t10 = l1 * static_cast(r0); + Th t11 = l1 * static_cast(r1); + Tl mid = (t01 & msk) + (static_cast(t10) & msk) + (t00 >> half_width); + Th o1 = t11 + static_cast(mid >> half_width) + (t10 >> half_width) + static_cast(t01 >> half_width); + Tl o0 = (t00 & msk) | ((mid & msk) << half_width); + return {o1, o0}; + } + + template + requires(std::numeric_limits::digits <= base_width) + [[nodiscard]] friend constexpr auto operator*(const double_width_int& lhs, Rhs rhs) + { + using RT = conditional, std::make_signed_t, Tl>; + auto lo_prod = double_width_int::wide_product_of(rhs, lhs.lo_); + // Normal C++ rules; with respect to signedness, the wider type always wins. + using ret_t = double_width_int; + return ret_t{static_cast(lo_prod.hi_) + lhs.hi_ * static_cast(rhs), lo_prod.lo_}; + } + + template + requires(std::numeric_limits::digits <= base_width) + [[nodiscard]] friend constexpr auto operator*(Lhs lhs, const double_width_int& rhs) + { + return rhs * lhs; + } + + // Schoolbook multiplication truncated to `width` bits: + // (lhs.hi*B + lhs.lo) * (rhs.hi*B + rhs.lo) where B = 2^base_width; + // the lhs.hi*rhs.hi*B^2 term overflows past `width` bits and is discarded. + // The two cross-products contribute only to the hi half, so modular arithmetic in Tl is enough. + [[nodiscard]] friend constexpr double_width_int operator*(const double_width_int& lhs, const double_width_int& rhs) + { + const auto lo_lo = double_width_int::wide_product_of(lhs.lo_, rhs.lo_); + const Tl mid = static_cast(lhs.hi_) * rhs.lo_ + lhs.lo_ * static_cast(rhs.hi_); + return {static_cast(lo_lo.hi_ + mid), lo_lo.lo_}; + } + + template + requires(std::numeric_limits::digits <= base_width) + [[nodiscard]] friend constexpr double_width_int operator/(const double_width_int& lhs, Rhs rhs) + { + // Normal C++ rules; with respect to signedness, the bigger type always wins. + using ret_t = double_width_int; + if constexpr (std::is_signed_v) { + if (rhs < 0) { + return (-lhs) / static_cast(-rhs); + } else { + return lhs / static_cast(rhs); + } + } else if constexpr (is_signed) { + if (lhs.hi_ < 0) { + return -((-lhs) / rhs); + } else { + using unsigned_t = double_width_int; + auto tmp = unsigned_t{static_cast(lhs.hi_), lhs.lo_} / rhs; + return ret_t{static_cast(tmp.hi_), tmp.lo_}; + } + } else { + Th res_hi = lhs.hi_ / rhs; + // unfortunately, wide division is hard: https://en.wikipedia.org/wiki/Division_algorithm. + // Here, we just provide a somewhat naive implementation of long division. + Tl rem_hi = lhs.hi_ % rhs; + Tl rem_lo = lhs.lo_; + Tl res_lo = 0; + for (std::size_t i = 0; i < base_width; ++i) { + // shift in one bit + rem_hi = (rem_hi << 1u) | (rem_lo >> (base_width - 1)); + rem_lo <<= 1u; + res_lo <<= 1u; + // perform one bit of long division + if (rem_hi >= rhs) { + rem_hi -= rhs; + res_lo |= 1u; + } + } + return ret_t{res_hi, res_lo}; + } + } + + // Division between two double-width values. Falls back to the narrow operator/ when rhs fits + // in T (with sign extension); otherwise reduces to an unsigned bit-by-bit long division. + [[nodiscard]] friend constexpr double_width_int operator/(const double_width_int& lhs, const double_width_int& rhs) + { + if constexpr (is_signed) { + // Fast path: rhs is representable in Th (positive or negative). + if (rhs.hi_ == Th{0} && static_cast(rhs.lo_) >= Th{0}) return lhs / static_cast(rhs.lo_); + if (rhs.hi_ == Th{-1} && static_cast(rhs.lo_) < Th{0}) return lhs / static_cast(rhs.lo_); + // General signed path: reduce to unsigned via absolute values, then re-sign. + using unsigned_t = double_width_int; + const bool neg_lhs = lhs.hi_ < 0; + const bool neg_rhs = rhs.hi_ < 0; + const auto abs_lhs = neg_lhs ? -lhs : lhs; + const auto abs_rhs = neg_rhs ? -rhs : rhs; + const auto uq = unsigned_t::from_hi_lo(static_cast(abs_lhs.hi_), abs_lhs.lo_) / + unsigned_t::from_hi_lo(static_cast(abs_rhs.hi_), abs_rhs.lo_); + const auto q = double_width_int::from_hi_lo(static_cast(uq.hi_), uq.lo_); + return (neg_lhs != neg_rhs) ? -q : q; + } else { + // Fast path: rhs is narrow. + if (rhs.hi_ == Th{0}) return lhs / rhs.lo_; + // Unsigned restoring long division. + double_width_int quotient{}; + double_width_int remainder{}; + for (std::size_t i = 0; i < width; ++i) { + const std::size_t bit_idx = width - 1 - i; + const Tl bit = + (bit_idx >= base_width) ? ((lhs.hi_ >> (bit_idx - base_width)) & Tl{1}) : ((lhs.lo_ >> bit_idx) & Tl{1}); + remainder = remainder << 1u; + remainder.lo_ |= bit; + quotient = quotient << 1u; + if (!(remainder < rhs)) { + remainder = remainder - rhs; + quotient.lo_ |= Tl{1}; + } + } + return quotient; + } + } + + template + requires(std::numeric_limits::digits <= base_width) + [[nodiscard]] friend constexpr double_width_int operator+(const double_width_int& lhs, Rhs rhs) + { + Th rhi = lhs.hi_; + Tl rlo = lhs.lo_; + if constexpr (std::is_signed_v) { + // sign extension; no matter if lhs is signed, negative rhs sign extend + if (rhs < 0) --rhi; + } + rlo += static_cast(rhs); + if (rlo < lhs.lo_) { + // carry bit + ++rhi; + } + return {rhi, rlo}; + } + + template + [[nodiscard]] friend constexpr double_width_int operator+(Lhs lhs, const double_width_int& rhs) + { + return rhs + lhs; + } + + [[nodiscard]] friend constexpr double_width_int operator+(const double_width_int& lhs, const double_width_int& rhs) + { + const Tl rlo = lhs.lo_ + rhs.lo_; + Th rhi = lhs.hi_ + rhs.hi_; + if (rlo < lhs.lo_) ++rhi; // carry + return {rhi, rlo}; + } + + template + requires(std::numeric_limits::digits <= base_width) + [[nodiscard]] friend constexpr double_width_int operator-(const double_width_int& lhs, Rhs rhs) + { + Th rhi = lhs.hi_; + Tl rlo = lhs.lo_; + if constexpr (std::is_signed_v) { + // sign extension; no matter if lhs is signed, negative rhs sign extend + if (rhs < 0) ++rhi; + } + rlo -= static_cast(rhs); + if (rlo > lhs.lo_) { + // carry bit + --rhi; + } + return {rhi, rlo}; + } + + template + [[nodiscard]] friend constexpr double_width_int operator-(Lhs lhs, const double_width_int& rhs) + { + Th rhi = 0; + Tl rlo = static_cast(lhs); + if constexpr (std::is_signed_v) { + // sign extension; no matter if rhs is signed, negative lhs sign extend + if (lhs < 0) --rhi; + } + rhi -= rhs.hi_; + if (rhs.lo_ > rlo) { + // carry bit + --rhi; + } + rlo -= rhs.lo_; + return {rhi, rlo}; + } + + [[nodiscard]] friend constexpr double_width_int operator-(const double_width_int& lhs, const double_width_int& rhs) + { + const Tl rlo = lhs.lo_ - rhs.lo_; + Th rhi = lhs.hi_ - rhs.hi_; + if (rlo > lhs.lo_) --rhi; // borrow + return {rhi, rlo}; + } + + [[nodiscard]] constexpr double_width_int operator-() const + { + return {(lo_ > 0 ? static_cast(-1) : Th{0}) - hi_, Tl{0} - lo_}; + } + + [[nodiscard]] constexpr double_width_int operator>>(unsigned n) const + { + if (n >= base_width) { + return {static_cast(hi_ < 0 ? -1 : 0), static_cast(hi_ >> (n - base_width))}; + } + return {hi_ >> n, (static_cast(hi_) << (base_width - n)) | (lo_ >> n)}; + } + + [[nodiscard]] constexpr double_width_int operator<<(unsigned n) const + { + if (n >= base_width) { + return {static_cast(lo_ << (n - base_width)), 0}; + } + return {(hi_ << n) + static_cast(lo_ >> (base_width - n)), lo_ << n}; + } + + // Bitwise operators (treat the value as a 2*base_width bit pattern). + [[nodiscard]] constexpr double_width_int operator~() const { return {static_cast(~hi_), static_cast(~lo_)}; } + + [[nodiscard]] friend constexpr double_width_int operator&(const double_width_int& lhs, const double_width_int& rhs) + { + return {static_cast(lhs.hi_ & rhs.hi_), static_cast(lhs.lo_ & rhs.lo_)}; + } + [[nodiscard]] friend constexpr double_width_int operator|(const double_width_int& lhs, const double_width_int& rhs) + { + return {static_cast(lhs.hi_ | rhs.hi_), static_cast(lhs.lo_ | rhs.lo_)}; + } + [[nodiscard]] friend constexpr double_width_int operator^(const double_width_int& lhs, const double_width_int& rhs) + { + return {static_cast(lhs.hi_ ^ rhs.hi_), static_cast(lhs.lo_ ^ rhs.lo_)}; + } + + // Compound assignment operators, defined in terms of the binary forms above. + constexpr double_width_int& operator+=(const double_width_int& rhs) { return *this = *this + rhs; } + constexpr double_width_int& operator-=(const double_width_int& rhs) { return *this = *this - rhs; } + constexpr double_width_int& operator*=(const double_width_int& rhs) { return *this = *this * rhs; } + constexpr double_width_int& operator/=(const double_width_int& rhs) { return *this = *this / rhs; } + constexpr double_width_int& operator&=(const double_width_int& rhs) { return *this = *this & rhs; } + constexpr double_width_int& operator|=(const double_width_int& rhs) { return *this = *this | rhs; } + constexpr double_width_int& operator^=(const double_width_int& rhs) { return *this = *this ^ rhs; } + constexpr double_width_int& operator<<=(unsigned n) { return *this = *this << n; } + constexpr double_width_int& operator>>=(unsigned n) { return *this = *this >> n; } + + // Pre/post increment and decrement. + constexpr double_width_int& operator++() + { + ++lo_; + if (lo_ == Tl{0}) ++hi_; + return *this; + } + constexpr double_width_int operator++(int) + { + auto tmp = *this; + ++*this; + return tmp; + } + constexpr double_width_int& operator--() + { + if (lo_ == Tl{0}) --hi_; + --lo_; + return *this; + } + constexpr double_width_int operator--(int) + { + auto tmp = *this; + --*this; + return tmp; + } + + static constexpr double_width_int max() { return {std::numeric_limits::max(), std::numeric_limits::max()}; } + static constexpr double_width_int min() + { + if constexpr (is_signed) + return {std::numeric_limits::min(), 0}; + else + return {0, 0}; + } + + Th hi_; + Tl lo_; +}; + +template +constexpr std::size_t integer_rep_width_v> = double_width_int::width; + +} // namespace mp_units::detail + +namespace std { + +// Specialize std::numeric_limits for the synthetic double-width integer. Without this, +// generic algorithms that probe ::max(), ::min(), ::digits, etc. (e.g. checked_int_pow, +// magnitude-folding code in unit_magnitude.h) silently see the default-unspecialized +// limits (all zeros) and break. On platforms with native __int128 the std-library +// specializations already exist; this specialization fills the gap on MSVC. +template + requires std::integral +class numeric_limits> { + using dwi = mp_units::detail::double_width_int; +public: + static constexpr bool is_specialized = true; + static constexpr bool is_signed = dwi::is_signed; + static constexpr bool is_integer = true; + static constexpr bool is_exact = true; + static constexpr bool has_infinity = false; + static constexpr bool has_quiet_NaN = false; + static constexpr bool has_signaling_NaN = false; + static constexpr float_round_style round_style = round_toward_zero; + static constexpr bool is_iec559 = false; + static constexpr bool is_bounded = true; + static constexpr bool is_modulo = !dwi::is_signed; + static constexpr int digits = static_cast(dwi::width) - (dwi::is_signed ? 1 : 0); + static constexpr int digits10 = static_cast(digits * 3010 / 10000); // floor(digits * log10(2)) + static constexpr int max_digits10 = 0; + static constexpr int radix = 2; + static constexpr int min_exponent = 0; + static constexpr int min_exponent10 = 0; + static constexpr int max_exponent = 0; + static constexpr int max_exponent10 = 0; + static constexpr bool traps = numeric_limits::traps; + static constexpr bool tinyness_before = false; + + static constexpr dwi min() noexcept { return dwi::min(); } + static constexpr dwi lowest() noexcept { return dwi::min(); } + static constexpr dwi max() noexcept { return dwi::max(); } + static constexpr dwi epsilon() noexcept { return dwi{}; } + static constexpr dwi round_error() noexcept { return dwi{}; } + static constexpr dwi infinity() noexcept { return dwi{}; } + static constexpr dwi quiet_NaN() noexcept { return dwi{}; } + static constexpr dwi signaling_NaN() noexcept { return dwi{}; } + static constexpr dwi denorm_min() noexcept { return dwi{}; } +}; + +} // namespace std diff --git a/src/core/include/mp-units/bits/fixed_point.h b/src/core/include/mp-units/bits/fixed_point.h index e5fa6a2b9..4376797ee 100644 --- a/src/core/include/mp-units/bits/fixed_point.h +++ b/src/core/include/mp-units/bits/fixed_point.h @@ -22,6 +22,7 @@ #pragma once +#include #include // IWYU pragma: keep #include // IWYU pragma: keep @@ -33,276 +34,14 @@ import std; #include #include #include -#include #include -#include #include #include -#include -// becomes partially freestanding in C++26 -// before that, there is no guarantee about this header even existing -// (GCC 14 has it, but actively #error's out) -#if __STDC_HOSTED__ -#include -#endif #endif #endif namespace mp_units::detail { -template -constexpr std::size_t integer_rep_width_v = std::numeric_limits>::digits; - -template -[[nodiscard]] consteval T int_power(T base, int exponent) -{ -#if __STDC_HOSTED__ && defined(__cpp_lib_constexpr_cmath) && __cpp_lib_constexpr_cmath >= 202202L - return std::ldexp(base, exponent); -#else - if (exponent < 0) { - base = T{1} / base; - exponent = -exponent; - } - T ret = 1; - while (exponent) { - if (exponent & 1) ret *= base; - exponent >>= 1; - base *= base; - } - return ret; - -#endif -} - -// this class synthesizes a double-width integer from two base-width integers. -template -struct double_width_int { - static constexpr bool is_signed = std::is_signed_v; - static constexpr std::size_t base_width = integer_rep_width_v; - static constexpr std::size_t width = 2 * base_width; - - using Th = T; - using Tl = std::make_unsigned_t; - - constexpr double_width_int() = default; - -#if !MP_UNITS_COMP_GCC || MP_UNITS_COMP_GCC > 12 -private: -#endif - constexpr double_width_int(Th hi, Tl lo) : hi_(hi), lo_(lo) {} - - friend struct double_width_int, std::make_signed_t>>; - -public: - static constexpr double_width_int from_hi_lo(Th hi, Tl lo) { return {hi, lo}; } - - explicit constexpr double_width_int(long double v) - { - constexpr auto scale = int_power(2, base_width); - constexpr auto iscale = 1.l / scale; - auto scaled = v * iscale; - hi_ = static_cast(scaled); - auto resid = (scaled - static_cast(hi_)); - if (resid < 0) { - --hi_; - resid += 1; - } - lo_ = static_cast(resid * scale); - } - template - requires(is_signed || !std::is_signed_v) - explicit(false) constexpr double_width_int(U v) - { - if constexpr (is_signed) { - hi_ = v < 0 ? Th{-1} : Th{0}; - } else { - hi_ = 0; - } - lo_ = static_cast(v); - } - - template - explicit constexpr operator U() const - { - if constexpr (integer_rep_width_v > base_width) { - return (static_cast(hi_) << base_width) + static_cast(lo_); - } else { - return static_cast(lo_); - } - } - - [[nodiscard]] constexpr auto operator<=>(const double_width_int&) const = default; - - // calculates the double-width product of two base-size integers; this implementation requires at least one of them to - // be unsigned - static constexpr double_width_int wide_product_of(Th lhs, Tl rhs) - { - constexpr std::size_t half_width = base_width / 2; - constexpr Tl msk = (Tl(1) << half_width) - 1u; - Th l1 = lhs >> half_width; - Tl l0 = static_cast(lhs) & msk; - Tl r1 = rhs >> half_width; - Tl r0 = rhs & msk; - Tl t00 = l0 * r0; - Tl t01 = l0 * r1; - Th t10 = l1 * static_cast(r0); - Th t11 = l1 * static_cast(r1); - Tl m = (t01 & msk) + (static_cast(t10) & msk) + (t00 >> half_width); - Th o1 = t11 + static_cast(m >> half_width) + (t10 >> half_width) + static_cast(t01 >> half_width); - Tl o0 = (t00 & msk) | ((m & msk) << half_width); - return {o1, o0}; - } - - template - requires(std::numeric_limits::digits <= base_width) - [[nodiscard]] friend constexpr auto operator*(const double_width_int& lhs, Rhs rhs) - { - using RT = conditional, std::make_signed_t, Tl>; - auto lo_prod = double_width_int::wide_product_of(rhs, lhs.lo_); - // Normal C++ rules; with respect to signedness, the wider type always wins. - using ret_t = double_width_int; - return ret_t{static_cast(lo_prod.hi_) + lhs.hi_ * static_cast(rhs), lo_prod.lo_}; - } - - template - requires(std::numeric_limits::digits <= base_width) - [[nodiscard]] friend constexpr auto operator*(Lhs lhs, const double_width_int& rhs) - { - return rhs * lhs; - } - - template - requires(std::numeric_limits::digits <= base_width) - [[nodiscard]] friend constexpr double_width_int operator/(const double_width_int& lhs, Rhs rhs) - { - // Normal C++ rules; with respect to signedness, the bigger type always wins. - using ret_t = double_width_int; - if constexpr (std::is_signed_v) { - if (rhs < 0) { - return (-lhs) / static_cast(-rhs); - } else { - return lhs / static_cast(rhs); - } - } else if constexpr (is_signed) { - if (lhs.hi_ < 0) { - return -((-lhs) / rhs); - } else { - using unsigned_t = double_width_int; - auto tmp = unsigned_t{static_cast(lhs.hi_), lhs.lo_} / rhs; - return ret_t{static_cast(tmp.hi_), tmp.lo_}; - } - } else { - Th res_hi = lhs.hi_ / rhs; - // unfortunately, wide division is hard: https://en.wikipedia.org/wiki/Division_algorithm. - // Here, we just provide a somewhat naive implementation of long division. - Tl rem_hi = lhs.hi_ % rhs; - Tl rem_lo = lhs.lo_; - Tl res_lo = 0; - for (std::size_t i = 0; i < base_width; ++i) { - // shift in one bit - rem_hi = (rem_hi << 1u) | (rem_lo >> (base_width - 1)); - rem_lo <<= 1u; - res_lo <<= 1u; - // perform one bit of long division - if (rem_hi >= rhs) { - rem_hi -= rhs; - res_lo |= 1u; - } - } - return ret_t{res_hi, res_lo}; - } - } - - template - requires(std::numeric_limits::digits <= base_width) - [[nodiscard]] friend constexpr double_width_int operator+(const double_width_int& lhs, Rhs rhs) - { - Th rhi = lhs.hi_; - Tl rlo = lhs.lo_; - if constexpr (std::is_signed_v) { - // sign extension; no matter if lhs is signed, negative rhs sign extend - if (rhs < 0) --rhi; - } - rlo += static_cast(rhs); - if (rlo < lhs.lo_) { - // carry bit - ++rhi; - } - return {rhi, rlo}; - } - - template - [[nodiscard]] friend constexpr double_width_int operator+(Lhs lhs, const double_width_int& rhs) - { - return rhs + lhs; - } - - template - requires(std::numeric_limits::digits <= base_width) - [[nodiscard]] friend constexpr double_width_int operator-(const double_width_int& lhs, Rhs rhs) - { - Th rhi = lhs.hi_; - Tl rlo = lhs.lo_; - if constexpr (std::is_signed_v) { - // sign extension; no matter if lhs is signed, negative rhs sign extend - if (rhs < 0) ++rhi; - } - rlo -= static_cast(rhs); - if (rlo > lhs.lo_) { - // carry bit - --rhi; - } - return {rhi, rlo}; - } - - template - [[nodiscard]] friend constexpr double_width_int operator-(Lhs lhs, const double_width_int& rhs) - { - Th rhi = 0; - Tl rlo = static_cast(lhs); - if constexpr (std::is_signed_v) { - // sign extension; no matter if rhs is signed, negative lhs sign extend - if (lhs < 0) --rhi; - } - rhi -= rhs.hi_; - if (rhs.lo_ > rlo) { - // carry bit - --rhi; - } - rlo -= rhs.lo_; - return {rhi, rlo}; - } - - [[nodiscard]] constexpr double_width_int operator-() const - { - return {(lo_ > 0 ? static_cast(-1) : Th{0}) - hi_, -lo_}; - } - - [[nodiscard]] constexpr double_width_int operator>>(unsigned n) const - { - if (n >= base_width) { - return {static_cast(hi_ < 0 ? -1 : 0), static_cast(hi_ >> (n - base_width))}; - } - return {hi_ >> n, (static_cast(hi_) << (base_width - n)) | (lo_ >> n)}; - } - - [[nodiscard]] constexpr double_width_int operator<<(unsigned n) const - { - if (n >= base_width) { - return {static_cast(lo_ << (n - base_width)), 0}; - } - return {(hi_ << n) + static_cast(lo_ >> (base_width - n)), lo_ << n}; - } - - static constexpr double_width_int max() { return {std::numeric_limits::max(), std::numeric_limits::max()}; } - -#if !MP_UNITS_COMP_GCC || MP_UNITS_COMP_GCC > 12 -private: -#endif - Th hi_; - Tl lo_; -}; - #if defined(__SIZEOF_INT128__) MP_UNITS_DIAGNOSTIC_PUSH @@ -332,9 +71,6 @@ constexpr std::size_t max_native_width = 64; #endif -template -constexpr std::size_t integer_rep_width_v> = double_width_int::width; - template constexpr bool is_signed_v = std::is_signed_v; diff --git a/test/static/CMakeLists.txt b/test/static/CMakeLists.txt index c46bcfd75..ff3e1ee09 100644 --- a/test/static/CMakeLists.txt +++ b/test/static/CMakeLists.txt @@ -43,6 +43,7 @@ add_library( custom_rep_test_min_impl.cpp dimension_test.cpp dimension_symbol_test.cpp + double_width_int_test.cpp fixed_point_test.cpp fixed_string_test.cpp hep_test.cpp diff --git a/test/static/double_width_int_test.cpp b/test/static/double_width_int_test.cpp new file mode 100644 index 000000000..27589a9f5 --- /dev/null +++ b/test/static/double_width_int_test.cpp @@ -0,0 +1,140 @@ +// The MIT License (MIT) +// +// Copyright (c) 2018 Mateusz Pusz +// +// Permission is hereby granted, free of charge, to any person obtaining a copy +// of this software and associated documentation files (the "Software"), to deal +// in the Software without restriction, including without limitation the rights +// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +// copies of the Software, and to permit persons to whom the Software is +// furnished to do so, subject to the following conditions: +// +// The above copyright notice and this permission notice shall be included in all +// copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +// SOFTWARE. + +#include +#ifdef MP_UNITS_IMPORT_STD +import std; +#else +#include +#include +#include +#endif + +using namespace mp_units; + +namespace { + +using i128 = detail::double_width_int; +using u128 = detail::double_width_int; +using i32x2 = detail::double_width_int; +using u32x2 = detail::double_width_int; + +// Shift/divide round-trip on i128 (exercises the narrow-rhs paths) +static_assert((((83 * 79 * 73) * (i128{97} << 64u) / 89) >> 64u) == (83 * 79 * 73 * 97) / 89); + +// ---- Arithmetic between two double-width operands ---- +// These compile-time checks pin down the operators that the synthetic int128 emulation must +// provide on platforms without native __int128 (issue #795). + +// + with carry from the low half to the high half +static_assert(u128{1, 2} + u128{3, 4} == u128{4, 6}); +static_assert(u128{0, ~std::uint64_t{0}} + u128{0, 1} == u128{1, 0}); +static_assert(i128{0, ~std::uint64_t{0}} + i128{0, 1} == i128{1, 0}); + +// - with borrow from the high half to the low half +static_assert(u128{4, 6} - u128{3, 4} == u128{1, 2}); +static_assert(u128{1, 0} - u128{0, 1} == u128{0, ~std::uint64_t{0}}); +static_assert(-i128{0, 1} == i128{-1, ~std::uint64_t{0}}); + +// * truncated to width bits; the high*high contribution overflows and is discarded +static_assert(u128{0, 5} * u128{0, 3} == u128{0, 15}); +static_assert(u128{0, ~std::uint64_t{0}} * u128{0, 2} == u128{1, ~std::uint64_t{0} - 1}); // 2 * (2^64 - 1) +static_assert(u128{0, 1u << 16u} * u128{0, 1u << 16u} == u128{0, std::uint64_t{1} << 32u}); +static_assert(i128{0, 5} * i128{-1, ~std::uint64_t{0} - 2} == i128{-1, ~std::uint64_t{0} - 14}); // 5 * -3 = -15 + +// / between two dwints — narrow fast path and wide long-division path +static_assert(u128{0, 100} / u128{0, 7} == u128{0, 14}); +static_assert(u128{1, 0} / u128{0, 2} == u128{0, std::uint64_t{1} << 63u}); // 2^64 / 2 == 2^63 +static_assert(i128{0, 100} / i128{-1, ~std::uint64_t{0} - 6} == i128{-1, ~std::uint64_t{0} - 13}); // 100 / -7 == -14 +static_assert(u128{1, 5} / u128{2, 0} == u128{0, 0}); // divisor wider than dividend → 0 + +// Bitwise operators +static_assert(~u128{0, 0} == u128{~std::uint64_t{0}, ~std::uint64_t{0}}); +static_assert(~u128{~std::uint64_t{0}, ~std::uint64_t{0}} == u128{0, 0}); +static_assert((u128{0xFF, 0xFF} & u128{0x0F, 0xF0}) == u128{0x0F, 0xF0}); +static_assert((u128{0xF0, 0x00} | u128{0x0F, 0xFF}) == u128{0xFF, 0xFF}); +static_assert((u128{0xFF, 0x00} ^ u128{0x0F, 0xFF}) == u128{0xF0, 0xFF}); + +// Compound assignment operators (each verified by reading back the value after one operation) +constexpr auto check_compound = [] { + u128 v{1, 2}; + v += u128{2, 3}; + if (!(v == u128{3, 5})) return false; + v -= u128{1, 1}; + if (!(v == u128{2, 4})) return false; + v *= u128{0, 3}; + if (!(v == u128{6, 12})) return false; + v /= u128{0, 2}; + if (!(v == u128{3, 6})) return false; + v &= u128{~std::uint64_t{0}, 0xF}; + if (!(v == u128{3, 6})) return false; + v |= u128{0, 0xF0}; + if (!(v == u128{3, 0xF6})) return false; + v ^= u128{0, 0xFF}; + if (!(v == u128{3, 9})) return false; + v <<= 4u; + if (!(v == u128{0x30, 0x90})) return false; + v >>= 4u; + if (!(v == u128{3, 9})) return false; + return true; +}; +static_assert(check_compound()); + +// Pre/post increment and decrement +constexpr auto check_inc_dec = [] { + u128 v{0, ~std::uint64_t{0}}; + ++v; + if (!(v == u128{1, 0})) return false; // carry across halves + --v; + if (!(v == u128{0, ~std::uint64_t{0}})) return false; + auto post_inc = v++; + if (!(post_inc == u128{0, ~std::uint64_t{0}} && v == u128{1, 0})) return false; + auto post_dec = v--; + if (!(post_dec == u128{1, 0} && v == u128{0, ~std::uint64_t{0}})) return false; + return true; +}; +static_assert(check_inc_dec()); + +// numeric_limits specialization (required by checked_int_pow etc.) +static_assert(std::numeric_limits::is_specialized); +static_assert(std::numeric_limits::is_integer); +static_assert(std::numeric_limits::is_exact); +static_assert(!std::numeric_limits::is_signed); +static_assert(std::numeric_limits::is_signed); +static_assert(std::numeric_limits::digits == 128); +static_assert(std::numeric_limits::digits == 127); +static_assert(std::numeric_limits::digits == 64); +static_assert(std::numeric_limits::digits == 63); +static_assert(std::numeric_limits::radix == 2); +static_assert(std::numeric_limits::min() == u128{0, 0}); +static_assert(std::numeric_limits::min() == i128{std::numeric_limits::min(), 0}); +static_assert(std::numeric_limits::max() == u128{~std::uint64_t{0}, ~std::uint64_t{0}}); +static_assert(std::numeric_limits::max() == i128{std::numeric_limits::max(), ~std::uint64_t{0}}); + +// Conversion to long double (round-trips small integers exactly; the signed branch goes +// through unary minus to avoid catastrophic cancellation on platforms where long double +// has fewer than `width` bits of mantissa, e.g. MSVC where long double == double). +static_assert(static_cast(u128{0, 12345}) == 12345.0L); +static_assert(static_cast(i128{0, 1}) == 1.0L); +static_assert(static_cast(-i128{0, 1}) == -1.0L); + +} // namespace diff --git a/test/static/fixed_point_test.cpp b/test/static/fixed_point_test.cpp index fa3a4f261..e31190e7d 100644 --- a/test/static/fixed_point_test.cpp +++ b/test/static/fixed_point_test.cpp @@ -35,6 +35,7 @@ using namespace mp_units; namespace { +// min_width_uint_t selects the narrowest standard unsigned type that holds N bits static_assert(std::is_same_v, std::uint8_t>); static_assert(std::is_same_v, std::uint8_t>); static_assert(std::is_same_v, std::uint8_t>); @@ -43,11 +44,6 @@ static_assert(std::is_same_v, std::uint32_t>); static_assert(std::is_same_v, std::uint32_t>); static_assert(std::is_same_v, std::uint64_t>); -using i128 = detail::double_width_int; -using u128 = detail::double_width_int; - -static_assert((((83 * 79 * 73) * (i128{97} << 64u) / 89) >> 64u) == (83 * 79 * 73 * 97) / 89); - // scale(M{}, value) — integer-to-integer path (exact arithmetic, no floating point) // integral factor: exact integer multiply