100万倍速いプログラムを書く

25
いいね
0
コメント

Akai_Banana

2017年07月07日に更新

この記事はなんなの

プログラミングを始めたばかりで高速化の大枠が全くわからず意味不明なことをしていた在学時、こんな資料があったら良かったのになあ、と思って書いたもの。
書いて、在学時研究室に押し付けた後紛失したと思われていたものが発掘されたもの。
SS検索の宣伝。とっても頑張って作ったので、せめて一目見てやってください。気に入ったら、使ってやってください。感想・バグ報告

要約

ライブラリがあるならそれを使う。

ライブラリが無ければ、ボトルネック部分を探してそこだけ高速な言語で書きなおすか、可能なら事前コンパイルする。

最初から全てを Low-Level な言語で書くと大変、でも結果のプログラムは速い。

以下の時間の計測ではインポートにかかる時間は除いています。

使用するもの

Python（3系）
Numba
Scipy
Line Profiler
Fortran（gfortran）
QUADPACK

QUADPACK以外の導入方法の説明は色んな所にあるので各自でお願いします。上3つに関しては、個人的にはAnacondaをおすすめします。QUADPACKの導入は当記事中で簡単に説明します。簡単に（嘘）。

プログラミングターゲット

任意の関数 $f(x)$ について、

$$ \left( \prod_{m=1}^{10} \left( \int_{0}^{1} \left( f \left( x \right) \right) ^{m} {\rm d}x \right) \right) ^{-1} $$

つまり、

$$
\frac{1}
{
\int_{0}^{1} \left( f \left( x \right) \right) ^{1} {\rm d}x \cdot
\int_{0}^{1} \left( f \left( x \right) \right) ^{2} {\rm d}x \cdots
\int_{0}^{1} \left( f \left( x \right) \right) ^{10} {\rm d}x
}
$$

の値を求めるプログラムを作りましょう。

自力でPythonで書く（8桁491ミリ秒）

式中で積分以外に悩む部分は無いので、どうやって積分を計算すればいいのかを考える事にします。

図に示すように、積分は複数の長方形の和で近似することができます。

Image from http://commons.wikimedia.org/wiki/File:Riemann_sum_convergence.png
Lisence: CC BY-SA http://creativecommons.org/licenses/by-sa/3.0/deed.en

高校で習う公式をそのまま使えば、

$$ \int_{0}^{1} f \left( x \right) {\rm d}x=\lim_{n\rightarrow \infty } \left( \sum_{k=0}^{n-1}{\frac {1}{n}f \left( {\frac {k}{n}} \right) } \right) $$

です（Σ の n-1 は図と辻褄を合わせるために導入したもの）。つまり、元の式中の積分は

$$ \int_{0}^{1}! \left( f \left( x \right) \right) ^{m}\,{\rm d}x \approx \sum_{k=0}^{10000-1}{\frac{1}{10000} \left( f \left( {\frac{k}{10000}} \right) \right)^{m}} $$

として近似することができます。ここで n=10000 は適当な大きな値です。これを元の式に代入すると、プログラムとしては

$$ \left( \prod_{m=1}^{10} \left( \sum_{k=0}^{10000-1}{\frac{1}{10000} \left( f \left( {\frac{k}{10000}} \right) \right)^{m}} \right) \right)^{-1} $$

という値を求めればいいことになります。

以下では、プログラムが正しく動いていることを確認するために、手計算でも値が求まる関数として $f(x)={\rm sin}(\pi x)$ を使います。このとき、これはウォリス積分というものになり、厳密な値は

$$ \left( \prod_{m=1}^{10} \left( \int_{0}^{1} \left( {\rm sin}\left(\pi x \right) \right)^{m} {\rm d}x \right) \right) ^{-1}=120 {\pi}^{5} $$

であり、これは 36722.362174233774391... となります。

つまり、出来る限り高速・正確に 36722.362174233774391... を数値積分で得るのが私達の目的になります。ここで、勿論、$ f(x) $ は任意の関数に変更できなければなりません。

さて、これを求めるプログラムを例えばPythonで書くと次のようになります：

# import * は一般的にはあまりいい習慣では無いけども…
from numpy import *

def f(x):
    return sin(pi*x)

n = 10000
# range(1, 11)は間違いではない
product = prod( [sum((1/n) * f(k/n)**m for k in range(n)) for m in range(1, 11)] )
result = 1/product
print(result)
# 36722.3624763

9行目について、このような1行に詰め込む書き方は普通しません、が、数式と比べやすいようにあえてこう書いています（後で書き換えないといけなくなるのですが…）。

次のように対応します:

このプログラムの実行にかかった時間は491ミリ秒で、結果は 36722.3624763 になります。

計算時間0.5秒で厳密値と比べて8桁まで正しい答えが得られました。やったね！

愚直に精度をあげてみる（11桁53秒）

精度を上げることを考えましょう。10桁を超える精度がほしいところです。

精度を上げるには、短冊形の分割数、つまりnを増やせばいいはずです。

n を1万から100万に変えてみましょう:

from numpy import *

def f(x):
    return sin(pi*x)

n = 1000000
product = prod( [sum((1/n) * f(k/n)**m for k in range(n)) for m in range(1, 11)] )
result = 1.0/product
print(result)
# 36722.3621743

結果は 36722.3621743 となり、11桁一致しています。

単純にnを増やしただけで、さらに高い精度を得ることが出来ました。やったね！

…しかしここで問題があります。

このプログラムの実行には私のパソコンだと53秒かかりました。

実は、このプログラムはまた別のプログラムに組み込まれ、$f(x)$ の中身を変えて色んな場面で何千回も使用されることになっています。

…いや、嘘ですが、このような状況は珍しくありません。

まあとにかく、すると、この部分だけで相当の時間がかかってしまうことになります。仮に5000回使用すると丸三日ほどです。

どうにかならないのでしょうか？

Pythonは遅すぎる、ので部分的にコンパイルしてみる（12桁1.26秒）

コード中で時間を食っている部分のことをボトルネック（又はホットスポット）と言います。「プログラムの処理にかかる時間の80%はコード全体の20%の部分が占める」という 80/20 rule (パレートの法則) というものがありますが、要するに、所謂律速段階であり、ここさえ早くなれば全体が急激に高速化されます。

逆に言うと、それ以外の部分が速くなったところであまり意味はありません。

今、どこがボトルネックとなっているかはほぼ自明であり、これは100万回の繰り返し部分（これはsin(πx)を含む）になります。

（なお、自明でない場合はプロファイリングという作業を行います。このページでも後のほうでちょっとだけ行います。）

従って、基本としては、次の項で述べるように、その部分だけを高速な言語で書き直せばいいのです。

