Skip to content

Commit

Permalink
wip
Browse files Browse the repository at this point in the history
  • Loading branch information
burnpanck committed Sep 15, 2024
1 parent 9ca56f6 commit cf9c05f
Show file tree
Hide file tree
Showing 3 changed files with 298 additions and 40 deletions.
224 changes: 224 additions & 0 deletions src/core/include/mp-units/bits/fixed_point.h
Original file line number Diff line number Diff line change
@@ -0,0 +1,224 @@
// 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

// IWYU pragma: private, include <mp-units/framework.h>
#include <mp-units/framework/magnitude.h>

#ifndef MP_UNITS_IN_MODULE_INTERFACE
#ifdef MP_UNITS_IMPORT_STD
import std;
#else
#include <bit>
#include <concepts>
#include <cstdint>
#include <cstdlib>
#include <limits>
#include <numbers>
#endif
#endif

namespace mp_units {
namespace detail {

// this class synthesizes a double-width integer from two base-width integers.
template<std::integral T>
struct double_width_int {
static constexpr bool is_signed = std::is_signed_v<T>;
static constexpr std::size_t base_width = std::numeric_limits<std::make_unsigned_t<T>>::digits;
static constexpr std::size_t width = 2 * base_width;

using Th = T;
using Tl = std::make_unsigned_t<T>;

constexpr double_width_int() = default;
constexpr double_width_int(const double_width_int&) = default;

constexpr double_width_int& operator=(const double_width_int&) = default;

constexpr double_width_int(Th hi, Tl lo) : hi(hi), lo(lo) {}

explicit(true) constexpr double_width_int(long double v)
{
constexpr auto scale = int_power<long double>(2, base_width);
constexpr auto iscale = 1.l / scale;
hi = static_cast<Th>(v * iscale);
lo = static_cast<Tl>(v - (hi * scale));
}
template<std::integral U>
requires(is_signed || !std::is_signed_v<U>)
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<Tl>(v);
}

explicit(true) constexpr operator Th() const { return static_cast<Th>(lo); }

constexpr auto operator<=>(const double_width_int&) const = default;

// calculates the double-width product of two base-size integers
static constexpr double_width_int wide_product_of(Th lhs, Tl rhs)
{
static constexpr std::size_t half_width = base_width / 2;
static constexpr Tl msk = Tl(1) << half_width;
Th l1 = lhs >> half_width;
Tl l0 = static_cast<Tl>(lhs) & msk;
Tl r1 = rhs >> half_width;
Tl r0 = rhs & msk;
Tl t00 = l0 * r0;
Tl t01 = l0 * r1;
Th t10 = l1 * r0;
Th t11 = l1 * r1;
Tl m = (t01 & msk) + (static_cast<Tl>(t10) & msk) + (t00 >> half_width);
Th o1 = t11 + (m >> half_width) + (t10 >> half_width) + static_cast<Th>(t01 >> half_width);
Tl o0 = (t00 & msk) | ((m & msk) << half_width);
return {o1, o0};
}

template<std::integral Lhs>
friend constexpr auto operator*(Lhs lhs, const double_width_int& rhs)
{
using ret_t = double_width_int<std::common_type_t<Lhs, Th>>;
auto ret = ret_t::wide_product_of(lhs, rhs.lo);
ret.hi += lhs * rhs.hi;
return ret;
};


constexpr double_width_int operator+(Tl rhs) const
{
auto ret = lo + rhs;
return {hi + (ret < lo ? 1 : 0), ret};
}

constexpr double_width_int operator>>(unsigned n) const
{
if (n >= base_width) {
return {0, hi >> (n - base_width)};
}
return {hi >> n, (static_cast<Tl>(hi) << (base_width - n)) | (lo >> n)};
}
constexpr double_width_int operator<<(unsigned n) const
{
if (n >= base_width) {
return {static_cast<Th>(lo >> (2 * base_width - n)), 0};
}
return {(hi << n) + static_cast<Th>(lo >> (base_width - n)), lo << n};
}

static constexpr double_width_int max() { return {std::numeric_limits<Th>::max(), std::numeric_limits<Tl>::max()}; }

Th hi;
Tl lo;
};

#ifdef __SIZEOF_INT128__
using int128_t = __int128;
using uint128_t = unsigned __int128;
inline constexpr std::size_t max_native_width = 128;
#else
using int128_t = double_width_int<std::int64_t>;
using uint128_t = double_width_int<std::uint64_t>;
inline constexpr std::size_t max_native_width = 64;
#endif

template<typename T>
inline constexpr std::size_t integer_rep_width_v = std::numeric_limits<std::make_unsigned_t<T>>::digits;
template<typename T>
inline constexpr std::size_t integer_rep_width_v<double_width_int<T>> = double_width_int<T>::width;

template<typename T>
inline constexpr bool is_signed_v = std::is_signed_v<T>;
template<typename T>
inline constexpr bool is_signed_v<double_width_int<T>> = double_width_int<T>::is_signed;

template<typename T>
using make_signed_t = std::make_signed_t<T>;

template<std::size_t N>
using min_width_uint_t =
std::tuple_element_t<std::bit_width(N) - 4 + (N > (1u << (std::bit_width(N) - 1)) ? 1 : 0),
std::tuple<std::uint8_t, std::uint16_t, std::uint32_t, std::uint64_t, uint128_t>>;

static_assert(std::is_same_v<min_width_uint_t<31>, std::uint32_t>);
static_assert(std::is_same_v<min_width_uint_t<32>, std::uint32_t>);
static_assert(std::is_same_v<min_width_uint_t<33>, std::uint64_t>);

template<std::size_t N>
using min_width_int_t = make_signed_t<min_width_uint_t<N>>;

template<typename T>
using double_width_int_for_t = std::conditional_t<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<typename Lhs, typename Rhs>
constexpr auto wide_product_of(Lhs lhs, Rhs rhs)
{
if constexpr (integer_rep_width_v<Lhs> + integer_rep_width_v<Rhs> <= max_native_width) {
using T = std::common_type_t<double_width_int_for_t<Lhs>, double_width_int_for_t<Rhs>>;
return static_cast<T>(lhs) * static_cast<T>(rhs);
} else {
using T = double_width_int<std::common_type_t<Lhs, Rhs>>;
return T::wide_product_of(lhs, rhs);
}
}

// This class represents rational numbers using a fixed-point representation, with a symmetric number of digits (bits)
// on either side of the decimal point. The template argument `T` specifies the range of the integral part,
// thus this class uses twice as many bits as the provided type, but is able to precisely store exactly all integers
// from the declared type, as well as efficiently describe all rational factors that can be applied to that type
// and neither always cause underflow or overflow.
template<std::integral T>
struct fixed_point {
using repr_t = double_width_int_for_t<T>;
static constexpr std::size_t fractional_bits = integer_rep_width_v<T>;

constexpr fixed_point() = default;
constexpr fixed_point(const fixed_point&) = default;

constexpr fixed_point& operator=(const fixed_point&) = default;

explicit constexpr fixed_point(long double v) :
int_repr_is_an_implementation_detail_(static_cast<repr_t>(v * int_power<long double>(2, fractional_bits)))
{
}

template<std::integral U>
requires(integer_rep_width_v<U> <= integer_rep_width_v<T>)
auto scale(U v) const
{
auto res = v * int_repr_is_an_implementation_detail_;
return static_cast<std::conditional_t<is_signed_v<decltype((res))>, std::make_signed_t<U>, U>>(res >>
fractional_bits);
}

repr_t int_repr_is_an_implementation_detail_;
};

} // namespace detail
} // namespace mp_units
109 changes: 69 additions & 40 deletions src/core/include/mp-units/bits/sudo_cast.h
Original file line number Diff line number Diff line change
Expand Up @@ -22,6 +22,7 @@

