Monte Carlo Integration

原始题目：LeetGPU - Monte Carlo Integration

题目描述

在 GPU 上实现蒙特卡洛积分。给定在区间 $[a, b]$ 上均匀分布的随机采样点 $x_i$ 处的函数值 $y_i = f(x_i)$，估计定积分：

$$\int_a^b f(x)\ dx \approx (b - a) \cdot \frac{1}{n}\sum_{i=1}^{n} y_i$$

蒙特卡洛方法通过计算函数值的平均值并乘以区间宽度来近似积分。

实现要求

• 不允许使用外部库。
• solve 函数签名必须保持不变。
• 最终结果必须存储在 result 变量中。
• 绝对容差 $10^{-2}$，相对容差 $10^{-2}$。

示例

Input:  a = 0, b = 2, n_samples = 8
        y_samples = [0.0625, 0.25, 0.5625, 1.0, 1.5625, 2.25, 3.0625, 4.0]
Output: result = 3.1875

约束条件

• $1 \le n\_samples \le 100{,}000{,}000$。
• $-1000.0 \le a < b \le 1000.0$。
• $-10000.0 \le$ 函数值 $\le 10000.0$。
• 性能测试在 $n\_samples = 10{,}000{,}000$ 的规模下进行。

解题思路

蒙特卡洛积分在 GPU 上是天然并行的：每个线程计算一个采样点，最后做平均值规约。核心是高效的并行规约求总和。由于各采样点完全独立，没有数据依赖，计算可以在全 GPU 上充分展开。当 $n$ 很大时，需要注意浮点累加的顺序误差——可以使用 Kahan 求和或分块累加来减少精度损失。

代码实现

CUDA

#include <cuda_runtime.h>
__global__ void mc_kernel(const float* y, float* res, int n, float a, float b) {
    __shared__ float s[256]; int tid=threadIdx.x; float sum=0.0f;
    for(int i=blockIdx.x*blockDim.x+tid;i<n;i+=gridDim.x*blockDim.x)sum+=y[i];
    s[tid]=sum;__syncthreads();
    for(int ss=blockDim.x/2;ss>0;ss>>=1){if(tid<ss)s[tid]+=s[tid+ss];__syncthreads();}
    if(tid==0)atomicAdd(res,s[0]/n*(b-a));
}
extern "C" void solve(const float* y, float* res, int n, float a, float b) {
    mc_kernel<<<min((n+255)/256,1024),256>>>(y,res,n,a,b);
    cudaDeviceSynchronize();
}

Triton

import triton, triton.language as tl
@triton.jit
def mc_kernel(y_ptr,res_ptr,n,a,b,BLOCK:tl.constexpr):
    idx=tl.program_id(0)*BLOCK+tl.arange(0,BLOCK); mask=idx<n
    y=tl.load(y_ptr+idx,mask=mask,other=0.0)
    tl.atomic_add(res_ptr,tl.sum(y,axis=0)/n*(b-a))