Modify ↓

Opened 16 years ago

Closed 12 years ago

## #62 closed Bugs (Rejected)

# Rational sum not optimal for overflows

Reported by: | nobody | Owned by: | Jonathan Turkanis |
---|---|---|---|

Milestone: | Component: | None | |

Version: | None | Severity: | |

Keywords: | Cc: |

### Description

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

### Attachments (0)

### Change History (5)

### comment:2 Changed 16 years ago by

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

### comment:3 Changed 16 years ago by

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

### comment:4 Changed 16 years ago by

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

### comment:5 Changed 12 years ago by

Status: | assigned → closed |
---|

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

**Note:**See TracTickets for help on using tickets.