Reliability of rational arithmetic

[khinsen]

Consider the following code that computes a cubic root using Newton's method from a given initial value:

function root(f, f′, x0, eps)
    x = x0
    fx = f(x)
    while abs(fx) > eps
        x = x - fx/f′(x)
        fx = f(x)
    end
    x
end

cube_root(x0, eps) = root(x->x^3-1//1, x->3*x^2, x0, eps)

Let's do this in the standard way, using floating-point arithmetic (Complex128):

cr_complex128 = cube_root(-149/100 - 3/2*im, 1/10)
println(cr_complex128)

The result is -0.5185349749124928 - 0.8466618927868804im, which is indeed one of the cubic roots of 1.

Since the root that the iteration converges to is hard to predict, let's try again using rational arithmetic:

cr_rational64 = cube_root(-149//100 - 3//2*im, 1//10)
println(cr_rational64)

Bad luck:

ERROR: integer division error
 in Rational at rational.jl:8
 in - at complex.jl:111
 in root at rational_root.jl:5

This can be tracked down to the usual issue with rational arithmetic: the denominators of all the intermediate results grow rapidly and overflow the Int64 that is available for storage. Don't even consider running this on a 32 bit machine!

So let's try BigInt:

cr_rational_big = cube_root(-149//big(100) - 3//2*im, 1//big(10))
println(cr_rational_big, " => ", convert(Complex128, cr_rational_big))

with the result:

-89570859265437960860752701388197033469884493568628250226810631152136502759702589412969256261337092079319544000799//172738317758708175397627917416341959088693237940646365283701126300442423569679831539574834489120285397196126445025 - 5850038042816378545405701188965568851268442922391815363213624265020034514028087001761881208976412103870377254856//6909532710348327015905116696653678363547729517625854611348045052017696942787193261582993379564811415887845057801*im => -0.5185349749124929 - 0.8466618927868804im

That's what I think everyone would have expected to get - but it takes considerable care and ugly syntax to get there.

Finally, let's see the remaining intermediate option, Int128:

cr_rational128 = cube_root(-149//convert(Int128, 100) - 3//2*im, 1//convert(Int128, 10))
println(cr_rational128)

The result is surprising: 1//1 + 0//1*im. I suspect this is due to denominator overflow as well, but it's much worse because it leads to a wrong result.

Question: are there any situations where arithmetic Rational{Int64} is a good choice?

In my personal experience, rational arithmetic is mainly used for infinite-precision computations, where correctness is the top priority. This implies that Rational{BigInt} should be the default, in particular for the result of `//``.

[ScottPJones]

@khinsen I did bring up the issue of how the arithmetic in julia is unchecked, but I just looked at the code for rational numbers, and they are supposed to be used checked arithmetic (i.e. checked_add, checked_mul, etc) so you should have gotten an error instead of an incorrect result.
This seems more like an outright bug to me...

[pao]

Int128 checked math may be broken due to upstream; I believe it relies on an LLVM intrinsic, and LLVM's support for 128-bit integers is far from perfect. Note that Int64 does throw an error as expected.

[ScottPJones]

@pao, yes, that's what seemed strange to me... @khinsen's test with Int64 got the error, and Int128 a bad result... I didn't know that was tied so much to LLVM though. Thanks for the info!

[pao]

A reproducible standalone checked-math failure for Int128 including output from @code_llvm would be worth opening a new issue for, @khinsen, if you're feeling up to that.

[mbauman]

It's issue #4905, and is fixed upstream.

[simonbyrne]

Question: are there any situations where arithmetic Rational{Int64} is a good choice?

See #11522.

[JeffBezanson]

Also #2960

[pao]

Thanks @mbauman; I didn't have the right search term for that.

[khinsen]

@pao: I am not familiar enough with Julia internals to prepare a bug report about checked Int128 arithmetic, sorry! If it's a fixed problem in LLVM, I guess we just have to wait anyway.

@simonbyrne: After reading through #11522, I understand that some people are interested in high-performance rational arithmetic within the limits of Int64 denominators, and that this somehow relates to interpolation, but I didn't find any more explicit reference to a use case.

[simonbyrne]

@khinsen That issue was prompted by JuliaMath/Interpolations.jl#37

[pao]

I am not familiar enough with Julia internals to prepare a bug report about checked Int128 arithmetic, sorry! If it's a fixed problem in LLVM, I guess we just have to wait anyway.

That was #4905 that @mbauman mentioned, so no need. It would have amounted to "here's a small reproducible case for a failed checked arithmetic call." Nothing scary!

[simonbyrne]

Given that we probably don't want to promote to BigInt without warning, it seems that the only other bug here is the Int128 checked arithmetic, which is covered by #4905. In that case, can we close this?

[simonbyrne]

Closed by #14362.

[tkelman]

Can we add this as a test?