How is factorial(BigInt) so fast?


#1

As an experiment, I copy/pasted Julia’s factorial() function and renamed it. But when I use it to evaluate a bigint, it is way slower than Julia’s factorial(BigInt). If the base function uses the same code, then why is the base function faster and is there a way I can achieve those speeds?

I’m guessing Julia’s either uses a dictionary lookup or somehow avoids most of the performance loss of bigints. The latter would be a useful trick to know.

I understand that creating a bigint longer than 9 digits is about 10x slower than creating a smaller bigint, and that another slowdown happens at 20 digits. Therefore, hypothetically if you evaluate BigInt(15)! manually, it’s faster to do

BigInt(15*14*13*12*11) * BigInt(10*9*8*7*6*5*4*3*2)

than it is to do

BigInt(15*14*13*12*11*10*9*8*7*6) * BigInt(5*4*3*2)

because (15 P 5) and 10! are both under 10 digits whereas (15 P 10) is not.

In fact, I know of no faster way to manually evaluate BigInt(15)! than to split it into two BigInts where each is smaller than 1 billion. It’s faster than BigInt(1514…*2), or BigInt(15)BigInt(14)…*BigInt(2), or any conceivable possibility.

However, even that is twice as slow as factorial(BigInt(15)) according to @btime. Julia has a trick up its sleeve and I want to know about it :slight_smile:

Thanks


#2

You copied the wrong function

julia> @which factorial(BigInt(1))
factorial(x::BigInt) in Base.GMP at gmp.jl:564

julia> @which factorial(1)
factorial(n::Union{Int64, UInt64}) in Base at combinatorics.jl:27

The correct one seems to be implemented by GMP.
You might also have luck with adapting the factorial code to be in-place (BigInts are mutable I think).


#3

Why are there multiple versions of factorial, and is this an ambiguity that happens more often, or does it only happen with factorial?


#4

It’s the strength of Julia by dispatching in type. The algorithm for factorial (or any other function) may be optimized for different types.


#5

Try methods(*) in the REPL.


#6

Right, I’m aware of the multiple dispatch of Julia, I use it all the time. But normally, if a function is extended, it was extended globally, but maybe the behavior is different in 1.0, so it should still have all the same methods available if it is reassigned to a new name, shouldnt it?


#7

I don’t know what extended globally means but factorial is no different from *

julia> methods(factorial)
# 7 methods for generic function "factorial":
[1] factorial(n::UInt128) in Base at combinatorics.jl:26
[2] factorial(n::Int128) in Base at combinatorics.jl:25
[3] factorial(x::BigFloat) in Base.MPFR at mpfr.jl:594
[4] factorial(x::BigInt) in Base.GMP at gmp.jl:564
[5] factorial(n::Union{Int64, UInt64}) in Base at combinatorics.jl:27
[6] factorial(n::Union{Int16, Int32, Int8, UInt16, UInt32, UInt8}) in Base at combinatorics.jl:33
[7] factorial(n::Integer) in Base at intfuncs.jl:820

#8

Well, if I do g = factorial and take methods of that, then all the definitions are retained. So, it behaves as I expected. I figured thats what OP meant by renaming the function.


#9

No, they copied the source code of one method.


#10

.


#11

Mystery solved, thanks!