蒙特卡罗方法估计π的值

这段代码的功能是使用蒙特卡罗方法估计π的值。它首先检查是否有足够的内存来存储所需的数组，然后在指定的设备（GPU或CPU）上生成随机点，并计算这些点落在单位圆内的比例。最后，根据这个比例估计π的值。这个过程会重复1000次，并将每次迭代的估计值写入一个名为"result.txt"的文件中。

注意升级路线：

1如果你想要提高估计π的准确性，可以尝试增加样本数量（num_samples）。这将使得更多的随机点被用于估计，从而提高准确性。但请注意，这可能会增加内存需求和计算时间。
2如果你的设备支持CUDA并且你想要利用GPU加速计算，确保你的PyTorch版本支持CUDA，并确保你的GPU有足够的内存来存储所需的数组。你可以调整num_samples的值以适应你的GPU内存。
3如果你想进一步提高性能，可以尝试使用更高效的数据类型（例如float16而不是float64），或者使用并行计算库（如cupy）来加速计算。
4如果你想要改进代码的可读性和可维护性，可以考虑将功能拆分成更小的函数，并为每个函数添加文档字符串以解释其作用和参数。

import torch
import cupy as cp
import numpy as np

def check_memory(num_samples, bytes_per_float):
    memory_per_array = num_samples * bytes_per_float
    return torch.cuda.get_device_properties(0).total_memory >= memory_per_array * 2

# 确定是否可以使用CUDA，如果可以，则设置为GPU
device = torch.device("cuda" if torch.cuda.is_available() else "cpu")

# 打印使用的设备
print(f"Using device: {device}")

# 假设每个float32数字占用4字节
bytes_per_float = 4

# 为x和y各分配约6GB内存
memory_per_array = 6 * 1024**3 # 6GB转换为字节

# 计算每个数组可以容纳的元素数量
num_samples_per_array = memory_per_array // bytes_per_float

# 因为我们需要为x和y两个数组分配内存，所以我们将样本数减半
num_samples = num_samples_per_array // 2

print(f"Using {num_samples} samples to estimate π.")

def estimate_pi(num_samples, device):
    # 在指定的设备上进行所有计算
    x = torch.rand(num_samples, device=device, dtype=torch.float64)
    y = torch.rand(num_samples, device=device, dtype=torch.float64)

    # 计算点是否在单位圆内
    inside_circle = (x**2 + y**2) <= 1.0
    num_inside_circle = inside_circle.sum().item()

    # 根据蒙特卡罗方法估计π的值
    pi_estimate = (4.0 * num_inside_circle) / num_samples
    return pi_estimate

# 根据你的GPU内存调整num_samples的值
num_samples = 100000000 # 示例值，根据你的GPU内存进行调整

# 检查是否有足够的内存来存储所需的数组
if not check_memory(num_samples, bytes_per_float):
    print("Not enough memory to run the simulation.")
else:
    # 运行10次蒙特卡罗方法估计π的值，并输出每个估计的实际数字
    with open("result.txt", "w") as f:
        for i in range(1000):
            pi_estimate = estimate_pi(num_samples, device)
            print(f"Iteration {i+1}: Estimated π： {pi_estimate:.20f}")
            f.write(str(pi_estimate) + "\n")