Line data Source code
1 : /** @file ratio.hpp
2 : *
3 : * Author: Roland Conybeare
4 : **/
5 :
6 : #pragma once
7 :
8 : #include "ratio_concept.hpp"
9 : #include "xo/flatstring/flatstring.hpp"
10 : #include <numeric>
11 : #include <compare>
12 : //#include <type_traits>
13 :
14 : namespace xo {
15 : namespace ratio {
16 : namespace detail {
17 : /** @brief converts ratio to lowest terms when feasible
18 : *
19 : * Falls back to identity function for non-totally-ordered Ratio::component_type
20 : */
21 : template <typename Ratio, bool EnabledFlag = std::totally_ordered<typename Ratio::component_type>>
22 : struct reducer_type;
23 :
24 : /** @brief promote value to ratio type **/
25 : template <typename Ratio, typename FromType, bool FromRatioFlag = ratio_concept<FromType>>
26 : struct promoter_type;
27 : }
28 :
29 : /** @class ratio
30 : * @brief represent a ratio of two Int values.
31 : **/
32 : template <typename Int> requires std::totally_ordered<Int>
33 : struct ratio
34 : {
35 : public:
36 : /** @defgroup ratio-types **/
37 : ///@{
38 : /** @brief representation for (numerator, denominator) **/
39 : using component_type = Int;
40 : ///@}
41 :
42 : public:
43 : /** @defgroup ratio-ctor **/
44 : ///@{
45 : /** @brief construct ratio 0/1 **/
46 : constexpr ratio() = default;
47 : /** @brief construct ratio with numerator @p n and denominator @p d.
48 : *
49 : * Ratio need not be normalized
50 : **/
51 23 : constexpr ratio(Int n, Int d) : num_{n}, den_{d} {}
52 : ///@}
53 :
54 : /** @defgroup ratio-static-methods **/
55 : ///@{
56 : /** @brief ratio in lowest commono terms
57 : *
58 : **/
59 : static constexpr ratio reduce(Int n, Int d) {
60 : return ratio(n, d).normalize();
61 : }
62 :
63 : /** @brief add ratios @p x and @p y
64 : *
65 : * @post result ratio is normalized
66 : **/
67 312 : static constexpr ratio add(const ratio & x,
68 : const ratio & y) {
69 : /* (a/b) + (c/d)
70 : * = a.d / (b.d) + b.c / (b.d)
71 : * = (a.d + b.c) / (b.d)
72 : */
73 :
74 312 : auto a = x.num();
75 312 : auto b = x.den();
76 312 : auto c = y.num();
77 312 : auto d = y.den();
78 :
79 312 : auto num = a*d + b*c;
80 312 : auto den = b*d;
81 :
82 624 : return ratio(num, den).maybe_reduce();
83 : }
84 :
85 : /** @brief subtract ratio @p y from ratio @p x
86 : *
87 : * @post result ratio is normalized
88 : **/
89 156 : static constexpr ratio subtract(const ratio & x,
90 : const ratio & y) {
91 312 : return add(x, y.negate());
92 : }
93 :
94 : /** @brief multiply ratios @p x and @p y
95 : *
96 : * @post result ratio is normalized
97 : **/
98 690 : static constexpr ratio multiply(const ratio & x,
99 : const ratio & y) {
100 : /* (a/b) * (c/d) = a.c / b.d */
101 :
102 : /* if x,y normalized,
103 : * opportunity to cancel common factor between (a, d) or (c, b)
104 : *
105 : * want to do this before multiplying to avoid overflow involving intermediate terms
106 : */
107 :
108 690 : auto a1 = x.num();
109 690 : auto b1 = x.den();
110 690 : auto c1 = y.num();
111 690 : auto d1 = y.den();
112 :
113 1380 : auto ad_gcf = std::gcd(a1, d1);
114 690 : auto bc_gcf = std::gcd(b1, c1);
115 :
116 690 : auto a = a1 / ad_gcf;
117 690 : auto b = b1 / bc_gcf;
118 690 : auto c = c1 / bc_gcf;
119 690 : auto d = d1 / ad_gcf;
120 :
121 690 : auto num = a*c;
122 690 : auto den = b*d;
123 :
124 690 : return ratio(num, den).maybe_reduce();
125 : }
126 :
127 : /** @brief divide ratio @p y into ratio @p x
128 : *
129 : * @post result ratio is normalized
130 : **/
131 : static constexpr ratio divide(const ratio & x,
132 : const ratio & y)
133 : {
134 : return multiply(x, y.reciprocal());
135 : }
136 :
137 : /** @brief 3-way compare two ratios **/
138 1092 : static constexpr auto compare(ratio x, ratio y) {
139 : /* ensure minus signs in numerators only */
140 1092 : if (x.den() < 0)
141 0 : return compare_aux(ratio(-x.num(), -x.den()), y);
142 1092 : if (y.den() < 0)
143 0 : return compare_aux(x, ratio(-y.num(), -y.den()));
144 :
145 1092 : return compare_aux(x, y);
146 : }
147 : ///@}
148 :
149 : /** @defgroup ratio-access **/
150 : ///@{
151 : /** @brief fetch ratio's numerator **/
152 3301 : constexpr Int num() const noexcept { return num_; }
153 : /** @brief fetch ratio's denominator **/
154 7289 : constexpr Int den() const noexcept { return den_; }
155 :
156 : /** @brief true if and only if ratio is equal to zero **/
157 : constexpr bool is_zero() const noexcept { return (num_ == 0) && (den_ != 0); }
158 :
159 : /** @brief true if and only if ratio is equal to one **/
160 : constexpr bool is_unity() const noexcept { return (num_ != 0) && (num_ == den_); }
161 :
162 : /** @brief true if and only if ratio represents an integer
163 : *
164 : * (denominator is +/- 1)
165 : **/
166 156 : constexpr bool is_integer() const noexcept { return den_ == 1 || den_ == -1; }
167 : ///@}
168 :
169 : /** @defgroup ratio-methods **/
170 : ///@{
171 : /** @brief r.negate() is the exact ratio representing @c -r **/
172 156 : constexpr ratio negate() const { return ratio(-num_, den_); }
173 : /** @brief r.reciprocal() is the eact ratio representing @c 1/r **/
174 62 : constexpr ratio reciprocal() const { return ratio(den_, num_); }
175 :
176 : /** @brief r.floor() is the largest integer x : x <= r **/
177 : constexpr Int floor() const { return (num_ / den_); }
178 :
179 : /** @brief r.ceil() is the smallest integer x : r < x. **/
180 : constexpr Int ceil() const { return floor() + 1; }
181 :
182 : /** @brief reduce to lowest terms
183 : *
184 : * @pre @c Int type must be totally ordered
185 : **/
186 1047 : constexpr ratio normalize() const requires std::totally_ordered<Int> {
187 1047 : if (den_ < 0)
188 23 : return ratio(-num_, -den_).normalize();
189 :
190 1024 : auto factor = std::gcd(num_, den_);
191 :
192 1024 : return ratio(num_ / factor,
193 1024 : den_ / factor);
194 : }
195 :
196 : /** @brief reduce to lowest terms, if Int representation admits
197 : *
198 : * Otherwise fallback to identity function
199 : **/
200 1002 : constexpr ratio maybe_reduce() const {
201 1002 : return detail::reducer_type<ratio>::attempt_reduce(*this);
202 : }
203 :
204 : /** @brief return fractional part of this ratio
205 : *
206 : * @pre @c Int type must be totally ordered
207 : **/
208 : constexpr ratio frac() const requires std::totally_ordered<Int> {
209 : return ratio::subtract(*this, ratio(this->floor(), 1));
210 : }
211 :
212 : /** @brief compute integer exponent @p p of this ratio
213 : *
214 : * @note time complexity is O(log p)
215 : **/
216 218 : constexpr ratio power(int p) const {
217 :
218 218 : constexpr ratio retval = ratio(1, 1);
219 :
220 218 : if (p == 0)
221 27 : return ratio(1, 1);
222 :
223 191 : if (p < 0)
224 62 : return this->reciprocal().power(-p);
225 :
226 : /* inv: x^p = aj.xj^pj */
227 129 : ratio aj = ratio(1, 1);
228 129 : ratio xj = *this;
229 129 : int pj = p;
230 :
231 507 : while (pj > 0) {
232 378 : if (pj % 2 == 0) {
233 : /* a.x^(2q) = a.(x^2)^q */
234 177 : xj = xj * xj;
235 177 : pj = pj / 2;
236 : } else {
237 : /* a.x^(2q+1) = (a.x).x^(2q) */
238 201 : aj = aj * xj;
239 201 : pj = (pj - 1);
240 : }
241 : }
242 :
243 : /* pj = 0, so: x^p = aj.xj^pj = aj.xj^0 = aj */
244 129 : return aj;
245 : }
246 :
247 : /** @brief negate operator **/
248 156 : constexpr ratio operator-() const { return ratio(-num_, den_); }
249 : ///@}
250 :
251 : /** @defgroup ratio-conversion **/
252 : ///@{
253 : /** @brief convert to non-ratio representation
254 : *
255 : * For example: to int or double
256 : **/
257 : template <typename Repr>
258 : constexpr Repr to() const noexcept { return num_ / static_cast<Repr>(den_); }
259 :
260 : /** @brief convert to short human-friendly flatstring representation
261 : *
262 : * Example:
263 : * @code
264 : * ratio(7,1).to_str<5>(); // "7"
265 : * ratio(1,7).to_str<5>(); // "(1/7)"
266 : * ratio(-1,7).to_str<10>(); // "(-1/7)"
267 : * ratio(-1,7).to_str<5>(); // "(-1/"
268 : * @endcode
269 : **/
270 : template <std::size_t N>
271 156 : constexpr flatstring<N> to_str() const noexcept {
272 156 : if (this->is_integer()) {
273 26 : return flatstring<N>::from_int(num_);
274 : } else {
275 130 : auto num_str = flatstring<N>::from_int(num_);
276 130 : auto den_str = flatstring<N>::from_int(den_);
277 :
278 : /* tmp capacity will be about 2N+3 */
279 260 : auto tmp = flatstring_concat(flatstring("("),
280 : num_str,
281 130 : flatstring("/"),
282 : den_str,
283 260 : flatstring(")"));
284 :
285 130 : flatstring<N> retval;
286 130 : retval.assign(tmp);
287 :
288 130 : return retval;
289 : }
290 130 : }
291 :
292 : /** @brief convert to representation using different integer types **/
293 : template <typename Ratio2>
294 : constexpr operator Ratio2 () const noexcept requires ratio_concept<Ratio2> {
295 : return Ratio2(num_, den_);
296 : }
297 : ///@}
298 :
299 : private:
300 : /** @brief 3-way compare auxiliary function.
301 : *
302 : * @pre @p x, @p y have non-negative denominator
303 : **/
304 1092 : static constexpr auto compare_aux(ratio x, ratio y) noexcept {
305 : /* control here: b>=0, d>=0 */
306 :
307 : /* (a/b) <=> (c/d)
308 : * (a.d/b) <=> c no sign change, since d >= 0
309 : * (a.d) <=> (b.c) no sign change, since b >= 0
310 : */
311 1092 : auto a = x.num();
312 1092 : auto b = x.den();
313 1092 : auto c = y.num();
314 1092 : auto d = y.den();
315 :
316 1092 : auto lhs = a*d;
317 1092 : auto rhs = b*c;
318 :
319 1092 : return lhs <=> rhs;
320 : }
321 :
322 : private:
323 : /** @defgroup ratio-instance-variables **/
324 : ///@{
325 : /** @brief numerator **/
326 : Int num_ = 0;
327 : /** @brief denominator **/
328 : Int den_ = 1;
329 : ///@}
330 : };
331 :
332 : namespace detail {
333 : template <typename Ratio, bool EnabledFlag>
334 : struct reducer_type {};
335 :
336 : template <typename Ratio>
337 : struct reducer_type<Ratio, true /*EnabledFlag*/> {
338 1002 : static constexpr Ratio attempt_reduce(Ratio x) { return x.normalize(); }
339 : };
340 :
341 : template <typename Ratio>
342 : struct reducer_type<Ratio, false /*!EnabledFlag*/> {
343 : static constexpr Ratio attempt_reduce(Ratio x) { return x; }
344 : };
345 : }
346 :
347 : namespace detail {
348 : template <typename Ratio, typename FromType, bool FromRatioFlag>
349 : struct promoter_type;
350 :
351 : template <typename Ratio, typename FromType>
352 : struct promoter_type<Ratio, FromType, true /*FromRatioFlag*/> {
353 : /* to 'promote' a ratio, rely on its conversion operator */
354 : static constexpr Ratio promote(FromType x) { return x; }
355 : };
356 :
357 : template <typename Ratio, typename FromType>
358 : struct promoter_type<Ratio, FromType, false /*!FromRatioFlag*/> {
359 : /* to 'promote' a non-ratio, use denominator=1 */
360 : static constexpr Ratio promote(FromType x) { return Ratio(x, 1); }
361 : };
362 : }
363 :
364 : /** @brief create a ratio in lowest terms from two integers **/
365 : template <typename Int1, typename Int2>
366 : constexpr auto
367 : make_ratio (Int1 n, Int2 d = 1) -> ratio<std::common_type_t<Int1, Int2>>
368 : {
369 : return ratio<std::common_type_t<Int1, Int2>>(n, d).maybe_reduce();
370 : }
371 :
372 : namespace detail {
373 : /** @brief auxiliary function for binary ratio operations
374 : *
375 : * Support binary ratio operations on combinations:
376 : * - (ratio<T>, ratio<U>)
377 : * - (ratio<T>, U) // where U is not a ratio
378 : * - (T, ratio(U)) // where T is not a ratio
379 : *
380 : * Goals:
381 : *
382 : * 1. Support expressions like
383 : *
384 : * @code
385 : * auto x = 1 + make_ratio(2,3);
386 : * @endcode
387 : *
388 : * 2. promote to wider types as needed
389 : *
390 : * @code
391 : * auto x = make_ratio(2,3) + make_ratio(1ul,2ul);
392 : * static_assert(std::same_as<x::component_type, unsigned long>);
393 : * @endcode
394 : *
395 : * 3. avoid interfering with other templates that may overload operator+
396 : *
397 : * @pre at least one of (Left,Right) must be known to be a ratio
398 : **/
399 : template <typename Left,
400 : typename Right,
401 : bool LeftIsRatio = ratio_concept<Left>,
402 : bool RightIsRatio = ratio_concept<Right>>
403 : struct op_aux_type;
404 :
405 : /** @brief specialization for two ratio types **/
406 : template <typename LeftRatio,
407 : typename RightRatio>
408 : requires (ratio_concept<LeftRatio> && ratio_concept<RightRatio>)
409 : struct op_aux_type<LeftRatio, RightRatio, true /*LeftIsRatio*/, true /*RightIsRatio*/> {
410 : using component_type = std::common_type_t<typename LeftRatio::component_type,
411 : typename RightRatio::component_type>;
412 :
413 : using ratio_type = ratio<component_type>;
414 :
415 156 : static constexpr ratio_type add (const LeftRatio & x,
416 : const RightRatio & y)
417 : {
418 156 : return ratio_type::add(x, y);
419 : }
420 :
421 156 : static constexpr ratio_type subtract (const LeftRatio & x,
422 : const RightRatio & y)
423 : {
424 156 : return ratio_type::subtract(x, y);
425 : }
426 :
427 690 : static constexpr ratio_type multiply (const LeftRatio & x,
428 : const RightRatio & y)
429 : {
430 690 : return ratio_type::multiply(x, y);
431 : }
432 :
433 : static constexpr ratio_type divide (const LeftRatio & x,
434 : const RightRatio & y)
435 : {
436 : return ratio_type::divide(x, y);
437 : }
438 :
439 1092 : static constexpr auto compare (const LeftRatio & x,
440 : const RightRatio & y)
441 : {
442 1092 : return ratio_type::compare(x, y);
443 : }
444 : };
445 :
446 : /** @brief specialization for left-hand ratio and right-hand integer value **/
447 : template <typename LeftRatio,
448 : typename Right>
449 : requires (ratio_concept<LeftRatio> && !ratio_concept<Right>)
450 : struct op_aux_type<LeftRatio, Right, true /*LeftIsRatio*/, false /*RightIsRatio*/> {
451 : using component_type = std::common_type_t<typename LeftRatio::component_type, Right>;
452 :
453 : using ratio_type = ratio<component_type>;
454 :
455 : static constexpr ratio_type add (const LeftRatio & x,
456 : const Right & y)
457 : {
458 : /* reminder: adding an integer can't introduce reduced terms */
459 : return ratio_type(x.num() + x.den() * y, x.den());
460 : }
461 :
462 : static constexpr ratio_type subtract (const LeftRatio & x,
463 : const Right & y)
464 : {
465 : /* reminder: subtracting an integer can't introduce reduced terms */
466 : return ratio_type(x.num() - x.den() * y, x.den());
467 : }
468 :
469 : static constexpr ratio_type multiply (const LeftRatio & x,
470 : const Right & yp)
471 : {
472 : auto gcf = std::gcd(x.den(), yp);
473 :
474 : auto a = x.num();
475 : auto b = x.den() / gcf;
476 : auto y = yp / gcf;
477 :
478 : return ratio_type(a*y, b);
479 : }
480 :
481 : static constexpr ratio_type divide (const LeftRatio & x,
482 : const Right & yp)
483 : {
484 : auto gcf = std::gcd(x.num(), yp);
485 :
486 : auto a = x.num() / gcf;
487 : auto b = x.den();
488 : auto y = yp / gcf;
489 :
490 : return ratio_type(a*y, b);
491 : }
492 :
493 : static constexpr auto compare (const LeftRatio & x,
494 : const Right & y)
495 : {
496 : /* note: in c++26 std::signof is constexpr, usable here */
497 : if (x.den() >= 0)
498 : return compare_aux(x, y);
499 : else
500 : return compare_aux(LeftRatio(-x.num(), -x.den()), y);
501 : }
502 :
503 : private:
504 : static constexpr auto compare_aux (const LeftRatio & x, const Right & y) {
505 : return (x.num() <=> x.den() * y);
506 : }
507 : };
508 :
509 : /** @brief specialization for left-hand integer value and right-hand ratio **/
510 : template <typename Left,
511 : typename RightRatio>
512 : requires (!ratio_concept<Left> && ratio_concept<RightRatio>)
513 : struct op_aux_type<Left, RightRatio, false /*LeftIsRatio*/, true /*RightIsRatio*/> {
514 : using component_type = std::common_type_t<Left, typename RightRatio::component_type>;
515 :
516 : using ratio_type = ratio<component_type>;
517 :
518 : static constexpr ratio_type add(const Left & x,
519 : const RightRatio & y)
520 : {
521 : /* reminder: adding an integer can't introduce reduced terms */
522 : return ratio_type(x * y.den() + y.num(), y.den());
523 : }
524 :
525 : static constexpr ratio_type subtract(const Left & x,
526 : const RightRatio & y)
527 : {
528 : /* reminder: subtracting an integer can't introduce reduced terms */
529 : return ratio_type(x * y.den() - y.num(), y.den());
530 : }
531 :
532 : static constexpr ratio_type multiply (const Left & xp,
533 : const RightRatio & y)
534 : {
535 : auto gcf = std::gcd(xp, y.den());
536 :
537 : auto x = xp / gcf;
538 : auto c = y.num();
539 : auto d = y.den() / gcf;
540 :
541 : return ratio_type(x*c, d);
542 : }
543 :
544 : static constexpr ratio_type divide (const Left & x,
545 : const RightRatio & y)
546 : {
547 : return multiply(x, y.reciprocal());
548 : }
549 :
550 : static constexpr auto compare(const Left & x,
551 : const RightRatio & y)
552 : {
553 : if (y.den() >= 0)
554 : return compare_aux(x, y);
555 : else
556 : return compare_aux(x, RightRatio(-y.num(), -y.den()));
557 : }
558 :
559 : private:
560 : static constexpr auto compare_aux (const Left & x,
561 : const RightRatio & y)
562 : {
563 : return (x * y.den() <=> y.num());
564 : };
565 : };
566 : } /*namespace detail*/
567 :
568 : /** @defgroup ratio-arithmetic **/
569 : ///@{
570 : /** @brief add two ratios.
571 : *
572 : * One argument may be a non-ratio type if it can be promoted to a ratio
573 : **/
574 : template <typename Ratio1, typename Ratio2>
575 : inline constexpr auto
576 156 : operator+ (const Ratio1 & x, const Ratio2 & y)
577 : requires (ratio_concept<Ratio1> || ratio_concept<Ratio2>)
578 : {
579 156 : return detail::op_aux_type<Ratio1, Ratio2>::add(x, y);
580 : }
581 :
582 : /** @brief subtract two ratios.
583 : *
584 : * One argument may be a non-ratio type if it can be promoted to a ratio
585 : **/
586 : template <typename Ratio1, typename Ratio2>
587 : inline constexpr auto
588 156 : operator- (const Ratio1 & x, const Ratio2 & y)
589 : requires (ratio_concept<Ratio1> || ratio_concept<Ratio2>)
590 : {
591 156 : return detail::op_aux_type<Ratio1, Ratio2>::subtract(x, y);
592 : }
593 :
594 : /** @brief multiply two ratios
595 : *
596 : * One argument may be a non-ratio type if it can be promoted to a ratio
597 : **/
598 : template <typename Ratio1, typename Ratio2>
599 : inline constexpr auto
600 690 : operator* (const Ratio1 & x, const Ratio2 & y)
601 : requires (ratio_concept<Ratio1> || ratio_concept<Ratio2>)
602 : {
603 690 : return detail::op_aux_type<Ratio1, Ratio2>::multiply(x, y);
604 : }
605 :
606 : /** @brief divide two ratios
607 : *
608 : * One argument may be a non-ratio type if it can be promoted to a ratio
609 : **/
610 : template <typename Ratio1, typename Ratio2>
611 : inline constexpr auto
612 : operator/ (const Ratio1 & x, const Ratio2 & y)
613 : requires (ratio_concept<Ratio1> || ratio_concept<Ratio2>)
614 : {
615 : return detail::op_aux_type<Ratio1, Ratio2>::divide(x, y);
616 : }
617 : ///@}
618 :
619 : /** @defgroup ratio-3way-compare 3way comparison **/
620 : ///@{
621 : /** @brief compare two ratios for equality
622 : *
623 : * One argument may be a non-ratio type if it can be promoted to a ratio
624 : **/
625 : template <typename Ratio1, typename Ratio2>
626 : inline constexpr bool
627 312 : operator== (const Ratio1 & x, const Ratio2 & y)
628 : requires (ratio_concept<Ratio1> || ratio_concept<Ratio2>)
629 : {
630 312 : return (detail::op_aux_type<Ratio1, Ratio2>::compare(x, y) == 0);
631 : }
632 :
633 : /** @brief compare two ratios
634 : *
635 : * One argument may be a non-ratio type if it can be promoted to a ratio
636 : **/
637 : template <typename Ratio1, typename Ratio2>
638 : inline constexpr auto
639 780 : operator<=> (const Ratio1 & x, const Ratio2 & y)
640 : requires (ratio_concept<Ratio1> || ratio_concept<Ratio2>)
641 : {
642 780 : return detail::op_aux_type<Ratio1, Ratio2>::compare(x, y);
643 : }
644 : ///@}
645 :
646 : } /*namespace ratio*/
647 : } /*namespace xo*/
648 :
649 : /** end ratio.hpp **/
|