Tonight I hated julia, REPL and overflows

I did a quick test of Int256. In this test it is slower that BigInt, probably because of divisions…

function test_hm(rtype,rank)
  hm(rank)=[rtype(1)//(n+m) for n in 1:rank, m in 1:rank]
  m=hm(rank)
  one(m)==m*inv(m)
end
julia> @btime test_hm(Int256,45)
  3.213 s (56720623 allocations: 1.12 GiB)
true

julia> @btime test_hm(BigInt,45)
  703.989 ms (10933914 allocations: 210.82 MiB)
true
1 Like

And that’s not what I was saying. 30 out of 10^9 doesn’t matter. The question is how many expressions in your code overflows. Only repeated calculation that can overflow matters.

Well, actually, it is a single one (complicated) expression which overflows (about 30 times in 10^9).
The values of the variables in the expression are different each time and I cannot predict which values
will make the expression overflow.

Errrr, a complicated expression is multiple expressions.

And in any case, what’s special about typemax(Int128) in your calculations for it to be a effective upper bound of those numbers or what’s the distribution of the numbers involved.

2^128 is just the largest number for which integer arithmetic is fast. The majority of numbers (perhaps 99%) actually fit in Int32. Perhaps 1% do not fit in Int32, and 0.001% do not fit in Int64. This kind of Poisson law is typical of what happens in my computations…

If this kind of distribution is what you have then it’s better to just detect the error and redo the ones with higher precision. In this case you still pretty much want to optimize only for the cases that fits instead of trying to make every operation safe and slower.

1 Like

detect the error and redo with higher precision is exactly what the union type we are discussing can do automatically for you, saving you a lot of difficult programming

1 Like

It would be great if someone who has a need for this would create a package for this kind of approach as an experiment.

Obviously there would be a performance penalty compared to Int, but I am not sure how big that would be with modern CPUs if the BigInt branch is actually rarely taken. It may just be a tolerable trade-off for some applications, especially compared to always using BigInt, but we won’t know until it is tried.

Look at the Nemo package, Nemo.fmpz is a wrapper over flint big integers, which uses GMP for actual big numbers, but otherwise stores an small integer inline. Unfortunately Nemo.fmpz is not currently a subtype of Integer, so this can limit its applicability as a replacement for BigInt.

2 Likes

No not for each operation, for the whole complex expression.

this situation actually begs for a julian Flint.jl which everybody would benefit from.
But I’m not going to open this can of worms :wink:

Indeed, (after 2 definitions which seem missing in Nemo):

julia> Base.one(::Type{fmpq})=fmpq(1)
julia> Base.:/(a::fmpq,b::fmpq)=a//b
julia> @btime test_hm(Nemo.fmpz,45)
  79.831 ms (466004 allocations: 15.86 MiB)
true

which is 10 times faster than BigInt!

3 Likes

Here’s a naïve proof of concept that someone might be interested in fleshing out (I’m only implementing multiplication here):

using Base: mul_with_overflow 

struct SometimesBigInt
    value::Union{Int, Base.RefValue{BigInt}}
end
SometimesBigInt(x::BigInt) = SometimesBigInt(Ref(x))

function Base.show(io::IO, x::SometimesBigInt)
    print(io, "SometimesBigInt(", value(x), ")")
end

function value(x::SometimesBigInt)::Union{Int, BigInt}
    x.value[]
end

function Base.:(*)(x::SometimesBigInt, y::SometimesBigInt)
    xv, yv = value(x), value(y)
    if xv isa Int && yv isa Int
        zv, overflowed = mul_with_overflow(xv, yv)
        if overflowed
            zv = big(xv) * big(yv)
        end
    else
        zv = xv * yv
    end
    SometimesBigInt(zv)
end

Okay, now let’s benchmark multiplication in the case where we have a 1% chance of overflow:

nums = (1:9..., 10^10) 
xs = rand(nums, 1000) # 10% chance of getting 10^10
ys = rand(nums, 1000) 
# for any element of xs .* ys we have 1% chance of overflow

function bench(xs, ys)
    sbigxs = SometimesBigInt.(xs)
    sbigys = SometimesBigInt.(ys)

    bigxs = big.(xs)
    bigys = big.(ys)
    print("Int mul:          ")  
    @btime $xs .* $ys
    print("SomeTimesBig mul: ")
    @btime $sbigxs .* $sbigys
    print("BigInt mul:       ")
    @btime $bigxs .* $bigys
    nothing
end


julia>  bench(xs, ys)
 Int mul:            515.380 ns (1 allocation: 7.94 KiB)
 SomeTimesBig mul:   22.607 μs (1274 allocations: 27.97 KiB)
 BigInt mul:         114.618 μs (3001 allocations: 54.81 KiB)

So we’re 42x slower than regular (overflowing) Int multiplication and 4x faster than making the entire array BigInt.

If we make it a 0% chance of overflow, we see the base overhead of SometimesBig relative to Int:

nums = (1:10...,) 
xs = rand(nums, 1000) 
ys = rand(nums, 1000) 

julia> bench(xs, ys)
 Int mul:            525.717 ns (1 allocation: 7.94 KiB)
 SomeTimesBig mul:   16.719 μs (1001 allocations: 23.56 KiB)
 BigInt mul:         113.401 μs (3001 allocations: 54.81 KiB)

and if we give it a 100% chance of overflowing, we see how much better off we would have been if we’d just used BigInt in the first place:

nums = ((10^10 .+ (1:10))...,) 
xs = rand(nums, 1000) 
ys = rand(nums, 1000) 

julia> bench(xs, ys)
 Int mul:            520.681 ns (1 allocation: 7.94 KiB)
 SomeTimesBig mul:   272.308 μs (9001 allocations: 164.19 KiB)
 BigInt mul:         116.686 μs (3001 allocations: 54.81 KiB)

I think given these benchmarks, there’s a case to be made for this sort of approach.

It would be better for performance to only store an Int and a Bool flag telling you whether to interpret that Int as a number or a pointer to a bigint, (that way SometimesBigInt would be an isbits type), but I don’t know enough about protecting pointers from GC to do that.

13 Likes

I’ve been thinking the same thing for some years, for integers and for decimal floating point - one thing would be to not use BigInt, but rather a Vector allocated with Base.StringVector (i.e. low memory overhead, only sizeof(Ptr).
Currently, because there’s no way of having Julia’s GC look at a value and decide what type (of a union type) it is (for example, using the bottom bits of a pointer to determine what type of pointer, or a to store an integer, or several characters, in the upper bits), I think one way would be to have a struct with two parts, one a Int, and the other a String that is set to a known value if the Int needs to be used, otherwise the String can hold the vector of chunks of the big int (or decimal)).

1 Like

I did a quick test of the approach I described, and it works out to be twice as fast as SomeTimesBigInt,
and only about a 12% penalty over using BigInt (instead of 2.34x)

Could you test the following code (inverting a Hilbert matrix)

function test_hm(rtype,rank)
  m=[rtype(1)//(n+m) for n in 1:rank, m in 1:rank]
  one(m)==m*inv(m)
end

For rank=45 and rtype=yourtype (you need enough definitions so that Rationals{yourtype} work but that should be mostly automatic).

The best I got so far is the type fmpz of Nemo, which is 10 times faster than rtype=BigInt.

2 Likes

I couldn’t get that to work, either with my type or yours - it gets the following error when I try it:

test_hm(SBig
ERROR: LoadError: OverflowError: 129997293855212511 * 3483 overflowed for type Int64
Stacktrace:
 [1] throw_overflowerr_binaryop(::Symbol, ::Int64, ::Int64) at ./checked.jl:154
 [2] checked_mul at ./checked.jl:288 [inlined]
 [3] - at ./rational.jl:264 [inlined]
 [4] generic_lufact!(::Array{Rational{Int64},2}, ::Val{true}; check::Bool) at /Users/scott/j/

It certainly works for the rtype=BigInt and rtype=fmpz (for the last one with the 2 additional definitions I give above). I thought the whole point of your type was to switch to BigInt when meeting overflow.
The error message suggests that you don’t handle an occuring overflow. Or perhaps you do not have the correct promotion rule. You should have

promote_type(yourtype, Int)=your_type 
promote_type(yourtype,BigInt)=BigInt # yourtype would also do

If you show your code I could try to debug it…

I think you need to explicitely overload Base.checked_mul for this to work.