
Boost Users : 
Subject: Re: [Boostusers] Rational: [PATCH] Eliminate excessive calls to gcd
From: Mario Lang (mlang_at_[hidden])
Date: 20120521 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 workhorse 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 mixedmode 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 mixedmode multiplication
>> involved) and the definition of operator%=(param_type) saves two
>> unnecessary calls to gcd since there is one mixedmode multiplication
>> and one mixedmode division involved. Obviously, all code that involves
>> mixedmode 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 Offload 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 @@
// Mixedmode 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, â¡â â —â Šâ •
Boostusers list run by williamkempf at hotmail.com, kalb at libertysoft.com, bjorn.karlsson at readsoft.com, gregod at cs.rpi.edu, wekempf at cox.net