C++ 如何在 Haskell 中提高这种数值计算的性能

我正把David Blei最初的C implementation个潜在的Dirichlet分配移植到Haskell，并且我试图决定是否在C.留下一些低级的东西.下面的函数是一个例子，它是lgamma的二阶导数的近似值:

double trigamma(double x)
{
    double p;
    int i;

    x=x+6;
    p=1/(x*x);
    p=(((((0.075757575757576*p-0.033333333333333)*p+0.0238095238095238)
         *p-0.033333333333333)*p+0.166666666666667)*p+1)/x+0.5*p;
    for (i=0; i<6 ;i++)
    {
        x=x-1;
        p=1/(x*x)+p;
    }
    return(p);
}

我把这句话翻译成或多或少习惯用的哈斯克尔，如下所示:

trigamma :: Double -> Double
trigamma x = snd $ last $ take 7 $ iterate next (x' - 1, p')
  where
    x' = x + 6
    p  = 1 / x' ^ 2
    p' = p / 2 + c / x'
    c  = foldr1 (\a b -> (a + b * p)) [1, 1/6, -1/30, 1/42, -1/30, 5/66]
    next (x, p) = (x - 1, 1 / x ^ 2 + p)

问题是，当我把这两个都运行到Criterion时，我的Haskell版本会慢六到七倍(我在GHC 6.12.1上使用-O2进行编译).有些类似的功能甚至更糟.

我对Haskell的性能几乎一无所知，我对digging through Core或类似的东西也不太感兴趣，因为我总是可以通过FFI调用少数数学密集型C函数.

但我很好奇是否有什么容易摘到的果子我错过了--某种扩展、库或注释，我可以用它们来加快数值运算的速度，同时又不会把它弄得太难看.

UPDATE:多亏了Don Stewart和Yitz，这里有两个更好的解决方案.我稍微修改了Yitz的答案，使用Data.Vector.

invSq x = 1 / (x * x)
computeP x = (((((5/66*p-1/30)*p+1/42)*p-1/30)*p+1/6)*p+1)/x+0.5*p
  where p = invSq x

trigamma_d :: Double -> Double
trigamma_d x = go 0 (x + 5) $ computeP $ x + 6
  where
    go :: Int -> Double -> Double -> Double
    go !i !x !p
        | i >= 6    = p
        | otherwise = go (i+1) (x-1) (1 / (x*x) + p)

trigamma_y :: Double -> Double
trigamma_y x = V.foldl' (+) (computeP $ x + 6) $ V.map invSq $ V.enumFromN x 6

两者的性能似乎几乎完全相同，一方或另一方以一两个百分点的优势获胜，这取决于编译器标志.

正如camccann说over at Reddit，这个故事的寓意是"为了获得最好的结果，请使用Don Stewart作为您的GHC后端代码生成器."除了该解决方案之外，最安全的 Select 似乎是直接将C控制 struct 转换为Haskell，尽管循环融合可以以更惯用的方式提供类似的性能.

我可能最终会在代码中使用Data.Vector方法.