Computing binomial(n,k) without overflow

[mike.ingold]

The function binomial(n,k) throws an OverflowError if the standard ::Int64 arguments are large enough that the result overflows an ::Int64. This can be modified to binomial(big(b),k) to promote the output to BigInt, e.g.:

julia> binomial(73, 24)
ERROR: OverflowError: binomial(73, 24) overflows
Stacktrace:
 [1] binomial(n::Int64, k::Int64)
   @ Base .\intfuncs.jl:1114
 [2] top-level scope
   @ REPL[32]:1

julia> binomial(big(73), 24)
11844267374132633700

I’m trying to write a function as generically as possible that internally computes a binomial, and occasionally the inputs are large enough to require a BigInt type. In those cases the output of the overall function winds up turning into BigFloat.

Does anybody have an elegant generic way to compute binomial where the types are promoted to BigInt i.f.f. a BigInt output would be required? A try/catch could work but seems probably suboptimal.

[mike.ingold]

Less abstractly, this is where the issue crops up for me.

mikeingold/LineIntegrals.jl/blob/ebb30955a9bea68360a83edc47afaa3aa6da1ba4/src/utils.jl#L14

In the short term I just forced everything up to BigInt, but the function outputs don’t require BigFloat level of precision even in those cases.

[jar1]

If your function’s return type depends on its argument values rather than just its inputs, then it isn’t “type stable”, which has some performance downsides and is more complex to use.

In the case where the result is small enough to fit in Int64, it’s faster to do try/catch than to always use big. In the overflow case, it’s faster to use big in the first place.

A third option is to convert to BigInt after the trycatch, so the function’s return value is statically known but it has internal dynamic types.

binomial_trycatch(n,k) =
    try
        binomial(n,k)
    catch e
        e isa OverflowError &&  binomial(big(n), k)
    end

binomial_trycatch_then_big(n,k)::BigInt =
    try
        binomial(n,k)
    catch e
        e isa OverflowError &&  binomial(big(n), k)
    end

binomial_big(n,k) = binomial(big(n),k)



using BenchmarkTools
# Fits in Int64.
@btime binomial_trycatch(73, 23);
@btime binomial_trycatch_then_big(73, 23);
@btime binomial_big(73, 23);
# Doesn't fit in Int64.
@btime binomial_trycatch(73, 24);
@btime binomial_trycatch_then_big(73, 24);
@btime binomial_big(73, 24);

julia> @btime binomial_trycatch(73, 23);
  177.726 ns (0 allocations: 0 bytes)

julia> @btime binomial_trycatch_then_big(73, 23);
  205.040 ns (2 allocations: 40 bytes)

julia> @btime binomial_big(73, 23);
  399.848 ns (12 allocations: 160 bytes)

julia> @btime binomial_trycatch(73, 24);
  37.892 μs (14 allocations: 248 bytes)

julia> @btime binomial_trycatch_then_big(73, 24);
  34.696 μs (14 allocations: 248 bytes)

julia> @btime binomial_big(73, 24);
  399.811 ns (11 allocations: 152 bytes)

[Benny]

Not sure from that exactly where the overhead of overflowing is, whether it’s recomputing binomial or catching an error, but I think the only way to cut down on that any further is to define a separate binomial function that breaks the while loop where Base.binomial would throw the error. Then, it resumes the algorithm with BigInt, possibly with a rewrite that picks up where the loop left off instead of starting over.

[mike.ingold]

If I set some threshold on n where this is likely to occur, is there any way to dispatch to a helper function based on that value?

In this situation the ‘n’ is essentially inherent to the BezierCurve type, but I don’t think it’s accessible as a type parameter.

[jar1]

Another option is to use Int128

julia> binomial_i128(n,k) = binomial(Int128(n),k)
binomial_i128 (generic function with 1 method)

julia> @btime binomial_i128(73, 23);
  50.938 ns (0 allocations: 0 bytes)

