[Boost-users] Rational: [PATCH] Eliminate excessive calls to gcd and add operator%

Subject: [Boost-users] Rational: [PATCH] Eliminate excessive calls to gcd and add operator%
From: Mario Lang (mlang_at_[hidden])
Date: 2012-05-20 14:58:13

Hi.

I am using boost/rational.hpp in my code to represent musical time with
rational numbers.  I have found some possible improvements.  These
improvements have accumulated into a small patch which I'd like to
present for review.  I think it is worth applying.  There doesn't appear
to be a maintainer for rational.hpp though, so I am writing to the list.

The story: In my application I need a remainder operation. Looking at
rational.hpp and its work-horse operators.hpp it became rather obvious
that rational.hpp can and (in my opinion) should implement operator%=
and friends.  operator%= can easily be defined as

return operator-= (other * rational_cast<int_type>(*this / other));

I have been using this definition for operator%= successfully in my code
for quite a while.  During a profiling session it showed up which made
me think a bit more.  Actually, the current way how rational.hpp
implements mixed-mode operators is a waste of performance.  It calls the
corresponding operator with a temporary rational where the denominator
is always equal to 1.  However, this means that there are excessive
(unnecessary) calls to boost::math::gcd performed by the operators.  For
instance, addition of two rational numbers needs two calls to gcd.  On
the other hand, addition of a rational number and an integer can be
performed without any gcd at all (n += d * i).  Similarily,
multiplication of two rational numbers needs two calls to gcd, while
multiplication of a rational with an integer only needs one call to gcd.
So it is definitely worth it to actually define the mixed mode operators
not in terms of the already existing operators but reimplement them with
the unnecessary calls to gcd removed.

After this, the definition of operator%=(const rational&) results in one
less call to gcd (because there is one mixed-mode multiplication
involved) and the definition of operator%=(param_type) saves two
unnecessary calls to gcd since there is one mixed-mode multiplication
and one mixed-mode division involved.  Obviously, all code that involves
mixed-mode operators benefits from these changes.

As a side effect of reading code I noticed that the definition of
operator=(param_type) (assign from an integer) also does unnecessary
gcd.  It currently calls assign(n, 1) which itself calls normalize()
which does gcd.  However, as above, we know that there is no need to
calculate gcd of (n, 1) since it will always be 1.  So I changed the
definition of operator=(param_type) to directly assign the parameter to
num and set den to 1, instead of calling assign().  Another useless call
of gcd squashed!

On top of these, we can observe that now that boost::rational provides
operator%= it is actually an instance of ordered_euclidean_ring_operators.

Comments and committers welcome.

Index: rational.hpp
===================================================================
--- rational.hpp	(Revision 78523)
+++ rational.hpp	(Arbeitskopie)
@@ -21,6 +21,7 @@
 //    Nickolay Mladenov, for the implementation of operator+=
 
 //  Revision History
+//  05 May 12  Add operator% and reduced use of implicit gcd (Mario Lang)
 //  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;
@@ -105,23 +106,10 @@
 
 template <typename IntType>
 class rational :
-    less_than_comparable < rational<IntType>,
-    equality_comparable < rational<IntType>,
-    less_than_comparable2 < rational<IntType>, IntType,
-    equality_comparable2 < rational<IntType>, IntType,
-    addable < rational<IntType>,
-    subtractable < rational<IntType>,
-    multipliable < rational<IntType>,
-    dividable < rational<IntType>,
-    addable2 < rational<IntType>, IntType,
-    subtractable2 < rational<IntType>, IntType,
-    subtractable2_left < rational<IntType>, IntType,
-    multipliable2 < rational<IntType>, IntType,
-    dividable2 < rational<IntType>, IntType,
-    dividable2_left < rational<IntType>, IntType,
-    incrementable < rational<IntType>,
-    decrementable < rational<IntType>
-    > > > > > > > > > > > > > > > >
+    ordered_euclidean_ring_operators < rational<IntType>,
+    ordered_euclidean_ring_operators2 < rational<IntType>, IntType,
+    unit_steppable < rational<IntType>
+    > > >
 {
     // Class-wide pre-conditions
     BOOST_STATIC_ASSERT( ::std::numeric_limits<IntType>::is_specialized );
@@ -141,7 +129,7 @@
     // Default copy constructor and assignment are fine
 
     // Add assignment from IntType
-    rational& operator=(param_type n) { return assign(n, 1); }
+    rational& operator=(param_type i) { num = i; den = 1; return *this; }
 
     // Assign in place
     rational& assign(param_type n, param_type d);
@@ -155,15 +143,17 @@
     rational& operator-= (const rational& r);
     rational& operator*= (const rational& r);
     rational& operator/= (const rational& r);
+    rational& operator%= (const rational& r);
 
-    rational& operator+= (param_type i);
-    rational& operator-= (param_type i);
+    rational& operator+= (param_type i) { num += i * den; return *this; }
+    rational& operator-= (param_type i) { num -= i * den; return *this; }
     rational& operator*= (param_type i);
     rational& operator/= (param_type i);
+    rational& operator%= (param_type i);
 
     // Increment and decrement
-    const rational& operator++();
-    const rational& operator--();
+    const rational& operator++() { num += den; return *this; }
+    const rational& operator--() { num -= den; return *this; }
 
     // Operator not
     bool operator!() const { return !num; }
@@ -205,6 +195,14 @@
     void normalize();
 };
 
