Re: [Boost-users] Rational: [PATCH] Eliminate excessive calls to gcd

21 May 2012

Steven Watanabe  writes:
...
On 05/20/2012 11:58 AM, Mario Lang wrote:
...
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(*this / other));
I have been using this definition for operator%= successfully in my code
for quite a while.
I'm not sure that this is a good idea.  For
integers, an important property is the relationship
between / and %:
(a / b) * a + a % b = a
This doesn't hold with your operator%.
In fact, your definition of % for rational
would also work for floating point, but
notice that C uses a named function, fmod
instead of an operator.  I'd prefer to
go that route with rational.
OK, I see this is a problem.  I guess what I am actually after is more
like a remainder operation then a modulo.
...
...
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.
This part looks good.
...
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!
This looks good.
Nice, thanks for your feedback.

Attached is another version of my proposed patch with all traces of
operator% removed, but with the elimination of useless gcd retained.

Index: rational.hpp
===================================================================

--- rational.hpp	(Revision 78525)
+++ rational.hpp	(Arbeitskopie)
@@ -21,6 +21,7 @@
 //    Nickolay Mladenov, for the implementation of operator+=
 
 //  Revision History
+//  05 May 12  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;
@@ -141,7 +142,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);
@@ -156,14 +157,14 @@
     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);
 
     // 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; }
@@ -330,46 +331,36 @@
 // Mixed-mode operators
 template <typename IntType>
 inline rational<IntType>&
-rational<IntType>::operator+= (param_type i)
+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;
 
-template <typename IntType>
-inline rational<IntType>&
-rational<IntType>::operator-= (param_type i)
-{
-    return operator-= (rational<IntType>(i));
+    return *this;
 }
 
 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>&
 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;
-    return *this;
-}
+    if (i == zero) throw bad_rational();
+    if (num == zero) return *this;
 
-template <typename IntType>
-inline const rational<IntType>& rational<IntType>::operator--()
-{
-    // This can never denormalise the fraction
-    num -= den;
+    // 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;
 }
 
@@ -477,12 +468,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>
@@ -597,10 +583,7 @@
 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,
  ⡍⠁⠗⠊⠕

    