速ければいいのだ！でも早すぎるのは嫌われるよね

Pythonのライブラリーのひとつにnumpyというのがある。数値計算に特化したライブラリである。高速に数値計算が行えるというが、残念ながらPython自体がコンパイラ言語ではないので、それほど速くなるというわけではない。やはりC/C++やFortranと較べると、速度的に遅いのは否めない。Pythonをコンパイルして実行することもでき、いろいろな方法があるようだが、あまり一般的ではないし、誰もが簡単にというわけにはいかない。

そこで考えられるのが、時間のかかる計算は、別のプログラミング言語にまかせようというやり方である。この方法は、Python以外のプログラミング言語を習得しなければならないが、Pythonでプログラミングを経験していれば、それほど難しくはないだろう。では、どのプログラミング言語を使うかというと、高速な数値計算において右に出るものはないといえるFortranである。レガシー言語でありながら、「まだまだ若いモンには負けられんわい」という気の強さがあり、元気ハツラツな言語だ。

方法としては、Fortranで時間のかかる部分のコードを書き、それをPythonで呼び出すというやり方である。「え〜、そんなの困っちゃう」と思われる方もおられるかもしれないが、実にレガシーさとモダンさがミックスしたやり方なのである。ワクワクするではないか。

さいわいにもPythonには、f2pyというライブラリが用意されており、そのライブラリでFortranのコードをコンパイル、その後、Pythonからimport文を使って呼び出すだけでいい。f2pyはnumpyをインストールすれば、一緒にインストールされるので、使わない手はないのである。

では実際にやってみよう。まず、numpyをpipを使ってインストールする（linuxでな）。

$ pip install numpy

f2pyが実行できるか試してみよう。おなじみのバージョン表示である。

$ f2py -v
1.21.4

検索すると、

sudo apt install python-numpy

と紹介されていることがあるが、この方法だと、コンパイル時にエラーを食らってしまうことがある。なので使わないほうがいいかもしれない。もし、上記の方法でインストールして、コンパイルしてしまったら、次のコマンドで削除しよう。Linuxのいいとこのひとつは、使った容量は、アンインストール時にちゃんと返してくれるところだ。

$ sudo apt purge python-numpy --autoremove

無事にインストールされたら、Fortranでコードを書いてみよう。コンパイルはf2pyを使うので、gfortranはインストールする必要がない。ただし、Fortranについて学習する場合は、インストールしておこう。Fortranについて知りたい方は、「Fortranで遊ぶ」に書いておいた。学習するサイトも案内してある。

Fortranのコードは、「Fortranで遊ぶ」に掲載されているjulian.f90を少し手直ししてPythonから呼び出してみよう。

! ユリウス日を求める
subroutine Julian(year, month, day, jd)
    implicit none

    ! 次の２行がポイント、入力と出力を指定する
    double precision, intent(In) :: year, month, day
    double precision, intent(Out) :: jd
    double precision :: yy, mm, dd
    
    yy = year
    mm = month
    dd = day

    if (mm < 3) then
        mm = mm + 12.0d0
        yy = yy - 1.0d0
    end if

    jd = int(yy * 365.25d0) + int(yy / 400.0d0) - int(yy / 100.0d0) + int(30.59d0 * (mm - 2.0d0)) + dd + 1721088.5d0
    
    return
end subroutine Julian

わざわざ引数で指定したのに、year、month、dayを別の変数に代入しているのかと思われる方もおられるかもしれない。これは決まりごとで、intent(In)で指定した変数は、サブルーチンの中で書き換えることができないことになっている。

さて、入力したコードに間違いがなければコンパイルする。

$ f2py -m flib -c julian.f90

オプションの-mに続いてモジュール名を指定し、-cオプションでソースファイルを指定する。なにやらずらずらと表示されるが、ErrorがでなければOKである。

作られたモジュール名は、次のようになっている。長ったらしいがこのままでも構わない。正常に呼び出される。長ったらしいのが嫌な人は名前を変更しておこう。

flib.cpython-38-x86_64-linux-gnu.so → flib.so

次はpythonで呼び出すコードjulianday.pyを書く。ファイル名は何でも構わない。

import flib  # 呼び出すモジュール名


year = 2021
month = 11
day = 15

jd = flib.julian(year, month, day)

print()
print('Julian day = %13.5f' % (jd))
print()

サブルーチンJulianは引数が４つなのに、３つでいいのかということが思うかもしれないが、最後の引数jdは戻り値に相当するので引数として書く必要はない。そのために、サブルーチンJulianの中で、intent(Out)と指定しているのである。

実行してみる。

Julian day = 2459533.50000

きちんと呼び出されているのがわかる。

まあ、この程度なら、わざわざf2pyを使うこともない。しかし、計算によってはPythonのみとPython+Fortranでは、計算時間に大きく差が出るのは確かなことなのである。500倍の違いが出るらしい。速いっていいよね。

最後に何がとは言わないが、世の中には早すぎて嫌われちゃうこともあるよね、という言葉を残して終わる。

老猫こてつ.f90