Boost C++ Libraries: Ticket #62: Rational sum not optimal for overflows

pmoore — Thu, 21 Feb 2002 15:56:33 GMT

Logged In: YES
user_id=113328
This algorithm calculates 3 GCDs, 3 divisions and 10
multiplications. This is far more expensive (particularly
for user-defined integer types) than the existing
algorithm, which uses 2 GCDs, 4 divisions and 3
multiplications.
The existing algorithm (which I did not design, BTW - see
the source for credits) is supposed to avoid overflow. If
you have a specific test case where the existing algorithm
overflows, and yours doesn't, please can you supply it.
I'm not guaranteeing to implement your algorithm even if
there is a case where overflow occurs - the difference in
cost between the two algorithms is likely to be
prohibitive. (I feel that efficiency when the underlying
integer type is a user-defined unlimited-size integer is
more important than avoiding overflow in all possible
borderline cases with limited-size integers).
If you can provide a redesigned version of your algorithm
which is more efficient, then I will consider incorporating
it.

nobody — Thu, 21 Feb 2002 19:57:56 GMT

Logged In: NO
The existing code fails if a rational is added to itself (r
+= r):
IntType g = gcd(den, r.den);
den /= g;
num = num * (r.den / g) + r.num * den;
             ^^^^^ Opps... value changed!
g = gcd(num, g);
num /= g;
den *= r.den / g;
       ^^^^^ Opps... value changed!
This requires a slight change:
IntType g = gcd(den, r.den);
// Must deal with r += r:
IntType rd = r.den;
den /= g;
num = num * (rd / g) + r.num * den;
g = gcd(num, g);
num /= g;
den *= rd / g;
This code requires 2 gcd(), 4 divisions and 3
multiplications.
I propose the following alternative:
IntType zero(0);
if(r.num == zero)
    return *this;
IntType gn = gcd(num, r.num);
IntType g = gcd(den, r.den);
// Must deal with r += r:
IntType rd = r.den;
den /= g;
num = num / gn * (rd / g) + r.num / gn * den;
g = gcd(num, g);
num = gn * (num / g);
den *= rd / g;
The "if" is necessary in order to prevent calling gcd with
both arguments zero (alternatively gcd might be extended to
deal with that case by returning 1).
This code requires 3 gcd() (one more), 6 divisions (two
more) and 4 multiplications (one more)
The existing code overflows more easily than the proposed
one.  Suppose *this is 13/21 and r is -13/66 (of course one
may find similar examples with larger numbers).
// num = 13, den = 21, r.num = -13, r.den = 66
IntType g = gcd(den, r.den); // g = 3
// Must deal with r += r:
IntType rd = r.den; // rd = 66
den /= g; // den = 7
num = num * (rd / g) + r.num * den; // num = 13 * (66 / 3) -
 13 * 7 = 13 * 15 = 195
g = gcd(num, g); // g = gcd(195, 3) = 3
num /= g; // num = 65
den *= rd / g; // den = 7 * 66 / 3 = 154
Compare with what happens in the alternative:
// num = 13, den = 21, r.num = -13, r.den = 66
IntType zero(0);
if(r.num == zero)
    return *this;
IntType gn = gcd(num, r.num); // gn = 13
IntType g = gcd(den, r.den); // g = 3
// Must deal with r += r:
IntType rd = r.den; // rd = 66
den /= g; // den = 7
num = num / gn * (rd / g) + r.num / gn * den; // num =
13/13 * (66/3) - 13/13 * 7 = 15
g = gcd(num, g); // g = gdc(15, 3) = 3
num = gn * (num / g); // num = 13 * (15 / 3) = 13 * 5 = 65
den *= rd / g; // den = 7 * 66 / 3 = 154
Notice that num reaches 195 in the original code but
remains at 15 in the version I propose.
Manuel.Sequeira@iscte.pt

nobody — Sat, 23 Feb 2002 22:17:28 GMT

Logged In: NO
(SF isn't letting me log in. Grr.)
First, the latest version (in CVS) has the self-assignment
bug in += (and the other assignment operators) fixed.
Thanks for pointing it out.
As to the overflow issue, I'm confused. Your example of
(13/21) + (-13/66) results in 65/154 in both versions,
which is correct. Can you provide an example of a
calculation which fails using rational<int> in the current
version, and succeeds with your version. I expect it to
require unusually complicated fractions (ie, very large
numerators and/or denominators).
Your version does 1 extra GCD, 2 extra divisions, and 1
extra multiplication. This could be a serious overhead for
the key case, rational<LongInt> where LongInt is a user-
defined unlimited-precision integer type. In this case,
overflow cannot happen. On the other hand, working with
bigger numbers *could* be more expensive (but I doubt that
it would be significantly so...)
In the case of rational<int>, the operations are cheap, and
overflow *is* possible. But dealing with the issues of
using rationals based on a limited-precision integer type
are complex, and less well-understood than the related
issues with floating point. Frankly, anyone using
rational<int> for numbers which get even close to the
overflow point, needs to know what they are doing. And for
that sort of person, boost::rational isn't meant to be
appropriate.
So overall, I'm not too keen on this change. It hurts the
case I care about, and doesn't add a lot to the (arguably
more common, but nevertheless dangerous) case of
rational<int>. Unless you come up with some more persuasive
arguments, I'm not likely to change the library. I may be
willing to make the documentation (section 9, "Design
Notes", "Limited-range integer types") stronger regarding
the danger, but I don't want to make it seem like
scaremongering...
Please feel free to provide further arguments, but make
sure you understand the design goals (unlimited-precision
integer is key, limited-precision is likely to be common,
but will never be 100% safe, and really isn't to be
encouraged...)
Paul Moore

nobody — Mon, 25 Feb 2002 11:15:54 GMT

Logged In: NO
> I expect it to require unusually complicated fractions
> (ie, very large numerators and/or denominators).
Indeed.  Here's an example:
920182/8915 + 918356/4005
I'm sure it is possible to find fractions with smaller
terms which also fail under the current algorithm but pass
under the one I propose.
> Your version does 1 extra GCD, 2 extra divisions, and 1
> extra multiplication. This could be a serious overhead
> for the key case, rational<LongInt> where LongInt is a
> user-defined unlimited-precision integer type. In this
> case, overflow cannot happen. On the other hand, working
> with bigger numbers *could* be more expensive (but I
> doubt that it would be significantly so...)
Yes.  But, being so, why not drop the attempt to avoid
overflow altogether?  In my opinion this should be dealt
with by adding a member to boost integer_traits which
allows us to check at compile time whether we are dealing
with limited- or unlimited-precision integers.  Two
specialized algorithms might then exist: one for unlimited-
precision integers, which does not attempt at all to avoid
overflow, and another for limited-precision integers, which
tries as hard as possible to avoid such overflows.  This
should be a relatively simple thing to implement.
Manuel Menezes de Sequeira
Manuel.Sequeira@iscte.pt

status changed

Markus Schöpflin — Tue, 22 Nov 2005 15:35:48 GMT

status assigned → closed

Logged In: YES
user_id=91733
Regarding the discussion below, I think this issue can
safely be declared as closed/rejected.