つれづれなるままに

_ [Vtuber][cover] アカシア (Acacia) - BUMP OF CHICKEN // covered by 松永依織

RIOT MUSIC の松永依織のカヴァー。

ポケモンの曲なんですか。

彼女も最近は歌は歌ってるのだけど、配信が多いものだから聞くことができないでいた人です。

この人の歌は結構好きなのですがねぇ。

_ [科学][PC] fortran 版 モンテカルロ円周率計算 と GNU Scientific Library

昨日のやつをFortranにしてみましたが、メルセンヌツイスターのライブラリの定番ってないんですよね。

昨日失敗していたのは Multiple stream Mersenne Twister PRNG にあったソースを元にしたもの。

ライブラリのコンパイルで少々てこずりましたが、gfortranだけでやりたかったのでMakefile で USE_NTL = no とし、FCFLAG のところに -If_jump_ahead_coeff -fallow-invalid-boz を加え、-traceback をはずしました。FC は gfortran にします。下のディレクトリのMakefile にも -I.. を入れる必要があったり。

できあがった *.o と *.mod をプログラムと同じディレクトリに入れてからコンパイルします。ソースはこちら。

program pi
 use mt_kind_defs
 use mt_stream
 implicit none
 real(kind(0d0))::r,x,y,c,ttl,p
 integer(int64)::count=0,i
 integer(int64)::total=10000000000
 integer(int32)::iseed
 type(mt_state)::mts

 call set_mt19937
 call new(mts)
 call init(mts,iseed)
 do i=1,total
   x=genrand_double1(mts)-0.5d0
   y=genrand_double1(mts)-0.5d0
   r=sqrt(x*x+y*y)
   if(r.le.0.5d0) then
     count=count+1
   end if
 end do
 c=real(count,kind(0d0))
 ttl=real(total,kind(0d0))
 p=c/ttl*4.d0
 print *,p
end program

gfortranのバグなのか、整数の範囲にひっかかってしまうので、このソースをコンパイルするときは -fno-range-check を入れる必要があります。

$ gfortran -o run -I. -fno-range-check pi.f90
$ time ./run.exe
  3.1415799220000000

real    5m2.022s
user    0m0.000s
sys     0m0.000s

あまり精度が良くないのはやっぱり変数の取り方が悪いせいか。

Cのライブラリになってしまいますが、GNU Scientific Library のFortranラッパーである fgsl を使ってみることにしました。

MSYS2 に入ってるので pacman でインストールしておきます。

GitHubの例や、日本語翻訳版マニュアルを参考にコードを組みます。

program pi
 use fgsl
 implicit none
 real(kind(0d0))::ra,x,y,c,ttl,p
 integer(selected_int_kind(16))::count=0,i
 integer(selected_int_kind(16))::total=10000000000
 type(fgsl_rng) :: r
 type(fgsl_rng_type) :: t

 t=fgsl_rng_env_setup()
 t=fgsl_rng_default
 r=fgsl_rng_alloc(t)

 do i=1,total
   x=fgsl_rng_uniform(r)-0.5d0
   y=fgsl_rng_uniform(r)-0.5d0
   ra=sqrt(x*x+y*y)
   if(ra.le.0.5d0) then
     count=count+1
   end if
 end do
 c=real(count,kind(0d0))
 ttl=real(total,kind(0d0))
 p=c/ttl*4.d0
 print *,p
end program

$ gfortran -o run $(pkg-config --cflags fgsl)  -I/mingw64/include/fgsl -fno-range-check pi.f90 $(pkg-config --libs fgsl)
$ ls
pi.f90  pi.f90~  pi.o  run.exe*
$ time ./run.exe
  3.1415942239999999

real    2m8.332s
user    0m0.000s
sys     0m0.000s