2026a

# 蒙特卡洛方法计算 π

以蒙特卡洛方法求 π 为示例讨论 Julia 中的串行和多种并行计算方式。

提示

这是经典范式中的第四个综合示例。它分别包含了多线程、多进程以及 CUDA 编程的示例，每个示例相对独立，你可以选择性地跳过你不感兴趣的内容。

当你在阅读具体示例之前，建议先阅读快速入门(多线程快速入门、多进程快速入门、或 GPU 计算快速入门 )来了解相关的内容背景。

多线程环境：以下示例需要配置多线程环境，演示在 4 线程环境下运行，注意需要通过 using Base.Threads 加载多线程模块。

本例中需要加载的包：

using Base.Threads
using BenchmarkTools
using Random
using Distributed
using SharedArrays

# 背景介绍

蒙特卡洛法是一种通过随机采样和统计方法来求解复杂问题的数学计算方法。它的名称来自于摩纳哥的蒙特卡罗赌场，它在统计学和计算机科学领域中被广泛应用于模拟随机事件。

蒙特卡洛法通常涉及随机模拟，以生成足够数量的数据样本，然后根据这些样本进行数学分析。它在估计各种量、解决优化问题、处理不确定性问题等方面有广泛的应用，特别是在统计学、金融学、物理学、工程学、计算机科学等领域。

例如，假设你要计算圆周率 π 的值。蒙特卡罗法可以通过生成大量的随机点，然后计算落在圆内的点数与总点数之比来估算 π 的值。这种方法的估计值的精度与模拟的随机样本数量有关，随着样本数量的增加，它可以逐渐接近真实值。

总之，蒙特卡洛法通过生成随机样本并使用统计学方法来求解复杂问题，是一种广泛应用于各种领域的计算方法。

而用 Monte-Carlo 方法计算 π 具体来说是，给定一个半径为的圆在一个边长为的正方形内，则圆的面积与正方形面积关系与 π 相关，

圆 的 面 积 正 方 形 面 积

随机生成在正方形内的均匀样本，并计算圆内有多少点。在圆中找到一个点的概率即是圆的面积除以正方形的面积的比值。要确定一个点是否在圆内，使用点的和坐标计算出离原点的距离即可，。如果小于圆半径，则该点在圆内。生成一个大样本，并计算圆内有多少点。则可以用以下公式来估计 π ，

圆 内 点 数 总 点 数

# Julia 串行实现

编写一个蒙特卡洛求 π 的串行程序：

function darts_in_circle(n, rng=Random.GLOBAL_RNG)
    inside = 0
    for i in 1:n
        # 随机生成点
        if rand(rng)^2 + rand(rng)^2 < 1
            inside += 1
        end
    end
    return inside
end
function pi_serial(n)
    return 4 * darts_in_circle(n) / n
end

如示意图，在的位置生成点，并计算在圆内比例来计算 π 。

现在测试串行计算的时间：

@btime pi_serial(10_000_000)

  23.492 ms (0 allocations: 0 bytes)
3.1417244

可以看到结果与预期一致。

# 多线程计算

考虑在单机上用共享内存多线程并行化这个计算。

查看线程数，此时 Julia 已经以多线程启动。

nthreads()

现在以多线程的方法改写这个程序并测试性能，多线程下需要不同线程独立正常地生成随机数。

但由于普通的rand()方式基于一个有状态的全局随机数生成器，每次生成一个新数字时它都会改变内部状态，这意味着当多个线程试图同时生成随机数时，会遇到类似多个线程试图更改内存中相同位置的问题。

即如果简单地改写串行算法，会得到错误的结果：

function darts_in_circle_wrong(n, rng=Random.GLOBAL_RNG)
    inside = 0
    @threads for i in 1:n
        # 随机生成点
        if rand(rng)^2 + rand(rng)^2 < 1
            inside += 1
        end
    end
    return inside
end
function pi_threads_wrong(n)
    return 4 * darts_in_circle_wrong(n) / n
end
@btime pi_threads_wrong(10_000_000)

  129.530 ms (7849468 allocations: 119.78 MiB)
0.8262664

这样的情况下，我们引入一个rnglist，来保证每个线程都有自己的随机数生成器。

在设置了随机数生成器后，编写函数以确保每个线程使用一个单独的 RNG(Random Number Generators)。通过使用线程 ID 索引到 RNG 列表中来确保这一点。

注意，由于计算中使用了线程 ID，故这里不推荐使用Threads.@spawn进行线程并行（可能会发生任务迁移而不能保证执行中的threadid如预期）。