#pragma once

#include <mp-units/bits/fixed_point.h>
#include <mp-units/ext/type_traits.h>
#include <mp-units/framework/magnitude.h>
#include <mp-units/framework/quantity_concepts.h>
Expand Down Expand Up @@ -53,14 +54,15 @@ template<Magnitude auto M, typename Rep1, typename Rep2>
struct conversion_type_traits {
using c_rep_type = maybe_common_type<Rep1, Rep2>;
using c_mag_type = common_magnitude_type<M>;
using multiplier_type = conditional<
treat_as_floating_point<c_rep_type>,
// ensure that the multiplier is also floating-point
conditional<std::is_arithmetic_v<value_type_t<c_rep_type>>,
// reuse user's type if possible
std::common_type_t<c_mag_type, value_type_t<c_rep_type>>, std::common_type_t<c_mag_type, double>>,
c_mag_type>;
using c_type = maybe_common_type<c_rep_type, multiplier_type>;
/* using multiplier_type = conditional<
treat_as_floating_point<c_rep_type>,
// ensure that the multiplier is also floating-point
conditional<std::is_arithmetic_v<value_type_t<c_rep_type>>,
// reuse user's type if possible
std::common_type_t<c_mag_type, value_type_t<c_rep_type>>, std::common_type_t<c_mag_type, double>>,
c_mag_type>;*/
using c_type = conditional<std::is_arithmetic_v<value_type_t<c_rep_type>>, value_type_t<c_rep_type>,
std::common_type_t<c_mag_type, double>>;
};