+// Type conversion
+template <typename T, typename IntType>
+inline T
+rational_cast(const rational<IntType>& r BOOST_APPEND_EXPLICIT_TEMPLATE_TYPE(T))
+{
+    return static_cast<T>(r.numerator()) / static_cast<T>(r.denominator());
+}
+
 // Assign in place
 template <typename IntType>
 inline rational<IntType>& rational<IntType>::assign(param_type n, param_type d)
@@ -327,50 +325,54 @@
     return *this;
 }
 
-// Mixed-mode operators
 template <typename IntType>
 inline rational<IntType>&
-rational<IntType>::operator+= (param_type i)
+rational<IntType>::operator%= (const rational<IntType>& r)
 {
-    return operator+= (rational<IntType>(i));
+    return operator-= (r * rational_cast<IntType>(*this / r));
 }
 
+// Mixed-mode operators
 template <typename IntType>
 inline rational<IntType>&
-rational<IntType>::operator-= (param_type i)
-{
-    return operator-= (rational<IntType>(i));
-}
-
-template <typename IntType>
-inline rational<IntType>&
 rational<IntType>::operator*= (param_type i)
 {
-    return operator*= (rational<IntType>(i));
+    // Avoid overflow and preserve normalization
+    IntType gcd = math::gcd(i, den);
+    num *= i / gcd;
+    den /= gcd;
+
+    return *this;
 }
 
 template <typename IntType>
-inline rational<IntType>&
+rational<IntType>&
 rational<IntType>::operator/= (param_type i)
 {
-    return operator/= (rational<IntType>(i));
-}
+    // Avoid repeated construction
+    IntType const zero(0);
 
-// Increment and decrement
-template <typename IntType>
-inline const rational<IntType>& rational<IntType>::operator++()
-{
-    // This can never denormalise the fraction
-    num += den;
+    if (i == zero) throw bad_rational();
+    if (num == zero) return *this;
+
+    // Avoid overflow and preserve normalization
+    IntType const gcd = math::gcd(num, i);
+    num /= gcd;
+    den *= i / gcd;
+
+    if (den < zero) {
+        num = -num;
+        den = -den;
+    }
+
     return *this;
 }
 
 template <typename IntType>
-inline const rational<IntType>& rational<IntType>::operator--()
+inline rational<IntType>&
+rational<IntType>::operator%= (param_type i)
 {
-    // This can never denormalise the fraction
-    num -= den;
-    return *this;
+    return operator-= (i * rational_cast<IntType>(*this / i));
 }
 
 // Comparison operators
@@ -477,12 +479,7 @@
 template <typename IntType>
 bool rational<IntType>::operator> (param_type i) const
 {
-    // Trap equality first
-    if (num == i && den == IntType(1))
-        return false;
-
-    // Otherwise, we can use operator<
-    return !operator<(i);
+    return operator==(i)? false: !operator<(i);
 }
 
 template <typename IntType>
@@ -583,24 +580,13 @@
     return os;
 }
 
-// Type conversion
-template <typename T, typename IntType>
-inline T rational_cast(
-    const rational<IntType>& src BOOST_APPEND_EXPLICIT_TEMPLATE_TYPE(T))
-{
-    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
 // difficulties involved (Koenig lookup required, there may not *be* an abs()
 // defined, etc etc).
 template <typename IntType>
 inline rational<IntType> abs(const rational<IntType>& r)
 {
-    if (r.numerator() >= IntType(0))
-        return r;
-
-    return rational<IntType>(-r.numerator(), r.denominator());
+    return r.numerator() >= IntType(0)? r: -r;
 }
 
 } // namespace boost

-- 
CYa,
  ⡍⠁⠗⠊⠕