Factorial function returning negative number for large input

Question

My factorial function seems to work for numbers between 1 and 6, but not for numbers much bigger than 6, for example starting with 21! the results are negative.

I cannot figure out why. Here's my function:

factorial :: Int -> Int
factorial 0 = 1
factorial 1 = 1
factorial num = num * factorial( num - 1)

And here's my binomial coefficient function that calls my factorial function (maybe the problem comes from this one ?):

binomialCoef :: Int -> Int -> Int
binomialCoef n 1 = n
binomialCoef n k = factorial n `div` 
                        ((factorial k) * factorial (n - k))

Answer 1

(…) realized my factorial function returns negative numbers starting at 21!, and I can't figure out why.

Because an Int has a fixed number of bits. An Int should at least represent all numbers between -2^-29 and 2²⁹-1, and on a 64-bit system, typically it will represent numbers between -2^-63 and 2⁶³-1, but regardless what bounds it represents, it will eventually run out of bits to represent such number.

You can work with Integer to represent arbitrary large numbers:

factorial :: Integer -> Integer
factorial 0 = 1
factorial 1 = 1
factorial num = num * factorial (num-1)

For example:

Prelude> factorial 21
51090942171709440000
Prelude> factorial 22
1124000727777607680000

Answer 2

The binomial coefficient is where ln(gamma) really shines:

Bi(n, k) = n!/(k!*(n-k)!)

Taking the natural log of both sides:

ln(Bi(n, k)) = ln(n!) - ln(k!) - ln((n-k)!)

But

gamma(n) = (n-1)!

Or

gamma(n+1) = n!

Substituting

ln(Bi(n, k)) = lngamma(n+1) - lngamma(k+1) -lngamma(n-k+1)

Taking the exponential of both sides gives the final result:

Bi(n, k) = exp(lngamma(n+1) - lngamma(k+1) - lngamma(n-k+1))

There's a Haskell implementation. I haven't looked at it, but it should return a Double instead of an Integer. You won't have overflow problems because of that fact. It'll also be better behaved because you will be subtracting logarithms instead of dividing a large numerator by a large product in the denominator.

Answer 3

Of course best way to avoid integer overflow and wrap-around while calculating a big factorial is not to calculate the factorial in the first place. Instead, since

factorial n = product [1..n]

keeping [1..n] as the representation of the factorial of n is as good -- or even much better -- as calculating the actual number. Postponing an action until absolutely unavoidable we get to pre-optimize it before post-calculating:

bincoef :: Int -> Int -> Int
bincoef n k = factorial n `div` 
                        ((factorial k) * factorial (n - k))
 = product [1 .. n] `div` 
        (product [1 .. k] * product [1 .. n-k])
 = product [n-k+1 .. n] `div` 
         product [1 .. k]
 = foldl' g 1 $ zip [n, n-1 .. n-k+1] [1 .. k]
     where g !acc (a,b) = (acc * a) `div` b

So now,

> mapM_ (\n -> print $ map (bincoef n) [5,10..n]) [20,30..60]
[15504,184756,15504,1]
[142506,30045015,155117520,30045015,142506,1]
[658008,847660528,40225345056,137846528820,40225345056,847660528,658008,1]
[2118760,10272278170,2250829575120,47129212243960,126410606437752,47129212243960,
2250829575120,10272278170,2118760,1]
[5461512,75394027566,53194089192720,4191844505805495,51915437974328292,1182645815
64861424,51915437974328292,4191844505805495,53194089192720,75394027566,5461512,1]

> head . filter (not . snd) $ map (\n -> (n, all (> 0) $ map (bincoef n) [1..n])) [1..]
(62,False)

the Int wrap-around error makes its first appearance at n=62. But it's still working at n=60, and we can see there are more than 16 digits in those numbers, so no Double-based calculation has a hope of working correctly, there.

To get into yet higher ranges still with the Int-based operations only, the next logical step is keeping the lists of integers as originally proposed, or better yet as their prime factorizations which are easy to multiply and divide; but at that point we'd be getting pretty close to re-implementing the bignum arithmetic ourselves, so might as well just use the simple Integer-based code,

bc :: Integer -> Integer -> Integer
bc n k = product [n-k+1 .. n] `div` product [1 .. k]

which "just works".

> bc 600 199
124988418115780688528958442419612410733294315465732363826979722360319899409241320138
666379143574138790334901309769571503484430553926248548697640619977793300443439200