const rnglist = [MersenneTwister() for i in 1:nthreads()];
function pi_threads(n, loops)
    inside = zeros(Int, loops)
    @threads for i in 1:loops
        rng = rnglist[threadid()]
        inside[threadid()] = darts_in_circle(n, rng)
    end
    return 4 * sum(inside) / (n * loops)
end
@btime pi_threads(2_500_000, 4)

  4.241 ms (22 allocations: 2.56 KiB)
3.1411724

在每个线程内部对整个问题的一部分进行串行计算，而不同线程之间是并行的，最后，将结果添加到inside数组的每个线程对应位置中。

在性能上也能看到这是在 4 核心以上的机器上运行 4 线程所期望的性能优化。

# 多进程计算

现在考虑用分布式内存多进程并行使用蒙特卡洛计算 π 。

在只含有主进程的 Julia 中使用 Distributed 库中的addprocs添加 4 个工作进程：

addprocs(4)
procs()

5-element Vector{Int64}:
 1
 2
 3
 4
 5

首先，广播给工作进程定义函数来计算单位圆内的命中次数, 然后通过pmap在多个节点上并行运行。

# 计算单位圆命中次数
@everywhere function darts_in_circle(n)
    inside = 0
    for i in 1:n
        # 随机生成点
        if rand()^2 + rand()^2 < 1
            inside += 1
        end
    end
    return inside
end
# pmap版本
function pi_pmap(N, loops)
    return 4 * sum(pmap((x) -> darts_in_circle(N), 1:loops)) / (loops * N)
end
@btime pi_pmap(1_000_000, 10)

  7.369 ms (642 allocations: 27.20 KiB)
3.1409112

由于是多进程，也不存在共享随机数生成器的问题。可以看到复杂计算的情况下多进程可以带来的性能提升。

# GPU 计算

提示

CUDA 支持：GPU 计算需要 Nvidia GPU 硬件设备支持，用户也需要提前安装 CUDA.jl 包并using CUDA

由于 GPU 更适合使用 Float32 类型计算，我们以 FLoat32 数据类型进行计算对比：

function pi_gpu_broadcast_f32(n)
    CUDA.@sync 4 * sum(CUDA.rand(Float32, n).^2 .+ CUDA.rand(Float32, n).^2 .<=1) / n
end
@btime pi_gpu_broadcast_f32(10_000_000)

  881.500 μs (124 allocations: 6.98 KiB)
3.1410404

可以根据这个风格重写一下串行的程序进行对比：

function pi_cpu_broadcast_f32(n)
    4 * sum(rand(Float32, n).^2 .+ rand(Float32, n).^2 .<=1) / n
end
@btime pi_cpu_broadcast_f32(10_000_000)

  14.121 ms (8 allocations: 77.49 MiB)
3.1419544

GPU 程序中只添加了CUDA.@sync进行同步和CUDA.rand生成 GPU 上的数据类型进行计算。可以看到成倍的性能提升。

同时，我们给出一个 Float64 数据类型下的基准测试结果便于横向对比，不过用户要注意，这并不是 GPU 计算推荐使用的类型。

function pi_gpu_broadcast_f64(n)
    CUDA.@sync 4 * sum(CUDA.rand(Float64, n).^2 .+ CUDA.rand(Float64, n).^2 .<=1) / n
end
@btime pi_gpu_broadcast_f64(10_000_000)

  1.684 ms (192 allocations: 10.36 KiB)
3.1423076

function pi_cpu_broadcast_f64(n)
    4 * sum(rand(Float64, n).^2 .+ rand(Float64, n).^2 .<=1) / n
end
@btime pi_cpu_broadcast_f64(10_000_000)

  32.448 ms (8 allocations: 153.78 MiB)
3.1428124

我们在此使用的 GPU 计算设备是 NVIDIA GeForce RTX 4060 Laptop GPU 。

# 性能对比

我们得到以下时间对比：

`Float64` 的 `@btime` 测试	时间
`pi_cpu_broadcast_f64` CPU 向量化计算	32.448 ms
`pi_gpu_broadcasts_f64` GPU 向量化计算	1.684 ms
`pi_serial`函数单线程 CPU 计算	23.492 ms
`pi_threads` 函数 4 线程 CPU 计算	4.241 ms
`pi_pmap`函数 4 进程 CPU 计算	7.369 ms

`Float32` 的 `@btime` 测试	时间
`pi_gpu_broadcasts_f32` GPU 向量化计算	881.500 μs
`pi_cpu_broadcast_f32` CPU 向量化计算	14.121 ms