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