/**
Expand All @@ -79,10 +81,56 @@ struct conversion_value_traits {
static constexpr Magnitude auto num = numerator(M);
static constexpr Magnitude auto den = denominator(M);
static constexpr Magnitude auto irr = M * (den / num);
static constexpr auto ratio = [] {
if constexpr (std::is_integral_v<T>) {
using U = long double;
return detail::fixed_point<T>{get_value<U>(num) / get_value<U>(den) * get_value<U>(irr)};
} else {
return get_value<T>(num) / get_value<T>(den) * get_value<T>(irr);
}
}();
static constexpr bool value_increases = ratio >= T{1};

template<typename V>
static constexpr auto scale(V value)
{
if constexpr (std::is_integral_v<T>) {
return ratio.scale(value);
} else {
return value * ratio;
}
}
};

template<Magnitude auto M, typename T>
requires(is_integral(M))
struct conversion_value_traits<M, T> {
static constexpr Magnitude auto num = numerator(M);
static constexpr T num_mult = get_value<T>(num);
static constexpr T den_mult = get_value<T>(den);
static constexpr T irr_mult = get_value<T>(irr);
static constexpr T ratio = num_mult / den_mult * irr_mult;
static constexpr bool value_increases = true;

static_assert(get_value<T>(denominator(M)) == 1);
static_assert(is_integral(M));

template<typename V>
static constexpr auto scale(V value)
{
return value * num_mult;
}
};

template<Magnitude auto M, typename T>
requires(is_integral(pow<-1>(M)) && !is_integral(M))
struct conversion_value_traits<M, T> {
static constexpr Magnitude auto den = denominator(M);
static constexpr T den_div = get_value<T>(den);
static constexpr bool value_increases = false;

template<typename V>
static constexpr auto scale(V value)
{
return value / den_div;
}
};


Expand Down Expand Up @@ -110,30 +158,11 @@ template<Quantity To, typename FwdFrom, Quantity From = std::remove_cvref_t<FwdF
} else {
constexpr Magnitude auto c_mag = get_canonical_unit(From::unit).mag / get_canonical_unit(To::unit).mag;
using type_traits = conversion_type_traits<c_mag, typename From::rep, typename To::rep>;
using multiplier_type = typename type_traits::multiplier_type;
// TODO the below crashed nearly every compiler I tried it on
// auto scale = [&](std::invocable<typename type_traits::c_type> auto func) {
auto scale = [&](auto func) {
auto res =
static_cast<To::rep>(func(static_cast<type_traits::c_type>(q.numerical_value_is_an_implementation_detail_)));
return To{res, To::reference};
};

// scale the number
if constexpr (is_integral(c_mag))
return scale([&](auto value) { return value * get_value<multiplier_type>(numerator(c_mag)); });
else if constexpr (is_integral(pow<-1>(c_mag)))
return scale([&](auto value) { return value / get_value<multiplier_type>(denominator(c_mag)); });
else {
using value_traits = conversion_value_traits<c_mag, multiplier_type>;
if constexpr (std::is_floating_point_v<multiplier_type>)
// this results in great assembly
return scale([](auto value) { return value * value_traits::ratio; });
else
// this is slower but allows conversions like 2000 m -> 2 km without loosing data
return scale(
[](auto value) { return value * value_traits::num_mult / value_traits::den_mult * value_traits::irr_mult; });
}
using value_traits = conversion_value_traits<c_mag, typename type_traits::c_type>;

auto res = static_cast<To::rep>(
value_traits::scale(static_cast<type_traits::c_rep_type>(q.numerical_value_is_an_implementation_detail_)));
return To{res, To::reference};
}
}

Expand Down Expand Up @@ -172,21 +201,21 @@ template<QuantityPoint ToQP, typename FwdFromQP, QuantityPoint FromQP = std::rem
// and target unit and representation.
constexpr Magnitude auto c_mag = get_canonical_unit(FromQP::unit).mag / get_canonical_unit(ToQP::unit).mag;
using type_traits = conversion_type_traits<c_mag, typename FromQP::rep, typename ToQP::rep>;
using value_traits = conversion_value_traits<c_mag, typename type_traits::multiplier_type>;
using c_rep_type = typename type_traits::c_rep_type;
if constexpr (value_traits::num_mult * value_traits::irr_mult > value_traits::den_mult) {
using c_type = type_traits::c_type;
using value_traits = conversion_value_traits<c_mag, c_type>;
if constexpr (value_traits::value_increases) {
// original unit had a larger unit magnitude; if we first convert to the common representation but retain the
// unit, we obtain the largest possible range while not causing truncation of fractional values. This is optimal
// for the offset computation.
return sudo_cast<ToQP>(
sudo_cast<quantity_point<FromQP::reference, FromQP::point_origin, c_rep_type>>(std::forward<FwdFromQP>(qp))
sudo_cast<quantity_point<FromQP::reference, FromQP::point_origin, c_type>>(std::forward<FwdFromQP>(qp))
.point_for(ToQP::point_origin));
} else {
// new unit may have a larger unit magnitude; we first need to convert to the new unit (potentially causing
// truncation, but no more than if we did the conversion later), but make sure we keep the larger of the two
// representation types. Then, we can perform the offset computation.
return sudo_cast<ToQP>(
sudo_cast<quantity_point<make_reference(FromQP::quantity_spec, ToQP::unit), FromQP::point_origin, c_rep_type>>(
sudo_cast<quantity_point<make_reference(FromQP::quantity_spec, ToQP::unit), FromQP::point_origin, c_type>>(
std::forward<FwdFromQP>(qp))
.point_for(ToQP::point_origin));
}
Expand Down
5 changes: 5 additions & 0 deletions test/static/international_test.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -37,6 +37,11 @@ using namespace mp_units::international;
using namespace mp_units::international::unit_symbols;

// Mass
constexpr Magnitude auto c_mag = get_canonical_unit(lb).mag / get_canonical_unit(si::kilogram).mag;
static_assert(get_value<int>(detail::numerator(c_mag)) == 45'359'237);
static_assert(get_value<int>(detail::denominator(c_mag)) == 100'000'000);
static_assert(!is_integral(c_mag));

static_assert(100'000'000 * isq::mass[lb] == 45'359'237 * isq::mass[si::kilogram]);
static_assert(1 * isq::mass[lb] == 16 * isq::mass[oz]);
static_assert(1 * isq::mass[oz] == 16 * isq::mass[dr]);
Expand Down

0 comments on commit cf9c05f

Please sign in to comment.