python、区分求積法を用いて円周率πを求めるプログラム

Question

python2.7で区分求積法でπを求めるプログラムを作成したのですが、精度と時間があまりよくないです。
改善するとしたらどうすればいいでしょうか？
（申し訳ありませんがインデントがうまくいかないため"_"を使っています。エディタ等で置換してください。）
回答よろしくお願いします。m(_ _)m
________________________________________________________________________

#-*- coding: utf-8 -*-
def cal_pi():
_r1 = 0
_r2 = 0
_for i in xrange(1,1000):
__r1 += mysqrt(1000000-i**2)
_for k in xrange(0,999):
__r2 += mysqrt(1000000-k**2)
_r = float((r1+r2)/500000)
_return r

def mysqrt(x):
_#初期値の設定
_t = 1
_i = -1
_while x > t:
__t = t*100
__i += 1
_r1 = 5*(10**i)
_#Newton_method
_prer = r1 
_r = float(r1 - ((r1**2)-x)/(2*r1))
_while round(r,10) != round(prer,10):
__prer = r
__r = float(r - ((r**2)-x)/(2*r))
_return r

cametan_42 · Accepted Answer

> とても速く正確で驚きました！

実は旧いタイプのプログラマ(例えばC/C++、Java、Pascal、VBA等を得意とする)なんかはリスト内包表記を「分かりづらい」と言うんで嫌う傾向があるんですね。
「リスト内包表記」はHaskellと言う、「関数型言語」と言われるパラダイムに属する言語から借りてきた機能です。

Haskell(リストとリスト内包表記):
https://ja.wikipedia.org/wiki/Haskell#.E3.83.AA.E3.82.B9.E3.83.88.E3.81.A8.E3.83.AA.E3.82.B9.E3.83.88.E5.86.85.E5.8C.85.E8.A1.A8.E8.A8.98

んで、確かに従来の手続き型/オブジェクト指向系の言語が備えてる機能とは見た目からして違う。
ただし、Pythonの場合、実はフツーにfor文でループを回すよりこっちの方が「実行速度が速い」んです。つまり、リスト内包表記には最適化が施されているんですね。
(これが意味する事は、Pythonの設計者がユーザーに「リスト内包表記をどんどん使って欲しい」って思ってる、って事です)
まあ、数値計算相手だけ、じゃなくって対象のリストは色々と考えられるんですが、特に今回のケースのように

* 数値の列をリストとして表現可能で、リストの各要素に「何らかの」計算を施さないといけない

場合は、取り敢えず「リスト内包表記の使用」ってのは一考の価値があります。
先に書いた通り、「フツーにforループを書いたコードより高速に動作する」んで、明らかに、「コードが短くなる」以上に実用的なメリットがある、って事だからです。

> sum関数は型変換も自動で行っているということでしょうか？

うーん、型変換・・・考えた事無かったですね(笑)。

>>> sum([1, 2, 3, 4, 5])
15
>>> sum([1.0, 2.0, 3.0, 4.0, 5.0])
15.0
>>>

まあそのまま、です(笑)。
ええと、元コード見る限り、恐らくC系とかの「静的型付け」言語の学習経験があるのかな、とか思ったんですが、Pythonの場合、「静的型付け」ではなくって、「動的型付け言語」って言われるモノに属しています。
これらの言語の特徴は、「実行時に型が決まる」って特徴があるんですね(他にはLisp、Ruby、JavaScript等がある)。
これは「型がない」わけではなくって、型はあるんですが、明示しなくても良い。コードが走る「その時点で」型が決まる、って特徴があるんです。
従って、sumにしても、エラーになるかならないか、「走ってみなけりゃ分からない」って事なんですね(笑)。いずれにせよ、コードを書いてる時点では、「型」は忘れて構わない・・・フツーは、ですね(笑)。

実は、Python3では改善されてるんですが、Python2.xだとちょっと困る事があるんですね。それはこれです。

>>> 1/2
0
>>>

