現在最速とされている、Chudnovskyの公式を用いて円周率を計算します。言語はRustです。実装にあたって、 こちらの記事を参考にしました。詳細な解説はリンク先をご覧ください。最終的な目標は分散システムでの並列計算ですが、本記事ではまず単一ノードでのスレッド並列を実装します。
実装
Binary Splitting Methodによる級数計算を、再帰関数を用いて実装します。$l$項目から$r$項目までの和を求める関数compute_xyzを、区間を2分割しながら再帰的に呼び出します。負の値を避けるため、区間の幅が2の場合の処理を別途用意しました。
fn compute_xk(k: u64) -> Integer {
if k == 0 {
return Integer::from(1);
}
else {
return Integer::from(k).pow(3) * C_24;
}
}
fn compute_yk(k: u64) -> Integer {
return A + B * Integer::from(k);
}
fn compute_zk(k: u64, n: u64) -> Integer {
if k == n - 1 {
return Integer::from(0);
}
else {
return Integer::from((6 * k + 1) * (2 * k + 1) * (6 * k + 5));
}
}
fn compute_xyz(l: u64, r: u64, n: u64) -> (Integer, Integer, Integer) {
if r - l == 1 {
return (compute_xk(l), compute_yk(l), -compute_zk(l, n));
}
else if r - l == 2 {
let _x1 = compute_xk(l);
let _x2 = compute_xk(l + 1);
let _y1 = compute_yk(l);
let _y2 = compute_yk(l + 1);
let _z1 = compute_zk(l, n);
let _z2 = compute_zk(l + 1, n);
return (_x1 * &_x2, &_x2 * _y1 - _y2 * &_z1, &_z1 * _z2);
}
else {
let m = (l + r) / 2;
let (x1, y1, z1) = compute_xyz(l, m, n);
let (x2, y2, z2) = compute_xyz(m, r, n);
return (x1 * &x2, &x2 * y1 + y2 * &z1, &z1 * z2);
}
}シングルスレッドの場合は$0$項目から$n$項目までを1スレッドで計算するため、compute_xyz(0, n, n)として計算します。マルチスレッドの場合はスレッドごとdisjointな区間を計算します。スレッドごとの計算終了後、計算結果を統合して円周率を計算します。x, y, zはそれぞれ(x1 * &x2, &x2 * y1 + y2 * &z1, &z1 * z2)にしたがって更新されるため、この計算式を用いることでスレッドごとの計算結果を統合することができます。マルチスレッドの計算にはstd::threadを使用しました。
use rug::{ops::Pow, Float, Integer};
use std::thread;
const A: u64 = 13591409;
const B: u64 = 545140134;
const C: u64 = 640320;
const C_CUBED: u64 = C.pow(3);
const C_24: u64 = C_CUBED / 24;
pub fn chudnovsky(n: u32, prec: u32, n_threads: u64) -> Float{
let __n: u64 = (n / 14).into();
let _n = __n + __n % 2;
if n_threads == 1 {
return chudnovsky_serial(n, prec);
}
let size = _n / n_threads;
let mut left = size + _n % n_threads;
let mut handles = Vec::new();
for _i in 1..n_threads {
let handle = thread::spawn(move || compute_xyz(left, left + size, _n));
handles.push(handle);
left += size;
}
let (xe, ye, ze) = compute_xyz(0, size + _n % n_threads, _n);
let (x, y, _z) = handles.into_iter().fold((xe, ye, ze), | (x, y, z), handle | {
let (xt, yt, zt) = handle.join().unwrap();
(x * &xt, &xt * y + yt * &z, &z * zt)
});
let pi = Float::with_val(prec, C_CUBED).sqrt() * x / 12 / y;
pi
}実験
検証に用いた計算機はM1 macbook air (8GB)、CPUコアはパフォーマンスコアが4コア、効率性コアが4コアの合計8コアです。コア数の制約もあってあまり桁数を増やせませんでしたが、並列化の恩恵がそれなりに感じられる結果となりました。
| 桁数 | 1スレッド | 2スレッド | 4スレッド | 8スレッド |
| $10^6$ | 9.61s | – | – | 9.53 s |
| $10^7$ | 13.41s | 11.27 s | 10.85 s | 11.54 s |
| $10^8$ | 73.38 s | 37.67 s | 31.60 s | 40.52 s |
この記事は役に立ちましたか？
もし参考になりましたら、下記のボタンで教えてください。
コメント