How to avoid spurious overflow in rational-number addition?

Using Swift
, ,
1

Suppose one wishes to implement a rational-number type. We can see examples of this in NumericAnnex, SwiftRational, NumberKit, and presumably more projects that I did not find in a brief search.

Suppose further that the Rational type uses Int to store its numerator and denominator, which are kept relatively prime with positive denominator, and that we have a gcd function for Int (eg. the one from Numerics).

As part of the implementation, one must write an addition operator. That is, given rational numbers u = a/b and v = c/d, we need to compute u + v = (a*d + b*c) / (b*d). And if the new numerator and denominator share any common factors, they must be canceled out.

The question is, how can we compute the result without spurious overflow? In other words, if the resulting numerator and denominator can both be represented as Int, how do we compute them without trapping?

Of course we can factor out g = gcd(b, d), and instead calculate (a*(d/g) + (b/g)*c) / (b*(d/g)). This greatly improves the situation, but spurious overflow may still occur.

To illustrate the problem, imagine for the moment that we used Int8 instead of Int, so the maximum representable value is 127. Now suppose we have the fractions: u = 16/35 and v = 39/85.

The sum u + v = 109/119 is indeed representable with Int8 numerator and denominator. But consider how to compute it:

a = 16, b = 35
c = 39, d = 85
g = gcd(35, 85) = 5
b/g = 35/5 = 7
d/g = 85/5 = 17
a*(d/g) = 16*17 = 272
(b/g)*c = 7*39 = 273
a*(d/g) + (b/g)*c = 272 + 273 = 545
b*(d/g) = 35*17 = 595

Thus we have a numerator of 545 and a denominator of 595. And only now, at the very end, does it become apparent that they share a common factor of 5, so the result reduces to (545/5) / (595/5) = 109/119.

But several of the intermediate steps have overflowed. Namely, 272, 273, 545, and 595 are all greater than 127. Indeed, they are even greater than 255, so working with magnitudes as UInt8 would not eliminate the difficulty.

We can compute the intermediate products and sums using a combination of multipliedFullWidth and addingReportingOverflow, which yields a pair (Int8, UInt8) for both numerator and denominator (rather, (Int, UInt) in a real implementation).

We could implement a gcd function on these pairs, then call dividingFullWidth to cancel the gcd out of the resulting numerator and denominator. This would work, though I’m not sure the best strategy for writing that gcd function.

(We might also have to preemptively check the size of the pair before dividing, because the documentation for dividingFullWidth merely says it “may” trap if the result overflows, but does not guarantee it.)

Is this the correct approach, or is there a better way?

3 Likes
2

…okay I think I’ve made progress.

Since gcd(a, b) = 1 and gcd(c, d) = 1, it follows that when we calculate the numerator as a*(d/g) + (b/g)*c, and the denominator as b*(d/g), any common factor they share must itself be a factor of g.

To see this, call the common factor f. If f does not divide g, then it has a prime factor which appears more times in one of b or d than the other. Without loss of generality suppose such a prime p appears in b/g and not d/g. Then p is not a factor of a because gcd(a, b) = 1. Hence a*(d/g) is not divisible by p, and so, because b/g is divisible by p, we find that the entire numerator is not divisible by p. But that means f cannot have p as a factor, because f is meant to be a common divisor of the numerator and denominator.

All of that is to say, we don’t need to compute the gcd of the full-width numerator and denominator. Instead, we only need to compute the gcd of g and the full-width numerator, which we can do by first full-width dividing and then taking the gcd of g and the remainder. In other words, we can write:

let wideNum = ...  // the full-width numerator as (Int, UInt)
let (_, remainder) = g.dividingFullWidth(wideNum)
let commonFactor = gcd(g, remainder)
let (numerator, _) = commonFactor.dividingFullWidth(wideNum)
let denominator = (b/commonFactor) * (d/g)

I did not include preemptive checking for dividingFullWidth potentially overflowing. I want it to trap on overflow. Do I actually need to check? If so, what’s the right way to do so?

Suppose instead that I want to report the overflow, as in a Rational.addingReportingOverflow method. How can I detect when Int.dividingFullWidth will overflow, without actually trapping?

2 Likes
3

You either need to have the not currently existing dividingFullWidthReportingOverflow or check all the conditions upfront manually to keep it from trapping.

