为什么这个蒙特卡罗的Haskell程序运行如此缓慢？

Question

更新：

我的编译命令是ghc -O2 Montecarlo.hs。我使用的随机版本是random-1.1，ghc版本是8.6.4，我的系统是macOS Big Sur 11.1（Intel芯片）。我用来测试速度的命令是time ./Montecarlo 10000000，它返回的结果是real 0m17.463s, user 0m17.176s, sys 0m0.162s。

---

以下是一个使用蒙特卡罗方法计算π的Haskell程序。然而，当输入为1000万时，该程序运行了20秒。在相同逻辑下编写的C程序只需要0.206秒。为什么会这样，并且如何加快速度？谢谢大家。

这是Haskell版本：

import System.Random
import Data.List
import System.Environment

montecarloCircle :: Int -> Double
montecarloCircle x
   = 4*fromIntegral
        (foldl' (\x y -> if y <= 1 then x+1 else x) 0
           $ zipWith (\x y -> (x**2 + y**2))
                     (take x $ randomRs (-1.0,1) (mkStdGen 1) :: [Double])
                     (take x $ randomRs (-1.0,1) (mkStdGen 2) :: [Double]) )
        / fromIntegral x

main = do
    num <- getArgs
    let n = read (num !! 0)
    print $ montecarloCircle n

以下是C语言版本：

#include <stdio.h>
#include <math.h>
#include <time.h>
#include <stdlib.h>

#define N       10000000

#define int_t   long // the type of N and M

// Rand from 0.0 to 1.0
double rand01()
{
    return rand()*1.0/RAND_MAX;
}

int main()
{
        srand((unsigned)time(NULL));

        double x,y;
        int_t M = 0;

        for (int_t i = 0;i < N;i++)
        {
                x = rand01();
                y = rand01();
                if (x*x+y*y<1) M++;
        }
        double pi = (double)4*M/N;
        printf("%lf\n", pi);
}

- XM Zg

---

无法重现。在我的机器上，使用gcc -O3编译的时间为0.11秒，使用ghc -O2编译的时间为0.69秒。请展示您如何编译和运行它们。如果性能关键，可能可以将它们的性能差距缩小到7倍以下--我的直觉是两者之间的伪随机数生成算法可能不同。但100倍绝对是错误的。 - Daniel Wagner

为什么你在 Haskell 中需要2个随机序列，而在 C 中只需要一个？顺便说一下，如果代码限制在80列左右（不需要水平标尺），事情会更容易，谢谢。 - jpmarinier

