并行化蒙特卡洛

测量 CPU 时间

上篇文章写了个简单的蒙特卡洛估计器，用于估计 的值，但是那个非常慢，可以通过 System.Diagnostics 提供的 Stopwatch() 方法算一下时间。

var watch = new Stopwatch();

stopwatch.Start();

for (ulong i = 0; i < N; i++)
{
    x = rnd.NextDouble();
    y = rnd.NextDouble();
    if (x * x + y * y <= 1)
    {
        circleCount++;
    }
}
double PI = 4 * (double)circleCount / N;

watch.Stop();

Console.WriteLine("Pi = {0}, time = {1}ms", PI, watch.ElapsedMilliseconds);

设置采样数为 5000000000，得到 Pi = 3.1415527664, time = 37026ms。为了计算这个稳定于 3.1415 以上的结果，花费了 37 秒。

多线程

现代 CPU 具有多个物理核心，可以把它们利用起来。用 ThreadLocal<Random>() 为各个线程单独生成随机数。

using System.Threading;
using System.Threading.Tasks;

void ParallelEst()
{
    // 全局计数器
    long totalCircleCount = 0;

    var watch = new Stopwatch();
    ThreadLocal<Random> threadRnd = new ThreadLocal<Random>(() => new Random(Guid.NewGuid().GetHashCode()));    //各线程通过唯一 GUID 生成哈希值，防止随机样本序列重复

    watch.Start();

    // 循环并行化
    Parallel.For(
        0,
        (long)N,
        () => 0UL,
        (i, loopState, localCount) =>   // i: 当前迭代，localCount: 当前线程的局部计数器
        {
            Random rnd = threadRnd.Value;
            double x = rnd.NextDouble();
            double y = rnd.NextDouble();
    
            if (x * x + y * y <= 1)
            {
                localCount++;
            }
            return localCount; 
        },
        localCount =>  // 累积局部计数器到全局计数器
        {
            Interlocked.Add(ref totalCircleCount, (long)localCount);
        }
    );

    threadRnd.Dispose();

    double PI = 4 * (double)totalCircleCount / N;

    watch.Stop();

    Console.WriteLine("Pi = {0}, time = {1}ms", PI, watch.ElapsedMilliseconds);
}

以下是运行时间：

1 线程	24 线程
36412 ms	6348 ms

多线程的确占用起来了：

taskmgr

踩坑

之前是这么写的，结果并行化反而慢了 10 倍：

void ParallelEst()
{
    ulong circleCount = 0;
    
    var watch = new Stopwatch();
    
    ThreadLocal<Random> threadRnd = new ThreadLocal<Random>(() => new Random(Guid.NewGuid().GetHashCode()));
    
    watch.Start();
    
    Parallel.For(0, (long)N, i =>
    {
        Random rnd = threadRnd.Value;
        double x = rnd.NextDouble();
        double y = rnd.NextDouble();

        if (x * x + y * y <= 1)
        {
            Interlocked.Increment(ref circleCount); // 锁在线程循环里面
        }
    });
    
    threadRnd.Dispose();
    
    double PI = 4 * (double)circleCount / N;
    
    watch.Stop();
    
    Console.WriteLine("Pi = {0}, time = {1}ms", PI, watch.ElapsedMilliseconds);
}

原因是每个线程每次累计一次 circleCount 都要调用一次 Interlocked 方法，导致线程之间产生锁竞争，大幅降低并发性能。

修复方案就是给单个线程分配一个局部的计数器 localCount，等各线程累计结束之后再调用原子操作。