代码之家 › 专栏 › 技术社区 › Travis Brown

如何提高哈斯克尔数值计算的性能?

haskell math performance c

Travis Brown · 技术社区 · 15 年前

我正在移植大卫布莱的原作 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);
}

我已经将其翻译成或多或少惯用的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 (\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版本慢了六到七倍(我正在编译 -O2 在GHC 6.12.1)上)。一些类似的功能甚至更糟。

我对哈斯克尔的表演几乎一无所知,而且我对 digging through Core 或者类似的事情,因为我可以通过ffi调用少量数学密集型C函数。

但我很好奇,是否有一些低挂水果,我错过了某种扩展、库或注释,我可以用它们来加快这些数字的速度,而不会使它变得太难看。

更新: 这有两个更好的解决方案,多亏了 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

两人的表现似乎几乎完全相同,其中一人或另一人赢了一个或两个百分点,这取决于编译器标志。

AS camccann 说 over at Reddit 这个故事的寓意是“为了达到最佳效果,使用Don Stewart作为您的GHC后端代码生成器。”除非有这个解决方案,否则最安全的选择似乎只是将C控制结构直接转换为Haskell,尽管循环融合可以以更惯用的方式提供类似的性能。

我可能最终会使用 矢量数据 我的代码中的方法。

2 回复 | 直到 15 年前

Bakuriu 11 年前

使用相同的控制和数据结构,产生:

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

我没有您的测试套件,但这会产生以下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或手动展开。如果你做了很多这样的事情,你可能会认为这是(一个th宏)。

实际上,ghc llvm后端确实展开了循环:—)

最后,如果你真的喜欢最初的haskell版本,用 stream fusion combinators, GHC将把它转换回循环。(为读者做练习)。

Yitz 15 年前

在优化工作之前,我不会说您的原始翻译是用haskell表达C代码所做工作的最惯用方法。

如果我们改为从以下步骤开始,优化过程将如何进行:

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