Data.Vector.Unboxed.Mutable.MVector 的索引真的这么慢吗？答案

【问题标题】：Is indexing of Data.Vector.Unboxed.Mutable.MVector really this slow?Data.Vector.Unboxed.Mutable.MVector 的索引真的这么慢吗？
【发布时间】：2012-03-13 15:45:41
【问题描述】：

我有一个应用程序花费大约 80% 的时间使用 Kahan summation algorithm 计算高维向量 (dim=100) 的大型列表 (10^7) 的质心。我已经尽我所能优化求和，但它仍然比等效的 C 实现慢 20 倍。分析表明罪魁祸首是来自Data.Vector.Unboxed.Mutable 的unsafeRead 和unsafeWrite 函数。我的问题是：这些功能真的这么慢还是我误解了分析统计数据？

这是两个实现。 Haskell 是使用 llvm 后端使用 ghc-7.0.3 编译的。 C 是用 llvm-gcc 编译的。

Haskell 中的 Kahan 求和：

{-# LANGUAGE BangPatterns #-}
module Test where

import Control.Monad ( mapM_ )
import Data.Vector.Unboxed ( Vector, Unbox )
import Data.Vector.Unboxed.Mutable ( MVector )
import qualified Data.Vector.Unboxed as U
import qualified Data.Vector.Unboxed.Mutable as UM
import Data.Word ( Word )
import Data.Bits ( shiftL, shiftR, xor )

prng :: Word -> Word
prng w = w' where
    !w1 = w  `xor` (w  `shiftL` 13)
    !w2 = w1 `xor` (w1 `shiftR` 7)
    !w' = w2 `xor` (w2 `shiftL` 17)

mkVect :: Word -> Vector Double
mkVect = U.force . U.map fromIntegral . U.fromList . take 100 . iterate prng

foldV :: (Unbox a, Unbox b) 
      => (a -> b -> a) -- componentwise function to fold
      -> Vector a      -- initial accumulator value
      -> [Vector b]    -- data vectors
      -> Vector a      -- final accumulator value
foldV fn accum vs = U.modify (\x -> mapM_ (liftV fn x) vs) accum where
    liftV f acc = fV where
        fV v = go 0 where
            n = min (U.length v) (UM.length acc)
            go i | i < n     = step >> go (i + 1)
                 | otherwise = return ()
                 where
                     step = {-# SCC "fV_step" #-} do
                         a <- {-# SCC "fV_read"  #-} UM.unsafeRead acc i
                         b <- {-# SCC "fV_index" #-} U.unsafeIndexM v i
                         {-# SCC "fV_write" #-} UM.unsafeWrite acc i $! {-# SCC "fV_apply" #-} f a b

kahan :: [Vector Double] -> Vector Double
kahan [] = U.singleton 0.0
kahan (v:vs) = fst . U.unzip $ foldV kahanStep acc vs where
    acc = U.map (\z -> (z, 0.0)) v

kahanStep :: (Double, Double) -> Double -> (Double, Double)
kahanStep (s, c) x = (s', c') where
    !y  = x - c
    !s' = s + y
    !c' = (s' - s) - y
{-# NOINLINE kahanStep #-}

zero :: U.Vector Double
zero = U.replicate 100 0.0

myLoop n = kahan $ map mkVect [1..n]

main = print $ myLoop 100000

使用 llvm 后端编译 ghc-7.0.3：

ghc -o Test_hs --make -fforce-recomp -O3 -fllvm -optlo-O3 -msse2 -main-is Test.main Test.hs

time ./Test_hs
real    0m1.948s
user    0m1.936s
sys     0m0.008s

分析信息：

16,710,594,992 bytes allocated in the heap
      33,047,064 bytes copied during GC
          35,464 bytes maximum residency (1 sample(s))
          23,888 bytes maximum slop
               1 MB total memory in use (0 MB lost due to fragmentation)

  Generation 0: 31907 collections,     0 parallel,  0.28s,  0.27s elapsed
  Generation 1:     1 collections,     0 parallel,  0.00s,  0.00s elapsed

  INIT  time    0.00s  (  0.00s elapsed)
  MUT   time   24.73s  ( 24.74s elapsed)
  GC    time    0.28s  (  0.27s elapsed)
  RP    time    0.00s  (  0.00s elapsed)
  PROF  time    0.00s  (  0.00s elapsed)
  EXIT  time    0.00s  (  0.00s elapsed)
  Total time   25.01s  ( 25.02s elapsed)

  %GC time       1.1%  (1.1% elapsed)

  Alloc rate    675,607,179 bytes per MUT second

  Productivity  98.9% of total user, 98.9% of total elapsed

    Thu Feb 23 02:42 2012 Time and Allocation Profiling Report  (Final)

       Test_hs +RTS -s -p -RTS

    total time  =       24.60 secs   (1230 ticks @ 20 ms)
    total alloc = 8,608,188,392 bytes  (excludes profiling overheads)

COST CENTRE                    MODULE               %time %alloc

fV_write                       Test                  31.1   26.0
fV_read                        Test                  27.2   23.2
mkVect                         Test                  12.3   27.2
fV_step                        Test                  11.7    0.0
foldV                          Test                   5.9    5.7
fV_index                       Test                   5.2    9.3
kahanStep                      Test                   3.3    6.5
prng                           Test                   2.2    1.8


                                                                                               individual    inherited
COST CENTRE              MODULE                                               no.    entries  %time %alloc   %time %alloc

MAIN                     MAIN                                                   1           0   0.0    0.0   100.0  100.0
 CAF:main1               Test                                                 339           1   0.0    0.0     0.0    0.0
  main                   Test                                                 346           1   0.0    0.0     0.0    0.0
 CAF:main2               Test                                                 338           1   0.0    0.0   100.0  100.0
  main                   Test                                                 347           0   0.0    0.0   100.0  100.0
   myLoop                Test                                                 348           1   0.2    0.2   100.0  100.0
    mkVect               Test                                                 350      400000  12.3   27.2    14.5   29.0
     prng                Test                                                 351     9900000   2.2    1.8     2.2    1.8
    kahan                Test                                                 349         102   0.0    0.0    85.4   70.7
     foldV               Test                                                 359           1   5.9    5.7    85.4   70.7
      fV_step            Test                                                 360     9999900  11.7    0.0    79.5   65.1
       fV_write          Test                                                 367    19999800  31.1   26.0    35.4   32.5
        fV_apply         Test                                                 368     9999900   1.0    0.0     4.3    6.5
         kahanStep       Test                                                 369     9999900   3.3    6.5     3.3    6.5
       fV_index          Test                                                 366     9999900   5.2    9.3     5.2    9.3
       fV_read           Test                                                 361     9999900  27.2   23.2    27.2   23.2
 CAF:lvl19_r3ei          Test                                                 337           1   0.0    0.0     0.0    0.0
  kahan                  Test                                                 358           0   0.0    0.0     0.0    0.0
 CAF:poly_$dPrimMonad3_r3eg Test                                                 336           1   0.0    0.0     0.0    0.0
  kahan                  Test                                                 357           0   0.0    0.0     0.0    0.0
 CAF:$dMVector2_r3ee     Test                                                 335           1   0.0    0.0     0.0    0.0
 CAF:$dVector1_r3ec      Test                                                 334           1   0.0    0.0     0.0    0.0
 CAF:poly_$dMonad_r3ea   Test                                                 333           1   0.0    0.0     0.0    0.0
 CAF:$dMVector1_r3e2     Test                                                 330           1   0.0    0.0     0.0    0.0
 CAF:poly_$dPrimMonad2_r3e0 Test                                                 328           1   0.0    0.0     0.0    0.0
  foldV                  Test                                                 365           0   0.0    0.0     0.0    0.0
 CAF:lvl11_r3dM          Test                                                 322           1   0.0    0.0     0.0    0.0
  kahan                  Test                                                 354           0   0.0    0.0     0.0    0.0
 CAF:lvl10_r3dK          Test                                                 321           1   0.0    0.0     0.0    0.0
  kahan                  Test                                                 355           0   0.0    0.0     0.0    0.0
 CAF:$dMVector_r3dI      Test                                                 320           1   0.0    0.0     0.0    0.0
  kahan                  Test                                                 356           0   0.0    0.0     0.0    0.0
 CAF                     GHC.Float                                            297           1   0.0    0.0     0.0    0.0
 CAF                     GHC.IO.Handle.FD                                     256           2   0.0    0.0     0.0    0.0
 CAF                     GHC.IO.Encoding.Iconv                                214           2   0.0    0.0     0.0    0.0
 CAF                     GHC.Conc.Signal                                      211           1   0.0    0.0     0.0    0.0
 CAF                     Data.Vector.Generic                                  182           1   0.0    0.0     0.0    0.0
 CAF                     Data.Vector.Unboxed                                  174           2   0.0    0.0     0.0    0.0

C 中的等效实现：

#include <stdint.h>
#include <stdio.h>


#define VDIM    100
#define VNUM    100000



uint64_t prng (uint64_t w) {
    w ^= w << 13;
    w ^= w >> 7;
    w ^= w << 17;
    return w;
};

void kahanStep (double *s, double *c, double x) {
    double y, t;
    y  = x - *c;
    t  = *s + y;
    *c = (t - *s) - y;
    *s = t;
}

void kahan(double s[], double c[]) {
    for (int i = 1; i <= VNUM; i++) {
        uint64_t w = i;
        for (int j = 0; j < VDIM; j++) {
                kahanStep(&s[j], &c[j], w);
                w = prng(w);
        }
    }
};


int main (int argc, char* argv[]) {
    double acc[VDIM], err[VDIM];
    for (int i = 0; i < VDIM; i++) {
        acc[i] = err[i] = 0.0;
    };
    kahan(acc, err);
    printf("[ ");
    for (int i = 0; i < VDIM; i++) {
        printf("%g ", acc[i]);
    };
    printf("]\n");
};

使用 llvm-gcc 编译：

>llvm-gcc -o Test_c -O3 -msse2 -std=c99 test.c

>time ./Test_c
real    0m0.096s
user    0m0.088s
sys     0m0.004s

更新 1： 我在 C 版本中取消了内联 kahanStep。它几乎没有对性能造成影响。我希望现在我们都能承认阿姆达尔定律并继续前进。作为 kahanStep 可能效率低下，unsafeRead 和 unsafeWrite 慢 9-10 倍。我希望有人可以阐明这一事实的可能原因。

另外，我应该说，因为我正在与一个使用 Data.Vector.Unboxed 的库进行交互，所以我在这一点上有点嫁给它，和它分开会很痛苦:-)

更新 2： 我想我最初的问题不够清楚。我不是在寻找加速这个微基准测试的方法。我正在寻找反直观分析统计信息的解释，因此我可以决定是否针对vector 提交错误报告。

【问题讨论】：

这两个实现绝不是等价的。另外，为什么kahanStep NOINLINE'd？内联 kahanStep、foldV 和 mkVect 弥补了一些差异，但对我来说它仍然比 C 版本慢很多。
SO 发布带有 GHC 分析图 (PNG) stackoverflow.com/questions/5939630/…
kahanStep 是 NOINLINED，否则它不会出现在我的分析信息中。如果我内联它，我会从比 C 慢 20 倍到 19.5 倍。那不是问题。如果查看分析信息，您可以看到 fV_read 和 fV_write 成本中心各占大约 30% 的时间，而 kahaneStep 仅占 3.3%。这就是问题所在：实际计算只是运行时间的一小部分。
kahanStep 本身不会占用太多时间，但是当它是 NOINLINE 时，元组和双精度会被装箱和拆箱，这会增加大量开销。我相信这个成本可能归因于读/写成本中心。 IME unsafeRead 等与 C 数组索引相当。
@John 是您在分析中看到的吗？因为那不是我所看到的。你能发布你的Test.prof吗？

标签： performance haskell vector profiling floating-point

【解决方案1】：

您的 C 版本不等同于您的 Haskell 实现。在 C 中，您自己内联了重要的 Kahan 求和步骤，在 Haskell 中，您创建了一个多态高阶函数，它可以做更多事情并将转换步骤作为参数。将 kahanStep 移动到 C 中的单独函数不是重点，编译器仍将内联它。即使您将其放入自己的源文件中，单独编译并链接而不进行链接时优化，您也只能解决部分差异。

我做了一个更接近 Haskell 版本的 C 版本，

kahan.h:

typedef struct DPair_st {
    double fst, snd;
    } DPair;

DPair kahanStep(DPair pr, double x);

kahanStep.c:

#include "kahan.h"

DPair kahanStep (DPair pr, double x) {
    double y, t;
    y  = x - pr.snd;
    t  = pr.fst + y;
    pr.snd = (t - pr.fst) - y;
    pr.fst = t;
    return pr;
}

main.c:

#include <stdint.h>
#include <stdio.h>
#include "kahan.h"


#define VDIM    100
#define VNUM    100000

uint64_t prng (uint64_t w) {
    w ^= w << 13;
    w ^= w >> 7;
    w ^= w << 17;
    return w;
};

void kahan(double s[], double c[], DPair (*fun)(DPair,double)) {
    for (int i = 1; i <= VNUM; i++) {
        uint64_t w = i;
        for (int j = 0; j < VDIM; j++) {
            DPair pr;
            pr.fst = s[j];
            pr.snd = c[j];
            pr = fun(pr,w);
            s[j] = pr.fst;
            c[j] = pr.snd;
            w = prng(w);
        }
    }
};


int main (int argc, char* argv[]) {
    double acc[VDIM], err[VDIM];
    for (int i = 0; i < VDIM; i++) {
        acc[i] = err[i] = 0.0;
    };
    kahan(acc, err,kahanStep);
    printf("[ ");
    for (int i = 0; i < VDIM; i++) {
        printf("%g ", acc[i]);
    };
    printf("]\n");
};

单独编译和链接，运行速度比这里的第一个 C 版本慢约 25%（0.1 秒对 0.079 秒）。

现在您在 C 中有了一个高阶函数，比原来的要慢得多，但仍然比 Haskell 代码快得多。一个重要的区别是 C 函数采用一对未装箱的 doubles 和一个未装箱的 double 作为参数，而 Haskell kahanStep 采用一对装箱的 Doubles 和一个装箱的 Double 并返回一对盒装的Doubles，需要在foldV 循环中进行昂贵的装箱和拆箱。这可以通过更多的内联来解决。使用 ghc-7.0.4 显式内联 foldV、kahanStep 和 step 使这里的时间从 0.90s 下降到 0.74s（它对 ghc-7.4.1 的输出影响较小，从 0.99s 下降到0.90 秒）。

但是，装箱和拆箱是差异的较小部分。 foldV 比 C 的 kahan 做得更多，它需要一个 向量列表 用于修改累加器。 C 代码中完全没有向量列表，这有很大的不同。所有这 100000 个向量都必须被分配、填充并放入一个列表中（由于懒惰，并非所有这些向量都同时处于活动状态，因此没有空间问题，但是它们以及列表单元格都必须被分配和垃圾收集，这需要相当长的时间）。在循环中，不是将Word# 传递到寄存器中，而是从向量中读取预先计算的值。

如果你使用 C 到 Haskell 的更直接的翻译，

{-# LANGUAGE CPP, BangPatterns #-}
module Main (main) where

#define VDIM 100
#define VNUM 100000

import Data.Array.Base
import Data.Array.ST
import Data.Array.Unboxed
import Control.Monad.ST
import GHC.Word
import Control.Monad
import Data.Bits

prng :: Word -> Word
prng w = w'
  where
    !w1 = w `xor` (w `shiftL` 13)
    !w2 = w1 `xor` (w1 `shiftR` 7)
    !w' = w2 `xor` (w2 `shiftL` 17)

type Vec s = STUArray s Int Double

kahan :: Vec s -> Vec s -> ST s ()
kahan s c = do
    let inner w j
            | j < VDIM  = do
                !cj <- unsafeRead c j
                !sj <- unsafeRead s j
                let !y = fromIntegral w - cj
                    !t = sj + y
                    !w' = prng w
                unsafeWrite c j ((t-sj)-y)
                unsafeWrite s j t
                inner w' (j+1)
            | otherwise = return ()
    forM_ [1 .. VNUM] $ \i -> inner (fromIntegral i) 0

calc :: ST s (Vec s)
calc = do
    s <- newArray (0,VDIM-1) 0
    c <- newArray (0,VDIM-1) 0
    kahan s c
    return s

main :: IO ()
main = print . elems $ runSTUArray calc

它要快得多。诚然，它仍然比 C 慢 3 倍左右，但原版在这里慢了 13 倍（而且我没有安装 llvm，所以我使用 vanilla gcc 和 GHC 的本机支持，使用 llvm 可能会产生略微不同的结果）。

我不认为索引真的是罪魁祸首。 vector 包严重依赖于编译器的魔法，但是为分析支持而编译会严重干扰这一点。对于像vector 或bytestring 这样使用自己的融合框架进行优化的包，分析干扰可能是相当灾难性的，分析结果完全无用。我倾向于相信我们这里有这样的案例。

在 Core 中，所有读取和写入都转换为 primops readDoubleArray#、indexDoubleArray# 和 writeDoubleArray#，它们很快。可能比 C 数组访问慢一点，但不是很慢。所以我相信这不是问题所在，也是造成巨大差异的原因。但是您已经在它们上添加了{-# SCC #-} 注释，因此禁用了涉及重新排列任何这些术语的任何优化。每次输入这些点之一时，都必须记录下来。我对分析器和优化器不够熟悉，无法知道到底发生了什么，但是，作为一个数据点，foldV、step 和 step 上的 {-# INLINE #-} 编译指示和 kahanStep 使用这些 SCC 进行了分析运行3.17 秒，并且 SCC fV_step、fV_read、fV_index、fV_write 和 fV_apply 被移除（没有其他改变），分析运行只需要 2.03 秒（两次都由 +RTS -P 报告，所以减去分析开销）。这种差异表明，具有廉价功能的 SCC 和过于细粒度的 SCC 会严重扭曲分析结果。现在，如果我们还在mkVect、kahan 和prng 上添加{-# INLINE #-} 编译指示，我们将得到一个完全没有信息的配置文件，但运行只需要1.23 秒。（但是，这些最后的内联对非分析运行没有影响，如果没有分析，它们会自动内联。）

因此，不要将分析结果视为无可置疑的事实。您的代码（直接或间接通过使用的库）对优化的依赖程度越高，就越容易受到禁用优化导致的误导性分析结果的影响。这也适用于堆分析以抑制空间泄漏，但程度要小得多。

当您有可疑的分析结果时，请检查删除某些 SCC 后会发生什么。如果这导致运行时间大幅下降，则该 SCC 不是您的主要问题（在修复其他问题后它可能再次成为问题）。

查看为您的程序生成的核心，跳出来的是您的 kahanStep - 顺便说一下，从中删除 {-# NOINLINE #-} 杂注，它会适得其反 - 产生了一对盒装的 Doubles在循环中，它立即被解构并且组件被拆箱。这种不必要的中间值装箱成本很高，并且会大大降低计算速度。

由于今天haskell-cafe 再次出现这种情况，有人使用 ghc-7.4.1 从上述代码中获得了糟糕的性能，tibbe 自行调查 GHC 生成的核心并发现 GHC 生成的代码不是最理想的从Word 到Double 的转换。将转换的fromIntegral 替换为仅使用（包装的）原语的自定义转换（并删除此处没有影响的爆炸模式，GHC 的严格分析器足以看穿算法，我应该学会信任更多 ;)，我们获得了一个与 gcc -O3 的原始 C 输出相当的版本：

{-# LANGUAGE CPP #-}
module Main (main) where

#define VDIM 100
#define VNUM 100000

import Data.Array.Base
import Data.Array.ST
import Data.Array.Unboxed
import Control.Monad.ST
import GHC.Word
import Control.Monad
import Data.Bits
import GHC.Float (int2Double)

prng :: Word -> Word
prng w = w'
  where
    w1 = w `xor` (w `shiftL` 13)
    w2 = w1 `xor` (w1 `shiftR` 7)
    w' = w2 `xor` (w2 `shiftL` 17)

type Vec s = STUArray s Int Double

kahan :: Vec s -> Vec s -> ST s ()
kahan s c = do
    let inner w j
            | j < VDIM  = do
                cj <- unsafeRead c j
                sj <- unsafeRead s j
                let y = word2Double w - cj
                    t = sj + y
                    w' = prng w
                unsafeWrite c j ((t-sj)-y)
                unsafeWrite s j t
                inner w' (j+1)
            | otherwise = return ()
    forM_ [1 .. VNUM] $ \i -> inner (fromIntegral i) 0

calc :: ST s (Vec s)
calc = do
    s <- newArray (0,VDIM-1) 0
    c <- newArray (0,VDIM-1) 0
    kahan s c
    return s

correction :: Double
correction = 2 * int2Double minBound

word2Double :: Word -> Double
word2Double w = case fromIntegral w of
                  i | i < 0 -> int2Double i - correction
                    | otherwise -> int2Double i

main :: IO ()
main = print . elems $ runSTUArray calc

【讨论】：

好吧，谢谢你的麻烦，但你回答了一个我没有问的问题（见我的更新）。时间是针对没有分析支持编译的程序，所以这不会是vector和分析器之间交互不良的问题
关于分析的部分的重点是你的分析结果可能毫无价值。请注意，您的分析时间大约是非分析时间的 12 倍，这是一个异常高的因素。在更新答案时再深入一点，请耐心等待，我是一个缓慢的作家。
所以你是说分析一般来说是没用的？因为我在问为什么 unsafeRead 和 unsafeWrite 在启用分析的情况下比 kahanStep 贵 10 倍，也启用了分析。如果您说由于分析开销而无法信任该度量，那么您通常是在声明分析器无用。
不，那未免言过其实。由于干扰优化，它可能无用，但不一定如此。它通常对不涉及高级和脆弱优化的代码很有用，例如vector 使用，但对于严重依赖它们的代码，它通常会完全破坏代码 - 但是请注意，ghc-7.4 中的新分析对优化的不良影响远小于早期的分析并产生更可靠的结果。了解何时以及在多大程度上信任分析器是一种魔法。
我现在已经在 GHC HEAD 中实现了 word2Double# 和 word2Float#。

【解决方案2】：

在所有这些看似Data.Vector 的代码中都有一个有趣的列表组合器混合。如果我做了第一个明显的修改，替换

mkVect = U.force . U.map fromIntegral . U.fromList . take 100 . iterate prng

正确使用Data.Vector.Unboxed:

mkVect = U.force . U.map fromIntegral . U.iterateN 100 prng

然后我的时间减少了三分之二——从real 0m1.306s 到real 0m0.429s 看起来除了prng 和zero 之外的所有顶级函数都有这个问题

【讨论】：

这对我来说时间减半（我不得不更新vector，Ubuntu中默认安装的版本没有iterateN）但这一切都归功于mkVect的改进我不感兴趣。它并不能解释或解释我的程序在 unsafeRead 和 unsafeWrite 上的花费是 kahanStep 的 10 倍以上。
你的想法是错误的。这只是一个独立的示例，它捕获了实际应用程序中出现的模式。我无法发布实际代码，因为它 a) 很大并且 b) 被 NDA 覆盖。
这个模块的核心是foldV，它是一个列表折叠。
几乎没有在分析统计中注册。并非世界上的所有事物都可以成为矢量。这是原始应用程序中的模式，所以我将其保留在独立示例中。
好地方，适用。 @Alinabi 不，您的程序在unsafeRead 和unsafeWrite 中的花费没有比在kahanStep 中的花费多10 倍（尽管那个也很便宜），它是一个分析工件。

【解决方案3】：

我知道您并没有要求改进此微基准测试的方法，但我会给您一个解释，在将来编写循环时可能会有所帮助：

一个未知的函数调用，例如对foldV 的高阶参数的调用，如果在循环中频繁执行，可能会很昂贵。特别是，它会抑制函数参数的拆箱，从而导致分配增加。它禁止参数拆箱的原因是我们不知道我们正在调用的函数在这些参数中是严格的，因此我们将参数传递为例如(Double, Double)，而不是Double# -> Double#。

如果循环（例如foldV）遇到循环体（例如kahanStep），编译器可以计算出严格性信息。出于这个原因，我建议人们INLINE 使用高阶函数。在这种情况下，内联 foldV 并删除 kahanStep 上的 NOINLINE 对我来说大大提高了运行时间。

在这种情况下，这并没有使性能与 C 相提并论，因为还有其他事情正在发生（正如其他人所评论的那样），但这是朝着正确方向迈出的一步（而且这是你可以不做的一步每个人都必须查看分析输出）。

【讨论】：

【解决方案4】：

这出现在邮件列表中，我发现 GHC 7.4.1 中的 Word->双重转换代码中存在错误（至少）。这个版本解决了这个 bug，和我机器上的 C 代码一样快：

{-# LANGUAGE CPP, BangPatterns, MagicHash #-}
module Main (main) where

#define VDIM 100
#define VNUM 100000

import Control.Monad.ST
import Data.Array.Base
import Data.Array.ST
import Data.Bits
import GHC.Word

import GHC.Exts

prng :: Word -> Word
prng w = w'
  where
    w1 = w `xor` (w `shiftL` 13)
    w2 = w1 `xor` (w1 `shiftR` 7)
    w' = w2 `xor` (w2 `shiftL` 17)

type Vec s = STUArray s Int Double

kahan :: Vec s -> Vec s -> ST s ()
kahan s c = do
    let inner !w j
            | j < VDIM  = do
                cj <- unsafeRead c j
                sj <- unsafeRead s j
                let y = word2Double w - cj
                    t = sj + y
                    w' = prng w
                unsafeWrite c j ((t-sj)-y)
                unsafeWrite s j t
                inner w' (j+1)
            | otherwise = return ()

        outer i | i <= VNUM = inner (fromIntegral i) 0 >> outer (i + 1)
                | otherwise = return ()
    outer (1 :: Int)

calc :: ST s (Vec s)
calc = do
    s <- newArray (0,VDIM-1) 0
    c <- newArray (0,VDIM-1) 0
    kahan s c
    return s

main :: IO ()
main = print . elems $ runSTUArray calc

{- I originally used this function, which isn't quite correct.
   We need a real bug fix in GHC.
word2Double :: Word -> Double
word2Double (W# w) = D# (int2Double# (word2Int# w))
-}

correction :: Double
correction = 2 * int2Double minBound

word2Double :: Word -> Double
word2Double w = case fromIntegral w of
                   i | i < 0 -> int2Double i - correction
                     | otherwise -> int2Double i

除了解决 Word->双重错误之外，我还删除了额外的列表以更好地匹配 C 版本。

【讨论】：