大输入的阶乘函数返回负数

Question

我的阶乘函数似乎适用于 1 到 6 之间的数字，但不适用于大于 6 的数字，例如从 21 开始！结果是否定的。

我不知道为什么。这是我的功能：

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

这是我的二项式系数函数，它调用了我的阶乘函数（也许问题来自这个？）：

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

Answer 1

(...) 意识到我的阶乘函数从 21 开始返回负数！，我不知道为什么。

因为Int 具有固定的位数。 Int 至少应该代表 -2^-29 和 2²⁹-1 之间的所有数字，并且在 64 位系统上，它通常代表介于 - 2^-63 和 2⁶³-1，但不管它代表什么边界，最终都会用完代表这个数字的位。

您可以使用Integer 来表示任意大数：

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

例如：

Prelude> factorial 21
51090942171709440000
Prelude> factorial 22
1124000727777607680000

Answer 2

二项式系数是ln(gamma)真正闪耀的地方：

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

取两边的自然对数：

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

但是

gamma(n) = (n-1)!

或者

gamma(n+1) = n!

替换

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

取两边的指数得到最终结果：

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

有一个Haskell implementation。我没有看过它，但它应该返回一个 Double 而不是 Integer。因为这个事实，你不会有溢出问题。它的表现也会更好，因为您将减去对数，而不是用大分子除以分母中的大乘积。

Answer 3

当然，在计算大阶乘时避免整数溢出和环绕的最佳方法是不首先计算阶乘。相反，因为

factorial n = product [1..n]

保持[1..n] 作为n 的阶乘表示与计算实际数字一样好——甚至更好。 推迟一个动作，直到绝对不可避免，我们可以在后计算之前对其进行预优化：

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

所以现在，

> 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)

Int 环绕错误首次出现在 n=62。但是它仍然在 n=60 工作，我们可以看到这些数字中有超过 16 个数字，所以没有基于 Double 的计算有希望正常工作，那里。

要进入更高的范围，仍然只使用基于Int 的操作，下一个合乎逻辑的步骤是保持整数的列表与最初提出的一样，或者更好的是作为它们的素数分解，它们是易于乘除；但那时我们已经非常接近于自己重新实现 bignum 算法了，所以不妨使用简单的基于Integer 的代码，

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

哪个“有效”。

> bc 600 199
124988418115780688528958442419612410733294315465732363826979722360319899409241320138
666379143574138790334901309769571503484430553926248548697640619977793300443439200