xo-unit: + Quantity division

This commit is contained in:
Roland Conybeare 2024-04-23 14:46:40 -04:00
commit 8548f26143
4 changed files with 353 additions and 25 deletions

View file

@ -1,4 +1,4 @@
/** @file Quantity2.hpp
/** @file Quantity.hpp
*
* Author: Roland Conybeare
**/
@ -34,7 +34,7 @@ namespace xo {
public:
constexpr Quantity(Repr scale,
const natural_unit<Int> & unit)
const natural_unit<Int> & unit)
: scale_{scale}, unit_{unit} {}
constexpr const repr_type & scale() const { return scale_; }
@ -63,6 +63,29 @@ namespace xo {
rr.natural_unit_);
}
template <typename Quantity2>
static constexpr
auto divide(const Quantity & x, const Quantity2 & y) {
using r_repr_type = std::common_type_t<typename Quantity::repr_type,
typename Quantity2::repr_type>;
using r_int_type = std::common_type_t<typename Quantity::ratio_int_type,
typename Quantity2::ratio_int_type>;
auto rr = detail::nu_ratio(x.unit(), y.unit());
/* note: nu_ratio() reports multiplicative outer scaling factors,
* so multiply is correct here
*/
r_repr_type r_scale = (::sqrt(rr.outer_scale_sq_)
* rr.outer_scale_exact_.template to<r_repr_type>()
* static_cast<r_repr_type>(x.scale())
/ static_cast<r_repr_type>(y.scale()));
return Quantity<r_repr_type, r_int_type>(r_scale,
rr.natural_unit_);
}
private:
/** @brief quantity represents this multiple of a unit amount **/
Repr scale_ = Repr{};
@ -89,8 +112,17 @@ namespace xo {
{
return Quantity::multiply(x, y);
}
/** note: won't have constexpr result until c++26 (when ::sqrt(), ::pow() are constexpr)
**/
template <typename Quantity, typename OtherQuantity>
constexpr auto
operator/ (const Quantity & x, const OtherQuantity & y)
{
return Quantity::divide(x, y);
}
} /*namespace qty*/
} /*namespace xo*/
/** end Quantity2.hpp **/
/** end Quantity.hpp **/

View file

@ -16,10 +16,11 @@ namespace xo {
/** @class natural_unit
* @brief an array representing the cartesian product of distinct basis-power-units
*
* 1. Each bpu in the array represents a power of a basis dimension, e.g. "meter" or "second^2".
* 2. Each bpu in an array has a different dimension id.
* 1. Quantities are represented as a multiple of a natural unit
* 2. Each bpu in the array represents a power of a basis dimension, e.g. "meter" or "second^2".
* 3. Each bpu in an array has a different dimension id.
* For example dim::time, if present, appears once.
* 3. Basis dimensions can appear in any order.
* 4. Basis dimensions can appear in any order.
* Order used for constructing abbreviations: will get @c "kg.m" or @c "m.kg"
* depending on the orderin of @c dim::distance and @c dim::mass in @c bpu_v_
**/
@ -46,6 +47,14 @@ namespace xo {
return retval;
}
/** @brief remove bpu at position @p p **/
constexpr void remove_bpu(size_t p) {
for (std::size_t i = p; i+1 < n_bpu_; ++i)
bpu_v_[i] = bpu_v_[i+1];
--n_bpu_;
}
constexpr void push_back(const bpu<Int> & bpu) {
if (n_bpu_ < n_dim)
bpu_v_[n_bpu_++] = bpu;
@ -187,8 +196,8 @@ namespace xo {
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());
p_target_bpu->scalefactor(),
p_target_bpu->power() + rhs_bpu_orig.power());
return outer_scalefactor_result<Int>(rhs_bpu_rr.outer_scale_exact_,
rhs_bpu_rr.outer_scale_sq_);
@ -197,21 +206,49 @@ namespace xo {
template <typename Int>
constexpr
outer_scalefactor_result<Int>
bpu_array_product_inplace(natural_unit<Int> * p_target,
const bpu<Int> & bpu)
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());
*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_);
}
template <typename Int>
constexpr
outer_scalefactor_result<Int>
nu_product_inplace(natural_unit<Int> * p_target,
const bpu<Int> & bpu)
{
std::size_t i = 0;
for (; i < p_target->n_bpu(); ++i) {
if ((*p_target)[i].native_dim() == bpu.native_dim()) {
outer_scalefactor_result<Int> retval = bpu_product_inplace(&((*p_target)[i]), bpu);
auto * p_target_bpu = &((*p_target)[i]);
/* TODO: strip 0 power */
if (p_target_bpu->native_dim() == bpu.native_dim()) {
outer_scalefactor_result<Int> 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 */
p_target->remove_bpu(i);
}
return retval;
}
}
/* control here: i=p_target->n_bpu() */
/* control here: i=p_target->n_bpu()
* Dimension represented by bpu does not already appear in *p_target.
* Adopt bpu's scalefactor
*/
p_target->push_back(bpu);
return outer_scalefactor_result<Int>
@ -219,6 +256,40 @@ namespace xo {
1.0 /*outer_scale_sq*/);
}
template <typename Int>
constexpr
outer_scalefactor_result<Int>
nu_ratio_inplace(natural_unit<Int> * p_target,
const bpu<Int> & bpu)
{
std::size_t i = 0;
for (; i < p_target->n_bpu(); ++i) {
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);
if (p_target_bpu->power().is_zero()) {
/* dimension assoc'd with *p_target_bpu has been cancelled */
p_target->remove_bpu(i);
}
return retval;
}
}
/* here: i=p_target->n_bpu()
* Dimension represented by bpu does not already appear in *p_target.
* Adopt bpu's scalefactor
*/
p_target->push_back(bpu.reciprocal());
return outer_scalefactor_result<Int>
(ratio::ratio<Int>(1, 1) /*outer_scale_exact*/,
1.0 /*outer_scale_sq*/);
}
template <typename Int>
constexpr natural_unit<Int>
nu_reciprocal(const natural_unit<Int> & nu)

View file

@ -67,15 +67,15 @@ namespace xo {
template <typename Int>
constexpr
scaled_unit<Int>
nu_product(const natural_unit<Int> & lhs_bpu_array,
const bpu<Int> & rhs_bpu)
nu_bpu_product(const natural_unit<Int> & lhs_bpu_array,
const bpu<Int> & rhs_bpu)
{
natural_unit<Int> prod = lhs_bpu_array;
auto rr = bpu_array_product_inplace(&prod, rhs_bpu);
auto rr = nu_product_inplace(&prod, rhs_bpu);
return scaled_unit<Int>(prod,
rr.outer_scale_exact_,
rr.outer_scale_sq_);
rr.outer_scale_exact_,
rr.outer_scale_sq_);
};
template <typename Int>
@ -95,7 +95,7 @@ namespace xo {
1.0 /*outer_scale_sq*/));
for (std::size_t i = 0; i < rhs_bpu_array.n_bpu(); ++i) {
auto sfr2 = bpu_array_product_inplace(&prod, rhs_bpu_array[i]);
auto sfr2 = nu_product_inplace(&prod, rhs_bpu_array[i]);
sfr.outer_scale_exact_ = sfr.outer_scale_exact_ * sfr2.outer_scale_exact_;
sfr.outer_scale_sq_ *= sfr2.outer_scale_sq_;
@ -106,6 +106,36 @@ namespace xo {
sfr.outer_scale_sq_);
}
template <typename Int>
constexpr
scaled_unit<Int>
nu_ratio(const natural_unit<Int> & nu_lhs,
const natural_unit<Int> & nu_rhs)
{
natural_unit<Int> ratio = nu_lhs;
/* accumulate product of scalefactors spun off by rescaling
* any basis-units in rhs_bpu_array that conflict with the same dimension
* in lh_bpu_array
*/
auto sfr = (detail::outer_scalefactor_result<Int>
(scalefactor_ratio_type(1, 1) /*outer_scale_exact*/,
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]);
/* 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_sq_ *= sfr2.outer_scale_sq_;
}
return scaled_unit<Int>(ratio,
sfr.outer_scale_exact_,
sfr.outer_scale_sq_);
}
}
template <typename Int>
@ -114,7 +144,21 @@ namespace xo {
const scaled_unit<Int> & y_unit)
{
auto rr = detail::nu_product(x_unit.natural_unit_,
y_unit.natural_unit_);
y_unit.natural_unit_);
return (scaled_unit<Int>
(rr.natural_unit_,
rr.outer_scale_exact_ * x_unit.outer_scale_exact_ * y_unit.outer_scale_exact_,
rr.outer_scale_sq_ * x_unit.outer_scale_sq_ * y_unit.outer_scale_sq_));
}
template <typename Int>
inline constexpr scaled_unit<Int>
operator/ (const scaled_unit<Int> & x_unit,
const scaled_unit<Int> & y_unit)
{
auto rr = detail::nu_ratio(x_unit.natural_unit_,
y_unit.natural_unit_);
return (scaled_unit<Int>
(rr.natural_unit_,

View file

@ -46,6 +46,9 @@ namespace xo {
using xo::qty::natural_unit;
using xo::qty::bpu_array_maker;
using xo::qty::detail::nu_product;
using xo::qty::detail::nu_bpu_product;
using xo::qty::detail::nu_ratio_inplace;
using xo::qty::detail::nu_ratio;
using xo::qty::unit_qty;
using xo::print::ccs;
@ -595,7 +598,7 @@ namespace xo {
static_assert(bpu.power() == power_ratio_type(-1, 2));
constexpr auto prod_rr = nu_product(v, bpu);
constexpr auto prod_rr = nu_bpu_product(v, bpu);
log && log(xtag("prod_rr.bpu_array", prod_rr.natural_unit_));
log && log(xtag("prod_rr.outer_scale_exact", prod_rr.outer_scale_exact_));
@ -642,7 +645,7 @@ namespace xo {
static_assert(bpu.power() == power_ratio_type(-1, 2));
constexpr auto prod_rr = nu_product(v, bpu);
constexpr auto prod_rr = nu_bpu_product(v, bpu);
log && log(xtag("prod_rr.bpu_array", prod_rr.natural_unit_));
log && log(xtag("prod_rr.outer_scale_exact", prod_rr.outer_scale_exact_));
@ -686,7 +689,7 @@ namespace xo {
static_assert(bpu.power() == power_ratio_type(-1, 1));
constexpr auto prod_rr = nu_product(v, bpu);
constexpr auto prod_rr = nu_bpu_product(v, bpu);
log && log(xtag("prod_rr.bpu_array", prod_rr.natural_unit_));
log && log(xtag("prod_rr.outer_scale_exact", prod_rr.outer_scale_exact_));
@ -832,6 +835,104 @@ namespace xo {
}
} /*TEST_CASE(natural_unit2)*/
TEST_CASE("natural_unit3", "[natural_unit]") {
constexpr bool c_debug_flag = true;
// can get bits from /dev/random by uncommenting the 2nd line below
//uint64_t seed = xxx;
//rng::Seed<xoshio256ss> seed;
//auto rng = xo::rng::xoshiro256ss(seed);
scope log(XO_DEBUG2(c_debug_flag, "TEST_CASE.natural_unit3"));
//log && log("(A)", xtag("foo", foo));
{
constexpr natural_unit<int64_t> v
= (bpu_array_maker<int64_t>::make_bpu_array
(bpu<int64_t>(dim::mass, scalefactor_ratio_type(1000, 1), power_ratio_type(1, 1)),
bpu<int64_t>(dim::distance, scalefactor_ratio_type(1, 1), power_ratio_type(1, 1))));
static_assert(v.n_bpu() == 2);
log && log(xtag("v.abbrev", v.abbrev()));
static_assert(v.abbrev().size() > 0);
static_assert(v.abbrev() == flatstring("kg.m"));
{
natural_unit<int64_t> w = v;
nu_ratio_inplace(&w,
bpu<int64_t>(dim::mass, scalefactor_ratio_type(1000, 1), power_ratio_type(1, 1)));
log && log(xtag("w.abbrev", w.abbrev()));
REQUIRE(w.n_bpu() == 1);
REQUIRE(w[0].native_dim() == dim::distance);
REQUIRE(w.abbrev() == flatstring("m"));
}
{
constexpr natural_unit<int64_t> w
= (bpu_array_maker<int64_t>::make_bpu_array
(bpu<int64_t>(dim::mass, scalefactor_ratio_type(1000, 1), power_ratio_type(1, 1))));
static_assert(w.n_bpu() == 1);
log && log(xtag("w.abbrev", w.abbrev()));
constexpr auto rr = nu_ratio(v, w);
log && log(xtag("rr", rr));
REQUIRE(rr.natural_unit_.n_bpu() == 1);
REQUIRE(rr.natural_unit_[0].native_dim() == dim::distance);
REQUIRE(rr.natural_unit_.abbrev() == flatstring("m"));
}
{
constexpr natural_unit<int64_t> w
= (bpu_array_maker<int64_t>::make_bpu_array
(bpu<int64_t>(dim::time, scalefactor_ratio_type(1, 1), power_ratio_type(1, 1))));
static_assert(w.n_bpu() == 1);
log && log(xtag("w.abbrev", w.abbrev()));
constexpr auto rr = nu_ratio(v, w);
log && log(xtag("rr", rr));
REQUIRE(rr.natural_unit_.n_bpu() == 3);
REQUIRE(rr.natural_unit_[0].native_dim() == dim::mass);
REQUIRE(rr.natural_unit_[1].native_dim() == dim::distance);
REQUIRE(rr.natural_unit_[2].native_dim() == dim::time);
REQUIRE(rr.natural_unit_.abbrev() == flatstring("kg.m.s^-1"));
}
{
natural_unit<int64_t> w = v;
REQUIRE(w.n_bpu() == 2);
REQUIRE(w[0].native_dim() == dim::mass);
nu_ratio_inplace(&w,
bpu<int64_t>(dim::time, scalefactor_ratio_type(1, 1), power_ratio_type(2, 1)));
REQUIRE(w.n_bpu() == 3);
REQUIRE(w[0].native_dim() == dim::mass);
REQUIRE(w[1].native_dim() == dim::distance);
REQUIRE(w[2].native_dim() == dim::time);
log && log(xtag("w.abbrev", w.abbrev()));
REQUIRE(w.n_bpu() == 3);
REQUIRE(w.abbrev() == flatstring("kg.m.s^-2"));
}
}
} /*TEST_CASE(natural_unit3)*/
TEST_CASE("scaled_unit0", "[scaled_unit0]") {
constexpr bool c_debug_flag = true;
@ -855,6 +956,29 @@ namespace xo {
static_assert(ng2.natural_unit_.n_bpu() == 1);
} /*TEST_CASE(scaled_unit0)*/
TEST_CASE("scaled_unit1", "[scaled_unit1]") {
constexpr bool c_debug_flag = true;
// can get bits from /dev/random by uncommenting the 2nd line below
//uint64_t seed = xxx;
//rng::Seed<xoshio256ss> seed;
//auto rng = xo::rng::xoshiro256ss(seed);
scope log(XO_DEBUG2(c_debug_flag, "TEST_CASE.scaled_unit1"));
//log && log("(A)", xtag("foo", foo));
constexpr auto ng = su2::nanogram;
constexpr auto ng2 = ng / ng;
log && log(xtag("ng", ng));
log && log(xtag("ng/ng", ng2));
//log && log(xtag("ng/ng",
static_assert(ng.natural_unit_.n_bpu() == 1);
static_assert(ng2.natural_unit_.n_bpu() == 0);
} /*TEST_CASE(scaled_unit1)*/
TEST_CASE("Quantity", "[Quantity]") {
constexpr bool c_debug_flag = true;
@ -908,6 +1032,27 @@ namespace xo {
scope log(XO_DEBUG2(c_debug_flag, "TEST_CASE.Quantity3"));
//log && log("(A)", xtag("foo", foo));
/* not constexpr until c++26 */
Quantity ng = unit_qty(su2::nanogram);
auto ng0 = ng / ng;
log && log(xtag("ng/ng", ng0));
REQUIRE(ng0.scale() == 1);
} /*TEST_CASE(Quantity3)*/
TEST_CASE("Quantity4", "[Quantity]") {
constexpr bool c_debug_flag = true;
// can get bits from /dev/random by uncommenting the 2nd line below
//uint64_t seed = xxx;
//rng::Seed<xoshio256ss> seed;
//auto rng = xo::rng::xoshiro256ss(seed);
scope log(XO_DEBUG2(c_debug_flag, "TEST_CASE.Quantity4"));
//log && log("(A)", xtag("foo", foo));
/* not constexpr until c++26 */
Quantity ng = unit_qty(su2::nanogram);
Quantity ug = unit_qty(su2::microgram);
@ -937,7 +1082,43 @@ namespace xo {
}
//REQUIRE(ng2.scale() == 1);
} /*TEST_CASE(Quantity3)*/
} /*TEST_CASE(Quantity4)*/
TEST_CASE("Quantity5", "[Quantity]") {
constexpr bool c_debug_flag = true;
// can get bits from /dev/random by uncommenting the 2nd line below
//uint64_t seed = xxx;
//rng::Seed<xoshio256ss> seed;
//auto rng = xo::rng::xoshiro256ss(seed);
scope log(XO_DEBUG2(c_debug_flag, "TEST_CASE.Quantity5"));
//log && log("(A)", xtag("foo", foo));
/* not constexpr until c++26 */
Quantity ng = unit_qty(su2::nanogram);
Quantity ug = unit_qty(su2::microgram);
{
auto ratio1 = ng / ug;
log && log(xtag("ng/ug", ratio1));
/* units will be nanograms, since that's on lhs */
REQUIRE(ratio1.unit().n_bpu() == 0);
REQUIRE(ratio1.scale() == 0.001);
}
{
auto ratio2 = ug / ng;
log && log(xtag("ug/ng", ratio2));
REQUIRE(ratio2.unit().n_bpu() == 0);
REQUIRE(ratio2.scale() == 1000.0);
}
//REQUIRE(ng2.scale() == 1);
} /*TEST_CASE(Quantity5)*/
} /*namespace ut*/
} /*namespace xo*/