Files
mp-units/test/runtime/fixed_point_test.cpp
T
Mateusz Pusz efbc844199 fix: fixed-point arithmetic for integer unit conversions (#580) (#764)
* Fix #580: use fixed-point arithmetic for integer unit conversions

Introduce a fixed-point implementation for unit conversions involving
integer representations, avoiding loss of significant digits that
previously occurred when the conversion factor was not a whole number.

New files:
- src/core/include/mp-units/bits/fixed_point.h: double_width_int<T> and
  fixed_point<T,n> types for exact rational scaling of integer values.
  Uses __int128 when available (__SIZEOF_INT128__) for 64-bit integers.
- src/core/include/mp-units/framework/scaling.h: public scaling_traits<>
  customization point and scale<To>(M, value) free function. Provides
  built-in specializations for floating-point and integer-like types.
- test/static/fixed_point_test.cpp: static assertions for the new types.
- test/runtime/fixed_point_test.cpp: runtime arithmetic edge-case tests.

Modified:
- sudo_cast.h: replace hand-rolled conversion_value_traits / sudo_cast_value
  machinery with a single scale<To::rep>(c_mag, ...) call.
- representation_concepts.h: add MagnitudeScalable concept; replace
  ComplexScalar with HasComplexOperations (which is its definition).
- customization_points.h: add unspecified_rep tag and declare the primary
  scaling_traits<> template.
- framework.h / CMakeLists.txt: wire in the new headers.
- hacks.h: add MP_UNITS_DIAGNOSTIC_IGNORE_PEDANTIC and
  MP_UNITS_DIAGNOSTIC_IGNORE_SIGN_CONVERSION macros.
- example/measurement.cpp: add scaling_traits specializations for
  measurement<T> to demonstrate the customization point.
- test/static/{international,usc}_test.cpp: disable two tests that are
  blocked on issue #614.

Co-authored-by: Tobias Hanhart <burnpanck@users.noreply.github.com>

* Fix value_Type typo in floating_point_scaling_factor_type specialization

The partial specialization for types with a nested value_type used
'value_Type' (capital T) instead of 'value_type', making the entire
specialization dead code as the requires-clause could never be satisfied.

Also fix 'mantiassa' -> 'mantissa' in the adjacent comment.

* Fix docstring typos in scaling_traits documentation

- 'quantitiy' -> 'quantity'
- 'dictatet' -> 'dictated'
- 'convetrible' -> 'convertible'
- 'implemenation' -> 'implementation'
- 'availabe' -> 'available'

* Fix conflict resolution error: keep ComplexScalar name from master

When resolving the merge conflict in representation_concepts.h, the
PR's renamed version of the concept ('HasComplexOperations') was used
instead of master's established name ('ComplexScalar'). The two concepts
are semantically equivalent — burnpanck simply renamed it in his branch.

Revert to the canonical 'ComplexScalar' name while retaining the new
'MagnitudeScalable' concept which was the actual addition from the PR.

* Fix measurement.cpp: remove duplicate class definition from merge

The PR branched from a version where measurement<T> was defined inline
in measurement.cpp. Master later moved the class to example/include/
measurement.h and changed measurement.cpp to #include that header.

The squash merge therefore introduced a duplicate definition: the class
from the header and the PR's inline class were both visible, causing
an 'ambiguous reference' error. Remove the now-redundant inline class;
the scaling_traits specializations added by the PR work correctly with
the class from measurement.h.

* style: pre-commit

* docs: chapters anchors improved in the "custom representation" chapter

* docs: value conversions chapter improved

* refactor: scaling support refactored

* fix: clang-16 crash fixed

* docs: `measurement` example documentation updated to match changes

* fix: use exact wide-integer arithmetic for rational unit conversions on all platforms

On ARM / Apple Silicon, long double == double (64-bit mantissa).  The old
fixed_point<T>(long double) initialiser lost ~12 bits of precision for 64-bit
integer types when representing the scaling ratio, producing an error of ~49
units for the 10/9 (degree → gradian) conversion with a 10^18 input value.

Fix by splitting the integer-path else-branch into two cases:

  • Pure rational M (is_integral(M * (denominator(M) / numerator(M))) == true):
    use (value * numerator) / denominator via double_width_int_for_t<> arithmetic.
    This is exact on every platform regardless of long double width.

  • Irrational M (involves π etc.): keep the long double fixed_point approximation.
    These conversions are inherently approximate; small values still produce correct
    truncated results on all platforms.

Update the test comment to reflect the new exact-arithmetic path.

Fixes CI failures on clang-18/ARM and apple-clang-16.

* fix: replace floating-point TeX-point test with exact integer equivalent

72.27 is not exactly representable as double (it rounds to 72.2699...96).
Multiplying by the conversion factor 100/7227 via long double gives a result
≥ 1.0 on x86 (80-bit long double, 64-bit mantissa) only by chance, but
0.99999...978 on ARM / Apple Silicon where long double == double (52-bit).

The correct mathematical statement is: 7227 tex_point = 100 inch (exact
rational relationship).  Use that integer form instead of the inexact 72.27
double literal so the test is correct and platform-independent.

---------

Co-authored-by: Tobias Hanhart <burnpanck@users.noreply.github.com>
2026-03-07 21:02:37 +01:00

239 lines
8.1 KiB
C++
Raw Blame History

This file contains invisible Unicode characters
This file contains invisible Unicode characters that are indistinguishable to humans but may be processed differently by a computer. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.
// 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 <catch2/catch_test_macros.hpp>
#include <catch2/matchers/catch_matchers.hpp>
#include <mp-units/bits/fixed_point.h>
#ifdef MP_UNITS_IMPORT_STD
import std;
#else
#include <array>
#include <cstdint>
#include <limits>
#include <tuple>
#include <vector>
#endif
using namespace mp_units;
using namespace mp_units::detail;
template<std::size_t N, typename... T, std::size_t... I>
requires(N == sizeof...(T) && N == sizeof...(I))
std::tuple<T...> at(const std::array<std::size_t, N>& idx, std::integer_sequence<std::size_t, I...>,
const std::vector<T>&... src)
{
return {src[idx[I]]...};
}
template<typename... T>
std::vector<std::tuple<T...>> cartesian_product(const std::vector<T>&... src)
{
std::vector<std::tuple<T...>> ret;
constexpr std::size_t N = sizeof...(src);
std::array<std::size_t, N> sizes;
{
std::size_t n = 1;
std::size_t k = 0;
for (std::size_t s : {src.size()...}) {
sizes[k++] = s;
n *= s;
}
ret.reserve(n);
}
std::array<std::size_t, N> idx = {};
bool done = false;
while (!done) {
ret.push_back(at(idx, std::make_index_sequence<N>{}, src...));
for (std::size_t k = 0; k < idx.size(); ++k) {
if (++idx[k] < sizes[k]) break;
if (k + 1 >= idx.size()) {
done = true;
break;
}
idx[k] = 0;
}
}
return ret;
}
template<std::integral T>
using half_width_int_for_t = std::conditional_t<std::is_signed_v<T>, min_width_int_t<integer_rep_width_v<T> / 2>,
min_width_uint_t<integer_rep_width_v<T> / 2>>;
template<std::integral Hi, std::unsigned_integral Lo>
requires(integer_rep_width_v<Hi> == integer_rep_width_v<Lo>)
auto combine_bits(Hi hi, Lo lo)
{
using ret_t = double_width_int_for_t<Hi>;
return (static_cast<ret_t>(hi) << integer_rep_width_v<Lo>)+static_cast<ret_t>(lo);
}
template<std::integral T, typename V>
void check(double_width_int<T> value, V&& visitor)
{
using DT = double_width_int_for_t<T>;
auto as_standard_int = static_cast<DT>(value);
auto expected = visitor(as_standard_int);
auto actual = visitor(value);
auto actual_as_standard = static_cast<DT>(actual);
REQUIRE(actual_as_standard == expected);
}
// Produce some test integers in the vicinity (~ +-1, modulo overflow)
// of those areas in their representation at risk of causing problems:
// intmin,intmin/2,zero,intmax/2,intmax...
template<std::integral T>
std::vector<T> test_values()
{
using U = std::make_unsigned_t<T>;
std::vector<T> ret;
for (int msb : {0, 1, 2, 3}) {
// vicinities reached from msb=:
// 0: signed: zero; unsigned: zero, intmin and intmax
// 1: signed: intmax/2
// 2: unsigned: intmax/2; signed: intmin and intmax
// 3: signed: intmin/2
auto ref = static_cast<U>(msb) << (integer_rep_width_v<U> - 2);
for (int lsb_corr : {-2, -1, 0, 1, 2}) {
auto corr = static_cast<U>(lsb_corr);
U value = ref + corr;
ret.push_back(static_cast<T>(value));
}
}
return ret;
}
using u32 = std::uint32_t;
using i32 = std::int32_t;
using u64 = std::uint64_t;
using i64 = std::int64_t;
using du32 = double_width_int<u32>;
using di32 = double_width_int<i32>;
MP_UNITS_DIAGNOSTIC_PUSH
// double_width_int implements the same sign-conversion rules as the standard int types, and we want to verify that;
// even if those sign-conversion rules are frowned upon.
MP_UNITS_DIAGNOSTIC_IGNORE_SIGN_CONVERSION
TEST_CASE("double_width_int addition and subtraction", "[double_width_int]")
{
SECTION("u32x2 +/- u32")
{
for (auto [lhi, llo, rhs] : cartesian_product(test_values<u32>(), test_values<u32>(), test_values<u32>())) {
CAPTURE(lhi, llo, rhs);
auto lhs = double_width_int<u32>::from_hi_lo(lhi, llo);
check(lhs, [r = rhs](auto v) { return v + r; });
check(lhs, [r = rhs](auto v) { return v - r; });
check(lhs, [r = rhs](auto v) { return r - v; });
}
}
SECTION("u32x2 +/- i32")
{
for (auto [lhi, llo, rhs] : cartesian_product(test_values<u32>(), test_values<u32>(), test_values<i32>())) {
CAPTURE(lhi, llo, rhs);
auto lhs = double_width_int<u32>::from_hi_lo(lhi, llo);
check(lhs, [r = rhs](auto v) { return v + r; });
check(lhs, [r = rhs](auto v) { return v - r; });
check(lhs, [r = rhs](auto v) { return r - v; });
}
}
SECTION("i32x2 +/- u32")
{
for (auto [lhi, llo, rhs] : cartesian_product(test_values<i32>(), test_values<u32>(), test_values<u32>())) {
CAPTURE(lhi, llo, rhs);
auto lhs = double_width_int<i32>::from_hi_lo(lhi, llo);
check(lhs, [r = rhs](auto v) { return v + r; });
check(lhs, [r = rhs](auto v) { return v - r; });
check(lhs, [r = rhs](auto v) { return r - v; });
}
}
SECTION("i32x2 +/- i32")
{
for (auto [lhi, llo, rhs] : cartesian_product(test_values<i32>(), test_values<u32>(), test_values<i32>())) {
CAPTURE(lhi, llo, rhs);
auto lhs = double_width_int<i32>::from_hi_lo(lhi, llo);
check(lhs, [r = rhs](auto v) { return v + r; });
check(lhs, [r = rhs](auto v) { return v - r; });
check(lhs, [r = rhs](auto v) { return r - v; });
}
}
}
TEST_CASE("double_width_int multiplication", "[double_width_int]")
{
SECTION("u32 * u32")
{
for (auto [lhs, rhs] : cartesian_product(test_values<u32>(), test_values<u32>())) {
CAPTURE(lhs, rhs);
u64 expected = u64{lhs} * u64{rhs};
auto actual = double_width_int<u32>::wide_product_of(lhs, rhs);
auto actual_as_std = static_cast<u64>(actual);
REQUIRE(actual_as_std == expected);
}
}
SECTION("i32 * u32")
{
for (auto [lhs, rhs] : cartesian_product(test_values<i32>(), test_values<u32>())) {
CAPTURE(lhs, rhs);
i64 expected = i64{lhs} * i64{rhs};
auto actual = double_width_int<i32>::wide_product_of(lhs, rhs);
auto actual_as_std = static_cast<i64>(actual);
REQUIRE(actual_as_std == expected);
}
}
SECTION("u32x2 * u32")
{
for (auto [lhi, llo, rhs] : cartesian_product(test_values<u32>(), test_values<u32>(), test_values<u32>())) {
CAPTURE(lhi, llo, rhs);
auto lhs = double_width_int<u32>::from_hi_lo(lhi, llo);
check(lhs, [r = rhs](auto v) { return v * r; });
}
}
SECTION("u32x2 * i32")
{
for (auto [lhi, llo, rhs] : cartesian_product(test_values<u32>(), test_values<u32>(), test_values<i32>())) {
CAPTURE(lhi, llo, rhs);
auto lhs = double_width_int<u32>::from_hi_lo(lhi, llo);
check(lhs, [r = rhs](auto v) { return v * r; });
}
}
SECTION("i32x2 * u32")
{
for (auto [lhi, llo, rhs] : cartesian_product(test_values<i32>(), test_values<u32>(), test_values<u32>())) {
CAPTURE(lhi, llo, rhs);
auto lhs = double_width_int<i32>::from_hi_lo(lhi, llo);
check(lhs, [r = rhs](auto v) { return v * r; });
}
}
SECTION("i32x2 * i32")
{
for (auto [lhi, llo, rhs] : cartesian_product(test_values<i32>(), test_values<u32>(), test_values<i32>())) {
CAPTURE(lhi, llo, rhs);
auto lhs = double_width_int<i32>::from_hi_lo(lhi, llo);
check(lhs, [r = rhs](auto v) { return v * r; });
}
}
}
MP_UNITS_DIAGNOSTIC_POP