…しかし、書きなおしたり、複数の言語を使うのが面倒な場合、部分的に事前/実行時コンパイルするという方法で比較的手軽に高速化が可能な場合があります。

Pythonはインタープリタ型の言語であることを思い出します。つまり、Pythonはコードを一行ずつ機械語(010011011...）に変換（コンパイル）して動作します。

最初からある程度コンパイルしておくことは出来ないのでしょうか。

何度も繰り返し実行される部分だけでも予めコンパイルしてしまえば、繰り返しのたびにその行を逐次コンパイル、という作業がなくなって速くなるのではないでしょうか。

実際にそれをやってくれるものは存在して、Pythonの場合はCythonやNumbaというものがあります。雑な言い方をすれば、Cythonは事前コンパイラ、Numbaは実行時コンパイラとして動作します。

個人的にはCythonの方が汎用的で好みですが、その分難しいので、ここではNumbaを使って100万回の繰り返し部分を実行時コンパイルさせます。

具体的には、実行時コンパイルしたい部分を関数として部品化し、定義の上に @jit と書くだけです:

from numpy import *
# numbaを使うために必要
from numba.decorators import jit

@jit
def f(x):
    return sin(pi*x)

# Σ部分を取り出したもの
@jit
def integrate(m):
    summation = 0
    for k in range(n):
        summation += (1/n) * f(k/n)**m
    return summation

n = 1000000
product = prod( [integrate(m) for m in range(1, 11)] )
result = 1.0/product
print(result)
# 36722.3621743

（ちなみに、@ はデコレーター構文であり、Javadoc等とは全く違います。jit は just in time (compiler) の略です。
また、Numbaが内包表記（for 文を1行に縮めた書き方）に対応していないので、Σ部分は普通の for に書き直しています）

Numbaを使っているのは sin(πx) の部分と総和の部分、つまり100万回繰り返される部分だけです。

先に述べたように、他の部分の処理は大したものではなく、@jit をつけても高速化は望めません。

結果として、このコードの実行時間はたった1.26秒になりました。数十倍の短縮です！

遅い部分だけ自力でLow-LevelなFortran言語で書く（11桁716ミリ秒）

Low-Level （低レベル）の言語というのは、程度の低い駄目な言語…という意味ではなく、ハードウェアに近い言語という意味です。

Python・Ruby・Java等が

ソースコード→バーチャルマシン向け中間コード→CPU向けネイティブコード

という過程を経るのに対して、C言語などはCPU向けのネイティブコードをソースコードからいきなり生成します。結果として、ソースコードを書く人間ははハードウェアレベルの挙動を理解していることが求められる上、バーチャルマシンによる安全装置が無いため、ちょっとのミスで致命的なバグを起こしがちです。

この手の言語で圧倒的に最も有名なものはC言語でしょうが、ここではFortranというものを使用します。

Fortranを一文で表すと「書くのは（とてつもなく）面倒臭い、でも計算実行速度は最速」な言語です。

（最も的を得ていると思うFortranの説明：http://vipprog.net/wiki/prog_lang_list.html#gb564d21 ）

（非常に参考になったFortranとその他の言語の比較（英語）：http://www.fortran.bcs.org/2012/nm_talk.pdf ）

今回の場合、プログラムは非常に短く、（次の項で実際やってみるように）全部Fortranで書いてしまってもいいのですが、一般的に言って、それなりの大きさを持つプログラム全体を全てFortranで書くのは至難です。

従って、ボトルネックな部分だけFortranで書きなおすことにしましょう。

ところで、複数の言語を組み合わせるにはそもそもどうすればいいのでしょうか。

2つの言語を実際に混ぜてしまうという方法は幾つか存在します。PythonとFortran(f2py)や、RubyとC(RubyのC拡張)などです。

これらを使うと、ライブラリのところで説明しますが、2つの言語を「本当に」混ぜてしまうことが出来ます。しかし、方法はそれぞれで異なっており、色々とややこしい点が多いです。

どのような組み合わせでも使える簡単な方法は、「片方の言語から、別の言語で書かれたプログラムを単に外部プログラムとして実行し、その出力を取ってくる」という方法です。

実際にやってみましょう。

元の問題に戻ると、高速化したいのは総和の部分です。これをFortranで書きます（後のコードと統一性を図るために、Fortran77、つまり固定形式で書いています）:

      IMPLICIT NONE
      INTEGER M,K,N
      DOUBLE PRECISION SUMMATION,DN,F
      CHARACTER(2) ARG
      N = 1000000
      DN = N
      CALL GETARG(1, ARG)
      READ (ARG, *) M
      SUMMATION = 0.0D0
      DO K = 0, N-1
      SUMMATION = SUMMATION + (1/DN) * F(K/DN)**M
      END DO
      PRINT *, SUMMATION
      STOP
      END
C
      DOUBLE PRECISION FUNCTION F(X)
      IMPLICIT NONE
      DOUBLE PRECISION X,PI
      PI = 4.0D0*ATAN(1.0D0)
      F = SIN(PI*X)
      END

総和を実際に計算しているのは8-11行目であり、その他はそんなに気にしなくても大丈夫です。強いて言うなら、外部から m の値を取ってこれるようにしていることがポイントです。

これを simple.f として保存し、例えば gfortran を使って、

gfortran -O3 simple.f -o simple

とコンパイルします。-O3 は最適化オプションです。-o simple は単に出力ファイル名を指定しているだけです。

例えば、コマンドライン上で

./simple 2

とすると、0.49999999999997380、つまりほぼ0.5の値が結果として得られます。これは$\int_{0}^{1}! \left( \sin \left( \pi x \right) \right) ^{2} \,{\rm d}x$ に相当します。

つまり、ただ単に、このプログラムを他の任意の言語から実行し、出力を回収すればいいのです。

このような機能は外部コマンドの呼び出しといい、（多分）どの言語にも存在します。Pythonだとcheck_output という関数を用います。なお、simple は以下のコードを実行する場所と同じフォルダにある必要があります：

from numpy import *
from subprocess import check_output
product = prod([float(check_output(["./simple", str(i)])) for i in range(1, 11)])
result = 1/product
print(result)
# 36722.3621743

結果は同じ 36722.3621743 でやはり11桁まで等しく、実行時間は716ミリ秒でした。

ボトルネック部分をFortranで書き直すことで、1分近くかかっていたものが1秒を切るようになりました！

全て自力でFortranで書く（11桁670ミリ秒）

前述したように、この程度のプログラムならば全てFortranで書いても問題ありません（というか、そのほうが楽です）。

実際にやってみると次のようになります:

      IMPLICIT NONE
      INTEGER M,K,N
      DOUBLE PRECISION SUMMATION,PRODUCT,DN,F
      N = 1000000
      DN = N
      PRODUCT = 1.0D0
      DO M = 1, 10
      SUMMATION = 0.0D0
      DO K = 0, N-1
      SUMMATION = SUMMATION + (1/DN) * F(K/DN, M)
      END DO
      PRODUCT = PRODUCT * SUMMATION
      END DO
      PRINT *, 1/PRODUCT
      STOP
      END
C
      DOUBLE PRECISION FUNCTION F(X, M)
      IMPLICIT NONE
      DOUBLE PRECISION X,PI
      INTEGER M
      PI = 4.0D0*ATAN(1.0D0)
      F = SIN(PI*X)**M
      END

コマンドライン上で実行すると、結果は 36722.362174284935 となって、やはり厳密解に11桁一致しています。

実行時間は670ミリ秒でした。またちょっと速くなりました！

ライブラリを使ってPythonで書く（15桁2ミリ秒）

ちょっと待ってください、数値積分なんて重要そうなプログラム、誰かが既に書いて公開してくれているのでは無いでしょうか？

…そうなんです。ここまで頑張って自分でプログラムを書いてきましたが、実はそんな必要は全くないのです。

…というか、よほど優秀でない限り自分で書いてはいけないのです！

数値積分などという作業は、科学計算においてメジャー中のメジャーのような計算であり、そのための手法は研究し尽くされており、限界まで速度と精度を高めたプログラムが既に存在するのです。

ただし、そのプログラムが自分の使っている言語で使用可能なのかというのは別の問題です。まあ、自分の使っている言語に無ければ、2つ前の項のように組み合わせて使えばいいのです。

外部から呼ばれるためのある種のプログラムをまとめたものをライブラリと呼びます。

自分のやりたい作業に合ったライブラリを探すというのはそれ自体が難しいことが多いのですが、というか見つけたらそれでクリアなことも多いのですが、科学計算の場合は LAPACK, Netlib, GSL, Scipy 等がよく聞く名前です（この並べ方はとても不適切）。

最近は機械学習が流行りなので、TensorFlow や Chainer などの名前は聞いたことがある人も多いのではないでしょうか？

さて、Pythonは科学計算ライブラリが異様に充実していることで知られています。ツールボックスなしのMatlabよりは確実に充実しているレベルです。

実際に使ってみましょう。既に Scipy ライブラリが入っている場合（Anacondaからインストールした場合）は簡単です。

最初に書いたものを積分プログラム quad を用いて書き換えます:

from numpy import *
from scipy.integrate import quad

def f(x, m):
    return sin(pi*x)**m

# ちなみに、quad という関数名は Gaussian Quadrature (ガウス求積) から来ている。
# 「積分範囲、等分割以外のいい感じの分割方法があるんじゃね？あったわ。」
# という数値積分アルゴリズム。
product = prod( [quad(f, 0, 1, args=m)[0] for m in range(1, 11)] )
result = 1/product
print("{:.12f}".format(result))
# 36722.362174233764

（注: f は引数を2つ持つように変えてあり、quad は args でfの第二引数mを渡しています。あと、Scipy知ってる方は「これFortranやん」ってツッコミを入れたくなると思いますが、後述するのでちょっと待ってください。）

結果は 36722.362174233742 であり、厳密値に15桁等しく、実行時間は1.97ミリ秒でした。

…そうです！私達が頑張って書いたFortranのプログラムよりも圧倒的に速くて、かつ正確なのです！

名前の付いている計算には全て既存のプログラムやライブラリがあると考えていいと思います。

例を上げれば、簡単なものには数値微分・積分や行列の対角化、もう少しマイナーなものには最近傍探索にグラフの同型判定などがあります。

結果として、自力のFortran実装よりも圧倒的に速くて正確なPythonのコードを書くことが出来ました！

遅い部分だけライブラリを使ってFortranで書く（15桁60ミリ秒）

現在のボトルネックはどこでしょうか。

ボトルネックがすぐにわからない場合は、前述したように、プロファイリングという方法を用います。

既にいい資料がたくさんあるのでプロファイリングの詳しい説明はしませんが、Pythonだと line_profiler が便利です。

細かい使用方法は説明しませんが、例えば次のコードを Jupyter Notebook 上で実行すると、依然として95%以上の時間は quad、つまり積分に使われていることがわかります:

from numpy import *
from scipy.integrate import quad

def f(x, m):
    return sin(pi*x)**m

# Notebook を使用しない場合はここにデコレータをつける
def main():
    product = 1
    for m in range(1, 11):
        summation = quad(f, 0, 1, args=m)[0]
        product *= summation
    result = 1/product
    print(result)

%load_ext line_profiler
%lprun -f main main()

結果を抜粋します:

% Time  Line Contents
=====================
        def main():
   0.1      product = 1
   0.7      for m in range(1, 11):
  96.0          summation = quad(f, 0, 1, args=m)[0]
!! ↑ main 中の 96.0% の時間が この行で使用されている

…ということは、ここをライブラリを用いたFortranで書き換えれば更に早くなりそうです。

Fortranの数値積分ライブラリには QUADPACK というデファクトスタンダードなものがあるので、これを使うことにしましょう。

使い方はすごくややこしいです。

なお、以下の説明では最も単純な方法を用います。.a や .lib のようなライブラリファイルは生成しません。

QUADPACKは http://www.netlib.org/quadpack/ に置いてあります。このページを開くだけで少しやる気が失せるのですが、頑張ってみましょう。

使い方の説明はページ中のdocというところに書いてあります。軽く概要をつかむだけならばWikipediaのQUADPACKのページの方が分かりやすいです。

たくさんあるファイルは、全てで一つのプログラムを成しているわけではなく、それぞれで違う積分手法を実装しています。

どれを使えばいいのかというのは難しい問題なのですが、というか私も分からないのですが、わからない場合は dqags.f あたりを使っておけばいいようです。

dqags.f ではなく、その横の dqags.f plus dependenciesからファイルを落とし解凍すると、.f ファイルがワラワラと出てきます。

使おうとしてみるとわかるのですが、どこにも書いてないくせに、これだけではファイルが足りません（ダウンロード画面にBLASが必要って書いてあるけど大嘘）。

d1mach.f と xerror.f というファイルが必要です。d1mach.f は同じくnetlibのBLASのページ http://www.netlib.no/netlib/blas/d1mach.f に置いてあります（でもBLASには含まれない、ファック）。

xerror.f もnetlibにありますが、さらにこいつが他のものを呼び出したりしていて鬱陶しいので、Scipyに使われているダミーを使います：https://raw.githubusercontent.com/scipy/scipy/master/scipy/integrate/mach/xerror.f

（場所が変わった場合はScipyのGithubから頑張って探してください…）

これらの .f ファイルを全て同じフォルダに入れると、ようやく準備が整います。

積分部分をQUADPACKを用いたFortranで書いたものは次のようになります:

      IMPLICIT NONE
      DOUBLE PRECISION A,ABSERR,B,EPSABS,EPSREL,F,RESULT,WORK,PI
      INTEGER IER,IWORK,LAST,LENW,LIMIT,NEVAL,M
      DIMENSION IWORK(50),WORK(200)
      EXTERNAL F
      COMMON M,PI
      CHARACTER(2) ARG
      PI = 4.0D0*ATAN(1.0D0)
      CALL GETARG(1, ARG)
      READ (ARG, *) M
      A = 0.0D0
      B = 1.0D0
      EPSABS = 1.49D-8
      EPSREL = 1.49D-8
      LIMIT = 50
      LENW = LIMIT*4
      CALL DQAGS(F,A,B,EPSABS,EPSREL,RESULT,ABSERR,NEVAL,IER,
     *  LIMIT,LENW,LAST,IWORK,WORK)
      PRINT *, RESULT
      STOP
      END
C
      DOUBLE PRECISION FUNCTION F(X)
      IMPLICIT NONE
      DOUBLE PRECISION X,PI
      INTEGER M
      COMMON M,PI
      F = SIN(PI*X)**M
      RETURN
      END

今回の場合は元のプログラムが非常に単純なために、ライブラリを使うことでむしろ長くなってしまいましたが、一般的には圧倒的に短くなります。

これを quad_a.f などとして先ほどと同じフォルダに保存し、

gfortran -O3 *.f -o quad_a

とすると、大量のWarningが出ますが無視して、dqags を用いたFortranプログラムが完成します。ちなみに、説明した以外の.fファイルがフォルダ中にあると多分失敗するので注意してください。しかも大量のWarningの中にErrorが紛れて分かりづらいです。

これを用いるPython側のコードは、先程のものの呼び出すプログラム名を変更しただけです:

from numpy import *
from subprocess import check_output

product = prod([float(check_output(["./quad_a", str(i)])) for i in range(1, 11)])
result = 1/product
print("{:.12f}".format(result))
# 36722.362174233749

結果は 36722.362174233749 でやはり15桁等しく、実行時間は…60ミリ秒でした。

ライブラリを用いたPythonだけで書いた場合に比べて数十倍遅くなってしまいました。苦労したのに！

原因は、「外部プログラムの呼び出し自体が quad_a の実行よりも遅いから」です。

中身が空のFortranプログラムを書いてみるとわかるのですが、外部プログラムの飛び出し（というかディスクI/O）というのは実はあまり早い作業ではありません。

なんと今回の場合、「Fortranプログラムが積分を計算する時間」よりも「そのプログラムをPythonから呼び出す時間」の方が長いのです。

私の環境の場合、一度の呼び出しに6ミリ秒ぐらいかかってしまいます。上のコードでは外部プログラムを10回呼び出しているので、実行時間の殆どはこの呼び出しのためだけに使われていることになります。

勿論、例えば一回の呼び出し中で（超）巨大な行列の対角化をする場合など、計算するべき部分が明らかに呼び出し時間よりも重たくなった場合は、（Python中でライブラリを用いた場合に比べても）遥かな高速化が望めます。

しかし、今回のように、外部で計算させるべき量が大したものではなく、またそれを繰り返し（今の場合は10回）呼び出す必要がある場合、元の言語中でのライブラリを使うか、下手をすれば普通に自分で書いたほうが速いことすらあります。

ちなみにこれを改善するのが、先に述べた、言語を「本当に」混ぜるという行為です。

ネタばらしをすると、実は、先ほど用いたPython（Scipy）のquadとQUADPACKは中身は全く同じものなのです。

QUADPACKのたくさんの手法の中からdqagsを選んだのは、それが無難だからというよりも、quadで呼ばれている手法がこの場合はdqagsだからです（より正確にはdqagseですが、これはあまりにも面倒なので避けました）。

Scipyの正体は「CやFortranで書かれた科学計算ライブラリに、便利なインターフェースを与えてPythonの一部として使用可能にしたもの」です。f2py等を用いてPythonのCインターフェースに繋げてあり、Pythonそのものとして使えるようにしてあるため、呼び出し時間は最小限に抑えられています。ちなみに、この作業を自力で行うのは非常に難しいので普通の人にはあまりおすすめしません、私も出来ません。

しかし、それでも呼び出し時間（と様々なオーバーヘッド）が0になるわけではなく、やはり100%Fortranで書かれたものには勝てません。

次の項でそれを見ます。

ライブラリを使って全てFortranで書く（15桁32マイクロ秒）

実行時間だけを考えるのならばこれが最適解です。実際にやってみましょう。

上で作ったフォルダから quad_a.f を取り除き、下記のプログラムを quad_b.f などとして保存し、同じように

gfortran -O3 *.f -o quad_b

とコンパイルします:

      IMPLICIT NONE
      DOUBLE PRECISION A,ABSERR,B,EPSABS,EPSREL,F,RESULT,WORK,
     *  PRODUCT, PI
      INTEGER IER,IWORK,LAST,LENW,LIMIT,NEVAL,M
      DIMENSION IWORK(50),WORK(200)
      EXTERNAL F
      COMMON PI,M
      PI = 4.0D0*ATAN(1.0D0)
      A = 0.0D0
      B = 1.0D0
      EPSABS = 1.49D-8
      EPSREL = 1.49D-8
      LIMIT = 50
      LENW = LIMIT*4
      PRODUCT = 1.0D0
      DO M = 1, 10
      CALL DQAGS(F,A,B,EPSABS,EPSREL,RESULT,ABSERR,NEVAL,IER,
     *  LIMIT,LENW,LAST,IWORK,WORK)
      PRODUCT = PRODUCT * RESULT
      END DO
      PRINT *, 1/PRODUCT
      STOP
      END
C
      DOUBLE PRECISION FUNCTION F(X)
      IMPLICIT NONE
      DOUBLE PRECISION X,PI
      INTEGER M
      COMMON PI,M
      F = SIN(PI*X)**M
      RETURN
      END

やってることは同じなのでやはり精度は15桁です。

実行時間は…なんと32マイクロ秒（0.032ミリ秒）でした！

ちなみに、あまりにも速く終わるので、実際には中身を100万回繰り返して得た時間をその100万で割っています。

（※: PRINT の出力はnullに向けています。回数を変えるとちゃんと時間は線形スケールになったので大丈夫だと思うのですが、この場合これを考慮するべきなのかはわかりません、ので誰か教えて下さい）。

Python+Scipyに比べても、さらに数十倍の速さになりました！！

精度を落とす（11桁19マイクロ秒）

最後に、精度と速度を交換してみましょう。

（なお、色々と試しましたが、DQAGS以外を用いてちょうどいい精度を得ることは不可能なように思えます。）

精度を指定しているパラメータは EPSABS と EPSREL です。適当に、両方1にしてみましょう：

      IMPLICIT NONE
      DOUBLE PRECISION A,ABSERR,B,EPSABS,EPSREL,F,RESULT,WORK,
     *  PRODUCT, PI
      INTEGER IER,IWORK,LAST,LENW,LIMIT,NEVAL,M
      DIMENSION IWORK(50),WORK(200)
      EXTERNAL F
      COMMON PI,M
      PI = 4.0D0*ATAN(1.0D0)
      A = 0.0D0
      B = 1.0D0
      EPSABS = 1.0D0
      EPSREL = 1.0D0
      LIMIT = 5
      LENW = LIMIT*4
      PRODUCT = 1.0D0
      DO M = 1, 10
      CALL DQAGS(F,A,B,EPSABS,EPSREL,RESULT,ABSERR,NEVAL,IER,
     *  LIMIT,LENW,LAST,IWORK,WORK)
      PRODUCT = PRODUCT * RESULT
      END DO
      PRINT *, 1.0D0/PRODUCT
      STOP
      END
C
      DOUBLE PRECISION FUNCTION F(X)
      IMPLICIT NONE
      DOUBLE PRECISION X,PI
      INTEGER M
      COMMON PI, M
      F = SIN(PI*X)**M
      RETURN
      END

結果は 36722.362174170281 で11桁等しく、実行時間は19マイクロ秒となりました。

精度を犠牲に、更に計算速度を1/3ほど削ることが出来ました。

ここで、最初に作ったものと比べてみましょう。

既に述べたように、最初に作ったもの（n=100万にして精度11桁にしたもの）は、5000回実行すると丸三日必要です。仮に100万回実行するのであれば10ヶ月以上かかる計算になります。

それに比べてこれはどうでしょうか。マイクロは100万分の1の単位です。つまり、同様に100万回実行したところで、実行時間は19秒になるのです。

結果として、プログラムの実行速度はもとの300万倍近くになりました。やったね！！

# この記事はなんなの

- プログラミングを始めたばかりで高速化の大枠が全くわからず意味不明なことをしていた在学時、こんな資料があったら良かったのになあ、と思って書いたもの。
- 書いて、在学時研究室に押し付けた後紛失したと思われていたものが発掘されたもの。
- [SS検索](https://ficsearch.net/ja/) の宣伝。とっても頑張って作ったので、せめて一目見てやってください。気に入ったら、使ってやってください。[感想・バグ報告](https://twitter.com/akai_banana)

# 要約

ライブラリがあるならそれを使う。

ライブラリが無ければ、ボトルネック部分を探してそこだけ高速な言語で書きなおすか、可能なら事前コンパイルする。

最初から全てを Low-Level な言語で書くと大変、でも結果のプログラムは速い。

-----

以下の時間の計測ではインポートにかかる時間は除いています。

-----

# 使用するもの

- Python（3系）
- Numba
- Scipy
- Line Profiler
- Fortran（gfortran）
- QUADPACK

QUADPACK以外の導入方法の説明は色んな所にあるので各自でお願いします。上3つに関しては、個人的には[Anaconda](https://www.continuum.io/downloads)をおすすめします。QUADPACKの導入は当記事中で簡単に説明します。簡単に（嘘）。

# プログラミングターゲット

任意の関数 $f(x)$ について、

$$ \left( \prod_{m=1}^{10} \left( \int_{0}^{1} \left( f \left( x \right)  \right) ^{m} {\rm d}x \right)  \right) ^{-1} $$

つまり、

$$
\frac{1}
{
\int_{0}^{1} \left( f \left( x \right)  \right) ^{1} {\rm d}x \cdot
\int_{0}^{1} \left( f \left( x \right)  \right) ^{2} {\rm d}x \cdots
\int_{0}^{1} \left( f \left( x \right)  \right) ^{10} {\rm d}x
}
$$

の値を求めるプログラムを作りましょう。

# 自力でPythonで書く（8桁491ミリ秒）

式中で積分以外に悩む部分は無いので、どうやって積分を計算すればいいのかを考える事にします。

図に示すように、積分は複数の長方形の和で近似することができます。

![480px-Riemann_sum_convergence.png](https://qiita-image-store.s3.amazonaws.com/0/81610/e69a9b43-7802-2006-4461-c1fb40522c05.png)

Image from http://commons.wikimedia.org/wiki/File:Riemann_sum_convergence.png Lisence: CC BY-SA http://creativecommons.org/licenses/by-sa/3.0/deed.en

高校で習う公式をそのまま使えば、

$$ \int_{0}^{1} f \left( x \right) {\rm d}x=\lim_{n\rightarrow  \infty } \left( \sum_{k=0}^{n-1}{\frac {1}{n}f \left( {\frac {k}{n}} \right) } \right) $$

です（Σ の n-1 は図と辻褄を合わせるために導入したもの）。つまり、元の式中の積分は

$$ \int_{0}^{1}\! \left( f \left( x \right)  \right) ^{m}\,{\rm d}x \approx \sum_{k=0}^{10000-1}{\frac{1}{10000} \left( f \left( {\frac{k}{10000}} \right)  \right)^{m}} $$

として近似することができます。ここで n=10000 は適当な大きな値です。これを元の式に代入すると、プログラムとしては

$$ \left( \prod_{m=1}^{10} \left( \sum_{k=0}^{10000-1}{\frac{1}{10000} \left( f \left( {\frac{k}{10000}} \right)  \right)^{m}} \right)  \right)^{-1} $$

という値を求めればいいことになります。

以下では、プログラムが正しく動いていることを確認するために、手計算でも値が求まる関数として $f(x)={\rm sin}(\pi x)$ を使います。このとき、これは[ウォリス積分](https://ja.wikipedia.org/wiki/%E3%82%A6%E3%82%A9%E3%83%AA%E3%82%B9%E7%A9%8D%E5%88%86)というものになり、厳密な値は

$$ \left( \prod_{m=1}^{10} \left( \int_{0}^{1} \left( {\rm sin}\left(\pi x \right)  \right)^{m} {\rm d}x \right)  \right) ^{-1}=120 {\pi}^{5} $$

であり、これは 36722.362174233774391... となります。

つまり、**出来る限り高速・正確に 36722.362174233774391... を数値積分で得るのが私達の目的**になります。ここで、勿論、$ f(x) $ は任意の関数に変更できなければなりません。

さて、これを求めるプログラムを例えばPythonで書くと次のようになります：

```py3
# import * は一般的にはあまりいい習慣では無いけども…
from numpy import *

def f(x):
    return sin(pi*x)

n = 10000
# range(1, 11)は間違いではない
product = prod( [sum((1/n) * f(k/n)**m for k in range(n)) for m in range(1, 11)] )
result = 1/product
print(result)
# 36722.3624763
```

次のように対応します:

![code_eq.png](https://qiita-image-store.s3.amazonaws.com/0/81610/bf8daf97-d1f2-42ab-cf08-48dd2b58b675.png)

このプログラムの実行にかかった時間は491ミリ秒で、結果は **36722.362**4763 になります。

計算時間0.5秒で厳密値と比べて8桁まで正しい答えが得られました。やったね！

# 愚直に精度をあげてみる（11桁53秒）

精度を上げることを考えましょう。10桁を超える精度がほしいところです。

精度を上げるには、短冊形の分割数、つまりnを増やせばいいはずです。

`n` を1万から100万に変えてみましょう:

```py3
from numpy import *

def f(x):
    return sin(pi*x)

n = 1000000
product = prod( [sum((1/n) * f(k/n)**m for k in range(n)) for m in range(1, 11)] )
result = 1.0/product
print(result)
# 36722.3621743
```

結果は **36722.362174**3 となり、11桁一致しています。

単純にnを増やしただけで、さらに高い精度を得ることが出来ました。やったね！

…しかしここで問題があります。

このプログラムの実行には私のパソコンだと53秒かかりました。

実は、このプログラムはまた別のプログラムに組み込まれ、$f(x)$ の中身を変えて色んな場面で何千回も使用されることになっています。

…いや、嘘ですが、このような状況は珍しくありません。

まあとにかく、すると、この部分だけで相当の時間がかかってしまうことになります。仮に5000回使用すると丸三日ほどです。

どうにかならないのでしょうか？

# Pythonは遅すぎる、ので部分的にコンパイルしてみる（12桁1.26秒）

コード中で時間を食っている部分のことをボトルネック（又はホットスポット）と言います。「プログラムの処理にかかる時間の80%はコード全体の20%の部分が占める」という [80/20 rule (パレートの法則)](https://ja.wikipedia.org/wiki/%E3%83%91%E3%83%AC%E3%83%BC%E3%83%88%E3%81%AE%E6%B3%95%E5%89%87) というものがありますが、要するに、所謂律速段階であり、ここさえ早くなれば全体が急激に高速化されます。

逆に言うと、それ以外の部分が速くなったところであまり意味はありません。

今、どこがボトルネックとなっているかはほぼ自明であり、これは100万回の繰り返し部分（これはsin(πx)を含む）になります。

（なお、自明でない場合はプロファイリングという作業を行います。このページでも後のほうでちょっとだけ行います。）

従って、基本としては、次の項で述べるように、その部分だけを高速な言語で書き直せばいいのです。

最初からある程度コンパイルしておくことは出来ないのでしょうか。

個人的にはCythonの方が汎用的で好みですが、その分難しいので、ここではNumbaを使って100万回の繰り返し部分を実行時コンパイルさせます。

具体的には、実行時コンパイルしたい部分を関数として部品化し、定義の上に `@jit` と書くだけです:

```py3
from numpy import *
# numbaを使うために必要
from numba.decorators import jit

@jit
def f(x):
    return sin(pi*x)

# Σ部分を取り出したもの
@jit
def integrate(m):
    summation = 0
    for k in range(n):
        summation += (1/n) * f(k/n)**m
    return summation

n = 1000000
product = prod( [integrate(m) for m in range(1, 11)] )
result = 1.0/product
print(result)
# 36722.3621743
```

（ちなみに、`@` はデコレーター構文であり、Javadoc等とは全く違います。`jit` は just in time (compiler) の略です。
また、Numbaが内包表記（`for` 文を1行に縮めた書き方）に対応していないので、Σ部分は普通の `for` に書き直しています）

Numbaを使っているのは sin(πx) の部分と総和の部分、つまり100万回繰り返される部分だけです。

先に述べたように、他の部分の処理は大したものではなく、`@jit` をつけても高速化は望めません。

結果として、このコードの実行時間はたった1.26秒になりました。数十倍の短縮です！

# 遅い部分だけ自力でLow-LevelなFortran言語で書く（11桁716ミリ秒）

Low-Level （低レベル） の言語というのは、程度の低い駄目な言語…という意味ではなく、ハードウェアに近い言語という意味です。

Python・Ruby・Java等が

ソースコード→バーチャルマシン向け中間コード→CPU向けネイティブコード

この手の言語で圧倒的に最も有名なものはC言語でしょうが、ここではFortranというものを使用します。

Fortranを一文で表すと「書くのは（とてつもなく）面倒臭い、でも計算実行速度は最速」な言語です。
 （最も的を得ていると思うFortranの説明：http://vipprog.net/wiki/prog_lang_list.html#gb564d21 ）
 （非常に参考になったFortranとその他の言語の比較（英語）：http://www.fortran.bcs.org/2012/nm_talk.pdf ）

従って、ボトルネックな部分だけFortranで書きなおすことにしましょう。

ところで、複数の言語を組み合わせるにはそもそもどうすればいいのでしょうか。

2つの言語を実際に混ぜてしまうという方法は幾つか存在します。PythonとFortran(f2py)や、RubyとC(RubyのC拡張)などです。

実際にやってみましょう。

```fortran
      IMPLICIT NONE
      INTEGER M,K,N
      DOUBLE PRECISION SUMMATION,DN,F
      CHARACTER(2) ARG
      N = 1000000
      DN = N
      CALL GETARG(1, ARG)
      READ (ARG, *) M
      SUMMATION = 0.0D0
      DO K = 0, N-1
      SUMMATION = SUMMATION + (1/DN) * F(K/DN)**M
      END DO
      PRINT *, SUMMATION
      STOP
      END
C
      DOUBLE PRECISION FUNCTION F(X)
      IMPLICIT NONE
      DOUBLE PRECISION X,PI
      PI = 4.0D0*ATAN(1.0D0)
      F = SIN(PI*X)
      END
```

総和を実際に計算しているのは8-11行目であり、その他はそんなに気にしなくても大丈夫です。強いて言うなら、外部から `m` の値を取ってこれるようにしていることがポイントです。

これを `simple.f` として保存し、例えば gfortran を使って、

`gfortran -O3 simple.f -o simple`

とコンパイルします。`-O3` は最適化オプションです。`-o simple` は単に出力ファイル名を指定しているだけです。

例えば、コマンドライン上で

`./simple 2`

とすると、0.49999999999997380、つまりほぼ0.5の値が結果として得られます。これは$\int_{0}^{1}\! \left( \sin \left( \pi x \right)  \right) ^{2} \,{\rm d}x$ に相当します。

つまり、ただ単に、このプログラムを他の任意の言語から実行し、出力を回収すればいいのです。

このような機能は外部コマンドの呼び出しといい、（多分）どの言語にも存在します。Pythonだと`check_output` という関数を用います。なお、`simple` は以下のコードを実行する場所と同じフォルダにある必要があります：

```py3
from numpy import *
from subprocess import check_output
product = prod([float(check_output(["./simple", str(i)])) for i in range(1, 11)])
result = 1/product
print(result)
# 36722.3621743
```

結果は同じ **36722.362174**3 でやはり11桁まで等しく、実行時間は716ミリ秒でした。

ボトルネック部分をFortranで書き直すことで、1分近くかかっていたものが1秒を切るようになりました！

# 全て自力でFortranで書く（11桁670ミリ秒）

前述したように、この程度のプログラムならば全てFortranで書いても問題ありません（というか、そのほうが楽です）。

実際にやってみると次のようになります:

```fortran
      IMPLICIT NONE
      INTEGER M,K,N
      DOUBLE PRECISION SUMMATION,PRODUCT,DN,F
      N = 1000000
      DN = N
      PRODUCT = 1.0D0
      DO M = 1, 10
      SUMMATION = 0.0D0
      DO K = 0, N-1
      SUMMATION = SUMMATION + (1/DN) * F(K/DN, M)
      END DO
      PRODUCT = PRODUCT * SUMMATION
      END DO
      PRINT *, 1/PRODUCT
      STOP
      END
C
      DOUBLE PRECISION FUNCTION F(X, M)
      IMPLICIT NONE
      DOUBLE PRECISION X,PI
      INTEGER M
      PI = 4.0D0*ATAN(1.0D0)
      F = SIN(PI*X)**M
      END
```

コマンドライン上で実行すると、結果は **36722.362174**284935 となって、やはり厳密解に11桁一致しています。

実行時間は670ミリ秒でした。またちょっと速くなりました！

# ライブラリを使ってPythonで書く（15桁2ミリ秒）

ちょっと待ってください、**数値積分なんて重要そうなプログラム、誰かが既に書いて公開してくれているのでは無いでしょうか？**

…そうなんです。ここまで頑張って自分でプログラムを書いてきましたが、実はそんな必要は全くないのです。

…というか、**よほど優秀でない限り自分で書いてはいけない**のです！

外部から呼ばれるためのある種のプログラムをまとめたものをライブラリと呼びます。

最近は機械学習が流行りなので、TensorFlow や Chainer などの名前は聞いたことがある人も多いのではないでしょうか？

さて、Pythonは科学計算ライブラリが異様に充実していることで知られています。ツールボックスなしのMatlabよりは確実に充実しているレベルです。

実際に使ってみましょう。既に Scipy ライブラリが入っている場合（Anacondaからインストールした場合）は簡単です。

最初に書いたものを積分プログラム `quad` を用いて書き換えます:

```py3
from numpy import *
from scipy.integrate import quad

def f(x, m):
    return sin(pi*x)**m

# ちなみに、quad という関数名は Gaussian Quadrature (ガウス求積) から来ている。
# 「積分範囲、等分割以外のいい感じの分割方法があるんじゃね？あったわ。」
# という数値積分アルゴリズム。
product = prod( [quad(f, 0, 1, args=m)[0] for m in range(1, 11)] )
result = 1/product
print("{:.12f}".format(result))
# 36722.362174233764
```

（注: `f` は引数を2つ持つように変えてあり、`quad` は `args` でfの第二引数mを渡しています。あと、Scipy知ってる方は「これFortranやん」ってツッコミを入れたくなると思いますが、後述するのでちょっと待ってください。）

結果は **36722.3621742337**42 であり、厳密値に15桁等しく、実行時間は1.97ミリ秒でした。

…そうです！私達が頑張って書いたFortranのプログラムよりも圧倒的に速くて、かつ正確なのです！

名前の付いている計算には全て既存のプログラムやライブラリがあると考えていいと思います。

例を上げれば、簡単なものには数値微分・積分や行列の対角化、もう少しマイナーなものには最近傍探索にグラフの同型判定などがあります。

結果として、自力のFortran実装よりも圧倒的に速くて正確なPythonのコードを書くことが出来ました！

# 遅い部分だけライブラリを使ってFortranで書く（15桁60ミリ秒）

現在のボトルネックはどこでしょうか。

ボトルネックがすぐにわからない場合は、前述したように、プロファイリングという方法を用います。

既にいい資料がたくさんあるのでプロファイリングの詳しい説明はしませんが、Pythonだと line_profiler が便利です。

細かい使用方法は説明しませんが、例えば次のコードを Jupyter Notebook 上で実行すると、依然として95%以上の時間は `quad`、つまり積分に使われていることがわかります:

```py3
from numpy import *
from scipy.integrate import quad

def f(x, m):
    return sin(pi*x)**m

# Notebook を使用しない場合はここにデコレータをつける
def main():
    product = 1
    for m in range(1, 11):
        summation = quad(f, 0, 1, args=m)[0]
        product *= summation
    result = 1/product
    print(result)
    
%load_ext line_profiler
%lprun -f main main()
```

結果を抜粋します:

```
% Time  Line Contents
=====================
        def main():
   0.1      product = 1
   0.7      for m in range(1, 11):
  96.0          summation = quad(f, 0, 1, args=m)[0]
!! ↑ main 中の 96.0% の時間が この行で使用されている
```

…ということは、ここをライブラリを用いたFortranで書き換えれば更に早くなりそうです。

Fortranの数値積分ライブラリには QUADPACK というデファクトスタンダードなものがあるので、これを使うことにしましょう。

使い方はすごくややこしいです。

なお、以下の説明では最も単純な方法を用います。`.a` や `.lib` のようなライブラリファイルは生成しません。

QUADPACKは http://www.netlib.org/quadpack/ に置いてあります。このページを開くだけで少しやる気が失せるのですが、頑張ってみましょう。

使い方の説明はページ中のdocというところに書いてあります。軽く概要をつかむだけならばWikipediaのQUADPACKのページの方が分かりやすいです。

たくさんあるファイルは、全てで一つのプログラムを成しているわけではなく、それぞれで違う積分手法を実装しています。

どれを使えばいいのかというのは難しい問題なのですが、**というか私も分からない**のですが、わからない場合は dqags.f あたりを使っておけばいいようです。

dqags.f **ではなく**、その横の dqags.f plus dependenciesからファイルを落とし解凍すると、`.f` ファイルがワラワラと出てきます。

xerror.f もnetlibにありますが、さらにこいつが他のものを呼び出したりしていて鬱陶しいので、Scipyに使われているダミーを使います：https://raw.githubusercontent.com/scipy/scipy/master/scipy/integrate/mach/xerror.f
 （場所が変わった場合はScipyのGithubから頑張って探してください…）

これらの `.f` ファイルを全て同じフォルダに入れると、ようやく準備が整います。

積分部分をQUADPACKを用いたFortranで書いたものは次のようになります:

```fortran
      IMPLICIT NONE
      DOUBLE PRECISION A,ABSERR,B,EPSABS,EPSREL,F,RESULT,WORK,PI
      INTEGER IER,IWORK,LAST,LENW,LIMIT,NEVAL,M
      DIMENSION IWORK(50),WORK(200)
      EXTERNAL F
      COMMON M,PI
      CHARACTER(2) ARG
      PI = 4.0D0*ATAN(1.0D0)
      CALL GETARG(1, ARG)
      READ (ARG, *) M
      A = 0.0D0
      B = 1.0D0
      EPSABS = 1.49D-8
      EPSREL = 1.49D-8
      LIMIT = 50
      LENW = LIMIT*4
      CALL DQAGS(F,A,B,EPSABS,EPSREL,RESULT,ABSERR,NEVAL,IER,
     *  LIMIT,LENW,LAST,IWORK,WORK)
      PRINT *, RESULT
      STOP
      END
C
      DOUBLE PRECISION FUNCTION F(X)
      IMPLICIT NONE
      DOUBLE PRECISION X,PI
      INTEGER M
      COMMON M,PI
      F = SIN(PI*X)**M
      RETURN
      END
```

これを `quad_a.f` などとして先ほどと同じフォルダに保存し、

`gfortran -O3 *.f -o quad_a`

これを用いるPython側のコードは、先程のものの呼び出すプログラム名を変更しただけです:

```
from numpy import *
from subprocess import check_output

product = prod([float(check_output(["./quad_a", str(i)])) for i in range(1, 11)])
result = 1/product
print("{:.12f}".format(result))
# 36722.362174233749
```

結果は **36722.3621742337**49 でやはり15桁等しく、実行時間は…60ミリ秒でした。

ライブラリを用いたPythonだけで書いた場合に比べて数十倍遅くなってしまいました。苦労したのに！

原因は、「外部プログラムの呼び出し自体が `quad_a` の実行よりも遅いから」です。

なんと今回の場合、「Fortranプログラムが積分を計算する時間」よりも「そのプログラムをPythonから呼び出す時間」の方が長いのです。

ちなみにこれを改善するのが、先に述べた、言語を「本当に」混ぜるという行為です。

ネタばらしをすると、実は、先ほど用いたPython（Scipy）のquadとQUADPACKは中身は全く同じものなのです。

しかし、それでも呼び出し時間（と様々なオーバーヘッド）が0になるわけではなく、やはり100%Fortranで書かれたものには勝てません。

次の項でそれを見ます。

# ライブラリを使って全てFortranで書く（15桁32マイクロ秒）

実行時間だけを考えるのならばこれが最適解です。実際にやってみましょう。

上で作ったフォルダから `quad_a.f` を取り除き、下記のプログラムを `quad_b.f` などとして保存し、同じように

`gfortran -O3 *.f -o quad_b`

とコンパイルします:

```fortran
      IMPLICIT NONE
      DOUBLE PRECISION A,ABSERR,B,EPSABS,EPSREL,F,RESULT,WORK,
     *  PRODUCT, PI
      INTEGER IER,IWORK,LAST,LENW,LIMIT,NEVAL,M
      DIMENSION IWORK(50),WORK(200)
      EXTERNAL F
      COMMON PI,M
      PI = 4.0D0*ATAN(1.0D0)
      A = 0.0D0
      B = 1.0D0
      EPSABS = 1.49D-8
      EPSREL = 1.49D-8
      LIMIT = 50
      LENW = LIMIT*4
      PRODUCT = 1.0D0
      DO M = 1, 10
      CALL DQAGS(F,A,B,EPSABS,EPSREL,RESULT,ABSERR,NEVAL,IER,
     *  LIMIT,LENW,LAST,IWORK,WORK)
      PRODUCT = PRODUCT * RESULT
      END DO
      PRINT *, 1/PRODUCT
      STOP
      END
C
      DOUBLE PRECISION FUNCTION F(X)
      IMPLICIT NONE
      DOUBLE PRECISION X,PI
      INTEGER M
      COMMON PI,M
      F = SIN(PI*X)**M
      RETURN
      END
```

やってることは同じなのでやはり精度は15桁です。

実行時間は…なんと32マイクロ秒（0.032ミリ秒）でした！

ちなみに、あまりにも速く終わるので、実際には中身を100万回繰り返して得た時間をその100万で割っています。
 （※: `PRINT` の出力はnullに向けています。回数を変えるとちゃんと時間は線形スケールになったので大丈夫だと思うのですが、この場合これを考慮するべきなのかはわかりません、ので誰か教えて下さい）。

Python+Scipyに比べても、さらに数十倍の速さになりました！！

# 精度を落とす（11桁19マイクロ秒）

最後に、精度と速度を交換してみましょう。

（なお、色々と試しましたが、DQAGS以外を用いてちょうどいい精度を得ることは不可能なように思えます。）

精度を指定しているパラメータは `EPSABS` と `EPSREL` です。適当に、両方1にしてみましょう：

```fortran
      IMPLICIT NONE
      DOUBLE PRECISION A,ABSERR,B,EPSABS,EPSREL,F,RESULT,WORK,
     *  PRODUCT, PI
      INTEGER IER,IWORK,LAST,LENW,LIMIT,NEVAL,M
      DIMENSION IWORK(50),WORK(200)
      EXTERNAL F
      COMMON PI,M
      PI = 4.0D0*ATAN(1.0D0)
      A = 0.0D0
      B = 1.0D0
      EPSABS = 1.0D0
      EPSREL = 1.0D0
      LIMIT = 5
      LENW = LIMIT*4
      PRODUCT = 1.0D0
      DO M = 1, 10
      CALL DQAGS(F,A,B,EPSABS,EPSREL,RESULT,ABSERR,NEVAL,IER,
     *  LIMIT,LENW,LAST,IWORK,WORK)
      PRODUCT = PRODUCT * RESULT
      END DO
      PRINT *, 1.0D0/PRODUCT
      STOP
      END
C
      DOUBLE PRECISION FUNCTION F(X)
      IMPLICIT NONE
      DOUBLE PRECISION X,PI
      INTEGER M
      COMMON PI, M
      F = SIN(PI*X)**M
      RETURN
      END
```

結果は **36722.362174**170281 で11桁等しく、実行時間は19マイクロ秒となりました。

精度を犠牲に、更に計算速度を1/3ほど削ることが出来ました。

-----

ここで、最初に作ったものと比べてみましょう。

それに比べてこれはどうでしょうか。マイクロは100万分の1の単位です。つまり、同様に100万回実行したところで、実行時間は19秒になるのです。

結果として、プログラムの実行速度はもとの300万倍近くになりました。やったね！！