在我的机器上（旧的Intel x86-64），你的C代码需要163毫秒，所以我们的机器看起来大致相似。使用GHC v8.8.4，即使关闭了优化，我也无法使Haskell代码变慢超过2秒。使用稍微更具惯用性的Haskell代码，例如insiderCount = length $ filter (<= 1.0) squaredNorms，并像您的C代码一样将x ** 2替换为x * x，我只需要455毫秒，因此比C代码慢不到3倍。为什么您看到100倍的速度差异仍然是个谜团：-( - jpmarinier

@XMZg 你使用的 Haskell random-1.1 很重要。我的结果是基于 random-1.2。而且版本1.2主要是为了回应抱怨，即在 v1.1 中 StdGen 太慢了（https://alexey.kuleshevi.ch/blog/2021/01/29/random-interface/）。顺便说一句，你可能想要编辑你的问题，包括你后来在评论中提供的所有性能相关信息，希望 SO 的“权力机构”能重新审视他们决定冻结你的问题。 - jpmarinier

@jpmarinier 是的！你说得对。当我按照你说的去做时，程序的运行时间急剧下降到了1.22秒。这让我意识到随机数算法中有很多奥秘。再次感谢你的建议和思考，对于像我这样的新手来说非常有帮助。 - XM Zg

Answer 1

在数值应用程序中，Haskell代码比使用传统命令式语言（如C/C++/Fortran）编写的等效代码要慢2倍到5倍左右。然而，Haskell代码比C代码慢100倍非常出乎意料！

检验结果的正确性

在一台配备Intel Celeron N2840 64位处理器，运行Linux内核5.3、GLIBC 2.28、GCC 8.3、GHC 8.2.2和GHC random-1.1的旧笔记本电脑上，我们得到以下计时结果：

C代码，gcc -O3:    950毫秒
Haskell代码，ghc -O3:    104秒

因此，在这些配置下，Haskell代码的运行速度大约是C代码的100倍。

为什么会出现这种情况？

首先，需要指出的是，利用伪随机数生成实现π的计算看起来非常简单，因此运行时间可能主要受到生成2000万个伪随机数的影响。此外，我们有充分的理由相信Haskell和C使用完全不同的伪随机数生成算法。

因此，我们可以比较这两种语言使用的随机数生成算法及其各自的实现，而不是将Haskell与C进行比较。

在C语言中，语言标准没有指定srand()函数应该使用哪个算法，但通常会使用旧的、简单和快速的算法。

另一方面，在Haskell中，我们传统上看到了很多关于StdGen效率低下的抱怨，就像2006年的这篇文章中所述：

https://gitlab.haskell.org/ghc/ghc/-/issues/427

一位领先的Haskell专家提到，StdGen可能会变得快30倍。

幸运的是，帮助已经在路上了。这篇最近的博客文章解释了新的Haskell random-1.2包如何使用完全不同的算法（称为splitmix）解决了StdGen速度不足的问题。

伪随机数生成（PRNG）领域是一个相当活跃的领域。算法通常会被更新和改进。为了更好地了解情况，这里有一篇相对较新的（2018年）关于该主题的综述论文。

转向更加现代、更好的Haskell组件：

使用另一台稍微更强大的机器，配备Intel Core i5-4440 3GHz 64位CPU，运行Linux内核5.10、GLIBC 2.32、GCC 10.2和至关重要的Haskell包random-1.2：

C代码，gcc -O3：164毫秒
Haskell代码，ghc -O3：985毫秒

因此，Haskell现在“只是”比C慢6倍而不是100倍。

我们仍然需要解决的问题是，Haskell必须使用x**2+y**2而C则获得x*x+y*y的不公平之处，在使用random-1.1时这并不重要。这给了我们379毫秒！所以我们回到了通常的2倍到5倍的Haskell比C速度比较范围。

请注意，如果我们要求Haskell可执行文件以统计信息运行，我们将得到以下输出：

$ time q66441802.x +RTS -s -RTS 10000000
3.1415616
     925,771,512 bytes allocated in the heap
     ...
  Alloc rate    2,488,684,937 bytes per MUT second
Productivity  99.3% of total user, 99.3% of total elapsed
real   0m0,379s
user   0m0,373s
sys    0m0,004s

因此，发现Haskell在运行过程中分配了接近1GB的内存，这有助于理解速度差异。

代码修复：

我们注意到C代码使用单个随机序列，而Haskell代码使用两个序列，其中包括两个调用mkStdGen，种子为1和2。这不仅不公平，而且是不正确的。设计PRNG算法的应用数学家非常注意确保任何单个序列具有正确的统计特性，但几乎不保证可能存在的不同序列之间的相关性。甚至可以将种子用作单个全局序列的偏移量。

以下代码对此进行了修复，其性能没有显著改变：

computeNorms :: [Double] -> [Double]
computeNorms []  = []
computeNorms [x] = []
computeNorms (x:y:xs2) = (x*x + y*y) : (computeNorms xs2)

monteCarloCircle :: Int -> Double
monteCarloCircle nn =
    let
         randomSeed   =  1
         gen0         =  mkStdGen randomSeed
                      -- use a single random serie:
         uxys         =  (randomRs  (-1.0, 1.0)  gen0) :: [Double]
         norms        =  take nn (computeNorms uxys)
         insiderCount =  length  $  filter (<= 1.0) norms
    in
         (4.0::Double) * ((fromIntegral insiderCount) / (fromIntegral nn))

附注：

新的Haskell random-1.2包已经在这个最近的SO问题中提到，尽管是在新的单子界面的上下文中。

总结一下：

假设练习的目标是衡量C和Haskell语言的相对运行时速度，那么基本上有两种可能性：

一种是完全避免使用PRNG，因为使用了不同的算法。

另一种方法是通过手动提供一个相同的算法来控制随机数生成。例如，可以使用由Pierre L'Ecuyer设计的公开可用的MRG32k3a算法。 此处提供了候选的Haskell MRG32k3a实现（截至2021年3月）。留给读者作为练习。

- jpmarinier