日記 『読書の日、数学的羽田』
昨日は読書の日だった。 電車で読書、羽田まで。 羽田を見物、それから読書。 何を読んだかって? 数学者、藤原正彦の「天才の栄光と挫折」 MASAさんのブログで、天才数学者のシュリニヴァーサ・ラマヌジャンを知る。 『MASAラボ 鸚鵡の会議は白昼夢』 http://plaza.rakuten.co.jp/nakabisya/diary/201111020000/#comment それをアマゾンで注文したら、注文したはずないのにもう一冊注文になっていた本も~。『催眠誘導の極意』ナンダコレって感じの。 数学の公式って・・・・覚えられないし苦手だけれど、一度そのラビリンスに紛れ込んだら、ハマりこんでしまう、宇宙迷路のようなものかしら? たまたま、本日二男が来たので、彼のミクシーのブログを見せてもらうと、二男もパイの公式を考えたのだとか? なんだかよくわからないけれど、分かる方? 観てね~。 天才とはほど遠い我が二男のバイを考えた数式だそう。 息子はIT系の会社員です~。 数学が得意な方しか関心は無いと思いますが・・・。 * 超益体もない日記(part 1)円周率に関する考察というかほぼ備忘録 昨日のことだけど、 酔っ払った帰りの電車の中で 何故か円周率の計算を考え始めてしまった。 なんか調べたら負けな気がしていて(酔っ払ってたからな!) いろいろ考えてしまって 家に帰ってから適当な紙に計算書いて 円周率計算するプログラム作ってしまったので 折角なので日記に書いておきます。 ━━━━━━VBScript━━━━━━━ dim dx dim x dim sum dx = 0.0000001 sum = 0 for x = 0 to 1 step dx sum = sum + sqr(1-x^2) * dx next sum = sum * 4 msgbox sum ━━━━━━━━━━━━━━━━━ こう書いて見るとかなり簡単なプログラムだよね。 まぁ1000万回計算してこの精度じゃあんま実用的じゃないけど。 ちなみに実行すると、「3.14159285415673」って結果になります。 原理としてはこんな感じだ! まず、πを求める式から考えるんだ。 たぶん皆良く知ってると思うけど、円の面積の式って「πr^2」じゃん? んで、単位円(半径1の円)の場合、体積をSとすると S = π * 1^2 S = π になるんだよね。 つまり、半径1の円の体積を求めればπが出るってことだ。 次に、πを使わないで単位円の面積の求め方を考える。 まず考えたのは、積分による求め方で、上のプログラムもそれで書かれてるんだけど 単位円の中心を原点としたグラフで x > 0 , y > 0 の範囲で考える。 この範囲の体積出たら最後に4倍すればいいからね。 んで、この扇形のグラフのxとyの関係を考えるんだけど、 原点~x,0と、原点~x,yとx,0~x,yの地点で線を引くと常に直角三角形が出来るよね。 この時、 原点~x,0の長さ = a x,0~x,yの長さ = b 原点~x,yの長さ = c とすると、三平方の定理から、 a^2 + b^2 = c^2 が成り立つはずだ。 グラフ上の数値と変数に置き換えるとこうなる。 x^2 + y^2 = 1^2 んで、積分で求める時に必要になるのはyなんで y^2 = 1 - x^2 y = √(1 - x^2) としておく。 で、yがxから出せるようになったから、後は適当に関係を纏めると 1 S = ∫(√(1 - x^2) * Δx) 0 な感じになるよね。 高等数学とかあんまわかんないから、書き方とか若干違うかもしらんけども。 んで、これは円の1/4の扇形の面積だから、最後に4倍して、 1 S = 4∫(√(1 - x^2) * Δx) 0 がπになるわけだ。 だからこれをプログラム的なアルゴリズムに直すと、最初に書いたようなプログラムになるんだ。 ちなみにプログラム上のdxはΔxのことね。 でもね、これじゃあまりに計算に時間かかり過ぎるんだ。 ってことで、ちょっと考え直して見た。 まず、円なんて同じ三角形の集合で近似値が算出できるはずなのに このやり方じゃ4分割までしかできない。 これを無限に小さい鋭角な三角形の集合として計算できないものかと考えたんだ。 まず、それに当たって、考える三角形の形を変える。 つまり以下のような形だ。 原点~1,0の長さ = a 1,0~x,yの長さ = b 原点~x,yの長さ = c 更にこの時、最終的に何個あれば円になるのか分からないといけないから 三角関数を使って角度から求められるように考える。(角度をθとする) これを考える上で、上の三角形を x,y~x,0の線で二つの三角形に分ける。 するとそれぞれの三角形の面積は S1 = (sinθ * cosθ) / 2 S2 = (sinθ * (1 - cosθ) / 2 になるよね。 これを加算して適当に纏めると以下のような感じ。 S = (sinθ * cosθ) / 2 + (sinθ * (1 - cosθ) / 2 S = (sinθ * (cosθ + (1 - cosθ)) / 2 S = (sinθ * 1) / 2 S = sinθ / 2 こうすると大分簡単な式になるよね。 これに気づいた時、 俺天才なんじゃね!?!?! とか思って舞い上がってしまったわ。 ほんと酔っ払いだった。 ともあれ、これで角度から円上の扇形の近似値を求める式が出来た訳だな。 んで、後は何個集めたら円になるかの必要個数を求める式だけど これは簡単だよな。 必要個数 = 360 / θ で簡単に出る。 だからこれを組み合わせて S = (sinθ / 2) * (360 / θ) S = 360sinθ / 2θ な感じになる。 更に、無限に小さい三角形として考えるから、 S = 360sinΔθ / 2Δθ とすれば完成だよな。 実際に計算するときはθを小さい値にすれば勝手に精度は上がる。 これはプログラムにするまでも無いんで、適当に関数電卓で計算してみると (360*sin(0.0000001))/(2*0.0000001) = 3.1415926535897936 って感じだな。 仮に無限小の値を最初のプログラムと同じにして計算してみたけど 確保できた正確な桁数はこれくらい違うんだ。 最初のプログラム ⇒ 3.141592 今の計算式 ⇒ 3.141592653589793 しかも、最初のプログラムは計算にうちのPCで5秒くらいかかったけど 今の計算式の計算はブラウザ上の関数電卓で一瞬だった。 これちょっと頑張ればギネス記録狙えるんじゃね!? と思った一瞬だった。 http://science.slashdot.jp/story/11/10/18/0135247/%E8%87%AA%E5%AE%85-PC-%E3%81%A7%E5%86%86%E5%91%A8%E7%8E%87%E8%A8%88%E7%AE%97%E3%81%AE%E3%82%AE%E3%83%8D%E3%82%B9%E8%A8%98%E9%8C%B2%E3%82%92%E5%8F%96%E5%BE%97%E3%81%97%E3%81%9F%E7%94%B7%E6%80%A7%E3%80%81%E8%A8%98%E9%8C%B2%E3%82%92%E6%9B%B4%E6%96%B0 ギネスは10兆桁らしい。 世界の暇人は文字通り桁が違った_| ̄|○前回の話の続きなんだけど よく考えたら、あのままじゃギネスなんて狙いようも無いよな。 ってことに翌日になって気づいた。 そもそも、プログラム上である程度以上の桁数を算出しようとしたら どうしても型毎の桁数制限に引っかかる。 Unsigned Doubleとか使っても、せいぜい20桁くらい(?)しか正確な値求められないよね。 正直何桁か忘れたけど、最近はもっと増えてたりするんだろうか・・・。 ともあれ1000桁いけるとかはないはずだ。 まぁそんなこんな考えつつ、 前の式は 「π=(360sin(Δθ))/(2Δθ)」 みたいなんだったと思うけど、 この式一発で出そうとするようじゃダメだよね。 ってことでもう、考えるだけじゃ答えは出ないんで調べてみた。 (πの計算アルゴリズムとか調べるとつまらないから、この式を発展させる方向で) そしたらテイラー展開とかってやつが見つかったんだけど これすごいね。 1700年代に既にこの計算をしていた人もいるんだとか。 天才すぎる・・・。 こういう式を作り出す人って 何か別の世界を生きてるよね。 ゲームの中に生きる廃人ニートみたいに 現実とはかけ離れた法則に従って頭が構成されてるんじゃないだろうか。 テイラー展開の原理とか、 説明されたら一応分かるけど 目から鱗どかろか網膜ごと落ちそうな発想力だわ。 ともあれ、sinはテイラー展開で近似式に置き換えられることが分かった。 ∞ sin(x) = Σ(((-1)^n/(2n+1)!)x^(2n+1)) n=0 みたいになる。 まぁこのままじゃ近似式じゃないけど ∞を適当な数字で止めれば近似式になるよね。 でもこの式複雑だし、もっと簡単にならんかと思ったんだけど、 cos使えばいいんじゃね!? って思って調べたら ∞ cos(x) = Σ(((-1)^n/(2n)!)x^2n) n=0 素晴らしかった。 FFTとかでも離散的コサイン変換してるって言うし、 プログラムで三角関数絡んだら とりあえずコサイン使えばいいな。 そんなこんなで、 とりあえずこれを使って最初の式を置き換えると ∞ π = 360/(2Δθ) * Σ(((-1)^n/(2n)!)(Δθ+90)^2n) n=0 な感じになるよね。 +90ってのはsinからcosに変えた際の位相調整です。 うん 美しいな まぁこれで、 桁数問題の解決は間近な気がする。 n=5 S1 = 360/(2Δθ) * Σ(((-1)^n/(2n)!)(Δθ+90)^2n) n=0 n=10 S2 = 360/(2Δθ) * Σ(((-1)^n/(2n)!)(Δθ+90)^2n) n=6 π = S1 + S2 みたいにできるからな。 と、ここまでやっといて気づいたんだけど、 Δθがある限り精度上げられねぇ_| ̄|○ またなんか考えるべ・・・ ギネス載ったる! ;;;;;;; ☆ スパコン京が、世界1位になったとか! デスネ!