蒙特卡洛计算圆周率 PI
蒙特卡洛计算圆周率 PI 其实是一个非常经典的案例,在 R 仿真里面,也都有。
基本上案例如下:
- 在一个边长为 1 的正方形里面,有一个半径为 1/2 的圆。
- 正方形和圆共一个中心点。圆内切正方形。
- 向这个正方形里面投点,计算点落在圆内的数量。落在圆内点的数量比上总点数就是等于 π / 4 \pi/4 π/4。
那么按照这个思想,只要是不断的生成若干个点,计算占比即可。
C 代码
这里先分享一个 C 语言的代码:
#include <stdlib.h>
#include <stdio.h>
#include <chrono>
#include <array>
#define N_POINTS 10000000
#define N_REPEATS 5
float estimate_pi(int n_points) {
double x, y, radius_squared, pi;
int within_circle=0;
for (int i=0; i < n_points; i++) {
x = (double)rand() / RAND_MAX;
y = (double)rand() / RAND_MAX;
radius_squared = x*x + y*y;
if (radius_squared <= 1) within_circle++;
}
pi=(double)within_circle/N_POINTS * 4;
return pi;
}
int main() {
double avg_time = 0;
srand(42);
for (int i=0; i < N_REPEATS; i++) {
auto begin = std::chrono::high_resolution_clock::now();
double pi = estimate_pi(N_POINTS);
auto end = std::chrono::high_resolution_clock::now();
auto elapsed = std::chrono::duration_cast<std::chrono::nanoseconds>(end - begin);
avg_time += elapsed.count() * 1e-9;
printf("Pi is approximately %g and took %.5f seconds to calculate.\n", pi, elapsed.count() * 1e-9);
}
printf("\nEach loop took on average %.5f seconds to calculate.\n", avg_time / N_REPEATS);
}
这里在 clion 里面运行一下,可以看到非常的快,基本上都是在 0.19 秒左右。
python 代码
既然 c 语言这么快,那 python 是什么情况,运行效率如何,这里给到 python 代码:
import random
import time
import argparse
def estimate_pi(
n_points: int,
show_estimate: bool,
) -> None:
"""
Simple Monte Carlo Pi estimation calculation.
Parameters
----------
n_points
number of random numbers used to for estimation.
show_estimate
if True, will show the estimation of Pi, otherwise
will not output anything.
"""
within_circle = 0
for _ in range(n_points):
x, y = (random.uniform(-1, 1) for v in range(2))
radius_squared = x**2 + y**2
if radius_squared <= 1:
within_circle += 1
pi_estimate = 4 * within_circle / n_points
if not show_estimate:
print("Final Estimation of Pi=", pi_estimate)
def run_test(
n_points: int,
n_repeats: int,
only_time: bool,
) -> None:
"""
Perform the tests and measure required time.
Parameters
----------
n_points
number of random numbers used to for estimation.
n_repeats
number of times the test is repeated.
only_time
if True will only print the time, otherwise
will also show the Pi estimate and a neat formatted
time.
"""
start_time = time.time()
for _ in range(n_repeats):
estimate_pi(n_points, only_time)
if only_time:
print(f"{(time.time() - start_time)/n_repeats:.4f}")
else:
print(
f"Estimating pi took {(time.time() - start_time)/n_repeats:.4f} seconds per run."
)
def positive_integer(value: str) -> int:
"""Check for positive integer in arg_parse."""
int_value = int(value)
if int_value <= 0:
raise argparse.ArgumentTypeError(f"{value} is an invalid positive int value")
return int_value
def main(arguments=None):
"""Main loop in arg parser."""
parser = argparse.ArgumentParser()
parser.add_argument(
"-p",
"--n_points",
help="Number of random points to use for estimating Pi.",
type=positive_integer,
default=1_000_000,
)
parser.add_argument(
"-r",
"--n_repeats",
help="Number of times to repeat the calculation.",
type=positive_integer,
default=5,
)
parser.add_argument(
"--only-time",
action='store_true',
default=False,
)
args = parser.parse_args()
run_test(
args.n_points,
args.n_repeats,
args.only_time,
)
if __name__ == "__main__":
main()
这里是 python 的运行结果,基本上在 1.5 秒左右
rust 代码
[package]
name = "code091501"
version = "0.1.0"
edition = "2021"
[dependencies]
rand = "0.8.5"
[profile.release]
opt-level = 3
use rand::Rng;
use std::time::Instant;
// use std::time::Duration;
// use std::thread::sleep;
fn calculate_pi(sample_n: usize) -> f32 {
let mut in_circile = 0;
let mut rng = rand::thread_rng();
for _ in 0..sample_n {
let x = rng.gen_range(0.0..1.0);
let y = rng.gen_range(0.0..1.0);
let r = x * x + y * y;
if r < 1.0 {
in_circile += 1;
}
}
let pi = in_circile as f32 / sample_n as f32 * 4.0;
pi
}
fn main() {
for _ in 0..10 {
let now = Instant::now();
let pi = calculate_pi(10000000);
// let pi = 2.00;
// sleep(Duration::new(2, 0));
let ela = now.elapsed();
let use_time = ela.subsec_nanos() as f64 * 1e-9 + ela.as_secs() as f64;
println!("the pi is : {} use time is: {:.3} s", pi, use_time);
}
}
rust 运行效果如下:
R 语言
我的第一个入门语言,那这里我也展示一下 R 的运行效率。
library(tidyverse)
calculate_pi <- function(sample_n) {
in_circle = 0
for (index in 1:sample_n) {
x = runif(1, min = 0, max = 1)
y = runif(1, min = 0, max = 1)
r = x*x + y*y
if (r <1) {
in_circle = in_circle+1
}
}
pi = in_circle / sample_n * 4
return(pi)
}
for (index in 1:10) {
start = Sys.time()
pi = calculate_pi(10000000)
end = Sys.time()
els = end - start
cat("the pi is :", pi, "use time is : ", els, '\n')
}
R语言的运行结果:
结论
- C 和 rust 运行时间基本在 0.15 秒左右。
- python 运行时间基本在 1.5 秒左右。
- R 语言运行时间基本在 43 秒左右。
这里要说明一下
- C是在window11系统的wsl2虚拟机下面运行的。cpu:17-8750H,内存16G
- python是在window11系统下的anaconda下运行的,python版本是3.8。cpu:17-8750H,内存16G
- rust是在mac系统下运行的。rust版本是1.62。机器版本是macbook 16寸2019款,cpu:2.6 GHz 六核Intel Core i7
- R是在mac系统下运行的。R版本是4.2.1。机器版本是macbook 16寸2019款,cpu:2.6 GHz 六核Intel Core i7
整理表格如下:
语言 | 运行时间| |
---|---|
C 语言 | 0.19s |
Rust | 0.15s |
python | 1.5s |
R | 43.0s |
阅读更多
list