Equivalent of "lchoose" function from R?

General Usage
1

In R there’s a function lchoose(a,b) which computes the logarithm of the binomial coefficient or choose(a,b) number (number of ways to choose b items from a set of a). This is the only way to work with really big numbers.

I’m trying to work out the probability of certain address collision situations in networks of devices using ipv6. There are 2^64 host suffixes in an ipv6 address so I need to calculate things like choose(2^64,100) and obviously I need to do this with the logarithm because that number is going to be astronomically ridiculous (lchoose(2^64,100) ~ 4072 according to R)

Is there a Julia way to do this. Yes I’m aware I could use stirling’s approximation and work directly with that, but figured there’s probably already something defined somewhere?

stirlingappx(x) = x*log(x) - x
lchoose(a,b) = stirlingappx(a) - stirlingappx(b) - stirlingappx(a-b)

Works ok for my problem, but maybe there’s already a lchoose like function?

2

There might be a better answer, but I would like to point out that this seems to work fine:

julia> log(binomial(BigInt(2)^64, 100))
4072.402580028086489857865546272512049756065408011813147713715631592357374294182

julia> log(binomial(big"2"^64, 100))
4072.402580028086489857865546272512049756065408011813147713715631592357374294182
1 Like
3

This will be very slow for large arguments. I’d probably recommend using lgamma from SpecialFunctions

4

Looking at SpecialFunctions.jl, there seems to be a logabsbinomial that looks appropriate.

https://specialfunctions.juliamath.org/stable/functions_list/#SpecialFunctions.logabsbinomial

julia> using SpecialFunctions

julia> logabsbinomial(big"2"^64, 100) |> first
4072.40258002808648985786554627251204975606540801181314771372293997172075937879
1 Like
5

Perfect.

lchoose(a,b) = lgamma(a+1) - lgamma(b+1) - lgamma(a-b+1)

Works fine!

or logabsbinomial maybe is more accurate? they diverge in the 12 or 13th decimal place. Doesn’t really matter, the important thing is any of these approximate calcs show my probability is very close to 1 :slight_smile:

6 
julia> using BenchmarkTools, SpecialFunctions

julia> sf(k) = logabsbinomial(big"2"^64, k) |> first
sf (generic function with 1 method)

julia> b(k) = log(binomial(big"2"^64, k))
b (generic function with 1 method)

julia> lab_128(k) = logabsbinomial(int128"2"^64, k) |> first
lab_128 (generic function with 1 method)

julia> @benchmark sf(100)
BenchmarkTools.Trial: 10000 samples with 1 evaluation.
 Range (min … max):  78.879 μs …  13.882 ms  ┊ GC (min … max): 0.00% … 61.89%
 Time  (median):     80.917 μs               ┊ GC (median):    0.00%
 Time  (mean ± σ):   82.970 μs ± 138.064 μs  ┊ GC (mean ± σ):  1.04% ±  0.62%

        ▂▄▆█▇▇▅▃▁                                               
  ▁▁▂▄▄▇█████████▇▅▄▃▃▂▂▂▃▃▃▃▃▃▃▃▃▃▂▂▂▂▂▂▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁ ▃
  78.9 μs         Histogram: frequency by time         89.2 μs <

 Memory estimate: 4.05 KiB, allocs estimate: 112.

julia> @benchmark b(100)
BenchmarkTools.Trial: 10000 samples with 1 evaluation.
 Range (min … max):  17.262 μs …  5.228 ms  ┊ GC (min … max): 0.00% … 60.46%
 Time  (median):     17.745 μs              ┊ GC (median):    0.00%
 Time  (mean ± σ):   18.364 μs ± 52.122 μs  ┊ GC (mean ± σ):  1.72% ±  0.60%

          ▁▁▃▄▆▆▅█▇▆▆▅▅▃▃▂▁                                    
  ▂▂▂▃▃▄▆▇██████████████████▇▇▆▆▅▅▄▃▃▃▃▃▃▂▂▂▂▂▂▂▂▁▂▂▂▂▂▂▂▂▂▁▂ ▄
  17.3 μs         Histogram: frequency by time        18.9 μs <

 Memory estimate: 4.54 KiB, allocs estimate: 44.

julia> @benchmark lab_128(100)
BenchmarkTools.Trial: 10000 samples with 864 evaluations.
 Range (min … max):  137.354 ns … 206.201 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     137.398 ns               ┊ GC (median):    0.00%
 Time  (mean ± σ):   139.310 ns ±   4.802 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

  █        ▁                ▁                                   ▁
  █▃▄▄▄▄▄▅▅█▇▇▆▆▆▅▆▅▆▇▇▇▇█████████████▇▆▆▆▆▇▇▆▆▇▅▆▅▅▄▄▅▄▅▄▄▄▄▅▃ █
  137 ns        Histogram: log(frequency) by time        158 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

julia> sf(100)
4072.40258002808648985786554627251204975606540801181314771372293997172075937879

julia> b(100)
4072.402580028086489857865546272512049756065408011813147713715631592357374294182

julia> lab_128(100)
4072.402580028087

From the above the fastest way I have found so far is first(SpecialFunctions.logabsbinomial(int128"2"^64, 100))

7

I would expect logabsbinomial to be more accurate since the lgamma variant can have some cancellation in different parts of the domain, but this is speculation rather than informed information.

8

That’s what I figured too. Looks like

lchoose(a,b) = logabsbinomial(a,b)[1]

is probably the best plan.