Rustで円周率を数値積分で求めてみた（C++とFortran付き）

はじめに

先日の「Rustで円周率をモンテカルロ法で求めた」のコメントなどのお陰で
Rustの文法に多少詳しくなったので、今度は数値積分で円周率を求めてみた。
やることは単純で、$\pi = 4\int_0^1 \sqrt{(1-x^2)}dx$なので、
この積分をなんの工夫もしない区分求積で求めるだけである¹,²。
同じくらいの抽象度と思われる、C++とmodern Fortranのコードも作成し、比較できるようにした。

---

ソースコード（Rust）

use std::env;

fn main() {
    let argc = env::args().len();

    if argc != 2 {
        println!("need exactly single argument");
        println!("usage: /path/to/get_pi2 num_division");
        std::process::exit(1);
    }
    let num_division = env::args()
        .nth(1)
        .map(|arg| arg.parse::<u32>())
        .unwrap()
        .unwrap();

    println!("num_division: {}", num_division);

    let mut xlist = Vec::<f64>::new();
    let dx = 1.0 / num_division as f64;

    for i in 0..num_division {
        xlist.push(i as f64 * dx)
    }
    let ylist = xlist
        .into_iter()
        .map(|x| (1.0 - x * x).sqrt())
        .collect::<Vec<f64>>();

    let area = ylist.iter().sum::<f64>() * dx;

    println!("pi is calclated to be {}", area * 4.0);
    println!("while the actual pi is {}", std::f64::consts::PI);
}

---

実行結果

cargo run --release -- 100000などと実行すれば、

num_division: 100000
pi is calclated to be 3.1416126164020075
while the actual pi is 3.141592653589793

と円周率の計算結果が示される。

---

ソースコード（C++）

#include<iostream>
#include<vector>
#include<algorithm>
#include<iomanip>
#include<numeric>
#include<cmath>

int main (int argc, char** argv){

    if (argc != 2){
        std::cout << "need exactly single argument" << std::endl; 
        std::cout << "usage: /path/to/get_pi2_cxx num_division" << std::endl; 
        exit(EXIT_FAILURE);
    }

    const int num_division = std::atoi(argv[1]);

    std::cout << "num_division: " << num_division << std::endl; 

    std::vector<double> xlist;
    const double dx = 1.0/num_division;

    for(int i=0; i<num_division; i++){
        xlist.push_back(i*dx);
    }

    std::vector<double> ylist(num_division);

    std::transform(xlist.begin(), xlist.end(), ylist.begin(), [](double x){return std::sqrt(1.0-x*x);}) ;

    const double area = std::accumulate(ylist.begin(), ylist.end(), 0.0) *dx;

    std::cout << std::setprecision(17); 
    std::cout << "pi is calculated to be " << area*4.0 <<std::endl;
    std::cout << "while the actual value is " << M_PI <<std::endl;
}

---

ソースコード（Modern Fortran）

program main
implicit none

integer :: i
integer :: num_division

real(8), allocatable :: xlist(:)
real(8), allocatable :: ylist(:)
real(8) :: area;
real(8) :: dx;

character(len=128) :: string_num_division
real(8), parameter :: PI = 3.1415926535897931

if( iargc() /= 1 )then
  write(*,*) "need exactly single argument"
  write(*,*) "usage: /path/to/get_pi2 num_division"
  stop
endif

call getarg(1,string_num_division)
read(string_num_division, * )  num_division
write(*,*) "num_division: ", num_division

dx = 1.0d0 / num_division

xlist  = (/ (i*dx,i=0,num_division-1) /)

ylist = sqrt(1.0d0 -xlist(:)*xlist(:))

area = sum(ylist(:))*dx

write(*,*) "pi is calclated to be ", area*4d0
write(*,*) "while the actual pi is  ",  PI

stop
end program main

• iargc()とかgetarg(N,STRING)は厳密には非標準手続き。³

1. 真剣に数値計算するときは、台形公式くらいは使いたい。 ↩

2. 円周率を求める数学的に高度な方法こことか見て下さい ↩

3. http://jjoo.sakura.ne.jp/tips/f90/getarg.html ↩