如何在Haskell中提高这个数值计算的性能？
我正在把David Blei的C潜在Dirichlet分配的C实现移植到Haskell中，我正在试着决定是否在C中留下一些低级的东西。下面的函数就是一个例子 – 它是一个近似值的二阶导数：
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); }
我已经把它翻译成或多或less俗语的Haskell，如下所示：
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 (\ab -> (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编译）。 一些类似的function更糟。
我对Haskell的性能几乎一无所知，而且我对通过Core或类似的东西进行挖掘并不感兴趣，因为我总是可以通过FFI调用一些math密集的C函数。
但是我很好奇我是否缺less一些低级的function – 某种扩展或图书馆或注释，我可以使用这些function来加速这些数字，而不会使它变得太难看。
更新：感谢Don Stewart和Yitz ，这里有两个更好的解决scheme。 我稍微修改了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
两者的performance似乎几乎完全相同，根据编译器标志，其中一个或另一个赢得一个或两个百分点。
正如Camccann 在Reddit上所说，故事的寓意是“为了获得最佳结果，请使用Don Stewart作为您的GHC后端代码生成器”。 除了这个解决scheme之外，最安全的方法似乎只是把C控制结构直接转换成Haskell，尽pipe循环融合可以以更习惯的方式提供类似的性能。
我可能会最终在我的代码中使用Data.Vector方法。
使用相同的控制和数据结构，产生：
{-# LANGUAGE BangPatterns #-} {-# OPTIONS_GHC -fvia-C -optc-O3 -fexcess-precision -optc-march=native #-} {-# INLINE trigamma #-} trigamma :: Double -> Double trigamma x = go 0 (x' - 1) p' where 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 go :: Int -> Double -> Double -> Double go !i !x !p | i >= 6 = p | otherwise = go (i+1) (x-1) (1 / (x*x) + p)
我没有你的testing套件，但这产生了以下asm：
A_zdwgo_info: cmpq $5, %r14 jg .L3 movsd .LC0(%rip), %xmm7 movapd %xmm5, %xmm8 movapd %xmm7, %xmm9 mulsd %xmm5, %xmm8 leaq 1(%r14), %r14 divsd %xmm8, %xmm9 subsd %xmm7, %xmm5 addsd %xmm9, %xmm6 jmp A_zdwgo_info
看起来没问题 这是-fllvm后端做得很好的代码。
虽然GCC展开循环，唯一的方法是通过模板Haskell或手动展开。 你可能会认为（一个THmacros）如果做了这么多。
实际上，GHC LLVM后端会展开循环:-)
最后，如果你真的喜欢原始的Haskell版本，可以使用stream融合组合器来编写它，而GHC会把它转换成循环。 （为读者锻炼）。
在优化工作之前，我不会说你的原始翻译是在Haskell中用C代码expression的最习惯的方式。
如果我们从下面开始，优化过程如何进行？
trigamma :: Double -> Double trigamma x = foldl' (+) p' . map invSq . take 6 . iterate (+ 1) $ x where invSq y = 1 / (y * y) x' = x + 6 p = invSq x' p' =(((((0.075757575757576*p-0.033333333333333)*p+0.0238095238095238) *p-0.033333333333333)*p+0.166666666666667)*p+1)/x'+0.5*p
