From 06ed8fec8b9066ce0497eb91bc75b4785e3814b9 Mon Sep 17 00:00:00 2001 From: Roland Conybeare Date: Sat, 27 Apr 2024 09:57:29 -0400 Subject: [PATCH] xo-unit: use Quantity.scale repr for outer scale to avoid overflow --- include/xo/unit/Quantity.hpp | 6 +-- include/xo/unit/natural_unit.hpp | 88 ++++++++++++++++++-------------- include/xo/unit/scaled_unit.hpp | 41 +++++++++------ utest/Quantity.test.cpp | 8 +-- 4 files changed, 83 insertions(+), 60 deletions(-) diff --git a/include/xo/unit/Quantity.hpp b/include/xo/unit/Quantity.hpp index 0e80c4cf..d3bed13c 100644 --- a/include/xo/unit/Quantity.hpp +++ b/include/xo/unit/Quantity.hpp @@ -62,7 +62,7 @@ namespace xo { */ repr_type r_scale = (::sqrt(rr.outer_scale_sq_) - * rr.outer_scale_factor_.template to() + * rr.outer_scale_factor_.template convert_to() * this->scale_); return Quantity(r_scale, unit2); } else { @@ -101,7 +101,7 @@ namespace xo { auto rr = detail::su_product(x.unit(), y.unit()); r_repr_type r_scale = (::sqrt(rr.outer_scale_sq_) - * rr.outer_scale_factor_.template to() + * rr.outer_scale_factor_.template convert_to() * static_cast(x.scale()) * static_cast(y.scale())); @@ -125,7 +125,7 @@ namespace xo { * so multiply is correct here */ r_repr_type r_scale = (::sqrt(rr.outer_scale_sq_) - * rr.outer_scale_factor_.template to() + * rr.outer_scale_factor_.template convert_to() * static_cast(x.scale()) / static_cast(y.scale())); diff --git a/include/xo/unit/natural_unit.hpp b/include/xo/unit/natural_unit.hpp index e51050f5..1eba006e 100644 --- a/include/xo/unit/natural_unit.hpp +++ b/include/xo/unit/natural_unit.hpp @@ -155,40 +155,43 @@ namespace xo { * - (b/b')^q inexactly (as a double) **/ - template + template > struct outer_scalefactor_result { - constexpr outer_scalefactor_result(const ratio::ratio & outer_scale_exact, + constexpr outer_scalefactor_result(const OuterScale & outer_scale_factor, double outer_scale_sq) - : outer_scale_exact_{outer_scale_exact}, + : outer_scale_factor_{outer_scale_factor}, outer_scale_sq_{outer_scale_sq} {} /* (b/b')^p0 */ - ratio::ratio outer_scale_exact_; + OuterScale outer_scale_factor_; /* (b/b')^q -- until c++26 only allow q=0 or q=1/2 */ double outer_scale_sq_; }; - template + template > struct bpu2_rescale_result { constexpr bpu2_rescale_result(const bpu & bpu_rescaled, - const ratio::ratio & outer_scale_exact, + const OuterScale & outer_scale_factor, double outer_scale_sq) : bpu_rescaled_{bpu_rescaled}, - outer_scale_exact_{outer_scale_exact}, + outer_scale_factor_{outer_scale_factor}, outer_scale_sq_{outer_scale_sq} {} /* (b'.u)^p */ bpu bpu_rescaled_; /* (b/b')^p0 */ - ratio::ratio outer_scale_exact_; + OuterScale outer_scale_factor_; /* [(b/b')^q]^2 -- until c++26 only allow q=0 or q=1/2 */ double outer_scale_sq_; }; - template + template < typename Int, + typename OuterScale = ratio::ratio > constexpr - bpu2_rescale_result + bpu2_rescale_result bpu2_rescale(const bpu & orig, const scalefactor_ratio_type & new_scalefactor) { @@ -203,59 +206,67 @@ namespace xo { if (p_frac.den() == 1) { mult_sq = 1.0; } else if(p_frac.den() == 2) { - mult_sq = mult.template to(); + mult_sq = mult.template convert_to(); } else { // remaining possibilities not supported until c++26 } - return bpu2_rescale_result(bpu(orig.native_dim(), - new_scalefactor, - orig.power()), - mult.power(orig.power().floor()), - mult_sq); + ratio::ratio mult_p = mult.power(orig.power().floor()); + + return bpu2_rescale_result(bpu(orig.native_dim(), + new_scalefactor, + orig.power()), + mult_p.template convert_to(), + mult_sq); } - template + template < typename Int, + typename OuterScale > constexpr - outer_scalefactor_result + outer_scalefactor_result bpu_product_inplace(bpu * p_target_bpu, const bpu & rhs_bpu_orig) { assert(rhs_bpu_orig.native_dim() == p_target_bpu->native_dim()); - bpu2_rescale_result rhs_bpu_rr = bpu2_rescale(rhs_bpu_orig, - p_target_bpu->scalefactor()); + bpu2_rescale_result rhs_bpu_rr + = bpu2_rescale(rhs_bpu_orig, + p_target_bpu->scalefactor().template convert_to()); *p_target_bpu = bpu(p_target_bpu->native_dim(), p_target_bpu->scalefactor(), p_target_bpu->power() + rhs_bpu_orig.power()); - return outer_scalefactor_result(rhs_bpu_rr.outer_scale_exact_, + return outer_scalefactor_result(rhs_bpu_rr.outer_scale_factor_, rhs_bpu_rr.outer_scale_sq_); } - template + template < typename Int, + typename OuterScale > constexpr - outer_scalefactor_result + outer_scalefactor_result bpu_ratio_inplace(bpu * p_target_bpu, const bpu & rhs_bpu_orig) { assert(rhs_bpu_orig.native_dim() == p_target_bpu->native_dim()); - bpu2_rescale_result rhs_bpu_rr = bpu2_rescale(rhs_bpu_orig, - p_target_bpu->scalefactor()); + bpu2_rescale_result rhs_bpu_rr + = bpu2_rescale(rhs_bpu_orig, + p_target_bpu->scalefactor()); *p_target_bpu = bpu(p_target_bpu->native_dim(), p_target_bpu->scalefactor(), p_target_bpu->power() - rhs_bpu_orig.power()); - return outer_scalefactor_result(power_ratio_type(1,1) / rhs_bpu_rr.outer_scale_exact_, - 1.0 / rhs_bpu_rr.outer_scale_sq_); + return outer_scalefactor_result + (OuterScale(1) / rhs_bpu_rr.outer_scale_factor_, + 1.0 / rhs_bpu_rr.outer_scale_sq_); } - template + template < typename Int, + typename OuterScale > constexpr - outer_scalefactor_result + outer_scalefactor_result nu_product_inplace(natural_unit * p_target, const bpu & bpu) { @@ -264,7 +275,8 @@ namespace xo { auto * p_target_bpu = &((*p_target)[i]); if (p_target_bpu->native_dim() == bpu.native_dim()) { - outer_scalefactor_result retval = bpu_product_inplace(p_target_bpu, bpu); + outer_scalefactor_result retval + = bpu_product_inplace(p_target_bpu, bpu); if (p_target_bpu->power().is_zero()) { /* dimension assoc'd with *p_target_bpu has been cancelled */ @@ -282,14 +294,15 @@ namespace xo { p_target->push_back(bpu); - return outer_scalefactor_result - (ratio::ratio(1, 1) /*outer_scale_exact*/, + return outer_scalefactor_result + (OuterScale(1) /*outer_scale_factor*/, 1.0 /*outer_scale_sq*/); } - template + template < typename Int, + typename OuterScale = ratio::ratio > constexpr - outer_scalefactor_result + outer_scalefactor_result nu_ratio_inplace(natural_unit * p_target, const bpu & bpu) { @@ -298,7 +311,8 @@ namespace xo { auto * p_target_bpu = &((*p_target)[i]); if (p_target_bpu->native_dim() == bpu.native_dim()) { - outer_scalefactor_result retval = bpu_ratio_inplace(p_target_bpu, bpu); + outer_scalefactor_result retval + = bpu_ratio_inplace(p_target_bpu, bpu); if (p_target_bpu->power().is_zero()) { /* dimension assoc'd with *p_target_bpu has been cancelled */ @@ -316,8 +330,8 @@ namespace xo { p_target->push_back(bpu.reciprocal()); - return outer_scalefactor_result - (ratio::ratio(1, 1) /*outer_scale_exact*/, + return outer_scalefactor_result + (OuterScale(1) /*outer_scale_factor*/, 1.0 /*outer_scale_sq*/); } diff --git a/include/xo/unit/scaled_unit.hpp b/include/xo/unit/scaled_unit.hpp index f96ffef8..b13adeac 100644 --- a/include/xo/unit/scaled_unit.hpp +++ b/include/xo/unit/scaled_unit.hpp @@ -100,9 +100,11 @@ namespace xo { rr.outer_scale_sq_); }; - template > + template , + typename OuterScale = ratio::ratio> constexpr - scaled_unit + scaled_unit su_product(const natural_unit & lhs_bpu_array, const natural_unit & rhs_bpu_array) { @@ -113,26 +115,29 @@ namespace xo { * in lh_bpu_array */ auto sfr = (detail::outer_scalefactor_result - (ratio::ratio(1, 1) /*outer_scale_exact*/, + (OuterScale(1) /*outer_scale_factor*/, 1.0 /*outer_scale_sq*/)); for (std::size_t i = 0; i < rhs_bpu_array.n_bpu(); ++i) { - auto sfr2 = nu_product_inplace(&prod, rhs_bpu_array[i].template to_repr()); + auto sfr2 = nu_product_inplace(&prod, + rhs_bpu_array[i].template to_repr()); - sfr.outer_scale_exact_ = sfr.outer_scale_exact_ * sfr2.outer_scale_exact_; + sfr.outer_scale_factor_ = sfr.outer_scale_factor_ * sfr2.outer_scale_factor_; sfr.outer_scale_sq_ *= sfr2.outer_scale_sq_; } - return scaled_unit(prod.template to_repr(), - sfr.outer_scale_exact_, - sfr.outer_scale_sq_); + return scaled_unit(prod.template to_repr(), + sfr.outer_scale_factor_, + sfr.outer_scale_sq_); } /* use Int2x to accumulate scalefactor */ - template > + template < typename Int, + typename Int2x = width2x, + typename OuterScale = ratio::ratio > constexpr - scaled_unit + scaled_unit su_ratio(const natural_unit & nu_lhs, const natural_unit & nu_rhs) { @@ -142,23 +147,25 @@ namespace xo { * any basis-units in rhs_bpu_array that conflict with the same dimension * in lh_bpu_array */ - auto sfr = (detail::outer_scalefactor_result - (ratio::ratio(1, 1) /*outer_scale_exact*/, + auto sfr = (detail::outer_scalefactor_result + (OuterScale(1) /*outer_scale_factor*/, 1.0 /*outer_scale_sq*/)); for (std::size_t i = 0; i < nu_rhs.n_bpu(); ++i) { - auto sfr2 = nu_ratio_inplace(&ratio, nu_rhs[i].template to_repr()); + auto sfr2 = nu_ratio_inplace(&ratio, + nu_rhs[i].template to_repr()); /* note: nu_ratio_inplace() reports multiplicative outer scaling factors, * so multiply is correct here */ - sfr.outer_scale_exact_ = sfr.outer_scale_exact_ * sfr2.outer_scale_exact_; + sfr.outer_scale_factor_ = (sfr.outer_scale_factor_ + * sfr2.outer_scale_factor_); sfr.outer_scale_sq_ *= sfr2.outer_scale_sq_; } - return scaled_unit(ratio.template to_repr(), - sfr.outer_scale_exact_, - sfr.outer_scale_sq_); + return scaled_unit(ratio.template to_repr(), + sfr.outer_scale_factor_, + sfr.outer_scale_sq_); } } diff --git a/utest/Quantity.test.cpp b/utest/Quantity.test.cpp index 5d6c78d6..511cdbe0 100644 --- a/utest/Quantity.test.cpp +++ b/utest/Quantity.test.cpp @@ -237,8 +237,8 @@ namespace xo { REQUIRE(nu2_j.n_bpu() == 1); - double rx = (nu1_j[0].scalefactor().template to() - / nu2_j[0].scalefactor().template to()); + double rx = (nu1_j[0].scalefactor().template convert_to() + / nu2_j[0].scalefactor().template convert_to()); if ((rx > max_magdiff_per_bu) || (rx < 1.0/max_magdiff_per_bu)) { log && log(xtag("nu_z", p_nu_v->size()), xtag("nu2_j_ix", nu2_j_ix)); @@ -295,7 +295,9 @@ namespace xo { * */ auto su = xo::qty::detail::su_ratio(q1.unit(), q2.unit()); + decltype(q1)::ratio_int2x_type, + double> + (q1.unit(), q2.unit()); INFO(xtag("su.natural_unit", su.natural_unit_)); INFO(xtag("su.outer_scale_exact", su.outer_scale_factor_));