xo-unit: use Quantity.scale repr for outer scale to avoid overflow
This commit is contained in:
parent
b52ab13495
commit
06ed8fec8b
4 changed files with 83 additions and 60 deletions
|
|
@ -62,7 +62,7 @@ namespace xo {
|
|||
*/
|
||||
|
||||
repr_type r_scale = (::sqrt(rr.outer_scale_sq_)
|
||||
* rr.outer_scale_factor_.template to<repr_type>()
|
||||
* rr.outer_scale_factor_.template convert_to<repr_type>()
|
||||
* this->scale_);
|
||||
return Quantity(r_scale, unit2);
|
||||
} else {
|
||||
|
|
@ -101,7 +101,7 @@ namespace xo {
|
|||
auto rr = detail::su_product<r_int_type, r_int2x_type>(x.unit(), y.unit());
|
||||
|
||||
r_repr_type r_scale = (::sqrt(rr.outer_scale_sq_)
|
||||
* rr.outer_scale_factor_.template to<r_repr_type>()
|
||||
* rr.outer_scale_factor_.template convert_to<r_repr_type>()
|
||||
* static_cast<r_repr_type>(x.scale())
|
||||
* static_cast<r_repr_type>(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<r_repr_type>()
|
||||
* rr.outer_scale_factor_.template convert_to<r_repr_type>()
|
||||
* static_cast<r_repr_type>(x.scale())
|
||||
/ static_cast<r_repr_type>(y.scale()));
|
||||
|
||||
|
|
|
|||
|
|
@ -155,40 +155,43 @@ namespace xo {
|
|||
* - (b/b')^q inexactly (as a double)
|
||||
**/
|
||||
|
||||
template <typename Int>
|
||||
template <typename Int,
|
||||
typename OuterScale = ratio::ratio<Int> >
|
||||
struct outer_scalefactor_result {
|
||||
constexpr outer_scalefactor_result(const ratio::ratio<Int> & 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<Int> 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 <typename Int>
|
||||
template <typename Int,
|
||||
typename OuterScale = ratio::ratio<Int> >
|
||||
struct bpu2_rescale_result {
|
||||
constexpr bpu2_rescale_result(const bpu<Int> & bpu_rescaled,
|
||||
const ratio::ratio<Int> & 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<Int> bpu_rescaled_;
|
||||
/* (b/b')^p0 */
|
||||
ratio::ratio<Int> 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 <typename Int>
|
||||
template < typename Int,
|
||||
typename OuterScale = ratio::ratio<Int> >
|
||||
constexpr
|
||||
bpu2_rescale_result<Int>
|
||||
bpu2_rescale_result<Int, OuterScale>
|
||||
bpu2_rescale(const bpu<Int> & 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<double>();
|
||||
mult_sq = mult.template convert_to<double>();
|
||||
} else {
|
||||
// remaining possibilities not supported until c++26
|
||||
}
|
||||
|
||||
return bpu2_rescale_result<Int>(bpu<Int>(orig.native_dim(),
|
||||
new_scalefactor,
|
||||
orig.power()),
|
||||
mult.power(orig.power().floor()),
|
||||
mult_sq);
|
||||
ratio::ratio<Int> mult_p = mult.power(orig.power().floor());
|
||||
|
||||
return bpu2_rescale_result<Int, OuterScale>(bpu<Int>(orig.native_dim(),
|
||||
new_scalefactor,
|
||||
orig.power()),
|
||||
mult_p.template convert_to<OuterScale>(),
|
||||
mult_sq);
|
||||
}
|
||||
|
||||
template <typename Int>
|
||||
template < typename Int,
|
||||
typename OuterScale >
|
||||
constexpr
|
||||
outer_scalefactor_result<Int>
|
||||
outer_scalefactor_result<Int, OuterScale>
|
||||
bpu_product_inplace(bpu<Int> * p_target_bpu,
|
||||
const bpu<Int> & rhs_bpu_orig)
|
||||
{
|
||||
assert(rhs_bpu_orig.native_dim() == p_target_bpu->native_dim());
|
||||
|
||||
bpu2_rescale_result<Int> rhs_bpu_rr = bpu2_rescale(rhs_bpu_orig,
|
||||
p_target_bpu->scalefactor());
|
||||
bpu2_rescale_result<Int, OuterScale> rhs_bpu_rr
|
||||
= bpu2_rescale<Int, OuterScale>(rhs_bpu_orig,
|
||||
p_target_bpu->scalefactor().template convert_to<OuterScale>());
|
||||
|
||||
*p_target_bpu = bpu<Int>(p_target_bpu->native_dim(),
|
||||
p_target_bpu->scalefactor(),
|
||||
p_target_bpu->power() + rhs_bpu_orig.power());
|
||||
|
||||
return outer_scalefactor_result<Int>(rhs_bpu_rr.outer_scale_exact_,
|
||||
return outer_scalefactor_result<Int>(rhs_bpu_rr.outer_scale_factor_,
|
||||
rhs_bpu_rr.outer_scale_sq_);
|
||||
}
|
||||
|
||||
template <typename Int>
|
||||
template < typename Int,
|
||||
typename OuterScale >
|
||||
constexpr
|
||||
outer_scalefactor_result<Int>
|
||||
outer_scalefactor_result<Int, OuterScale>
|
||||
bpu_ratio_inplace(bpu<Int> * p_target_bpu,
|
||||
const bpu<Int> & rhs_bpu_orig)
|
||||
{
|
||||
assert(rhs_bpu_orig.native_dim() == p_target_bpu->native_dim());
|
||||
|
||||
bpu2_rescale_result<Int> rhs_bpu_rr = bpu2_rescale(rhs_bpu_orig,
|
||||
p_target_bpu->scalefactor());
|
||||
bpu2_rescale_result<Int, OuterScale> rhs_bpu_rr
|
||||
= bpu2_rescale<Int, OuterScale>(rhs_bpu_orig,
|
||||
p_target_bpu->scalefactor());
|
||||
|
||||
*p_target_bpu = bpu<Int>(p_target_bpu->native_dim(),
|
||||
p_target_bpu->scalefactor(),
|
||||
p_target_bpu->power() - rhs_bpu_orig.power());
|
||||
|
||||
return outer_scalefactor_result<Int>(power_ratio_type(1,1) / rhs_bpu_rr.outer_scale_exact_,
|
||||
1.0 / rhs_bpu_rr.outer_scale_sq_);
|
||||
return outer_scalefactor_result<Int, OuterScale>
|
||||
(OuterScale(1) / rhs_bpu_rr.outer_scale_factor_,
|
||||
1.0 / rhs_bpu_rr.outer_scale_sq_);
|
||||
}
|
||||
|
||||
template <typename Int>
|
||||
template < typename Int,
|
||||
typename OuterScale >
|
||||
constexpr
|
||||
outer_scalefactor_result<Int>
|
||||
outer_scalefactor_result<Int, OuterScale>
|
||||
nu_product_inplace(natural_unit<Int> * p_target,
|
||||
const bpu<Int> & 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<Int> retval = bpu_product_inplace(p_target_bpu, bpu);
|
||||
outer_scalefactor_result<Int, OuterScale> retval
|
||||
= bpu_product_inplace<Int, OuterScale>(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<Int>
|
||||
(ratio::ratio<Int>(1, 1) /*outer_scale_exact*/,
|
||||
return outer_scalefactor_result<Int, OuterScale>
|
||||
(OuterScale(1) /*outer_scale_factor*/,
|
||||
1.0 /*outer_scale_sq*/);
|
||||
}
|
||||
|
||||
template <typename Int>
|
||||
template < typename Int,
|
||||
typename OuterScale = ratio::ratio<Int> >
|
||||
constexpr
|
||||
outer_scalefactor_result<Int>
|
||||
outer_scalefactor_result<Int, OuterScale>
|
||||
nu_ratio_inplace(natural_unit<Int> * p_target,
|
||||
const bpu<Int> & 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<Int> retval = bpu_ratio_inplace(p_target_bpu, bpu);
|
||||
outer_scalefactor_result<Int, OuterScale> retval
|
||||
= bpu_ratio_inplace<Int, OuterScale>(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<Int>
|
||||
(ratio::ratio<Int>(1, 1) /*outer_scale_exact*/,
|
||||
return outer_scalefactor_result<Int, OuterScale>
|
||||
(OuterScale(1) /*outer_scale_factor*/,
|
||||
1.0 /*outer_scale_sq*/);
|
||||
}
|
||||
|
||||
|
|
|
|||
|
|
@ -100,9 +100,11 @@ namespace xo {
|
|||
rr.outer_scale_sq_);
|
||||
};
|
||||
|
||||
template <typename Int, typename Int2x = width2x<Int>>
|
||||
template <typename Int,
|
||||
typename Int2x = width2x<Int>,
|
||||
typename OuterScale = ratio::ratio<Int2x>>
|
||||
constexpr
|
||||
scaled_unit<Int>
|
||||
scaled_unit<Int, OuterScale>
|
||||
su_product(const natural_unit<Int> & lhs_bpu_array,
|
||||
const natural_unit<Int> & rhs_bpu_array)
|
||||
{
|
||||
|
|
@ -113,26 +115,29 @@ namespace xo {
|
|||
* in lh_bpu_array
|
||||
*/
|
||||
auto sfr = (detail::outer_scalefactor_result<Int2x>
|
||||
(ratio::ratio<Int2x>(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<Int2x>());
|
||||
auto sfr2 = nu_product_inplace<Int2x, OuterScale>(&prod,
|
||||
rhs_bpu_array[i].template to_repr<Int2x>());
|
||||
|
||||
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<Int>(prod.template to_repr<Int>(),
|
||||
sfr.outer_scale_exact_,
|
||||
sfr.outer_scale_sq_);
|
||||
return scaled_unit<Int, OuterScale>(prod.template to_repr<Int>(),
|
||||
sfr.outer_scale_factor_,
|
||||
sfr.outer_scale_sq_);
|
||||
}
|
||||
|
||||
/* use Int2x to accumulate scalefactor
|
||||
*/
|
||||
template <typename Int, typename Int2x = width2x<Int>>
|
||||
template < typename Int,
|
||||
typename Int2x = width2x<Int>,
|
||||
typename OuterScale = ratio::ratio<Int2x> >
|
||||
constexpr
|
||||
scaled_unit<Int>
|
||||
scaled_unit<Int, OuterScale>
|
||||
su_ratio(const natural_unit<Int> & nu_lhs,
|
||||
const natural_unit<Int> & 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<Int2x>
|
||||
(ratio::ratio<Int2x>(1, 1) /*outer_scale_exact*/,
|
||||
auto sfr = (detail::outer_scalefactor_result<Int2x, OuterScale>
|
||||
(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<Int2x>());
|
||||
auto sfr2 = nu_ratio_inplace<Int2x, OuterScale>(&ratio,
|
||||
nu_rhs[i].template to_repr<Int2x>());
|
||||
|
||||
/* 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<Int>(ratio.template to_repr<Int>(),
|
||||
sfr.outer_scale_exact_,
|
||||
sfr.outer_scale_sq_);
|
||||
return scaled_unit<Int, OuterScale>(ratio.template to_repr<Int>(),
|
||||
sfr.outer_scale_factor_,
|
||||
sfr.outer_scale_sq_);
|
||||
}
|
||||
}
|
||||
|
||||
|
|
|
|||
|
|
@ -237,8 +237,8 @@ namespace xo {
|
|||
|
||||
REQUIRE(nu2_j.n_bpu() == 1);
|
||||
|
||||
double rx = (nu1_j[0].scalefactor().template to<double>()
|
||||
/ nu2_j[0].scalefactor().template to<double>());
|
||||
double rx = (nu1_j[0].scalefactor().template convert_to<double>()
|
||||
/ nu2_j[0].scalefactor().template convert_to<double>());
|
||||
|
||||
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<decltype(q1)::ratio_int_type,
|
||||
decltype(q1)::ratio_int2x_type>(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_));
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue