Line data Source code
1 : /** @file natural_unit.hpp
2 : *
3 : * Author: Roland Conybeare
4 : **/
5 :
6 : #pragma once
7 :
8 : #include "bpu.hpp"
9 : #include <cmath>
10 : #include <cassert>
11 :
12 : namespace xo {
13 : namespace qty {
14 : using nu_abbrev_type = flatstring<32>;
15 :
16 : template <typename Int>
17 : class natural_unit;
18 :
19 : /** @class natural_unit
20 : * @brief an array representing the cartesian product of distinct basis-power-units
21 : *
22 : * 1. Quantities are represented as a multiple of a natural unit
23 : * 2. Each bpu in the array represents a power of a basis dimension, e.g. "meter" or "second^2".
24 : * 3. Each bpu in an array has a different dimension id.
25 : * For example dim::time, if present, appears once.
26 : * 4. Basis dimensions can appear in any order.
27 : * Order used for constructing abbreviations: will get @c "kg.m" or @c "m.kg"
28 : * depending on the orderin of @c dim::distance and @c dim::mass in @c bpu_v_
29 : **/
30 : template <typename Int>
31 : class natural_unit {
32 : public:
33 : using ratio_int_type = Int;
34 :
35 : public:
36 : constexpr natural_unit() : n_bpu_{0} {}
37 :
38 102 : constexpr std::size_t n_bpu() const { return n_bpu_; }
39 16 : constexpr bool is_dimensionless() const { return n_bpu_ == 0; }
40 :
41 : constexpr bpu<Int> * bpu_v() const { return bpu_v_; }
42 :
43 : constexpr natural_unit reciprocal() const {
44 : natural_unit retval;
45 :
46 : for (std::size_t i = 0; i < this->n_bpu(); ++i)
47 : retval.push_back((*this)[i].reciprocal());
48 :
49 : return retval;
50 : }
51 :
52 12 : constexpr nu_abbrev_type abbrev() const {
53 12 : nu_abbrev_type retval;
54 :
55 34 : for (std::size_t i = 0; i < n_bpu_; ++i) {
56 22 : if (i > 0)
57 10 : retval.append(".");
58 22 : retval.append(bpu_v_[i].abbrev(), 0, -1);
59 : }
60 :
61 12 : return retval;
62 : }
63 :
64 : /** @brief remove bpu at position @p p **/
65 20 : constexpr void remove_bpu(size_t p) {
66 21 : for (std::size_t i = p; i+1 < n_bpu_; ++i)
67 1 : bpu_v_[i] = bpu_v_[i+1];
68 :
69 20 : --n_bpu_;
70 20 : }
71 :
72 5 : constexpr void push_back(const bpu<Int> & bpu) {
73 4 : if (n_bpu_ < n_dim)
74 5 : bpu_v_[n_bpu_++] = bpu;
75 : }
76 :
77 25 : constexpr bpu<Int> & operator[](std::size_t i) { return bpu_v_[i]; }
78 47 : constexpr const bpu<Int> & operator[](std::size_t i) const { return bpu_v_[i]; }
79 :
80 : template <typename Int2>
81 4 : constexpr natural_unit<Int2> to_repr() const {
82 16 : natural_unit<Int2> retval;
83 :
84 : std::size_t i = 0;
85 8 : for (; i < n_bpu_; ++i)
86 8 : retval.push_back(bpu_v_[i].template to_repr<Int2>());
87 :
88 4 : return retval;
89 : }
90 :
91 : private:
92 : /** @brief the number of occupied slots in @c bpu_v_ **/
93 : std::size_t n_bpu_;
94 :
95 : /** @brief storage for basis power units **/
96 : bpu<Int> bpu_v_[n_dim];
97 : };
98 :
99 : namespace detail {
100 : template <typename Int, typename... Ts>
101 : constexpr void
102 : push_bpu_array(natural_unit<Int> * p_target, Ts... args);
103 :
104 : template <typename Int>
105 : constexpr void
106 : push_bpu_array(natural_unit<Int> * p_target) {}
107 :
108 : template <typename Int, typename T0, typename... Ts>
109 : constexpr void
110 : push_bpu_array(natural_unit<Int> * p_target, T0 && bpu0, Ts... args) {
111 : p_target->push_back(bpu0);
112 : push_bpu_array(p_target, args...);
113 : }
114 : }
115 :
116 : template <typename Int>
117 : struct bpu_array_maker {
118 : template <typename... Ts>
119 : static constexpr natural_unit<Int>
120 : make_bpu_array(Ts... args) {
121 : natural_unit<Int> bpu_array;
122 : detail::push_bpu_array(&bpu_array, args...);
123 : return bpu_array;
124 : }
125 : };
126 :
127 : namespace detail {
128 : /**
129 : * Given bpu ~ (b.u)^p:
130 : * - b = bpu.scalefactor
131 : * - u = bpu.native_dim
132 : * - p = bpu.power
133 : *
134 : * want to rewrite in the form a'.(b'.u)^p
135 : *
136 : * Can compute a' exactly iff p is integral.
137 : * In that case:
138 : * (b.u)^p = ((b/b').b'.u)^p
139 : * = (b/b')^p.(b'.u)^p
140 : * = a'.(b'.u)^p with a' = (b/b')^p
141 : *
142 : * Can write p = p0 + q, with p0 = floor(p) integral, q = frac(p) in [0,1)
143 : *
144 : * Then
145 : * (b/b')^p = (b/b')^p0 * (b/b')^q
146 : *
147 : * we'll compute:
148 : * - (b/b')^p0 exactly (as a ratio)
149 : * - (b/b')^q inexactly (as a double)
150 : **/
151 :
152 : template <typename Int>
153 : struct outer_scalefactor_result {
154 46 : constexpr outer_scalefactor_result(const ratio::ratio<Int> & outer_scale_exact,
155 : double outer_scale_sq)
156 46 : : outer_scale_exact_{outer_scale_exact},
157 22 : outer_scale_sq_{outer_scale_sq} {}
158 :
159 : /* (b/b')^p0 */
160 : ratio::ratio<Int> outer_scale_exact_;
161 : /* (b/b')^q -- until c++26 only allow q=0 or q=1/2 */
162 : double outer_scale_sq_;
163 : };
164 :
165 : template <typename Int>
166 : struct bpu2_rescale_result {
167 23 : constexpr bpu2_rescale_result(const bpu<Int> & bpu_rescaled,
168 : const ratio::ratio<Int> & outer_scale_exact,
169 : double outer_scale_sq)
170 23 : : bpu_rescaled_{bpu_rescaled},
171 23 : outer_scale_exact_{outer_scale_exact},
172 23 : outer_scale_sq_{outer_scale_sq}
173 : {}
174 :
175 : /* (b'.u)^p */
176 : bpu<Int> bpu_rescaled_;
177 : /* (b/b')^p0 */
178 : ratio::ratio<Int> outer_scale_exact_;
179 : /* [(b/b')^q]^2 -- until c++26 only allow q=0 or q=1/2 */
180 : double outer_scale_sq_;
181 : };
182 :
183 : template <typename Int>
184 : constexpr
185 : bpu2_rescale_result<Int>
186 23 : bpu2_rescale(const bpu<Int> & orig,
187 : const scalefactor_ratio_type & new_scalefactor)
188 : {
189 23 : ratio::ratio<Int> mult = (orig.scalefactor() / new_scalefactor);
190 :
191 : /* inv: p_frac in [0, 1) */
192 23 : auto p_frac = orig.power().frac();
193 :
194 : /* asof c++26: replace mult_sq with ::pow(mult, p_frac) */
195 23 : double mult_sq = std::numeric_limits<double>::quiet_NaN();
196 :
197 23 : if (p_frac.den() == 1) {
198 : mult_sq = 1.0;
199 0 : } else if(p_frac.den() == 2) {
200 0 : mult_sq = mult.template to<double>();
201 : } else {
202 : // remaining possibilities not supported until c++26
203 : }
204 :
205 23 : return bpu2_rescale_result<Int>(bpu<Int>(orig.native_dim(),
206 : new_scalefactor,
207 : orig.power()),
208 23 : mult.power(orig.power().floor()),
209 23 : mult_sq);
210 : }
211 :
212 : template <typename Int>
213 : constexpr
214 : outer_scalefactor_result<Int>
215 3 : bpu_product_inplace(bpu<Int> * p_target_bpu,
216 : const bpu<Int> & rhs_bpu_orig)
217 : {
218 3 : assert(rhs_bpu_orig.native_dim() == p_target_bpu->native_dim());
219 :
220 3 : bpu2_rescale_result<Int> rhs_bpu_rr = bpu2_rescale(rhs_bpu_orig,
221 : p_target_bpu->scalefactor());
222 :
223 3 : *p_target_bpu = bpu<Int>(p_target_bpu->native_dim(),
224 : p_target_bpu->scalefactor(),
225 3 : p_target_bpu->power() + rhs_bpu_orig.power());
226 :
227 : return outer_scalefactor_result<Int>(rhs_bpu_rr.outer_scale_exact_,
228 3 : rhs_bpu_rr.outer_scale_sq_);
229 : }
230 :
231 : template <typename Int>
232 : constexpr
233 : outer_scalefactor_result<Int>
234 20 : bpu_ratio_inplace(bpu<Int> * p_target_bpu,
235 : const bpu<Int> & rhs_bpu_orig)
236 : {
237 20 : assert(rhs_bpu_orig.native_dim() == p_target_bpu->native_dim());
238 :
239 20 : bpu2_rescale_result<Int> rhs_bpu_rr = bpu2_rescale(rhs_bpu_orig,
240 : p_target_bpu->scalefactor());
241 :
242 20 : *p_target_bpu = bpu<Int>(p_target_bpu->native_dim(),
243 : p_target_bpu->scalefactor(),
244 20 : p_target_bpu->power() - rhs_bpu_orig.power());
245 :
246 20 : return outer_scalefactor_result<Int>(power_ratio_type(1,1) / rhs_bpu_rr.outer_scale_exact_,
247 20 : 1.0 / rhs_bpu_rr.outer_scale_sq_);
248 : }
249 :
250 : template <typename Int>
251 : constexpr
252 : outer_scalefactor_result<Int>
253 3 : nu_product_inplace(natural_unit<Int> * p_target,
254 : const bpu<Int> & bpu)
255 : {
256 3 : std::size_t i = 0;
257 3 : for (; i < p_target->n_bpu(); ++i) {
258 3 : auto * p_target_bpu = &((*p_target)[i]);
259 :
260 3 : if (p_target_bpu->native_dim() == bpu.native_dim()) {
261 3 : outer_scalefactor_result<Int> retval = bpu_product_inplace(p_target_bpu, bpu);
262 :
263 3 : if (p_target_bpu->power().is_zero()) {
264 : /* dimension assoc'd with *p_target_bpu has been cancelled */
265 0 : p_target->remove_bpu(i);
266 : }
267 :
268 3 : return retval;
269 : }
270 : }
271 :
272 : /* control here: i=p_target->n_bpu()
273 : * Dimension represented by bpu does not already appear in *p_target.
274 : * Adopt bpu's scalefactor
275 : */
276 :
277 0 : p_target->push_back(bpu);
278 :
279 0 : return outer_scalefactor_result<Int>
280 0 : (ratio::ratio<Int>(1, 1) /*outer_scale_exact*/,
281 0 : 1.0 /*outer_scale_sq*/);
282 : }
283 :
284 : template <typename Int>
285 : constexpr
286 : outer_scalefactor_result<Int>
287 21 : nu_ratio_inplace(natural_unit<Int> * p_target,
288 : const bpu<Int> & bpu)
289 : {
290 21 : std::size_t i = 0;
291 23 : for (; i < p_target->n_bpu(); ++i) {
292 22 : auto * p_target_bpu = &((*p_target)[i]);
293 :
294 22 : if (p_target_bpu->native_dim() == bpu.native_dim()) {
295 20 : outer_scalefactor_result<Int> retval = bpu_ratio_inplace(p_target_bpu, bpu);
296 :
297 20 : if (p_target_bpu->power().is_zero()) {
298 : /* dimension assoc'd with *p_target_bpu has been cancelled */
299 20 : p_target->remove_bpu(i);
300 : }
301 :
302 20 : return retval;
303 : }
304 : }
305 :
306 : /* here: i=p_target->n_bpu()
307 : * Dimension represented by bpu does not already appear in *p_target.
308 : * Adopt bpu's scalefactor
309 : */
310 :
311 2 : p_target->push_back(bpu.reciprocal());
312 :
313 1 : return outer_scalefactor_result<Int>
314 1 : (ratio::ratio<Int>(1, 1) /*outer_scale_exact*/,
315 1 : 1.0 /*outer_scale_sq*/);
316 : }
317 :
318 : } /*namespace detail*/
319 :
320 : namespace nu2 {
321 : constexpr auto nanogram = bpu_array_maker<std::int64_t>::make_bpu_array(make_unit_power<std::int64_t>(bu::nanogram));
322 : constexpr auto microgram = bpu_array_maker<std::int64_t>::make_bpu_array(make_unit_power<std::int64_t>(bu::microgram));
323 : }
324 : } /*namespace qty*/
325 : } /*namespace xo*/
326 :
327 : /** end natural_unit.hpp **/
|
