Rational sum not optimal for overflows

[nobody]

The sum of n1/d1 + n2/d2 should be calculated as:

L*((n1'*d2'+n2'*d1')/m)
-----------------------
     (k/m)*d1'*d2'

where

m = gcd(k, n1'*d2'+n2'*d1')
k = gcd(d1, d2)
L = gcd(n1, n2)

k*d1'=d1
k*d2'=d2
L*n1'=n1
L*n2'=n2

This way overflows are avoided while calculatins the 
resulting numerator (since division by m is performed 
before multiplication by L).

Manuel.Sequeira@iscte.pt

[pmoore]

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]

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]

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]

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

[Markus Schöpflin]

Logged In: YES 
user_id=91733

Regarding the discussion below, I think this issue can
safely be declared as closed/rejected.