Qiita Teams that are logged in
You are not logged in to any team

Log in to Qiita Team
Community
OrganizationAdvent CalendarQiitadon (β)
Service
Qiita JobsQiita ZineQiita Blog
0
Help us understand the problem. What is going on with this article?
@tnagai-github

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

はじめに

先日の「Rustで円周率をモンテカルロ法で求めた」のコメントなどのお陰で
Rustの文法に多少詳しくなったので、今度は数値積分で円周率を求めてみた。
やることは単純で、$\pi = 4\int_0^1 \sqrt{(1-x^2)}dx$なので、
この積分をなんの工夫もしない区分求積で求めるだけである1,2
同じくらいの抽象度と思われる、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)は厳密には非標準手続き。3

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

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

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

0
Help us understand the problem. What is going on with this article?
Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away

Comments

No comments
Sign up for free and join this conversation.
Sign Up
If you already have a Qiita account Login
0
Help us understand the problem. What is going on with this article?