julia> @btime binomial_i128(73, 24);
  49.638 ns (0 allocations: 0 bytes)

[mike.ingold]

Ooh, good point. Looks like converting n to Int128 extends the non-erroring input range into sufficiently large numbers with good performance and prevents the need for the wrapping function to return a BigFloat.

@jar1 I marked your earlier post as the solution since it actually answered the original question, but this Int128 tip was the solution I actually needed. Thanks!

[Dan]

The accuracy of Float64 should be enough for most purposes. So perhaps the following will be enough:

using SpecialFunctions

logbinomial(n,k) = logfactorial(n)-logfactorial(k)-logfactorial(n-k)
B(i,n) = t -> exp(logbinomial(n,i) + i*log(t) + (n-i)*log1p(-t))

As an example:

julia> using UnicodePlots

julia> lineplot(0,1,B(3,6))
            ┌────────────────────────────────────────┐        
        0.4 │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│ #198(x)
            │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│        
            │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│        
            │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣀⠔⠒⠒⠢⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│        
            │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⠎⠁⠀⠀⠀⠀⠀⠑⡄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│        
            │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡰⠃⠀⠀⠀⠀⠀⠀⠀⠀⠘⢆⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│        
            │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡸⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⢣⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│        
   f(x)     │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡜⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⢣⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│        
            │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡜⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢣⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│        
            │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡜⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢣⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│        
            │⠀⠀⠀⠀⠀⠀⠀⠀⠀⡜⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢣⠀⠀⠀⠀⠀⠀⠀⠀⠀│        
            │⠀⠀⠀⠀⠀⠀⠀⢀⡜⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢣⡀⠀⠀⠀⠀⠀⠀⠀│        
            │⠀⠀⠀⠀⠀⠀⢀⠎⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠱⡄⠀⠀⠀⠀⠀⠀│        
            │⠀⠀⠀⠀⠀⡠⠃⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⢆⠀⠀⠀⠀⠀│        
          0 │⣀⣀⣀⠤⠊⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⠦⣄⣀⣀│        
            └────────────────────────────────────────┘        
            ⠀0⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀1⠀        

If performance becomes an issue, there might be faster ways to get at the results without the anonymous functions.

[PeterSimon]

The Int128 solution of @Jar1 always returns a Float64 for integer inputs. As pointed out by @Dan this is perfectly fine for @mike.ingold 's specific application because he’s using these coefficients in a floating point calculation. So another, faster way to get a floating point approximation of the binomial is

binomial_float(n, k) = prod((n+1-i)/i for i in 1:min(k, n-k); init=1.0)

julia> @btime binomial_i128(73, 23)
  51.368 ns (0 allocations: 0 bytes)
5.685248339583665e18

julia> @btime binomial_float(73, 23)
  19.138 ns (0 allocations: 0 bytes)
5.685248339583665e18

Edit: After looking over the source of binomial, when called with a non-integer first argument, I see it’s almost the same as above. So a better definition for binomial_float would be

binomial_float(n, k) = binomial(float(n), k)

[Dan]

@PeterSimon

19.138 ns (0 allocations: 0 bytes)

I think the timings here are due to compile-time optimization.

[PeterSimon]

julia> n = 73; k = 23; @btime binomial_float($n, $k)
  20.060 ns (0 allocations: 0 bytes)
5.685248339583665e18

[Dan]

Okay… still surprises me how fast computers have become over the years.

[PeterSimon]

Maybe SIMD?

[Dan]

My version also in the same ballpark:

julia> n = 73; k = 23; @btime exp(logbinomial($n, $k))
  37.290 ns (0 allocations: 0 bytes)
5.685248339583862e18

[PeterSimon]

Edited: After re-starting Julia I’m getting results consistent with you:

julia> n=73; k=23; @btime exp(logbinomial($n, $k))
  57.259 ns (0 allocations: 0 bytes)
5.685248339583862e18