Rational: [PATCH] Eliminate excessive calls to gcd and add operator%
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, ⡍⠁⠗⠊⠕
AMDG 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.
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.
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.
In Christ, Steven Watanabe
Steven Watanabe <watanabesj@gmail.com> 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, ⡍⠁⠗⠊⠕
AMDG On 05/21/2012 03:11 AM, Mario Lang wrote:
Attached is another version of my proposed patch with all traces of operator% removed, but with the elimination of useless gcd retained.
<snip>
This version looks good. The tests pass with msvc and gcc. I've also verified by inspection that this patch can't cause integer overflow. Committed in r78537. In Christ, Steven Watanabe
participants (2)
-
Mario Lang -
Steven Watanabe