// in all copies. This software is provided "as is" without express or
// implied warranty, and with no claim as to its suitability for any purpose.
+// boostinspect:nolicense (don't complain about the lack of a Boost license)
+// (Paul Moore hasn't been in contact for years, so there's no way to change the
+// license.)
+
// See http://www.boost.org/libs/rational for documentation.
// Credits:
// Nickolay Mladenov, for the implementation of operator+=
// Revision History
+// 05 Nov 06 Change rational_cast to not depend on division between different
+// types (Daryle Walker)
+// 04 Nov 06 Off-load GCD and LCM to Boost.Math; add some invariant checks;
+// add std::numeric_limits<> requirement to help GCD (Daryle Walker)
+// 31 Oct 06 Recoded both operator< to use round-to-negative-infinity
+// divisions; the rational-value version now uses continued fraction
+// expansion to avoid overflows, for bug #798357 (Daryle Walker)
+// 20 Oct 06 Fix operator bool_type for CW 8.3 (Joaquín M López Muñoz)
+// 18 Oct 06 Use EXPLICIT_TEMPLATE_TYPE helper macros from Boost.Config
+// (Joaquín M López Muñoz)
// 27 Dec 05 Add Boolean conversion operator (Daryle Walker)
// 28 Sep 02 Use _left versions of operators from operators.hpp
// 05 Jul 01 Recode gcd(), avoiding std::swap (Helmut Zeisel)
#define BOOST_RATIONAL_HPP
#include <iostream> // for std::istream and std::ostream
-#include <iomanip> // for std::noskipws
+#include <ios> // for std::noskipws
#include <stdexcept> // for std::domain_error
#include <string> // for std::string implicit constructor
#include <boost/operators.hpp> // for boost::addable etc
#include <cstdlib> // for std::abs
#include <boost/call_traits.hpp> // for boost::call_traits
#include <boost/config.hpp> // for BOOST_NO_STDC_NAMESPACE, BOOST_MSVC
+#include <boost/detail/workaround.hpp> // for BOOST_WORKAROUND
+#include <boost/assert.hpp> // for BOOST_ASSERT
+#include <boost/math/common_factor_rt.hpp> // for boost::math::gcd, lcm
+#include <limits> // for std::numeric_limits
+#include <boost/static_assert.hpp> // for BOOST_STATIC_ASSERT
+
+// Control whether depreciated GCD and LCM functions are included (default: yes)
+#ifndef BOOST_CONTROL_RATIONAL_HAS_GCD
+#define BOOST_CONTROL_RATIONAL_HAS_GCD 1
+#endif
namespace boost {
-// Note: We use n and m as temporaries in this function, so there is no value
-// in using const IntType& as we would only need to make a copy anyway...
+#if BOOST_CONTROL_RATIONAL_HAS_GCD
template <typename IntType>
IntType gcd(IntType n, IntType m)
{
- // Avoid repeated construction
- IntType zero(0);
-
- // This is abs() - given the existence of broken compilers with Koenig
- // lookup issues and other problems, I code this explicitly. (Remember,
- // IntType may be a user-defined type).
- if (n < zero)
- n = -n;
- if (m < zero)
- m = -m;
-
- // As n and m are now positive, we can be sure that %= returns a
- // positive value (the standard guarantees this for built-in types,
- // and we require it of user-defined types).
- for(;;) {
- if(m == zero)
- return n;
- n %= m;
- if(n == zero)
- return m;
- m %= n;
- }
+ // Defer to the version in Boost.Math
+ return math::gcd( n, m );
}
template <typename IntType>
IntType lcm(IntType n, IntType m)
{
- // Avoid repeated construction
- IntType zero(0);
-
- if (n == zero || m == zero)
- return zero;
-
- n /= gcd(n, m);
- n *= m;
-
- if (n < zero)
- n = -n;
- return n;
+ // Defer to the version in Boost.Math
+ return math::lcm( n, m );
}
+#endif // BOOST_CONTROL_RATIONAL_HAS_GCD
class bad_rational : public std::domain_error
{
decrementable < rational<IntType>
> > > > > > > > > > > > > > > >
{
+ // Class-wide pre-conditions
+ BOOST_STATIC_ASSERT( ::std::numeric_limits<IntType>::is_specialized );
+
+ // Helper types
typedef typename boost::call_traits<IntType>::param_type param_type;
struct helper { IntType parts[2]; };
bool operator!() const { return !num; }
// Boolean conversion
+
+#if BOOST_WORKAROUND(__MWERKS__,<=0x3003)
+ // The "ISO C++ Template Parser" option in CW 8.3 chokes on the
+ // following, hence we selectively disable that option for the
+ // offending memfun.
+#pragma parse_mfunc_templ off
+#endif
+
operator bool_type() const { return operator !() ? 0 : &helper::parts; }
+#if BOOST_WORKAROUND(__MWERKS__,<=0x3003)
+#pragma parse_mfunc_templ reset
+#endif
+
// Comparison operators
bool operator< (const rational& r) const;
bool operator== (const rational& r) const;
// times. normalized form is defined as gcd(num,den) == 1 and den > 0.
// In particular, note that the implementation of abs() below relies
// on den always being positive.
+ bool test_invariant() const;
void normalize();
};
IntType r_num = r.num;
IntType r_den = r.den;
- IntType g = gcd(den, r_den);
+ IntType g = math::gcd(den, r_den);
den /= g; // = b1 from the calculations above
num = num * (r_den / g) + r_num * den;
- g = gcd(num, g);
+ g = math::gcd(num, g);
num /= g;
den *= r_den/g;
// This calculation avoids overflow, and minimises the number of expensive
// calculations. It corresponds exactly to the += case above
- IntType g = gcd(den, r_den);
+ IntType g = math::gcd(den, r_den);
den /= g;
num = num * (r_den / g) - r_num * den;
- g = gcd(num, g);
+ g = math::gcd(num, g);
num /= g;
den *= r_den/g;
IntType r_den = r.den;
// Avoid overflow and preserve normalization
- IntType gcd1 = gcd<IntType>(num, r_den);
- IntType gcd2 = gcd<IntType>(r_num, den);
+ IntType gcd1 = math::gcd(num, r_den);
+ IntType gcd2 = math::gcd(r_num, den);
num = (num/gcd1) * (r_num/gcd2);
den = (den/gcd2) * (r_den/gcd1);
return *this;
return *this;
// Avoid overflow and preserve normalization
- IntType gcd1 = gcd<IntType>(num, r_num);
- IntType gcd2 = gcd<IntType>(r_den, den);
+ IntType gcd1 = math::gcd(num, r_num);
+ IntType gcd2 = math::gcd(r_den, den);
num = (num/gcd1) * (r_den/gcd2);
den = (den/gcd2) * (r_num/gcd1);
bool rational<IntType>::operator< (const rational<IntType>& r) const
{
// Avoid repeated construction
- IntType zero(0);
+ int_type const zero( 0 );
+
+ // This should really be a class-wide invariant. The reason for these
+ // checks is that for 2's complement systems, INT_MIN has no corresponding
+ // positive, so negating it during normalization keeps it INT_MIN, which
+ // is bad for later calculations that assume a positive denominator.
+ BOOST_ASSERT( this->den > zero );
+ BOOST_ASSERT( r.den > zero );
+
+ // Determine relative order by expanding each value to its simple continued
+ // fraction representation using the Euclidian GCD algorithm.
+ struct { int_type n, d, q, r; } ts = { this->num, this->den, this->num /
+ this->den, this->num % this->den }, rs = { r.num, r.den, r.num / r.den,
+ r.num % r.den };
+ unsigned reverse = 0u;
+
+ // Normalize negative moduli by repeatedly adding the (positive) denominator
+ // and decrementing the quotient. Later cycles should have all positive
+ // values, so this only has to be done for the first cycle. (The rules of
+ // C++ require a nonnegative quotient & remainder for a nonnegative dividend
+ // & positive divisor.)
+ while ( ts.r < zero ) { ts.r += ts.d; --ts.q; }
+ while ( rs.r < zero ) { rs.r += rs.d; --rs.q; }
+
+ // Loop through and compare each variable's continued-fraction components
+ while ( true )
+ {
+ // The quotients of the current cycle are the continued-fraction
+ // components. Comparing two c.f. is comparing their sequences,
+ // stopping at the first difference.
+ if ( ts.q != rs.q )
+ {
+ // Since reciprocation changes the relative order of two variables,
+ // and c.f. use reciprocals, the less/greater-than test reverses
+ // after each index. (Start w/ non-reversed @ whole-number place.)
+ return reverse ? ts.q > rs.q : ts.q < rs.q;
+ }
+
+ // Prepare the next cycle
+ reverse ^= 1u;
+
+ if ( (ts.r == zero) || (rs.r == zero) )
+ {
+ // At least one variable's c.f. expansion has ended
+ break;
+ }
+
+ ts.n = ts.d; ts.d = ts.r;
+ ts.q = ts.n / ts.d; ts.r = ts.n % ts.d;
+ rs.n = rs.d; rs.d = rs.r;
+ rs.q = rs.n / rs.d; rs.r = rs.n % rs.d;
+ }
- // If the two values have different signs, we don't need to do the
- // expensive calculations below. We take advantage here of the fact
- // that the denominator is always positive.
- if (num < zero && r.num >= zero) // -ve < +ve
- return true;
- if (num >= zero && r.num <= zero) // +ve or zero is not < -ve or zero
+ // Compare infinity-valued components for otherwise equal sequences
+ if ( ts.r == rs.r )
+ {
+ // Both remainders are zero, so the next (and subsequent) c.f.
+ // components for both sequences are infinity. Therefore, the sequences
+ // and their corresponding values are equal.
return false;
-
- // Avoid overflow
- IntType gcd1 = gcd<IntType>(num, r.num);
- IntType gcd2 = gcd<IntType>(r.den, den);
- return (num/gcd1) * (r.den/gcd2) < (den/gcd2) * (r.num/gcd1);
+ }
+ else
+ {
+#ifdef BOOST_MSVC
+#pragma warning(push)
+#pragma warning(disable:4800)
+#endif
+ // Exactly one of the remainders is zero, so all following c.f.
+ // components of that variable are infinity, while the other variable
+ // has a finite next c.f. component. So that other variable has the
+ // lesser value (modulo the reversal flag!).
+ return ( ts.r != zero ) != static_cast<bool>( reverse );
+#ifdef BOOST_MSVC
+#pragma warning(pop)
+#endif
+ }
}
template <typename IntType>
bool rational<IntType>::operator< (param_type i) const
{
// Avoid repeated construction
- IntType zero(0);
-
- // If the two values have different signs, we don't need to do the
- // expensive calculations below. We take advantage here of the fact
- // that the denominator is always positive.
- if (num < zero && i >= zero) // -ve < +ve
- return true;
- if (num >= zero && i <= zero) // +ve or zero is not < -ve or zero
- return false;
-
- // Now, use the fact that n/d truncates towards zero as long as n and d
- // are both positive.
- // Divide instead of multiplying to avoid overflow issues. Of course,
- // division may be slower, but accuracy is more important than speed...
- if (num > zero)
- return (num/den) < i;
- else
- return -i < (-num/den);
+ int_type const zero( 0 );
+
+ // Break value into mixed-fraction form, w/ always-nonnegative remainder
+ BOOST_ASSERT( this->den > zero );
+ int_type q = this->num / this->den, r = this->num % this->den;
+ while ( r < zero ) { r += this->den; --q; }
+
+ // Compare with just the quotient, since the remainder always bumps the
+ // value up. [Since q = floor(n/d), and if n/d < i then q < i, if n/d == i
+ // then q == i, if n/d == i + r/d then q == i, and if n/d >= i + 1 then
+ // q >= i + 1 > i; therefore n/d < i iff q < i.]
+ return q < i;
}
template <typename IntType>
return ((den == IntType(1)) && (num == i));
}
+// Invariant check
+template <typename IntType>
+inline bool rational<IntType>::test_invariant() const
+{
+ return ( this->den > int_type(0) ) && ( math::gcd(this->num, this->den) ==
+ int_type(1) );
+}
+
// Normalisation
template <typename IntType>
void rational<IntType>::normalize()
return;
}
- IntType g = gcd<IntType>(num, den);
+ IntType g = math::gcd(num, den);
num /= g;
den /= g;
num = -num;
den = -den;
}
+
+ BOOST_ASSERT( this->test_invariant() );
}
namespace detail {
// Type conversion
template <typename T, typename IntType>
-inline T rational_cast(const rational<IntType>& src)
+inline T rational_cast(
+ const rational<IntType>& src BOOST_APPEND_EXPLICIT_TEMPLATE_TYPE(T))
{
- return static_cast<T>(src.numerator())/src.denominator();
+ return static_cast<T>(src.numerator())/static_cast<T>(src.denominator());
}
// Do not use any abs() defined on IntType - it isn't worth it, given the
projects / lyx.git / blobdiff