factorial and Int64

[my-little-repository]

Factorial returns a wrong result for all Int64 integers bigger than 20, due to overflow.

 julia> factorial(21)
-4249290049419214848

I suggest converting the argument to BigInt if it is an Int64 bigger than 20, or at least
returning a warning that the result is wrong.

It seems that factorial convert other number types like Int32 or Int8 to Int64 if there is an overflow for these types, so that handling Int64 should be sufficient to handle all non-Bool Integer types.

[tknopp]

converting to BigInt is problematic as this would make the function type instable.

That integers are automatically converted to Int64 is because all integer arithmetic is done in the native integer type in Julia. This is something not everybody is happy with.

Whether to raise an error depends on if this would impact the performance or not. For instance 2^60 * 2^60 will not raise an error as this would have a major performance impact.

[my-little-repository]

Really, this should just amount to add a test at start of the function definition, e.g.if arg > 21 warn("you may want bigint for arguments bigger than 21 in your factorial") end and this is unlikely to impact the performance.

Actually everything would be better than returning a negative integer as the value of factorial(21).
Also, if the goal is to return a result as fast as possible, a hash of the factorial of the 21 first INT64 integers values, and 0 + warning otherwise, would probably be faster and more correct than the current implementation.

[tknopp]

In this specific case the warning/error might be really the right thing. My personal opinion is that its not worth to fix some individual cases while the majority of code will not be integer overflow secure anyway. It might be a pretty extreme example but factorial(2^64) would not warn you even when factorial contains a warning.

[ViralBShah]

Good idea to have the warning as well as precomputing the results.

[tknopp]

Note that both cases Int==Int32 and Int==Int64 have to handled when tackling this.

[StefanKarpinski]

Marked as up for grabs if someone want to take a crack at this. It's nicely isolated bit of work.

[hisq]

Maybe one should also change the definition for negative integer arguments, see http://en.wikipedia.org/wiki/Factorial#Non-extendability_to_negative_integers

The correct value for negative integers is complex infinity, for example 1/factorial(-1) should give zero. If one cannot do that, then probably factorial should give an error for negative integers.

[JeffBezanson]

There's no reason to give a warning instead of an overflow error; we know the exact argument for which it will start to overflow.

[my-little-repository]

@JeffBezanson
Then I guess this is just a matter of adding the following line at start of the definition of factorial on top of the combinatorics.jl file.

 if n>20 throw(OverFlowError) end

I am not yet skilled enough with github to submit a pull request though.

Note that this does not seem to lead to any significant performance penalty on my machine.

julia> @time curated_factorial(20)
elapsed time: 5.448e-6 seconds (64 bytes allocated)
2432902008176640000

julia> @time factorial(20)
elapsed time: 5.168e-6 seconds (64 bytes allocated)
2432902008176640000

Also precomputing the result does not seem to improve speed.

julia> let h = [20 => 2432902008176640000] ; @time h[20] ; end
elapsed time: 5.238e-6 seconds (64 bytes allocated)
2432902008176640000

[nalimilan]

@my-little-repository You need to run the function may times in a loop to get reasonable timings.

[my-little-repository]

@nalimilan I can run it as long as you want.

julia> @time for i = 1:100000000 curated_factorial(20) end
elapsed time: 2.595464475 seconds (0 bytes allocated)

julia> @time for i = 1:100000000 factorial(20) end
elapsed time: 2.594375864 seconds (0 bytes allocated)

julia> let h = Dict{Int64,Int64}([20 => 2432902008176640000])
         @time for i = 1:100000000 h[20] end
       end
elapsed time: 3.034445719 seconds (0 bytes allocated)

Sincerely,

[tknopp]

The first two timings would be more accurate if you call the function once before. In this way compilation time will not be included. But in this case the ifshould not really matter.

The last should also be put in a function. A Dict seems unnecessary here. An array will be faster.

But again, the 20 will only be for 64 bit systems. On 32 bit Julia it has to be 12 when Integer!=Int64 .

[my-little-repository]

@tknopp I assure you that the two timings were indeed preceded by a call to the function before, in order to insure compilation. I don't think that the last snippet should be put in a function; the let form instantiates a local environment, hence there is no call to global variables there. You may want to do your own benchmark to convince yourself on that point. You are indeed right that an array is faster, but the cost of the computation of the factorial is negligible with respect to the cost of the function call, and the function call is simulated by the dictionnary.

julia> let h = Int64[factorial(big(i)) for i = 1:20]
         function curated_factorial_bis (n::Int64)
            n<21? h[n]: error(OverflowError)
         end
       end
curated_factorial_bis (generic function with 1 method)

julia> let g = Int64[factorial(big(i)) for i = 1:20]
         @time for i=1:100000000  curated_factorial_bis(20) end
         @time for i=1:100000000  g[20] end
       end  
elapsed time: 4.004516008 seconds (1600000000 bytes allocated)
elapsed time: 0.041102083 seconds (0 bytes allocated)

As you see, comparing once again to the previous results, buffering values does not speed up the computation.

You have a good point concerning 32bit julia. Unfortunately my knowledge of Julia is not deep enough to suggest a fix here.

[tknopp]

The issue with curated_factorial_bis is that it is not type stable. It has either return type Int64 or Void.
The (1600000000 bytes allocated) is an indicator for these type of issues.

With

function curated_factorial_bis (n::Int64)
  if n>20 error(OverflowError) end; h[n]
end

I get about 0.5s vs 0.05s for the direct access to h. However,

function curated_factorial_bis (n::Int64)
  h[n]
end

gives me also 0.05s, which is what I expected as such a short function should be inlined.
So you can see how the performance impact can be when introducing an if statement in such a short function. I would not have thought that it is so big in this case.

Note that implementing the proposed error is not as simple as introducing the if statement in the factorial function in combinatorics.jl. It is defined for type Integer which can be any integer such as also a BigInt or an Int128. Thus, one would have to specialize the function for several integer types.