1///////////////////////////////////////////////////////////////
2// Copyright 2020 John Maddock. Distributed under the Boost
3// Software License, Version 1.0. (See accompanying file
4// LICENSE_1_0.txt or copy at https://www.boost.org/LICENSE_1_0.txt
5
6#ifndef BOOST_MP_RATIONAL_ADAPTOR_HPP
7#define BOOST_MP_RATIONAL_ADAPTOR_HPP
8
9#include <boost/multiprecision/number.hpp>
10#include <boost/multiprecision/detail/hash.hpp>
11#include <boost/multiprecision/detail/float128_functions.hpp>
12#include <boost/multiprecision/detail/no_exceptions_support.hpp>
13
14namespace boost {
15namespace multiprecision {
16namespace backends {
17
18template <class Backend>
19struct rational_adaptor
20{
21 //
22 // Each backend need to declare 3 type lists which declare the types
23 // with which this can interoperate. These lists must at least contain
24 // the widest type in each category - so "long long" must be the final
25 // type in the signed_types list for example. Any narrower types if not
26 // present in the list will get promoted to the next wider type that is
27 // in the list whenever mixed arithmetic involving that type is encountered.
28 //
29 typedef typename Backend::signed_types signed_types;
30 typedef typename Backend::unsigned_types unsigned_types;
31 typedef typename Backend::float_types float_types;
32
33 typedef typename std::tuple_element<0, unsigned_types>::type ui_type;
34
35 static Backend get_one()
36 {
37 Backend t;
38 t = static_cast<ui_type>(1);
39 return t;
40 }
41 static Backend get_zero()
42 {
43 Backend t;
44 t = static_cast<ui_type>(0);
45 return t;
46 }
47
48 static const Backend& one()
49 {
50 static const Backend result(get_one());
51 return result;
52 }
53 static const Backend& zero()
54 {
55 static const Backend result(get_zero());
56 return result;
57 }
58
59 void normalize()
60 {
61 using default_ops::eval_gcd;
62 using default_ops::eval_eq;
63 using default_ops::eval_divide;
64 using default_ops::eval_get_sign;
65
66 int s = eval_get_sign(m_denom);
67
68 if(s == 0)
69 {
70 BOOST_MP_THROW_EXCEPTION(std::overflow_error("Integer division by zero"));
71 }
72 else if (s < 0)
73 {
74 m_num.negate();
75 m_denom.negate();
76 }
77
78 Backend g, t;
79 eval_gcd(g, m_num, m_denom);
80 if (!eval_eq(g, one()))
81 {
82 eval_divide(t, m_num, g);
83 m_num.swap(t);
84 eval_divide(t, m_denom, g);
85 m_denom = std::move(t);
86 }
87 }
88
89 // We must have a default constructor:
90 rational_adaptor()
91 : m_num(zero()), m_denom(one()) {}
92
93 rational_adaptor(const rational_adaptor& o) : m_num(o.m_num), m_denom(o.m_denom) {}
94 rational_adaptor(rational_adaptor&& o) = default;
95
96 // Optional constructors, we can make this type slightly more efficient
97 // by providing constructors from any type we can handle natively.
98 // These will also cause number<> to be implicitly constructible
99 // from these types unless we make such constructors explicit.
100 //
101 template <class Arithmetic>
102 rational_adaptor(const Arithmetic& val, typename std::enable_if<std::is_constructible<Backend, Arithmetic>::value && !std::is_floating_point<Arithmetic>::value>::type const* = nullptr)
103 : m_num(val), m_denom(one()) {}
104
105 //
106 // Pass-through 2-arg construction of components:
107 //
108 template <class T, class U>
109 rational_adaptor(const T& a, const U& b, typename std::enable_if<std::is_constructible<Backend, T const&>::value && std::is_constructible<Backend, U const&>::value>::type const* = nullptr)
110 : m_num(a), m_denom(b)
111 {
112 normalize();
113 }
114 template <class T, class U>
115 rational_adaptor(T&& a, const U& b, typename std::enable_if<std::is_constructible<Backend, T>::value && std::is_constructible<Backend, U>::value>::type const* = nullptr)
116 : m_num(static_cast<T&&>(a)), m_denom(b)
117 {
118 normalize();
119 }
120 template <class T, class U>
121 rational_adaptor(T&& a, U&& b, typename std::enable_if<std::is_constructible<Backend, T>::value && std::is_constructible<Backend, U>::value>::type const* = nullptr)
122 : m_num(static_cast<T&&>(a)), m_denom(static_cast<U&&>(b))
123 {
124 normalize();
125 }
126 template <class T, class U>
127 rational_adaptor(const T& a, U&& b, typename std::enable_if<std::is_constructible<Backend, T>::value && std::is_constructible<Backend, U>::value>::type const* = nullptr)
128 : m_num(a), m_denom(static_cast<U&&>(b))
129 {
130 normalize();
131 }
132 //
133 // In the absense of converting constructors, operator= takes the strain.
134 // In addition to the usual suspects, there must be one operator= for each type
135 // listed in signed_types, unsigned_types, and float_types plus a string constructor.
136 //
137 rational_adaptor& operator=(const rational_adaptor& o) = default;
138 rational_adaptor& operator=(rational_adaptor&& o) = default;
139 template <class Arithmetic>
140 inline typename std::enable_if<!std::is_floating_point<Arithmetic>::value, rational_adaptor&>::type operator=(const Arithmetic& i)
141 {
142 m_num = i;
143 m_denom = one();
144 return *this;
145 }
146 rational_adaptor& operator=(const char* s)
147 {
148 using default_ops::eval_eq;
149
150 std::string s1;
151 multiprecision::number<Backend> v1, v2;
152 char c;
153 bool have_hex = false;
154 const char* p = s; // saved for later
155
156 while ((0 != (c = *s)) && (c == 'x' || c == 'X' || c == '-' || c == '+' || (c >= '0' && c <= '9') || (have_hex && (c >= 'a' && c <= 'f')) || (have_hex && (c >= 'A' && c <= 'F'))))
157 {
158 if (c == 'x' || c == 'X')
159 have_hex = true;
160 s1.append(n: 1, c: c);
161 ++s;
162 }
163 v1.assign(s1);
164 s1.erase();
165 if (c == '/')
166 {
167 ++s;
168 while ((0 != (c = *s)) && (c == 'x' || c == 'X' || c == '-' || c == '+' || (c >= '0' && c <= '9') || (have_hex && (c >= 'a' && c <= 'f')) || (have_hex && (c >= 'A' && c <= 'F'))))
169 {
170 if (c == 'x' || c == 'X')
171 have_hex = true;
172 s1.append(n: 1, c: c);
173 ++s;
174 }
175 v2.assign(s1);
176 }
177 else
178 v2 = 1;
179 if (*s)
180 {
181 BOOST_MP_THROW_EXCEPTION(std::runtime_error(std::string("Could not parse the string \"") + p + std::string("\" as a valid rational number.")));
182 }
183 multiprecision::number<Backend> gcd;
184 eval_gcd(gcd.backend(), v1.backend(), v2.backend());
185 if (!eval_eq(gcd.backend(), one()))
186 {
187 v1 /= gcd;
188 v2 /= gcd;
189 }
190 num() = std::move(std::move(v1).backend());
191 denom() = std::move(std::move(v2).backend());
192 return *this;
193 }
194 template <class Float>
195 typename std::enable_if<std::is_floating_point<Float>::value, rational_adaptor&>::type operator=(Float i)
196 {
197 using default_ops::eval_eq;
198 BOOST_MP_FLOAT128_USING using std::floor; using std::frexp; using std::ldexp;
199
200 int e;
201 Float f = frexp(i, &e);
202#ifdef BOOST_HAS_FLOAT128
203 f = ldexp(f, std::is_same<float128_type, Float>::value ? 113 : std::numeric_limits<Float>::digits);
204 e -= std::is_same<float128_type, Float>::value ? 113 : std::numeric_limits<Float>::digits;
205#else
206 f = ldexp(f, std::numeric_limits<Float>::digits);
207 e -= std::numeric_limits<Float>::digits;
208#endif
209 number<Backend> num(f);
210 number<Backend> denom(1u);
211 if (e > 0)
212 {
213 num <<= e;
214 }
215 else if (e < 0)
216 {
217 denom <<= -e;
218 }
219 number<Backend> gcd;
220 eval_gcd(gcd.backend(), num.backend(), denom.backend());
221 if (!eval_eq(gcd.backend(), one()))
222 {
223 num /= gcd;
224 denom /= gcd;
225 }
226 this->num() = std::move(std::move(num).backend());
227 this->denom() = std::move(std::move(denom).backend());
228 return *this;
229 }
230
231 void swap(rational_adaptor& o)
232 {
233 m_num.swap(o.m_num);
234 m_denom.swap(o.m_denom);
235 }
236 std::string str(std::streamsize digits, std::ios_base::fmtflags f) const
237 {
238 using default_ops::eval_eq;
239 //
240 // We format the string ourselves so we can match what GMP's mpq type does:
241 //
242 std::string result = num().str(digits, f);
243 if (!eval_eq(denom(), one()))
244 {
245 result.append(n: 1, c: '/');
246 result.append(denom().str(digits, f));
247 }
248 return result;
249 }
250 void negate()
251 {
252 m_num.negate();
253 }
254 int compare(const rational_adaptor& o) const
255 {
256 std::ptrdiff_t s1 = eval_get_sign(*this);
257 std::ptrdiff_t s2 = eval_get_sign(o);
258 if (s1 != s2)
259 {
260 return s1 < s2 ? -1 : 1;
261 }
262 else if (s1 == 0)
263 return 0; // both zero.
264
265 bool neg = false;
266 if (s1 >= 0)
267 {
268 s1 = eval_msb(num()) + eval_msb(o.denom());
269 s2 = eval_msb(o.num()) + eval_msb(denom());
270 }
271 else
272 {
273 Backend t(num());
274 t.negate();
275 s1 = eval_msb(t) + eval_msb(o.denom());
276 t = o.num();
277 t.negate();
278 s2 = eval_msb(t) + eval_msb(denom());
279 neg = true;
280 }
281 s1 -= s2;
282 if (s1 < -1)
283 return neg ? 1 : -1;
284 else if (s1 > 1)
285 return neg ? -1 : 1;
286
287 Backend t1, t2;
288 eval_multiply(t1, num(), o.denom());
289 eval_multiply(t2, o.num(), denom());
290 return t1.compare(t2);
291 }
292 //
293 // Comparison with arithmetic types, default just constructs a temporary:
294 //
295 template <class A>
296 typename std::enable_if<boost::multiprecision::detail::is_arithmetic<A>::value, int>::type compare(A i) const
297 {
298 rational_adaptor t;
299 t = i; // Note: construct directly from i if supported.
300 return compare(t);
301 }
302
303 Backend& num() { return m_num; }
304 const Backend& num()const { return m_num; }
305 Backend& denom() { return m_denom; }
306 const Backend& denom()const { return m_denom; }
307
308 #ifndef BOOST_MP_STANDALONE
309 template <class Archive>
310 void serialize(Archive& ar, const std::integral_constant<bool, true>&)
311 {
312 // Saving
313 number<Backend> n(num()), d(denom());
314 ar& boost::make_nvp("numerator", n);
315 ar& boost::make_nvp("denominator", d);
316 }
317 template <class Archive>
318 void serialize(Archive& ar, const std::integral_constant<bool, false>&)
319 {
320 // Loading
321 number<Backend> n, d;
322 ar& boost::make_nvp("numerator", n);
323 ar& boost::make_nvp("denominator", d);
324 num() = n.backend();
325 denom() = d.backend();
326 }
327 template <class Archive>
328 void serialize(Archive& ar, const unsigned int /*version*/)
329 {
330 using tag = typename Archive::is_saving;
331 using saving_tag = std::integral_constant<bool, tag::value>;
332 serialize(ar, saving_tag());
333 }
334 #endif // BOOST_MP_STANDALONE
335
336 private:
337 Backend m_num, m_denom;
338};
339
340//
341// Helpers:
342//
343template <class T>
344inline constexpr typename std::enable_if<std::numeric_limits<T>::is_specialized && !std::numeric_limits<T>::is_signed, bool>::type
345is_minus_one(const T&)
346{
347 return false;
348}
349template <class T>
350inline constexpr typename std::enable_if<!std::numeric_limits<T>::is_specialized || std::numeric_limits<T>::is_signed, bool>::type
351is_minus_one(const T& val)
352{
353 return val == -1;
354}
355
356//
357// Required non-members:
358//
359template <class Backend>
360inline void eval_add(rational_adaptor<Backend>& a, const rational_adaptor<Backend>& b)
361{
362 eval_add_subtract_imp(a, a, b, true);
363}
364template <class Backend>
365inline void eval_subtract(rational_adaptor<Backend>& a, const rational_adaptor<Backend>& b)
366{
367 eval_add_subtract_imp(a, a, b, false);
368}
369
370template <class Backend>
371inline void eval_multiply(rational_adaptor<Backend>& a, const rational_adaptor<Backend>& b)
372{
373 eval_multiply_imp(a, a, b.num(), b.denom());
374}
375
376template <class Backend>
377void eval_divide(rational_adaptor<Backend>& a, const rational_adaptor<Backend>& b)
378{
379 using default_ops::eval_divide;
380 rational_adaptor<Backend> t;
381 eval_divide(t, a, b);
382 a = std::move(t);
383}
384//
385// Conversions:
386//
387template <class R, class IntBackend>
388inline typename std::enable_if<number_category<R>::value == number_kind_floating_point>::type eval_convert_to(R* result, const rational_adaptor<IntBackend>& backend)
389{
390 //
391 // The generic conversion is as good as anything we can write here:
392 //
393 ::boost::multiprecision::detail::generic_convert_rational_to_float(*result, backend);
394}
395
396template <class R, class IntBackend>
397inline typename std::enable_if<(number_category<R>::value != number_kind_integer) && (number_category<R>::value != number_kind_floating_point) && !std::is_enum<R>::value>::type eval_convert_to(R* result, const rational_adaptor<IntBackend>& backend)
398{
399 using default_ops::eval_convert_to;
400 R d;
401 eval_convert_to(result, backend.num());
402 eval_convert_to(&d, backend.denom());
403 *result /= d;
404}
405
406template <class R, class Backend>
407inline typename std::enable_if<number_category<R>::value == number_kind_integer>::type eval_convert_to(R* result, const rational_adaptor<Backend>& backend)
408{
409 using default_ops::eval_divide;
410 using default_ops::eval_convert_to;
411 Backend t;
412 eval_divide(t, backend.num(), backend.denom());
413 eval_convert_to(result, t);
414}
415
416//
417// Hashing support, not strictly required, but it is used in our tests:
418//
419template <class Backend>
420inline std::size_t hash_value(const rational_adaptor<Backend>& arg)
421{
422 std::size_t result = hash_value(arg.num());
423 std::size_t result2 = hash_value(arg.denom());
424 boost::multiprecision::detail::hash_combine(seed&: result, v: result2);
425 return result;
426}
427//
428// assign_components:
429//
430template <class Backend>
431void assign_components(rational_adaptor<Backend>& result, Backend const& a, Backend const& b)
432{
433 using default_ops::eval_gcd;
434 using default_ops::eval_divide;
435 using default_ops::eval_eq;
436 using default_ops::eval_is_zero;
437 using default_ops::eval_get_sign;
438
439 if (eval_is_zero(b))
440 {
441 BOOST_MP_THROW_EXCEPTION(std::overflow_error("Integer division by zero"));
442 }
443 Backend g;
444 eval_gcd(g, a, b);
445 if (eval_eq(g, rational_adaptor<Backend>::one()))
446 {
447 result.num() = a;
448 result.denom() = b;
449 }
450 else
451 {
452 eval_divide(result.num(), a, g);
453 eval_divide(result.denom(), b, g);
454 }
455 if (eval_get_sign(result.denom()) < 0)
456 {
457 result.num().negate();
458 result.denom().negate();
459 }
460}
461//
462// Again for arithmetic types, overload for whatever arithmetic types are directly supported:
463//
464template <class Backend, class Arithmetic1, class Arithmetic2>
465inline void assign_components(rational_adaptor<Backend>& result, const Arithmetic1& a, typename std::enable_if<std::is_arithmetic<Arithmetic1>::value && std::is_arithmetic<Arithmetic2>::value, const Arithmetic2&>::type b)
466{
467 using default_ops::eval_gcd;
468 using default_ops::eval_divide;
469 using default_ops::eval_eq;
470
471 if (b == 0)
472 {
473 BOOST_MP_THROW_EXCEPTION(std::overflow_error("Integer division by zero"));
474 }
475
476 Backend g;
477 result.num() = a;
478 eval_gcd(g, result.num(), b);
479 if (eval_eq(g, rational_adaptor<Backend>::one()))
480 {
481 result.denom() = b;
482 }
483 else
484 {
485 eval_divide(result.num(), g);
486 eval_divide(result.denom(), b, g);
487 }
488 if (eval_get_sign(result.denom()) < 0)
489 {
490 result.num().negate();
491 result.denom().negate();
492 }
493}
494template <class Backend, class Arithmetic1, class Arithmetic2>
495inline void assign_components(rational_adaptor<Backend>& result, const Arithmetic1& a, typename std::enable_if<!std::is_arithmetic<Arithmetic1>::value || !std::is_arithmetic<Arithmetic2>::value, const Arithmetic2&>::type b)
496{
497 using default_ops::eval_gcd;
498 using default_ops::eval_divide;
499 using default_ops::eval_eq;
500
501 Backend g;
502 result.num() = a;
503 result.denom() = b;
504
505 if (eval_get_sign(result.denom()) == 0)
506 {
507 BOOST_MP_THROW_EXCEPTION(std::overflow_error("Integer division by zero"));
508 }
509
510 eval_gcd(g, result.num(), result.denom());
511 if (!eval_eq(g, rational_adaptor<Backend>::one()))
512 {
513 eval_divide(result.num(), g);
514 eval_divide(result.denom(), g);
515 }
516 if (eval_get_sign(result.denom()) < 0)
517 {
518 result.num().negate();
519 result.denom().negate();
520 }
521}
522//
523// Optional comparison operators:
524//
525template <class Backend>
526inline bool eval_is_zero(const rational_adaptor<Backend>& arg)
527{
528 using default_ops::eval_is_zero;
529 return eval_is_zero(arg.num());
530}
531
532template <class Backend>
533inline int eval_get_sign(const rational_adaptor<Backend>& arg)
534{
535 using default_ops::eval_get_sign;
536 return eval_get_sign(arg.num());
537}
538
539template <class Backend>
540inline bool eval_eq(const rational_adaptor<Backend>& a, const rational_adaptor<Backend>& b)
541{
542 using default_ops::eval_eq;
543 return eval_eq(a.num(), b.num()) && eval_eq(a.denom(), b.denom());
544}
545
546template <class Backend, class Arithmetic>
547inline typename std::enable_if<std::is_convertible<Arithmetic, Backend>::value&& std::is_integral<Arithmetic>::value, bool>::type
548 eval_eq(const rational_adaptor<Backend>& a, Arithmetic b)
549{
550 using default_ops::eval_eq;
551 return eval_eq(a.denom(), rational_adaptor<Backend>::one()) && eval_eq(a.num(), b);
552}
553
554template <class Backend, class Arithmetic>
555inline typename std::enable_if<std::is_convertible<Arithmetic, Backend>::value&& std::is_integral<Arithmetic>::value, bool>::type
556 eval_eq(Arithmetic b, const rational_adaptor<Backend>& a)
557{
558 using default_ops::eval_eq;
559 return eval_eq(a.denom(), rational_adaptor<Backend>::one()) && eval_eq(a.num(), b);
560}
561
562//
563// Arithmetic operations, starting with addition:
564//
565template <class Backend, class Arithmetic>
566void eval_add_subtract_imp(rational_adaptor<Backend>& result, const Arithmetic& arg, bool isaddition)
567{
568 using default_ops::eval_multiply;
569 using default_ops::eval_divide;
570 using default_ops::eval_add;
571 using default_ops::eval_gcd;
572 Backend t;
573 eval_multiply(t, result.denom(), arg);
574 if (isaddition)
575 eval_add(result.num(), t);
576 else
577 eval_subtract(result.num(), t);
578 //
579 // There is no need to re-normalize here, we have
580 // (a + bm) / b
581 // and gcd(a + bm, b) = gcd(a, b) = 1
582 //
583 /*
584 eval_gcd(t, result.num(), result.denom());
585 if (!eval_eq(t, rational_adaptor<Backend>::one()) != 0)
586 {
587 Backend t2;
588 eval_divide(t2, result.num(), t);
589 t2.swap(result.num());
590 eval_divide(t2, result.denom(), t);
591 t2.swap(result.denom());
592 }
593 */
594}
595
596template <class Backend, class Arithmetic>
597inline typename std::enable_if<std::is_convertible<Arithmetic, Backend>::value && (std::is_integral<Arithmetic>::value || std::is_same<Arithmetic, Backend>::value)>::type
598 eval_add(rational_adaptor<Backend>& result, const Arithmetic& arg)
599{
600 eval_add_subtract_imp(result, arg, true);
601}
602
603template <class Backend, class Arithmetic>
604inline typename std::enable_if<std::is_convertible<Arithmetic, Backend>::value && (std::is_integral<Arithmetic>::value || std::is_same<Arithmetic, Backend>::value)>::type
605 eval_subtract(rational_adaptor<Backend>& result, const Arithmetic& arg)
606{
607 eval_add_subtract_imp(result, arg, false);
608}
609
610template <class Backend>
611void eval_add_subtract_imp(rational_adaptor<Backend>& result, const rational_adaptor<Backend>& a, const rational_adaptor<Backend>& b, bool isaddition)
612{
613 using default_ops::eval_eq;
614 using default_ops::eval_multiply;
615 using default_ops::eval_divide;
616 using default_ops::eval_add;
617 using default_ops::eval_subtract;
618 //
619 // Let a = an/ad
620 // b = bn/bd
621 // g = gcd(ad, bd)
622 // result = rn/rd
623 //
624 // Then:
625 // rn = an * (bd/g) + bn * (ad/g)
626 // rd = ad * (bd/g)
627 // = (ad/g) * (bd/g) * g
628 //
629 // And the whole thing can then be rescaled by
630 // gcd(rn, g)
631 //
632 Backend gcd, t1, t2, t3, t4;
633 //
634 // Begin by getting the gcd of the 2 denominators:
635 //
636 eval_gcd(gcd, a.denom(), b.denom());
637 //
638 // Do we have gcd > 1:
639 //
640 if (!eval_eq(gcd, rational_adaptor<Backend>::one()))
641 {
642 //
643 // Scale the denominators by gcd, and put the results in t1 and t2:
644 //
645 eval_divide(t1, b.denom(), gcd);
646 eval_divide(t2, a.denom(), gcd);
647 //
648 // multiply the numerators by the scale denominators and put the results in t3, t4:
649 //
650 eval_multiply(t3, a.num(), t1);
651 eval_multiply(t4, b.num(), t2);
652 //
653 // Add them up:
654 //
655 if (isaddition)
656 eval_add(t3, t4);
657 else
658 eval_subtract(t3, t4);
659 //
660 // Get the gcd of gcd and our numerator (t3):
661 //
662 eval_gcd(t4, t3, gcd);
663 if (eval_eq(t4, rational_adaptor<Backend>::one()))
664 {
665 result.num() = t3;
666 eval_multiply(result.denom(), t1, a.denom());
667 }
668 else
669 {
670 //
671 // Uncommon case where gcd is not 1, divide the numerator
672 // and the denominator terms by the new gcd. Note we perform division
673 // on the existing gcd value as this is the smallest of the 3 denominator
674 // terms we'll be multiplying together, so there's a good chance it's a
675 // single limb value already:
676 //
677 eval_divide(result.num(), t3, t4);
678 eval_divide(t3, gcd, t4);
679 eval_multiply(t4, t1, t2);
680 eval_multiply(result.denom(), t4, t3);
681 }
682 }
683 else
684 {
685 //
686 // Most common case (approx 60%) where gcd is one:
687 //
688 eval_multiply(t1, a.num(), b.denom());
689 eval_multiply(t2, a.denom(), b.num());
690 if (isaddition)
691 eval_add(result.num(), t1, t2);
692 else
693 eval_subtract(result.num(), t1, t2);
694 eval_multiply(result.denom(), a.denom(), b.denom());
695 }
696}
697
698
699template <class Backend>
700inline void eval_add(rational_adaptor<Backend>& result, const rational_adaptor<Backend>& a, const rational_adaptor<Backend>& b)
701{
702 eval_add_subtract_imp(result, a, b, true);
703}
704template <class Backend>
705inline void eval_subtract(rational_adaptor<Backend>& result, const rational_adaptor<Backend>& a, const rational_adaptor<Backend>& b)
706{
707 eval_add_subtract_imp(result, a, b, false);
708}
709
710template <class Backend, class Arithmetic>
711void eval_add_subtract_imp(rational_adaptor<Backend>& result, const rational_adaptor<Backend>& a, const Arithmetic& b, bool isaddition)
712{
713 using default_ops::eval_add;
714 using default_ops::eval_subtract;
715 using default_ops::eval_multiply;
716
717 if (&result == &a)
718 return eval_add_subtract_imp(result, b, isaddition);
719
720 eval_multiply(result.num(), a.denom(), b);
721 if (isaddition)
722 eval_add(result.num(), a.num());
723 else
724 BOOST_IF_CONSTEXPR(std::numeric_limits<Backend>::is_signed == false)
725 {
726 Backend t;
727 eval_subtract(t, a.num(), result.num());
728 result.num() = std::move(t);
729 }
730 else
731 {
732 eval_subtract(result.num(), a.num());
733 result.negate();
734 }
735 result.denom() = a.denom();
736 //
737 // There is no need to re-normalize here, we have
738 // (a + bm) / b
739 // and gcd(a + bm, b) = gcd(a, b) = 1
740 //
741}
742template <class Backend, class Arithmetic>
743inline typename std::enable_if<std::is_convertible<Arithmetic, Backend>::value && (std::is_integral<Arithmetic>::value || std::is_same<Arithmetic, Backend>::value)>::type
744 eval_add(rational_adaptor<Backend>& result, const rational_adaptor<Backend>& a, const Arithmetic& b)
745{
746 eval_add_subtract_imp(result, a, b, true);
747}
748template <class Backend, class Arithmetic>
749inline typename std::enable_if<std::is_convertible<Arithmetic, Backend>::value && (std::is_integral<Arithmetic>::value || std::is_same<Arithmetic, Backend>::value)>::type
750 eval_subtract(rational_adaptor<Backend>& result, const rational_adaptor<Backend>& a, const Arithmetic& b)
751{
752 eval_add_subtract_imp(result, a, b, false);
753}
754
755//
756// Multiplication:
757//
758template <class Backend>
759void eval_multiply_imp(rational_adaptor<Backend>& result, const rational_adaptor<Backend>& a, const Backend& b_num, const Backend& b_denom)
760{
761 using default_ops::eval_multiply;
762 using default_ops::eval_divide;
763 using default_ops::eval_gcd;
764 using default_ops::eval_get_sign;
765 using default_ops::eval_eq;
766
767 Backend gcd_left, gcd_right, t1, t2;
768 eval_gcd(gcd_left, a.num(), b_denom);
769 eval_gcd(gcd_right, b_num, a.denom());
770 //
771 // Unit gcd's are the most likely case:
772 //
773 bool b_left = eval_eq(gcd_left, rational_adaptor<Backend>::one());
774 bool b_right = eval_eq(gcd_right, rational_adaptor<Backend>::one());
775
776 if (b_left && b_right)
777 {
778 eval_multiply(result.num(), a.num(), b_num);
779 eval_multiply(result.denom(), a.denom(), b_denom);
780 }
781 else if (b_left)
782 {
783 eval_divide(t2, b_num, gcd_right);
784 eval_multiply(result.num(), a.num(), t2);
785 eval_divide(t1, a.denom(), gcd_right);
786 eval_multiply(result.denom(), t1, b_denom);
787 }
788 else if (b_right)
789 {
790 eval_divide(t1, a.num(), gcd_left);
791 eval_multiply(result.num(), t1, b_num);
792 eval_divide(t2, b_denom, gcd_left);
793 eval_multiply(result.denom(), a.denom(), t2);
794 }
795 else
796 {
797 eval_divide(t1, a.num(), gcd_left);
798 eval_divide(t2, b_num, gcd_right);
799 eval_multiply(result.num(), t1, t2);
800 eval_divide(t1, a.denom(), gcd_right);
801 eval_divide(t2, b_denom, gcd_left);
802 eval_multiply(result.denom(), t1, t2);
803 }
804 //
805 // We may have b_denom negative if this is actually division, if so just correct things now:
806 //
807 if (eval_get_sign(b_denom) < 0)
808 {
809 result.num().negate();
810 result.denom().negate();
811 }
812}
813
814template <class Backend>
815void eval_multiply(rational_adaptor<Backend>& result, const rational_adaptor<Backend>& a, const rational_adaptor<Backend>& b)
816{
817 using default_ops::eval_multiply;
818
819 if (&a == &b)
820 {
821 // squaring, gcd's are 1:
822 eval_multiply(result.num(), a.num(), b.num());
823 eval_multiply(result.denom(), a.denom(), b.denom());
824 return;
825 }
826 eval_multiply_imp(result, a, b.num(), b.denom());
827}
828
829template <class Backend, class Arithmetic>
830void eval_multiply_imp(Backend& result_num, Backend& result_denom, Arithmetic arg)
831{
832 if (arg == 0)
833 {
834 result_num = rational_adaptor<Backend>::zero();
835 result_denom = rational_adaptor<Backend>::one();
836 return;
837 }
838 else if (arg == 1)
839 return;
840
841 using default_ops::eval_multiply;
842 using default_ops::eval_divide;
843 using default_ops::eval_gcd;
844 using default_ops::eval_convert_to;
845
846 Backend gcd, t;
847 Arithmetic integer_gcd;
848 eval_gcd(gcd, result_denom, arg);
849 eval_convert_to(&integer_gcd, gcd);
850 arg /= integer_gcd;
851 if (boost::multiprecision::detail::unsigned_abs(arg) > 1)
852 {
853 eval_multiply(t, result_num, arg);
854 result_num = std::move(t);
855 }
856 else if (is_minus_one(arg))
857 result_num.negate();
858 if (integer_gcd > 1)
859 {
860 eval_divide(t, result_denom, integer_gcd);
861 result_denom = std::move(t);
862 }
863}
864template <class Backend>
865void eval_multiply_imp(Backend& result_num, Backend& result_denom, Backend arg)
866{
867 using default_ops::eval_multiply;
868 using default_ops::eval_divide;
869 using default_ops::eval_gcd;
870 using default_ops::eval_convert_to;
871 using default_ops::eval_is_zero;
872 using default_ops::eval_eq;
873 using default_ops::eval_get_sign;
874
875 if (eval_is_zero(arg))
876 {
877 result_num = rational_adaptor<Backend>::zero();
878 result_denom = rational_adaptor<Backend>::one();
879 return;
880 }
881 else if (eval_eq(arg, rational_adaptor<Backend>::one()))
882 return;
883
884 Backend gcd, t;
885 eval_gcd(gcd, result_denom, arg);
886 if (!eval_eq(gcd, rational_adaptor<Backend>::one()))
887 {
888 eval_divide(t, arg, gcd);
889 arg = t;
890 }
891 else
892 t = arg;
893 if (eval_get_sign(arg) < 0)
894 t.negate();
895
896 if (!eval_eq(t, rational_adaptor<Backend>::one()))
897 {
898 eval_multiply(t, result_num, arg);
899 result_num = std::move(t);
900 }
901 else if (eval_get_sign(arg) < 0)
902 result_num.negate();
903 if (!eval_eq(gcd, rational_adaptor<Backend>::one()))
904 {
905 eval_divide(t, result_denom, gcd);
906 result_denom = std::move(t);
907 }
908}
909
910template <class Backend, class Arithmetic>
911inline typename std::enable_if<std::is_convertible<Arithmetic, Backend>::value && (std::is_integral<Arithmetic>::value || std::is_same<Arithmetic, Backend>::value)>::type
912 eval_multiply(rational_adaptor<Backend>& result, const Arithmetic& arg)
913{
914 eval_multiply_imp(result.num(), result.denom(), arg);
915}
916
917template <class Backend, class Arithmetic>
918typename std::enable_if<std::is_convertible<Arithmetic, Backend>::value && std::is_integral<Arithmetic>::value>::type
919 eval_multiply_imp(rational_adaptor<Backend>& result, const Backend& a_num, const Backend& a_denom, Arithmetic b)
920{
921 if (b == 0)
922 {
923 result.num() = rational_adaptor<Backend>::zero();
924 result.denom() = rational_adaptor<Backend>::one();
925 return;
926 }
927 else if (b == 1)
928 {
929 result.num() = a_num;
930 result.denom() = a_denom;
931 return;
932 }
933
934 using default_ops::eval_multiply;
935 using default_ops::eval_divide;
936 using default_ops::eval_gcd;
937 using default_ops::eval_convert_to;
938
939 Backend gcd;
940 Arithmetic integer_gcd;
941 eval_gcd(gcd, a_denom, b);
942 eval_convert_to(&integer_gcd, gcd);
943 b /= integer_gcd;
944 if (boost::multiprecision::detail::unsigned_abs(b) > 1)
945 eval_multiply(result.num(), a_num, b);
946 else if (is_minus_one(b))
947 {
948 result.num() = a_num;
949 result.num().negate();
950 }
951 else
952 result.num() = a_num;
953 if (integer_gcd > 1)
954 eval_divide(result.denom(), a_denom, integer_gcd);
955 else
956 result.denom() = a_denom;
957}
958template <class Backend>
959inline void eval_multiply_imp(rational_adaptor<Backend>& result, const Backend& a_num, const Backend& a_denom, const Backend& b)
960{
961 result.num() = a_num;
962 result.denom() = a_denom;
963 eval_multiply_imp(result.num(), result.denom(), b);
964}
965
966template <class Backend, class Arithmetic>
967inline typename std::enable_if<std::is_convertible<Arithmetic, Backend>::value && (std::is_integral<Arithmetic>::value || std::is_same<Arithmetic, Backend>::value)>::type
968 eval_multiply(rational_adaptor<Backend>& result, const rational_adaptor<Backend>& a, const Arithmetic& b)
969{
970 if (&result == &a)
971 return eval_multiply(result, b);
972
973 eval_multiply_imp(result, a.num(), a.denom(), b);
974}
975
976template <class Backend, class Arithmetic>
977inline typename std::enable_if<std::is_convertible<Arithmetic, Backend>::value && (std::is_integral<Arithmetic>::value || std::is_same<Arithmetic, Backend>::value)>::type
978 eval_multiply(rational_adaptor<Backend>& result, const Arithmetic& b, const rational_adaptor<Backend>& a)
979{
980 return eval_multiply(result, a, b);
981}
982
983//
984// Division:
985//
986template <class Backend>
987inline void eval_divide(rational_adaptor<Backend>& result, const rational_adaptor<Backend>& a, const rational_adaptor<Backend>& b)
988{
989 using default_ops::eval_multiply;
990 using default_ops::eval_get_sign;
991
992 if (eval_get_sign(b.num()) == 0)
993 {
994 BOOST_MP_THROW_EXCEPTION(std::overflow_error("Integer division by zero"));
995 return;
996 }
997 if (&a == &b)
998 {
999 // Huh? Really?
1000 result.num() = result.denom() = rational_adaptor<Backend>::one();
1001 return;
1002 }
1003 if (&result == &b)
1004 {
1005 rational_adaptor<Backend> t(b);
1006 return eval_divide(result, a, t);
1007 }
1008 eval_multiply_imp(result, a, b.denom(), b.num());
1009}
1010
1011template <class Backend, class Arithmetic>
1012inline typename std::enable_if<std::is_convertible<Arithmetic, Backend>::value && (std::is_integral<Arithmetic>::value || std::is_same<Arithmetic, Backend>::value)>::type
1013 eval_divide(rational_adaptor<Backend>& result, const Arithmetic& b, const rational_adaptor<Backend>& a)
1014{
1015 using default_ops::eval_get_sign;
1016
1017 if (eval_get_sign(a.num()) == 0)
1018 {
1019 BOOST_MP_THROW_EXCEPTION(std::overflow_error("Integer division by zero"));
1020 return;
1021 }
1022 if (&a == &result)
1023 {
1024 eval_multiply_imp(result.denom(), result.num(), b);
1025 result.num().swap(result.denom());
1026 }
1027 else
1028 eval_multiply_imp(result, a.denom(), a.num(), b);
1029
1030 if (eval_get_sign(result.denom()) < 0)
1031 {
1032 result.num().negate();
1033 result.denom().negate();
1034 }
1035}
1036
1037template <class Backend, class Arithmetic>
1038typename std::enable_if<std::is_convertible<Arithmetic, Backend>::value && std::is_integral<Arithmetic>::value>::type
1039eval_divide(rational_adaptor<Backend>& result, Arithmetic arg)
1040{
1041 if (arg == 0)
1042 {
1043 BOOST_MP_THROW_EXCEPTION(std::overflow_error("Integer division by zero"));
1044 return;
1045 }
1046 else if (arg == 1)
1047 return;
1048 else if (is_minus_one(arg))
1049 {
1050 result.negate();
1051 return;
1052 }
1053 if (eval_get_sign(result) == 0)
1054 {
1055 return;
1056 }
1057
1058
1059 using default_ops::eval_multiply;
1060 using default_ops::eval_gcd;
1061 using default_ops::eval_convert_to;
1062 using default_ops::eval_divide;
1063
1064 Backend gcd, t;
1065 Arithmetic integer_gcd;
1066 eval_gcd(gcd, result.num(), arg);
1067 eval_convert_to(&integer_gcd, gcd);
1068 arg /= integer_gcd;
1069
1070 eval_multiply(t, result.denom(), boost::multiprecision::detail::unsigned_abs(arg));
1071 result.denom() = std::move(t);
1072 if (arg < 0)
1073 {
1074 result.num().negate();
1075 }
1076 if (integer_gcd > 1)
1077 {
1078 eval_divide(t, result.num(), integer_gcd);
1079 result.num() = std::move(t);
1080 }
1081}
1082template <class Backend>
1083void eval_divide(rational_adaptor<Backend>& result, const rational_adaptor<Backend>& a, Backend arg)
1084{
1085 using default_ops::eval_multiply;
1086 using default_ops::eval_gcd;
1087 using default_ops::eval_convert_to;
1088 using default_ops::eval_divide;
1089 using default_ops::eval_is_zero;
1090 using default_ops::eval_eq;
1091 using default_ops::eval_get_sign;
1092
1093 if (eval_is_zero(arg))
1094 {
1095 BOOST_MP_THROW_EXCEPTION(std::overflow_error("Integer division by zero"));
1096 return;
1097 }
1098 else if (eval_eq(a, rational_adaptor<Backend>::one()) || (eval_get_sign(a) == 0))
1099 {
1100 if (&result != &a)
1101 result = a;
1102 result.denom() = arg;
1103 return;
1104 }
1105
1106 Backend gcd, u_arg, t;
1107 eval_gcd(gcd, a.num(), arg);
1108 bool has_unit_gcd = eval_eq(gcd, rational_adaptor<Backend>::one());
1109 if (!has_unit_gcd)
1110 {
1111 eval_divide(u_arg, arg, gcd);
1112 arg = u_arg;
1113 }
1114 else
1115 u_arg = arg;
1116 if (eval_get_sign(u_arg) < 0)
1117 u_arg.negate();
1118
1119 eval_multiply(t, a.denom(), u_arg);
1120 result.denom() = std::move(t);
1121
1122 if (!has_unit_gcd)
1123 {
1124 eval_divide(t, a.num(), gcd);
1125 result.num() = std::move(t);
1126 }
1127 else if (&result != &a)
1128 result.num() = a.num();
1129
1130 if (eval_get_sign(arg) < 0)
1131 {
1132 result.num().negate();
1133 }
1134}
1135template <class Backend>
1136void eval_divide(rational_adaptor<Backend>& result, const Backend& arg)
1137{
1138 eval_divide(result, result, arg);
1139}
1140
1141template <class Backend, class Arithmetic>
1142typename std::enable_if<std::is_convertible<Arithmetic, Backend>::value && std::is_integral<Arithmetic>::value>::type
1143 eval_divide(rational_adaptor<Backend>& result, const rational_adaptor<Backend>& a, Arithmetic arg)
1144{
1145 if (&result == &a)
1146 return eval_divide(result, arg);
1147 if (arg == 0)
1148 {
1149 BOOST_MP_THROW_EXCEPTION(std::overflow_error("Integer division by zero"));
1150 return;
1151 }
1152 else if (arg == 1)
1153 {
1154 result = a;
1155 return;
1156 }
1157 else if (is_minus_one(arg))
1158 {
1159 result = a;
1160 result.num().negate();
1161 return;
1162 }
1163
1164 if (eval_get_sign(a) == 0)
1165 {
1166 result = a;
1167 return;
1168 }
1169
1170 using default_ops::eval_multiply;
1171 using default_ops::eval_divide;
1172 using default_ops::eval_gcd;
1173 using default_ops::eval_convert_to;
1174
1175 Backend gcd;
1176 Arithmetic integer_gcd;
1177 eval_gcd(gcd, a.num(), arg);
1178 eval_convert_to(&integer_gcd, gcd);
1179 arg /= integer_gcd;
1180 eval_multiply(result.denom(), a.denom(), boost::multiprecision::detail::unsigned_abs(arg));
1181
1182 if (integer_gcd > 1)
1183 {
1184 eval_divide(result.num(), a.num(), integer_gcd);
1185 }
1186 else
1187 result.num() = a.num();
1188 if (arg < 0)
1189 {
1190 result.num().negate();
1191 }
1192}
1193
1194//
1195// Increment and decrement:
1196//
1197template <class Backend>
1198inline void eval_increment(rational_adaptor<Backend>& arg)
1199{
1200 using default_ops::eval_add;
1201 eval_add(arg.num(), arg.denom());
1202}
1203template <class Backend>
1204inline void eval_decrement(rational_adaptor<Backend>& arg)
1205{
1206 using default_ops::eval_subtract;
1207 eval_subtract(arg.num(), arg.denom());
1208}
1209
1210//
1211// abs:
1212//
1213template <class Backend>
1214inline void eval_abs(rational_adaptor<Backend>& result, const rational_adaptor<Backend>& arg)
1215{
1216 using default_ops::eval_abs;
1217 eval_abs(result.num(), arg.num());
1218 result.denom() = arg.denom();
1219}
1220
1221} // namespace backends
1222
1223//
1224// Define a category for this number type, one of:
1225//
1226// number_kind_integer
1227// number_kind_floating_point
1228// number_kind_rational
1229// number_kind_fixed_point
1230// number_kind_complex
1231//
1232template<class Backend>
1233struct number_category<rational_adaptor<Backend> > : public std::integral_constant<int, number_kind_rational>
1234{};
1235
1236template <class Backend, expression_template_option ExpressionTemplates>
1237struct component_type<number<rational_adaptor<Backend>, ExpressionTemplates> >
1238{
1239 typedef number<Backend, ExpressionTemplates> type;
1240};
1241
1242template <class IntBackend, expression_template_option ET>
1243inline number<IntBackend, ET> numerator(const number<rational_adaptor<IntBackend>, ET>& val)
1244{
1245 return val.backend().num();
1246}
1247template <class IntBackend, expression_template_option ET>
1248inline number<IntBackend, ET> denominator(const number<rational_adaptor<IntBackend>, ET>& val)
1249{
1250 return val.backend().denom();
1251}
1252
1253template <class Backend>
1254struct is_unsigned_number<rational_adaptor<Backend> > : public is_unsigned_number<Backend>
1255{};
1256
1257
1258}} // namespace boost::multiprecision
1259
1260namespace std {
1261
1262 template <class IntBackend, boost::multiprecision::expression_template_option ExpressionTemplates>
1263 class numeric_limits<boost::multiprecision::number<boost::multiprecision::rational_adaptor<IntBackend>, ExpressionTemplates> > : public std::numeric_limits<boost::multiprecision::number<IntBackend, ExpressionTemplates> >
1264 {
1265 using base_type = std::numeric_limits<boost::multiprecision::number<IntBackend> >;
1266 using number_type = boost::multiprecision::number<boost::multiprecision::rational_adaptor<IntBackend> >;
1267
1268 public:
1269 static constexpr bool is_integer = false;
1270 static constexpr bool is_exact = true;
1271 static constexpr number_type(min)() { return (base_type::min)(); }
1272 static constexpr number_type(max)() { return (base_type::max)(); }
1273 static constexpr number_type lowest() { return -(max)(); }
1274 static constexpr number_type epsilon() { return base_type::epsilon(); }
1275 static constexpr number_type round_error() { return epsilon() / 2; }
1276 static constexpr number_type infinity() { return base_type::infinity(); }
1277 static constexpr number_type quiet_NaN() { return base_type::quiet_NaN(); }
1278 static constexpr number_type signaling_NaN() { return base_type::signaling_NaN(); }
1279 static constexpr number_type denorm_min() { return base_type::denorm_min(); }
1280 };
1281
1282 template <class IntBackend, boost::multiprecision::expression_template_option ExpressionTemplates>
1283 constexpr bool numeric_limits<boost::multiprecision::number<boost::multiprecision::rational_adaptor<IntBackend>, ExpressionTemplates> >::is_integer;
1284 template <class IntBackend, boost::multiprecision::expression_template_option ExpressionTemplates>
1285 constexpr bool numeric_limits<boost::multiprecision::number<boost::multiprecision::rational_adaptor<IntBackend>, ExpressionTemplates> >::is_exact;
1286
1287} // namespace std
1288
1289#endif
1290

source code of boost/libs/multiprecision/include/boost/multiprecision/rational_adaptor.hpp