1を2で割ると0になる・・・。これがPython2.xだと困るトコで、Python3以降だと改善されてるトコなんですね。これは答えだと0.5とかそれに類する値が欲しいトコなんですが、これがちょっと「出来の悪いC言語」みたいな値を返してくる(笑)。
Python3以降だと次のようになります。

>>> 1/2
0.5
>>>

Scheme(Lisp言語の一種)だと次のようになります。

> (/ 1 2)
1/2
>

いずれにせよ、「型を意識しなくて良い筈の」動的型付け言語で、「割り算」の時だけ型を意識しなくちゃいけない、ってのは結構不都合なんですよね。
Python2.xの場合、この手の問題を避ける時にやるのが(動的型付けなのに!・苦笑)、C言語で言う「キャスト」ですね。型変換。floatの出番、なんつーのはこの時くらいしかない。

>>> 1/float(2)
0.5

あるいは、僕自身は、参考コードとかであんまfloat使うの見たことないんですが、良くあるケースでは次のように書いて「お茶を濁す」ケースが殆どです。

>>> 1/(2 * 1.0)
0.5

つまり、ワザと「割る数に1.0(float型)をかけてfloatに直して」割る、と言う、結構みっともない(笑)ハックですよね。
こーゆーのがPython2.x系の「欠点」と言えば「欠点」なんです。

さて、そうなると、「どこでキャストするのか」ってのが問題になってきます。
例えば1/2を計算したいとして、

>>> 1/2
0

になっちゃう、ってのは既に見ました。
問題は、です。

>>> float(1/2)
0.0

これをついついやっちゃうんじゃないか、って事ですよね。1/2は元情報がInt(整数型)なんで当然1/2も整数型にならざるを得ないんで0を返します。
この「0をfloatにキャストしても」結果はfloatになるだけで計算結果としては役立たずなんですよ。
つまり、安全策としては「割る数をキャストしないと」意味を成さない、って事です。

>>> 1/float(2)
0.5

「割った後にキャストしても」ダメなんですよね。
さぁ、自分のコードを見なおしてみて下さい。「割った後に」floatに型変換してませんか?「割る数を」キャストしないと、「既に手遅れ」なんです。

Python2.x系だと、この「割る数の」型だけ気をつければ、通常は「型」の事なんざ忘れててもオーケーです。
ドツボになるのは「割り算の時だけ」です。
(しかもPython3.x系だともはやこの辺さえも忘れて大丈夫です)

cametan_42 · Answer

本当に「区分求積法」でイイんですか?

いや、別に区分求積法は区分求積法でイイんですが、要するに精度そのものは良くないんで、期待に沿わないかもしれません。

一応定義を確認してみますか。
次のページの通りでイイんですね?

区分求積法の基本式:
http://w3e.kanazawa-it.ac.jp/math/category/sekibun/henkan-tex.cgi?target=/math/category/sekibun/kubun-kyuuseki-hou.html

と言うのも、この辺の「数値計算による求積」ってのは色々なテクニックっつーかアルゴリズムが提案されてるんで、そういう意味だと区分求積法自体は一番原始的でそんなに精度が良くないから、です。

提示されたコードcal_pyだと結果はこうなりますよね。

>>> cal_pi()
3.664502
>>>

さて。ところで。実は。
「区分求積法」自体のアルゴリズムはPythonだと「たった二行で」実装出来ちゃいます(笑)。
ウソッ、って思うかもしれませんが本当です。Pythonと言うのはそういうプログラミング言語なんです。
次のコードが・・・まあ、簡単の為に求積区間を0〜1としますが、Pythonで書かれた区分求積法のアルゴリズムです。

def quadrature(func, n):
　　delta_x = 1 /(1.0 * n)
　　return sum([func(x)*delta_x for x in [k/(1.0*n) for k in range(n)]])

以上です。
まあ、3行ありますが、本体はたった二行ですね。これで区間0〜1のどんな関数でも区分求積してくれます。
何故なら、このコードは、基本的には、「数学的式」をそのままPython表現に置き換えただけ、ですから。

