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

Subject: Re: [Boost-users] Rational: [PATCH] Eliminate excessive calls to gcd
From: Mario Lang (mlang_at_[hidden])
Date: 2012-05-21 06:11:01

Steven Watanabe <watanabesj_at_[hidden]> 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<int_type>(*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,
  ⡍⠁⠗⠊⠕