1 // Boost rational.hpp header file ------------------------------------------//
3 // (C) Copyright Paul Moore 1999. Permission to copy, use, modify, sell and
4 // distribute this software is granted provided this copyright notice appears
5 // in all copies. This software is provided "as is" without express or
6 // implied warranty, and with no claim as to its suitability for any purpose.
8 // boostinspect:nolicense (don't complain about the lack of a Boost license)
9 // (Paul Moore hasn't been in contact for years, so there's no way to change the
12 // See http://www.boost.org/libs/rational for documentation.
15 // Thanks to the boost mailing list in general for useful comments.
16 // Particular contributions included:
17 // Andrew D Jewell, for reminding me to take care to avoid overflow
18 // Ed Brey, for many comments, including picking up on some dreadful typos
19 // Stephen Silver contributed the test suite and comments on user-defined
21 // Nickolay Mladenov, for the implementation of operator+=
24 // 05 Nov 06 Change rational_cast to not depend on division between different
25 // types (Daryle Walker)
26 // 04 Nov 06 Off-load GCD and LCM to Boost.Math; add some invariant checks;
27 // add std::numeric_limits<> requirement to help GCD (Daryle Walker)
28 // 31 Oct 06 Recoded both operator< to use round-to-negative-infinity
29 // divisions; the rational-value version now uses continued fraction
30 // expansion to avoid overflows, for bug #798357 (Daryle Walker)
31 // 20 Oct 06 Fix operator bool_type for CW 8.3 (Joaquín M López Muñoz)
32 // 18 Oct 06 Use EXPLICIT_TEMPLATE_TYPE helper macros from Boost.Config
33 // (Joaquín M López Muñoz)
34 // 27 Dec 05 Add Boolean conversion operator (Daryle Walker)
35 // 28 Sep 02 Use _left versions of operators from operators.hpp
36 // 05 Jul 01 Recode gcd(), avoiding std::swap (Helmut Zeisel)
37 // 03 Mar 01 Workarounds for Intel C++ 5.0 (David Abrahams)
38 // 05 Feb 01 Update operator>> to tighten up input syntax
39 // 05 Feb 01 Final tidy up of gcd code prior to the new release
40 // 27 Jan 01 Recode abs() without relying on abs(IntType)
41 // 21 Jan 01 Include Nickolay Mladenov's operator+= algorithm,
42 // tidy up a number of areas, use newer features of operators.hpp
43 // (reduces space overhead to zero), add operator!,
44 // introduce explicit mixed-mode arithmetic operations
45 // 12 Jan 01 Include fixes to handle a user-defined IntType better
46 // 19 Nov 00 Throw on divide by zero in operator /= (John (EBo) David)
47 // 23 Jun 00 Incorporate changes from Mark Rodgers for Borland C++
48 // 22 Jun 00 Change _MSC_VER to BOOST_MSVC so other compilers are not
49 // affected (Beman Dawes)
50 // 6 Mar 00 Fix operator-= normalization, #include <string> (Jens Maurer)
51 // 14 Dec 99 Modifications based on comments from the boost list
52 // 09 Dec 99 Initial Version (Paul Moore)
54 #ifndef BOOST_RATIONAL_HPP
55 #define BOOST_RATIONAL_HPP
57 #include <iostream> // for std::istream and std::ostream
58 #include <ios> // for std::noskipws
59 #include <stdexcept> // for std::domain_error
60 #include <string> // for std::string implicit constructor
61 #include <boost/operators.hpp> // for boost::addable etc
62 #include <cstdlib> // for std::abs
63 #include <boost/call_traits.hpp> // for boost::call_traits
64 #include <boost/config.hpp> // for BOOST_NO_STDC_NAMESPACE, BOOST_MSVC
65 #include <boost/detail/workaround.hpp> // for BOOST_WORKAROUND
66 #include <boost/assert.hpp> // for BOOST_ASSERT
67 #include <boost/math/common_factor_rt.hpp> // for boost::math::gcd, lcm
68 #include <limits> // for std::numeric_limits
69 #include <boost/static_assert.hpp> // for BOOST_STATIC_ASSERT
71 // Control whether depreciated GCD and LCM functions are included (default: yes)
72 #ifndef BOOST_CONTROL_RATIONAL_HAS_GCD
73 #define BOOST_CONTROL_RATIONAL_HAS_GCD 1
78 #if BOOST_CONTROL_RATIONAL_HAS_GCD
79 template <typename IntType
>
80 IntType
gcd(IntType n
, IntType m
)
82 // Defer to the version in Boost.Math
83 return math::gcd( n
, m
);
86 template <typename IntType
>
87 IntType
lcm(IntType n
, IntType m
)
89 // Defer to the version in Boost.Math
90 return math::lcm( n
, m
);
92 #endif // BOOST_CONTROL_RATIONAL_HAS_GCD
94 class bad_rational
: public std::domain_error
97 explicit bad_rational() : std::domain_error("bad rational: zero denominator") {}
100 template <typename IntType
>
103 template <typename IntType
>
104 rational
<IntType
> abs(const rational
<IntType
>& r
);
106 template <typename IntType
>
108 less_than_comparable
< rational
<IntType
>,
109 equality_comparable
< rational
<IntType
>,
110 less_than_comparable2
< rational
<IntType
>, IntType
,
111 equality_comparable2
< rational
<IntType
>, IntType
,
112 addable
< rational
<IntType
>,
113 subtractable
< rational
<IntType
>,
114 multipliable
< rational
<IntType
>,
115 dividable
< rational
<IntType
>,
116 addable2
< rational
<IntType
>, IntType
,
117 subtractable2
< rational
<IntType
>, IntType
,
118 subtractable2_left
< rational
<IntType
>, IntType
,
119 multipliable2
< rational
<IntType
>, IntType
,
120 dividable2
< rational
<IntType
>, IntType
,
121 dividable2_left
< rational
<IntType
>, IntType
,
122 incrementable
< rational
<IntType
>,
123 decrementable
< rational
<IntType
>
124 > > > > > > > > > > > > > > > >
126 // Class-wide pre-conditions
127 BOOST_STATIC_ASSERT( ::std::numeric_limits
<IntType
>::is_specialized
);
130 typedef typename
boost::call_traits
<IntType
>::param_type param_type
;
132 struct helper
{ IntType parts
[2]; };
133 typedef IntType (helper::* bool_type
)[2];
136 typedef IntType int_type
;
137 rational() : num(0), den(1) {}
138 rational(param_type n
) : num(n
), den(1) {}
139 rational(param_type n
, param_type d
) : num(n
), den(d
) { normalize(); }
141 // Default copy constructor and assignment are fine
143 // Add assignment from IntType
144 rational
& operator=(param_type n
) { return assign(n
, 1); }
147 rational
& assign(param_type n
, param_type d
);
149 // Access to representation
150 IntType
numerator() const { return num
; }
151 IntType
denominator() const { return den
; }
153 // Arithmetic assignment operators
154 rational
& operator+= (const rational
& r
);
155 rational
& operator-= (const rational
& r
);
156 rational
& operator*= (const rational
& r
);
157 rational
& operator/= (const rational
& r
);
159 rational
& operator+= (param_type i
);
160 rational
& operator-= (param_type i
);
161 rational
& operator*= (param_type i
);
162 rational
& operator/= (param_type i
);
164 // Increment and decrement
165 const rational
& operator++();
166 const rational
& operator--();
169 bool operator!() const { return !num
; }
171 // Boolean conversion
173 #if BOOST_WORKAROUND(__MWERKS__,<=0x3003)
174 // The "ISO C++ Template Parser" option in CW 8.3 chokes on the
175 // following, hence we selectively disable that option for the
177 #pragma parse_mfunc_templ off
180 operator bool_type() const { return operator !() ? 0 : &helper::parts
; }
182 #if BOOST_WORKAROUND(__MWERKS__,<=0x3003)
183 #pragma parse_mfunc_templ reset
186 // Comparison operators
187 bool operator< (const rational
& r
) const;
188 bool operator== (const rational
& r
) const;
190 bool operator< (param_type i
) const;
191 bool operator> (param_type i
) const;
192 bool operator== (param_type i
) const;
195 // Implementation - numerator and denominator (normalized).
196 // Other possibilities - separate whole-part, or sign, fields?
200 // Representation note: Fractions are kept in normalized form at all
201 // times. normalized form is defined as gcd(num,den) == 1 and den > 0.
202 // In particular, note that the implementation of abs() below relies
203 // on den always being positive.
204 bool test_invariant() const;
209 template <typename IntType
>
210 inline rational
<IntType
>& rational
<IntType
>::assign(param_type n
, param_type d
)
218 // Unary plus and minus
219 template <typename IntType
>
220 inline rational
<IntType
> operator+ (const rational
<IntType
>& r
)
225 template <typename IntType
>
226 inline rational
<IntType
> operator- (const rational
<IntType
>& r
)
228 return rational
<IntType
>(-r
.numerator(), r
.denominator());
231 // Arithmetic assignment operators
232 template <typename IntType
>
233 rational
<IntType
>& rational
<IntType
>::operator+= (const rational
<IntType
>& r
)
235 // This calculation avoids overflow, and minimises the number of expensive
236 // calculations. Thanks to Nickolay Mladenov for this algorithm.
239 // We have to compute a/b + c/d, where gcd(a,b)=1 and gcd(b,c)=1.
240 // Let g = gcd(b,d), and b = b1*g, d=d1*g. Then gcd(b1,d1)=1
242 // The result is (a*d1 + c*b1) / (b1*d1*g).
243 // Now we have to normalize this ratio.
244 // Let's assume h | gcd((a*d1 + c*b1), (b1*d1*g)), and h > 1
245 // If h | b1 then gcd(h,d1)=1 and hence h|(a*d1+c*b1) => h|a.
246 // But since gcd(a,b1)=1 we have h=1.
247 // Similarly h|d1 leads to h=1.
248 // So we have that h | gcd((a*d1 + c*b1) , (b1*d1*g)) => h|g
249 // Finally we have gcd((a*d1 + c*b1), (b1*d1*g)) = gcd((a*d1 + c*b1), g)
250 // Which proves that instead of normalizing the result, it is better to
251 // divide num and den by gcd((a*d1 + c*b1), g)
253 // Protect against self-modification
254 IntType r_num
= r
.num
;
255 IntType r_den
= r
.den
;
257 IntType g
= math::gcd(den
, r_den
);
258 den
/= g
; // = b1 from the calculations above
259 num
= num
* (r_den
/ g
) + r_num
* den
;
260 g
= math::gcd(num
, g
);
267 template <typename IntType
>
268 rational
<IntType
>& rational
<IntType
>::operator-= (const rational
<IntType
>& r
)
270 // Protect against self-modification
271 IntType r_num
= r
.num
;
272 IntType r_den
= r
.den
;
274 // This calculation avoids overflow, and minimises the number of expensive
275 // calculations. It corresponds exactly to the += case above
276 IntType g
= math::gcd(den
, r_den
);
278 num
= num
* (r_den
/ g
) - r_num
* den
;
279 g
= math::gcd(num
, g
);
286 template <typename IntType
>
287 rational
<IntType
>& rational
<IntType
>::operator*= (const rational
<IntType
>& r
)
289 // Protect against self-modification
290 IntType r_num
= r
.num
;
291 IntType r_den
= r
.den
;
293 // Avoid overflow and preserve normalization
294 IntType gcd1
= math::gcd(num
, r_den
);
295 IntType gcd2
= math::gcd(r_num
, den
);
296 num
= (num
/gcd1
) * (r_num
/gcd2
);
297 den
= (den
/gcd2
) * (r_den
/gcd1
);
301 template <typename IntType
>
302 rational
<IntType
>& rational
<IntType
>::operator/= (const rational
<IntType
>& r
)
304 // Protect against self-modification
305 IntType r_num
= r
.num
;
306 IntType r_den
= r
.den
;
308 // Avoid repeated construction
311 // Trap division by zero
313 throw bad_rational();
317 // Avoid overflow and preserve normalization
318 IntType gcd1
= math::gcd(num
, r_num
);
319 IntType gcd2
= math::gcd(r_den
, den
);
320 num
= (num
/gcd1
) * (r_den
/gcd2
);
321 den
= (den
/gcd2
) * (r_num
/gcd1
);
330 // Mixed-mode operators
331 template <typename IntType
>
332 inline rational
<IntType
>&
333 rational
<IntType
>::operator+= (param_type i
)
335 return operator+= (rational
<IntType
>(i
));
338 template <typename IntType
>
339 inline rational
<IntType
>&
340 rational
<IntType
>::operator-= (param_type i
)
342 return operator-= (rational
<IntType
>(i
));
345 template <typename IntType
>
346 inline rational
<IntType
>&
347 rational
<IntType
>::operator*= (param_type i
)
349 return operator*= (rational
<IntType
>(i
));
352 template <typename IntType
>
353 inline rational
<IntType
>&
354 rational
<IntType
>::operator/= (param_type i
)
356 return operator/= (rational
<IntType
>(i
));
359 // Increment and decrement
360 template <typename IntType
>
361 inline const rational
<IntType
>& rational
<IntType
>::operator++()
363 // This can never denormalise the fraction
368 template <typename IntType
>
369 inline const rational
<IntType
>& rational
<IntType
>::operator--()
371 // This can never denormalise the fraction
376 // Comparison operators
377 template <typename IntType
>
378 bool rational
<IntType
>::operator< (const rational
<IntType
>& r
) const
380 // Avoid repeated construction
381 int_type
const zero( 0 );
383 // This should really be a class-wide invariant. The reason for these
384 // checks is that for 2's complement systems, INT_MIN has no corresponding
385 // positive, so negating it during normalization keeps it INT_MIN, which
386 // is bad for later calculations that assume a positive denominator.
387 BOOST_ASSERT( this->den
> zero
);
388 BOOST_ASSERT( r
.den
> zero
);
390 // Determine relative order by expanding each value to its simple continued
391 // fraction representation using the Euclidian GCD algorithm.
392 struct { int_type n
, d
, q
, r
; } ts
= { this->num
, this->den
, this->num
/
393 this->den
, this->num
% this->den
}, rs
= { r
.num
, r
.den
, r
.num
/ r
.den
,
395 unsigned reverse
= 0u;
397 // Normalize negative moduli by repeatedly adding the (positive) denominator
398 // and decrementing the quotient. Later cycles should have all positive
399 // values, so this only has to be done for the first cycle. (The rules of
400 // C++ require a nonnegative quotient & remainder for a nonnegative dividend
401 // & positive divisor.)
402 while ( ts
.r
< zero
) { ts
.r
+= ts
.d
; --ts
.q
; }
403 while ( rs
.r
< zero
) { rs
.r
+= rs
.d
; --rs
.q
; }
405 // Loop through and compare each variable's continued-fraction components
408 // The quotients of the current cycle are the continued-fraction
409 // components. Comparing two c.f. is comparing their sequences,
410 // stopping at the first difference.
413 // Since reciprocation changes the relative order of two variables,
414 // and c.f. use reciprocals, the less/greater-than test reverses
415 // after each index. (Start w/ non-reversed @ whole-number place.)
416 return reverse
? ts
.q
> rs
.q
: ts
.q
< rs
.q
;
419 // Prepare the next cycle
422 if ( (ts
.r
== zero
) || (rs
.r
== zero
) )
424 // At least one variable's c.f. expansion has ended
428 ts
.n
= ts
.d
; ts
.d
= ts
.r
;
429 ts
.q
= ts
.n
/ ts
.d
; ts
.r
= ts
.n
% ts
.d
;
430 rs
.n
= rs
.d
; rs
.d
= rs
.r
;
431 rs
.q
= rs
.n
/ rs
.d
; rs
.r
= rs
.n
% rs
.d
;
434 // Compare infinity-valued components for otherwise equal sequences
437 // Both remainders are zero, so the next (and subsequent) c.f.
438 // components for both sequences are infinity. Therefore, the sequences
439 // and their corresponding values are equal.
445 #pragma warning(push)
446 #pragma warning(disable:4800)
448 // Exactly one of the remainders is zero, so all following c.f.
449 // components of that variable are infinity, while the other variable
450 // has a finite next c.f. component. So that other variable has the
451 // lesser value (modulo the reversal flag!).
452 return ( ts
.r
!= zero
) != static_cast<bool>( reverse
);
459 template <typename IntType
>
460 bool rational
<IntType
>::operator< (param_type i
) const
462 // Avoid repeated construction
463 int_type
const zero( 0 );
465 // Break value into mixed-fraction form, w/ always-nonnegative remainder
466 BOOST_ASSERT( this->den
> zero
);
467 int_type q
= this->num
/ this->den
, r
= this->num
% this->den
;
468 while ( r
< zero
) { r
+= this->den
; --q
; }
470 // Compare with just the quotient, since the remainder always bumps the
471 // value up. [Since q = floor(n/d), and if n/d < i then q < i, if n/d == i
472 // then q == i, if n/d == i + r/d then q == i, and if n/d >= i + 1 then
473 // q >= i + 1 > i; therefore n/d < i iff q < i.]
477 template <typename IntType
>
478 bool rational
<IntType
>::operator> (param_type i
) const
480 // Trap equality first
481 if (num
== i
&& den
== IntType(1))
484 // Otherwise, we can use operator<
485 return !operator<(i
);
488 template <typename IntType
>
489 inline bool rational
<IntType
>::operator== (const rational
<IntType
>& r
) const
491 return ((num
== r
.num
) && (den
== r
.den
));
494 template <typename IntType
>
495 inline bool rational
<IntType
>::operator== (param_type i
) const
497 return ((den
== IntType(1)) && (num
== i
));
501 template <typename IntType
>
502 inline bool rational
<IntType
>::test_invariant() const
504 return ( this->den
> int_type(0) ) && ( math::gcd(this->num
, this->den
) ==
509 template <typename IntType
>
510 void rational
<IntType
>::normalize()
512 // Avoid repeated construction
516 throw bad_rational();
518 // Handle the case of zero separately, to avoid division by zero
524 IntType g
= math::gcd(num
, den
);
529 // Ensure that the denominator is positive
535 BOOST_ASSERT( this->test_invariant() );
540 // A utility class to reset the format flags for an istream at end
541 // of scope, even in case of exceptions
543 resetter(std::istream
& is
) : is_(is
), f_(is
.flags()) {}
544 ~resetter() { is_
.flags(f_
); }
546 std::istream::fmtflags f_
; // old GNU c++ lib has no ios_base
552 template <typename IntType
>
553 std::istream
& operator>> (std::istream
& is
, rational
<IntType
>& r
)
555 IntType n
= IntType(0), d
= IntType(1);
557 detail::resetter
sentry(is
);
563 is
.clear(std::istream::badbit
); // old GNU c++ lib has no ios_base
565 #if !defined(__GNUC__) || (defined(__GNUC__) && (__GNUC__ >= 3)) || defined __SGI_STL_PORT
568 is
.unsetf(ios::skipws
); // compiles, but seems to have no effect.
578 // Add manipulators for output format?
579 template <typename IntType
>
580 std::ostream
& operator<< (std::ostream
& os
, const rational
<IntType
>& r
)
582 os
<< r
.numerator() << '/' << r
.denominator();
587 template <typename T
, typename IntType
>
588 inline T
rational_cast(
589 const rational
<IntType
>& src
BOOST_APPEND_EXPLICIT_TEMPLATE_TYPE(T
))
591 return static_cast<T
>(src
.numerator())/static_cast<T
>(src
.denominator());
594 // Do not use any abs() defined on IntType - it isn't worth it, given the
595 // difficulties involved (Koenig lookup required, there may not *be* an abs()
596 // defined, etc etc).
597 template <typename IntType
>
598 inline rational
<IntType
> abs(const rational
<IntType
>& r
)
600 if (r
.numerator() >= IntType(0))
603 return rational
<IntType
>(-r
.numerator(), r
.denominator());
608 #endif // BOOST_RATIONAL_HPP