元のコードでxrange関数は使ってるみたいなんで、range関数の説明は省きますが、ここでは「リスト内包表記」ってカタチで次のようなトリックを使ってます。

リスト内包表記:
http://docs.python.jp/2/tutorial/datastructures.html

例えば n = 10の時、

[k/(1.0*10) for k in range(10)]

ってのは

[0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9]

と言うリストを返してくれます。
これはさっきの参照サイトで言うと、k_xって言う「0〜1をn(この場合は10)分割した際のx座標」のリストですね。
これを次々と与えられたf(x)と言う関数にぶち込んで行けばいいわけです。
そこでこのk_xリストをまたもやリスト内包表記を使って次のように操作します。

[func(x) * delta_x for x in [k_x]]

ここでk_xはさっきのリスト(と言うかリスト内包表記)です。delta_xは分割区間の大きさ(つまり面積で言うところの高さ)なんで、例えば今回みたいに10分割すれば1/10で0.1ですよね。
つまりk_xってリストの要素に順次func(x) * delta_xと言う要素が「適用」されていきます。

10分割の場合:
func(x) * 0.1と言う計算が[0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9]の各要素に適用される(今のトコ、func(x)が何の計算をするのかは問わない)

そして最後にsum、です。
これもPythonでは基本中の基本な関数でリストの要素を全部足し合わせます。

例:  sum([1, 2, 3, 4, 5]) -> 15

つまり数学的に言うとまさしくΣで、結局この上の段にあったf(k_x)*Δxと言う「計算結果」を全部足し合わせるので、区分求積法は一丁終了、って事ですね。
数式をそのまま意図するトコに置き換えただけなんで簡単でしょ(笑)?殆ど「プログラミング」なんざしてねぇのも同然ですね(笑)。

さて、定義したquadratureってコードは引数に関数funcと分割数nを取ります。分割数nはいいとして、関数funcを取る、ってのがこのコードのポイントで、こういう「関数を引数として取る関数」を一般に「高階関数」と言います。
こういうのはちょっと古い言語じゃ出来ないんですが、Pythonを初めとして、最近の新しい言語では割に「当たり前に」付いてる機能ですね。
こういうイカした機能は使わないと損です。そして、そのお陰で一般性を獲得出来るわけです。

現時点、要するに単位円の面積を求めればいいわけですよね。
単位円は x^2 + y^2 = 1^2なわけですけど、今f(x)のカタチに直すのなら、必要なのは次のような関数になります。

f(x) = √(1- x^2) ... ただし、f(x) ≧ 0 となる

これをPythonでプログラムすると、

import math
def circle(x):
　　return math.sqrt(1-x**2)

になりますね。・・・オリジナルのコードだとルート取るのも自分でやってんのかしら?
まあ、ここではメンド臭いんで(笑)、取り敢えずPython組み込みのmathライブラリを使います(自作が良いのなら、あとで自分で試験してみて下さい)。

さて、単位円の関数もある、区分求積法の関数もある、ってんで、ためしに使ってみましょう。
分割数nは・・・取り敢えず1000000で行ってみますか。

>>> quadrature(circle, 1000000)
0.7853986631034552
>>>

あれれれ、0.78とか出てきちゃったよ・・・当然ですね。何せ「単位円の1/4の面積」しか求めてないんで。当たり前です。
つまり、単位円の面積である円周率π、を求めるにはこれを4倍せなアカン、って事ですね。

>>> 4 * quadrature(circle, 1000000)
3.1415946524138207
>>>

まあまあ、円周率らしき数が出てきましたね。
Pythonの組み込みだと

>>> math.pi
3.141592653589793
>>>

なんで近いような遠いような・・・・・・。
まあ、この辺が区分求積法の限界ですかねぇ。

python、区分求積法を用いて円周率πを求めるプログラム

> とても速く正確で驚きました！

本当に「区分求積法」でイイんですか?

似たような質問が見つかりました

関連するカテゴリからQ&Aを探す

デイリーランキングこのカテゴリの人気デイリーQ&Aランキング

マンスリーランキングこのカテゴリの人気マンスリーQ&Aランキング