1 Like
4

The unsigned case is easier to analyze: dividingFullWidth can only overflow if the high word of the dividend is greater than or equal to the divisor. The signed case is about the same, just fussier because of the asymmetry in representable values. (I would probably recommend working with unsigned numerator and denominator plus a sign for this and other reasons, but it's not really any harder, just fussier...)

You can write an easy-and-obviously-correct-but-not-as-efficient-as-it-could-be check for the signed case by doing divisor.multipliedFullWidth(by: T.min), and also T.max, and then doing a double-word compare. The bounds this gives you plus one bigger and one smaller are the obvious test cases for the more efficient implementation.

(I haven't read it too closely, but I think that in your use, if the divide overflows then the actual result is not representable, so this might all be a moot point for you?)

3 Likes
5

…well, it’s moot if-and-only-if a trap is guaranteed to occur. But dividingFullWidth does not appear to require that behavior. Its documentation just says it “may” trap.

So if some implementation somewhere decides to implement that function by, say, silently wrapping on overflow, that would technically be valid. Which is not what I want.

(On the other hand, it might also be nice to fail gracefully. For example, when expanding a continued fraction as far as possible within the bounds of the type, we’d like to return the last convergent that is representable. But that requires attempting to compute the next one and seeing that it overflows.)

6

Oh, that's definitely not allowed (I think there's a open bug where this can happen in the REPL, but that's a bug). The documentation should just be clarified there.

1 Like
7

Should I open a PR to fix the documentation?

8

Sure!

I think that the current documentation was written to accurately reflect the current state of affairs, while allowing the desired implementation (trapping). So we should also fix the implementation, which has been somewhere on my list of low-priority stuff-that-should-be-fixed for a while, but without actual users who need the desired behavior, it was hard to justify bumping up the priority.

1 Like
9

I think those are not actually the correct bounds for signed integers.

Given that divisor is positive, the maximum possible result of the division is .max with a remainder of divisor - 1, which occurs when the full-width starting value is .max * divisor + divisor - 1. So the first positive input that overflows is .max * divisor + divisor, equivalent to (divisor >> 1, divisor.magnitude << (.bitWidth - 1)).

Similarly, the minimum possible result is .min with a remainder of -(divisor - 1), which occurs for starting value .min * divisor - (divisor - 1). So the first negative input that overflows is .min * divisor - divisor, and that is slightly more tricky to compute:

var high = (-divisor) >> 1
var low = divisor.magnitude << (.bitWidth - 1)
if low == 0 { high -= 1 }
low &-= divisor.magnitude
return (high, low)

I’m moderately confident I did that correctly, and your point about it being simpler in unsigned stands.

…of course, in the specific case of this thread, the final division by commonFactor is guaranteed to have a remainder of zero, which can perhaps simplify things.

1 Like
10

Yeah, sorry for being unclear, I was looking at your use site specifically. You have the right bounds for the general problem.

1 Like
11

If, instead of an Int, a list of prime factors were used, would the problem become easier to tackle? :slight_smile:

12

Easier in some ways, harder in others. It's a perfectly reasonable number representation, but not used very often for reasons you will discover if you implement it... =)

1 Like
13

I’m rather skeptical that prime factorization would make addition easier. Multiplication and division perhaps, but addition? Unlikely.

That said, for the specific things I’m doing, the denominators tend to have many small-prime factors, so I did consider a hybrid approach that tracks powers of the first few primes, alongside a non-factored fraction for the rest.

But that doesn’t really end up helping me, because the numerators I’m working with tend to have large prime factors.

1 Like
14

Why? It is simpler to start by implementing an arbitrary precision integer type (such as in some of the packages you link to) and use that for the numerator and denominator. Or does that defeat the purpose of whatever you are really trying to solve?

15

It's not how I read the headline post that presumably wanted a fixed width rational number. Using an arbitrary precision integer type in the implementation would be a major overkill (due to dynamic memory / ARC overhead).

16

FYI: it turns out that some implementations--in particular small integer types--did silently wrap on overflow! I've gone ahead and put together a PR to fix that and the documentation as well.

3 Likes
17

Ooo, thanks for getting to that! I'd noticed it as a problem back in 2018 but then promptly forgot.

1 Like