Haskell中Fisher精确检验的高效实现方案求助

嘿，我懂你在Haskell里实现Fisher精确检验时想要更高效计算的需求——直接算大阶乘不仅慢，还容易碰到数值溢出的问题。结合你给出的公式，我来分享几个更优的实现思路，避开直接硬算阶乘的坑：

直接计算阶乘的乘积很容易因为数值过大导致性能下降甚至溢出，我们可以用对数把乘法转换成加法，最后再通过指数运算还原结果。这种方法计算速度快，而且数值稳定性更好，适合不需要精确整数结果的场景：

import Data.List (foldl')
-- 可以使用number-theory库中的logFactorial，或者自己实现一个
import Math.NumberTheory.Logarithms (logFactorial)

fisherExactLog :: Integer -> Integer -> Integer -> Integer -> Double
fisherExactLog a b c d = exp $ 
  logFactorial (a+b) + logFactorial (a+c) + logFactorial (b+d) + logFactorial (c+d)
  - logFactorial a - logFactorial b - logFactorial c - logFactorial d
  - logFactorial (a+b+c+d)

如果不想依赖外部库，也可以自己实现对数阶乘的计算（用伽马函数或者递推求和），不过number-theory库的实现已经足够高效了。

如果需要精确的整数结果，而且输入的自然数可能很大，质因数分解是个好选择。我们先计算每个阶乘的质因数分解，然后把分子的质因数指数相加、分母的指数相加，再用分子的指数减去分母的指数，最后计算剩余质因数的乘积。这样完全避免了大整数的直接相乘，效率更高：

import Data.List (groupBy, sort)
import Data.Function (on)

-- 计算n!的质因数分解，返回[(质数, 指数)]
factorialPrimeFactors :: Integer -> [(Integer, Integer)]
factorialPrimeFactors n = concatMap countPrimeExponents $ primesUpTo n
  where
    countPrimeExponents p = 
      let total = sum [n `div` p^k | k <- [1..], p^k <= n]
      in (p, total)
    -- 埃氏筛法生成n以内的质数
    primesUpTo m = sieve [2..m]
    sieve [] = []
    sieve (x:xs) = x : sieve [y | y <- xs, y `mod` x /= 0]

-- 合并两个质因数分解列表，指数相加
mergeAdd :: [(Integer, Integer)] -> [(Integer, Integer)] -> [(Integer, Integer)]
mergeAdd a b = merge (sort a) (sort b)
  where
    merge [] ys = ys
    merge xs [] = xs
    merge ((p1,e1):xs) ((p2,e2):ys)
      | p1 == p2 = (p1, e1+e2) : merge xs ys
      | p1 < p2 = (p1,e1) : merge xs ((p2,e2):ys)
      | otherwise = (p2,e2) : merge ((p1,e1):xs) ys

-- 从第一个质因数分解中减去第二个的指数（只保留正指数）
mergeSubtract :: [(Integer, Integer)] -> [(Integer, Integer)] -> [(Integer, Integer)]
mergeSubtract a b = subtractFactors (sort a) (sort b)
  where
    subtractFactors [] _ = []
    subtractFactors xs [] = xs
    subtractFactors ((p1,e1):xs) ((p2,e2):ys)
      | p1 == p2 = if e1 > e2 then (p1, e1-e2) : subtractFactors xs ys else subtractFactors xs ys
      | p1 < p2 = (p1,e1) : subtractFactors xs ((p2,e2):ys)
      | otherwise = subtractFactors ((p1,e1):xs) ys

-- 计算质因数分解的乘积
computeProduct :: [(Integer, Integer)] -> Integer
computeProduct = foldl' (\acc (p,e) -> acc * p^e) 1

fisherExactPrimeFactor :: Integer -> Integer -> Integer -> Integer -> Integer
fisherExactPrimeFactor a b c d = computeProduct $ 
  numeratorFactors `mergeSubtract` denominatorFactors
  where
    numeratorFactors = foldl' mergeAdd [] 
      [factorialPrimeFactors (a+b), factorialPrimeFactors (a+c), 
       factorialPrimeFactors (b+d), factorialPrimeFactors (c+d)]
    denominatorFactors = foldl' mergeAdd []
      [factorialPrimeFactors a, factorialPrimeFactors b, factorialPrimeFactors c, 
       factorialPrimeFactors d, factorialPrimeFactors (a+b+c+d)]

观察你的公式可以发现，它可以转换成组合数的比值：

p = C(a+b, a) × C(c+d, c) / C(a+b+c+d, a+c)

其中C(n,k)是组合数“n选k”。组合数可以用递推乘法的方式高效计算，每一步都保持整数，而且避免了大阶乘的直接计算：

-- 高效计算组合数C(n,k)，优先计算较小的k减少运算量
comb :: Integer -> Integer -> Integer
comb n k
  | k < 0 || k > n = 0
  | k == 0 || k == n = 1
  | k > n - k = comb n (n - k) -- 取较小的k，减少计算步骤
  | otherwise = foldl' (\acc i -> acc * (n - k + i) `div` i) 1 [1..k]

-- 返回Rational类型，保持精确的分数结果
fisherExactComb :: Integer -> Integer -> Integer -> Integer -> Rational
fisherExactComb a b c d = 
  let total = a + b + c + d
      c1 = comb (a+b) a
      c2 = comb (c+d) c
      c3 = comb total (a+c)
  in fromIntegral (c1 * c2) / fromIntegral c3

这个方案代码最简洁，而且计算效率很高，返回的Rational类型可以完全保留精确值，适合大多数常规场景。

内容的提问来源于stack exchange，提问作者Sebi Danila