さてさて、面白いページを見つけた。
ここにこんな事が書いてある。
他の事はどーでもよい。
問題はこれをどうやってプログラミングするか、である。
なかなかの難問じゃないか?
もちろん、件のページにCで書いたソースコードは紹介されている。しかし、アレを読んだトコで上の「計算アルゴリズム」とは繋がらないだろう。何故なら人間コンパイラが数式をグチャグチャにしてしまったから、だ。
例えると、機械語の数値羅列を読んで、数学的に記述された式を想像しなさい、と言う無茶振りと事実上あんま変わらん、って事だ。いくら説明されても数式と「コード」を関連付けるのは難しいだろう。
Racket等のLisp処理系でこれに手を付けるにはどうすれば良いのだろう?
まずは手始めに、各項の一般的な記述を考えてみる。
なお、上の数式だとk = 1から、と言う前提らしいが、Racketだと0から数えるんで、k = 0から、と項をズラす。
いずれにせよ、項の基本は
(+ 2 (* (/ (+ k 1) (+ (* 2 k) 3))
になるのは想像に難しくない。
問題は、だ。上の第二項が掛け算の対象として自己再帰的な項を利用する、と言う事だ・・・・・・。この辺の計算が面倒なのである。
取り敢えず計算対象を何であれxとしてみよう。
そうすれば、上の式は次のようになる。
(+ 2 (* (/ (+ k 1) (+ (* 2 k) 3)) x))
これで、必要最小限の「項の計算」は行える・・・xがまだ未定義だけどな。
xを自己再帰的に引数に取ろうとすると、取り敢えずはxが引数になる「関数」になるだろ、ってのは想像が難くない。
もう一つ未定義のkがある。つまり、ここで上の「項の計算」を用いて関数を定義する際、kとxの引数二つの関数になる筈・・・・・・だろうか?
kは「何番目の項なのか」を決定するが、xは単に次の項を引数として取ってるだけ、である。そしてこの内部ではkがxを決定出来ないのだ。
そして元の数式を見ても「連鎖」が重要なのだが、「何番目の項」が決定してもその他は全く決定しない。言い換えるとkとxの決定のタイミングが違う、と言う話である。
一体こういう場合はどうすれば良いのだろうか?
その答えが「高階関数」である。kは随時与えて「関数のカタチ」を決定するが、引数のxを受け取ってないので、「xを受け取る関数」を返すようにすれば良い。これでkとxの決定タイミングをズラす事が出来る。
具体的には次のように書く。
(define (term-atan k)
(lambda (x)
(+ 2 (* (/ (+ k 1) (+ (* 2 k) 3)) x))))
つまり、第0項目はkに0を与えて返ってくる関数、第1項目はkに1を与えて返ってくる関数、となる(両者ともxを引数として取る)。
確認してみよう。
> ((term-atan 0) 1)
2 1/3
> ((term-atan 1) 1)
2 2/5
> ((term-atan 2) 1)
2 3/7
>
取り敢えず簡単の為、x = 1として与えている。
確かに元々の数式の第1項、第2項、第3項が返ってきてるのが分かるだろう。
これを使えばπの公式を再帰的に書くことが可能・・・な筈である。多分。
(define (π k)
(let loop ((i 0) (ini (term-atan 0)))
(if (= i k)
ini(loop (+ i 1) (??????)))))
さてさて。再帰に慣れてても上のコードの??????の部分はどう書けば良いのか、ちと悩むのではないだろうか。
と言うのも、この再帰、外側の関数が決定しても内側の関数は決定されない、って状態の連鎖なのである。これが困ったちゃんなのだ。
(決定 (決定 (決定 (決定 (......... 非決定)......))))
ところが、フツーのプログラミング言語のルールだと引数が決定してから関数が適用される。故に一番内側が決定してない状態はマズいし、翻るとk = 1の状態で、k = 0の関数がk = 1の関数に「だけ」適用されても、その時点ではk = 1の関数が何を返すのかは決定していない。やっぱりおかしな事になるのだ。
加えると、関数term-atanが返す関数の引数は数を仮定し、返り値も数である。当然「何かが決定しない限り」その関数は数を返さず、関数として返る。従ってxに関数が受け渡された場合、引数と返り値の間で型の整合性が崩れるのでエラーとなるのだ。
;; ini も関数なんで関数に関数を適用して・・・ってのはこの時点では上手く行かない(define (π k)
(let loop ((i 0) (ini (term-atan 0)))
(if (= i k)
ini
(loop (+ i 1) (ini (term-atan (+ i 1)))))))
一体どうすんべぇ、って、こんな時に使うべき関数がcomposeである。
関数composeも高階関数で、こいつは合成関数を返す関数だ。
Racketでは組み込みだが、他のScheme処理系でも、仕様にはないが、組み込みの可能性が高い。
また、仮になくてもポール・グレアムのOn Lispにその実装法が載っている。
;; Schemeでの実装例(require srfi/1) ;; ここは実装依存だがsrfi-1くらい使える筈
(define (compose . fns)
(if fns
(let ((reversed (reverse fns)))
(let ((fn1 (car reversed))
(fns (reverse (cdr reversed))))
(lambda args
(fold-right (lambda (x y) (x y)) (apply fn1 args) fns))))
(lambda (x) x)))
composeを使えば合成関数を再帰的に作り続けてくれる。
従って、実際に
と言う状態を作り出してくれるのだ。
よって、全体として作成する再帰関数は、
(define (π k)
(let loop ((i 0) (ini (term-atan 0)))
(if (= i k)
ini
(loop (+ i 1) (compose ini (term-atan (+ i 1)))))))
となる。
実際、またもやxとして1を与えて、項の累積がどうなってるのか見てみよう。
> ((π 0) 1)
2 1/3
> ((π 1) 1)
2 4/5
> ((π 2) 1)
2 104/105
>
確かに(2 + 1/3)、(2 + 1/3*(2 + 1/5))、そして(2 + 1/3*(2 + 1/5*(2 + 1/7)))の答えになっている。
どうやら成功したようだ。
ところですっかり忘れてたのかもしれないが、僕らはまたもや円周率πを計算している。
じゃあ、2と104/105はいくらくらいなんだろう。
> (exact->inexact ((π 2) 1))
2.9904761904761905
>
実は3にも届いていない。
一つの理由として、分数の状態を簡単に分かりやすくするため、xに1を与え続けてきたわけだが。
実際のトコ、
の部分はk -> ∞の時には1/2に収束する。そしてこれが「次の項」の係数になるのだが、この項を捨てて2のみ使うか、あるいはx = 2 + 1/2 = 5/2にするかは好みだろう。
いずれにせよ、その辺だと初期値周辺でもそこそこ円周率に近くなる。
> (exact->inexact ((π 2) 2))
3.0476190476190474
> (exact->inexact ((π 2) 5/2))
3.0761904761904764
>
kを一つ上げてみよう。
> (exact->inexact ((π 3) 2))
3.0984126984126985> (exact->inexact ((π 3) 5/2))
3.111111111111111
だんだん円周率っぽくなってきた。
そして、どうやら見ていくと、以降、kの値が4つ増える毎に桁数の精度が1上がる模様である。
> (exact->inexact ((π 7) 2))
3.139469680646151
> (exact->inexact ((π 7) 5/2))
3.1400547165253045
> (exact->inexact ((π 11) 2))
3.1414796489611403
> (exact->inexact ((π 11) 5/2))
3.1415099430713913
>
ちなみに、Schemeのexact->inexactと言う関数は分数等の「正確な数」を「不正確な数」である小数に直すわけだが。
そう、実は、最初見た通り、実はこれらの計算は元々分数で行われてるのである。
> ((π 11) 5/2)
3 47354789/334639305
>
つまり、47354789/334639305の部分が小数であり、ここの分子を例えば1000倍すれば整数部分が「円周率の小数点以下での正確な確定数」って事になってる。
> (* 1000 47354789/334639305)
141 34129399/66927861
>
141が「3.141...」の小数点以下3桁での決定部分だ、って事だよな。
考えてみると、「分数を使えないプログラミング言語」の場合、一々継続的に桁数繰り上げて出力、を繰り返さなアカンのだけど、「分数を使える」のなら一気に分子に必要なだけ10の累乗を掛けて・・・例えば1万桁表示させたいのならk = 40,000で計算して3引いて分子に10 ^ 10,000かけて出てきた整数読めばエエんちゃうの?って事になる。
このブログで20,000桁表示と言うのは無理な為(そもそも、記事の文字上限がたった30,000文字)取り敢えずスクショで。件の計算の分数部分を10^10000倍して冒頭60桁は取り敢えず合ってる。
「ふむ。やっぱり分数を扱えるプログラミング言語は超強力だな。」
と言う事は分かったが、同時に「小数点以下☓桁まで出力せよ」と言う題意は、分数を扱える以上、理論的にクソ簡単な事が発覚し、計算時間がかかって出力にも手間取る割にはあまり面白くない議題だ、って事がバレたので、急激に興味を失ってしまった・・・・・・そう、Cなら難しくてもLispなら単に面白くないだけの問題なのだ。
よって唐突にここで終わろうと思う。
