(1/3) へ (3/3) へ

C# で解く「なぜ関数プログラミングは重要か」(2/3)

2015-01-26 (鈴)

3. 「4. プログラムの貼り合せ」について

「なぜ関数プログラミングは重要か」§4「プログラムの貼り合せ」は,遅延評価のもとで

    g (f input)

の f, g が input に対してパイプライン的な動作をすることを利用し,プログラムを生成器と選択器のモジュールに分割して構成する方法を説明している。 生成器の出力は必ずしも停止しなくてよい。 結果を得るために必要なだけの出力が使われる。

もしここで1回の演算の対象がデータ列の1個の要素であり,それをプログラム中で1回しか参照しないと限定できるならば,そうしたプログラムは LINQ で自然に表現できる。 例えば,上記の f, g がデータ列に対してそれぞれなんらかの変換と選択をする関数であると仮定したとき,f1 を1要素についてその変換をする関数,g1 を1要素についてその選択をする述語関数とすれば LINQ で次のように表現できるだろう。

    input.Select(f1).Where(g1)

しかし,そのように限定できない場合は LINQ の背後にある IEnumerator<T> に立ち戻って補助的な関数やクラスを定義しなければならない。 本章では,そのような定義を導入しながら「なぜ関数プログラミングは重要か」の §4 の内容を C# で再現する。

3.1 「4.1 ニュートン-ラプソン法による平方根」と FindPair<T>

「なぜ関数プログラミングは重要か」§4.1「ニュートン-ラプソン法による平方根」の repeat 関数は次のように無限ループを使って書ける。 これは任意の型の初期値 a に対する a, f(a), f(f(a)), f(f(f(a))), ... の無限列を生成する。

    static IEnumerable<T> repeat<T>(Func<T, T> f, T a) {
        for (;;) {
            yield return a;
            a = f(a);
        }
    }

しかし within 関数は列の隣り合う2個の要素について述語関数を適用する必要があるから,既成の LINQ メソッドではうまく表現できない。 そこでまず,列の1個1個の要素を取り出す IEnumerator<T>MoveNext()Current をじかに使って次のような FindPair<T> メソッドを定義する。

    // seq から pred(A, B) が成り立つ隣り合った項 A, B を見つける。
    static Tuple<T, T> FindPair<T>(this IEnumerable<T> seq,
                                 Func<T, T, bool> pred) {
        using (var e = seq.GetEnumerator()) {
            if (! e.MoveNext())
                return null;
            T a = e.Current;
            if (! e.MoveNext())
                return null;
            T b = e.Current;
            for (;;) {
                if (pred(a, b))
                    return new Tuple<T, T>(a, b);
                a = b;
                if (! e.MoveNext())
                    return null;
                b = e.Current;
            }
        }
    }
今まで抽象度の高いコードばかりだったのが,ここに来て突然,手続き型まみれのコードが出現して不愉快に感じるかもしれません。 とはいっても C# のプログラムとしては平凡であり,特に複雑ではありません。

このとき within は次のように書ける。 FindPair は引数の述語関数を満たす a, b の二項組を返し,.Item2 はそのうちの b だけを選択する。

    static double within(double eps, IEnumerable<double> seq) {
        return seq.FindPair((a, b)=> Math.Abs(a - b) <= eps).Item2;
    }

平方根を求める sqrt の within を使ったバージョンは次のように書ける。

    static double sqrt_v1(double a0, double eps, double N) {
        return within(eps, repeat((x)=> (x + N/x) / 2.0, a0));
    }

relative は次のように書ける。

    static double relative(double eps, IEnumerable<double> seq) {
        var pair= seq.FindPair((a, b)=> Math.Abs(a - b) <= eps * Math.Abs(b));
        return pair.Item2;
    }

sqrt の relative を使ったバージョンは次のように書ける。

    static double sqrt_v2(double a0, double eps, double N) {
        return relative(eps, repeat((x)=> (x + N/x) / 2.0, a0));
    }

近似値の列を生成する部分を書き換える必要はない。

0.3 の平方根 0.54772… を 1 から開始して許容誤差 0.1 で求めてみよう。

        Console.WriteLine(sqrt_v1(1, 0.1, 0.3));
        Console.WriteLine(sqrt_v2(1, 0.1, 0.3));

実行結果の例を示す。

0.555769230769231
0.547780809156242

3.2 「4.2 数値微分」と SelectWithNeighbor<T>

微分係数を近似的に得る easydiff 関数は次のように書ける。 デバッグ用に途中結果を見るため Console.WriteLine の呼び出しを入れておく。

    static double easydiff(Func<double, double> f, double x, double h) {
        double diff = (f(x + h) - f(x)) / h;
        Console.WriteLine(" -- h={0}, diff={1}", h, diff); // DEBUG
        return diff;
    }

h の値を繰り返し半分にして近似値の精度をあげる素朴な数値微分 differentiate は次のように書ける。

    static IEnumerable<double>
    differentiate(double h0, Func<double, double> f, double x) {
        return repeat(h => h / 2.0, h0).Select(h => easydiff(f, x, h));
    }

許容誤差 10, h の初期値 10 として x3 の x = 100 での微分係数 3x2 = 30000 を求めてみよう。

        Console.WriteLine
            (within(10, differentiate(10, x => x*x*x, 100)));
 -- h=10, diff=33100
 -- h=5, diff=31525
 -- h=2.5, diff=30756.25
 -- h=1.25, diff=30376.5625
 -- h=0.625, diff=30187.890625
 -- h=0.3125, diff=30093.84765625
 -- h=0.15625, diff=30046.8994140625
 -- h=0.078125, diff=30023.4436035156
 -- h=0.0390625, diff=30011.7202758789
 -- h=0.01953125, diff=30005.8597564697
30005.8597564697

確かに誤差の範囲で結果が得られた。 このとき easydiff でちょうど必要なだけの回数の演算が行われていることがデバッグ用の出力から確認できる。

この easydiff が生み出す近似値の列の収束を速くしていこう。

elimerror 関数を C# で実現するためには,任意の列の各要素 ai を任意の関数 f を使って f(ai, ai+1) に置き換えた列を生成する LINQ メソッドがあればよい。 次のような SelectWithNeighbor<T> メソッドを定義する。

    // seq の各要素 a(i) を f(a(i), a(i+1)) に換えた列を返す。
    static IEnumerable<R> SelectWithNeighbor<T, R>(this IEnumerable<T> seq,
                                                   Func<T, T, R> f) {
        using (var e = seq.GetEnumerator()) {
            if (! e.MoveNext())
                yield break;
            T a = e.Current;
            if (! e.MoveNext())
                yield break;
            T b = e.Current;
            for (;;) {
                yield return f(a, b);
                a = b;
                if (! e.MoveNext())
                    yield break;
                b = e.Current;
            }
        }
    }

これを使って elimerror 関数を次のように書く。 ここでもデバッグ用に途中結果を見るため Console.WriteLine の呼び出しを入れておく。

    static IEnumerable<double> elimerror(double n, IEnumerable<double> seq) {
        double p = Math.Pow(2, n);
        if (Double.IsNaN(p) || p == 1) {
            Console.WriteLine(" -- p={0}", p);
            return seq;
        } else {
            return seq.SelectWithNeighbor((double a, double b) => {
                    double r =  (b * p - a) / (p - 1);
                    Console.WriteLine(" -- r={0}, b={1}, n={2}", r, b, n);
                    return r;
                });
        }
    }
ここでもし値が得られない場合は,さしあたり return seq で元の列をそのまま返すことにします。

列の最初の三つの近似値から値を計算する order 関数は次のように書ける。 foreach ループを三回で打ち切って列の最初の3要素の値を得る。

    static double order(IEnumerable<double> seq) {
        var a = new double[3];
        int i = 0;
        foreach (double e in seq) {
            a[i] = e;
            if (++i == 3)
                break;
        }
        double x = (a[0] - a[2])/(a[1] - a[2]) - 1;
        double n = Math.Round(Math.Log(x, 2.0));
        return n;
    }

しかし,この elimerrororder を組み合わせて次のように improve 関数を書くことは適切ではない。 IEnumerable<double> オブジェクト s が2回出現しており,それぞれ別々に列の先頭から要素の計算を繰り返すからである。 order 関数で最初の3要素を計算した後,elimerror 関数でまた最初の要素から計算するため,最初の3要素の計算を無駄に繰り返すことになる。

    static IEnumerable<double> improve(IEnumerable<double> s) {
        return elimerror(order(s), s);
    }
ここに遅延評価の call by need と LINQ の call by name の違いが現れています。

実際,同じように実験すると

        Console.WriteLine
            (within(10, improve(differentiate(10, x => x*x*x, 100))));
 -- h=10, diff=33100
 -- h=5, diff=31525
 -- h=2.5, diff=30756.25
 -- h=10, diff=33100
 -- h=5, diff=31525
 -- r=29950, b=31525, n=1
 -- h=2.5, diff=30756.25
 -- r=29987.5, b=30756.25, n=1
 -- h=1.25, diff=30376.5625
 -- r=29996.875, b=30376.5625, n=1
29996.875

予想どおり,残念ながら最初の h=10, h=5, h=2.5 について同じ計算を繰り返した。 この問題を解決するため,次節のような Cached<T> を定義しよう。

3.3 Cached<T> による遅延リスト

Cached<T> はコンストラクタ引数で与えられた列の透過的なラッパーである。 ただし,一度評価した要素の値をキャッシュし,次に同じ要素が求められたときはキャッシュした値を返す。 Cached<T> はマルチスレッドで安全に使うことができる。 ただし,その代償として,元の列の要素を評価するとき,この Cached<T> インスタンス自身に lock がかかる。

これを代償と呼ぶのは,元の列の要素を計算しているときに何かを lock をした場合,デッドロックが発生する潜在的な原因となるからです。 マルチスレッドでない限り無関係です。 そもそも,ことのついでに実装しただけで,このスレッドセーフ性そのものが「なぜ関数プログラミングは重要か」の範囲では不要です……つまり,正しく蛇足です。 しかし,本格的な実装のヒントになると思いますから,このまま使います。
    // 一度評価した値をキャッシュする列 (≒ lazy list)
    class Cached<T>: IEnumerable<T> {
        List<T> cache = new List<T>();
        IEnumerator<T> generator;
        int index = -1;         // generator で生成した値の添え字位置

        public Cached(IEnumerable<T> sequence) {
            generator = sequence.GetEnumerator();
        }

        IEnumerator IEnumerable.GetEnumerator() {
            return GetEnumerator();
        }

        // 注意: 要素の評価時に this を lock する。
        public IEnumerator<T> GetEnumerator() {
            for (int i = 0;; i++) {
                T e;
                lock (this) {
                    if (i <= index) {
                        e = cache[i];
                    } else {
                        if (! generator.MoveNext())
                            break;
                        e = generator.Current;
                        cache.Add(e);
                        index++;
                    }
                }
                yield return e;
            }
        }
    }

これを使って improve 関数を次のように書くことができる。

    static IEnumerable<double> improve(IEnumerable<double> seq) {
        var s = new Cached<double>(seq);
        return elimerror(order(s), s);
    }

同じように実験すると,今度は同じ計算を繰り返さないことが分かる。

        Console.WriteLine
            (within(10, improve(differentiate(10, x => x*x*x, 100))));
 -- h=10, diff=33100
 -- h=5, diff=31525
 -- h=2.5, diff=30756.25
 -- r=29950, b=31525, n=1
 -- r=29987.5, b=30756.25, n=1
 -- h=1.25, diff=30376.5625
 -- r=29996.875, b=30376.5625, n=1
29996.875

improve を二重に適用するとさらに速く収束する。 許容誤差を 10 としても 1 の桁の誤差を残すよりも速く収束し,double 型での正しい値 30000 が得られる。

        Console.WriteLine
            (within(10, improve(improve(differentiate(10, x => x*x*x, 100)))));
 -- h=10, diff=33100
 -- h=5, diff=31525
 -- h=2.5, diff=30756.25
 -- r=29950, b=31525, n=1
 -- r=29987.5, b=30756.25, n=1
 -- h=1.25, diff=30376.5625
 -- r=29996.875, b=30376.5625, n=1
 -- r=30000, b=29987.5, n=2
 -- r=30000, b=29996.875, n=2
30000

super 関数は,近似列に対し improve 関数を repeat 関数で繰り返し適用することによって,次第に改良されていく近似列の列を作り,それぞれの近似列の二番目を選り抜いて新しい近似列を構成する。 この関数では各近似列が繰り返し評価されるから,初期の近似列,および次第に改良されていくそれぞれの近似列を Cached<double> でラップする。

    static IEnumerable<double> super(IEnumerable<double> seq) {
        return repeat(s => new Cached<double>(improve(s)),
                      new Cached<double>(seq))
            .Select(s => s.Skip(1).First());
    }

同じように実験してみよう。

       Console.WriteLine
           (within(10, super(differentiate(10, x => x*x*x, 100))));
 -- h=10, diff=33100
 -- h=5, diff=31525
 -- h=2.5, diff=30756.25
 -- r=29950, b=31525, n=1
 -- r=29987.5, b=30756.25, n=1
 -- h=1.25, diff=30376.5625
 -- r=29996.875, b=30376.5625, n=1
 -- r=30000, b=29987.5, n=2
 -- r=30000, b=29996.875, n=2
 -- h=0.625, diff=30187.890625
 -- r=29999.21875, b=30187.890625, n=1
 -- r=30000, b=29999.21875, n=2
 -- p=NaN (非数値)
30000
elimerror の計算で使う数値が最後には NaN になっています。 そうした場合,ここでの elimerror の定義は元の列をそのまま返すようにしていますから,この例は無事終わります。

いわゆる car と cdr の cons セルからなる普通の遅延リストならば,car を利用した後,その参照を消せばやがてガベージコレクションで cdr だけが残る。 しかし,この Cached<T> の実装ではいったん作られた要素は Cached<T> インスタンス自身が消えるまでいつまでも消えない。 ここでの例では生成されるそれぞれの近似列もその寿命も短いから問題ないが,場合によっては問題になるおそれがある。 そのときは一定時間経つなどしたらキャッシュを消すなどの対策が必要になるだろう。

3.4 「4.3 数値積分」について

関数 f の二端点 a, b 間の積分を台形の面積で近似する easyintegrate 関数は次のように書ける。 途中結果を見るために Console.WriteLine の呼び出しを入れておく。

    static double easyintegrate(Func<double, double> f, double a, double b) {
        double integ = (f(a) + f(b)) * (b - a) / 2;
        Console.WriteLine(" -- a={0}, b={1}, integ={2}", a, b, integ);
        return integ;
    }

この easyintegrate から積分の近似値の列を得る integrate 関数は,やや驚くべきことに C# で次のようにそのまま再帰定義できる。

    static IEnumerable<double> integrate_v1
    (Func<double, double> f, double a, double b) {
        yield return easyintegrate(f, a, b);
        double mid = (a + b) / 2;
        var xx = integrate_v1(f, a, mid);
        var yy = integrate_v1(f, mid, b);
        foreach (double e in xx.Zip(yy, (x, y)=> x + y)) {
            yield return e;
        }
    }

ここで関数本体の最初の行の yield return が,ある意味,再帰の底になっており,列の最初の要素を返す。 それ以降は再帰的に二端点間 a, b の前半と後半の近似列 xx, yy を得て順次その要素 x, y を足し合わせ,自らの要素の値として yield return する。

x3 の 1 から 2 までの積分 [x4/4]21 = (16 - 1)/ 4 = 3.75 に対する近似列の最初の4要素を見よう。

       Console.WriteLine
           (string.Join(", ", integrate_v1(x => x*x*x, 1, 2).Take(4)));
 -- a=1, b=2, integ=4.5
 -- a=1, b=1.5, integ=1.09375
 -- a=1.5, b=2, integ=2.84375
 -- a=1, b=1.25, integ=0.369140625
 -- a=1.25, b=1.5, integ=0.666015625
 -- a=1.5, b=1.75, integ=1.091796875
 -- a=1.75, b=2, integ=1.669921875
 -- a=1, b=1.125, integ=0.1514892578125
 -- a=1.125, b=1.25, integ=0.2110595703125
 -- a=1.25, b=1.375, integ=0.2845458984375
 -- a=1.375, b=1.5, integ=0.3734130859375
 -- a=1.5, b=1.625, integ=0.4791259765625
 -- a=1.625, b=1.75, integ=0.6031494140625
 -- a=1.75, b=1.875, integ=0.7469482421875
 -- a=1.875, b=2, integ=0.9119873046875
4.5, 3.9375, 3.796875, 3.76171875

「なぜ関数プログラミングは重要か」の議論と同じように f の値を再計算しないバージョンを次のように書くことができる。

    static IEnumerable<double> integrate
    (Func<double, double> f, double a, double b) {
        return integ(f, a, b, f(a), f(b));
    }
    static IEnumerable<double> integ
    (Func<double, double> f, double a, double b, double fa, double fb) {
        yield return (fa + fb) * (b - a) / 2;
        double m = (a + b) / 2;
        double fm = f(m);
        var xx = integ(f, a, m, fa, fm);
        var yy = integ(f, m, b, fm, fb);
        foreach (double e in xx.Zip(yy, (x, y)=> x + y)) {
            yield return e;
        }
    }

f の値を何度計算したか確認できるように f の定義に Console.WriteLine の呼び出しを入れて, f(x) = x3 の 1 から 2 までの積分 3.75 に対する近似列の最初の4要素を見よう。

       Console.WriteLine
           (string.Join(", ", integrate(x => {
                       Console.WriteLine(" -- x={0}", x);
                       return x*x*x;
                   }, 1, 2).Take(4)));
 -- x=1
 -- x=2
 -- x=1.5
 -- x=1.25
 -- x=1.75
 -- x=1.125
 -- x=1.375
 -- x=1.625
 -- x=1.875
4.5, 3.9375, 3.796875, 3.76171875
この再帰では重複なく手戻りなく区間を細分化していますから,計算がなにも無駄に繰り返されていないことが分かります。 キャッシュは不要です。

許容誤差 0.005 で値を求めてみよう。

       Console.WriteLine
           (within(0.005, integrate(x => {
                       Console.WriteLine(" -- x={0}", x);
                       return x*x*x;
                   }, 1, 2)));
 -- x=1
 -- x=2
 -- x=1.5
 -- x=1.25
 -- x=1.75
 -- x=1.125
 -- x=1.375
 -- x=1.625
 -- x=1.875
 -- x=1.0625
 -- x=1.1875
 -- x=1.3125
 -- x=1.4375
 -- x=1.5625
 -- x=1.6875
 -- x=1.8125
 -- x=1.9375
 -- x=1.03125
 -- x=1.09375
 -- x=1.15625
 -- x=1.21875
 -- x=1.28125
 -- x=1.34375
 -- x=1.40625
 -- x=1.46875
 -- x=1.53125
 -- x=1.59375
 -- x=1.65625
 -- x=1.71875
 -- x=1.78125
 -- x=1.84375
 -- x=1.90625
 -- x=1.96875
3.750732421875

残念ながら収束が遅いことが分かる。 ここで「なぜ関数プログラミングは重要か」の議論にしたがって前節の improve を利用すると

       Console.WriteLine
           (within(0.005, improve(integrate(x => {
                           Console.WriteLine(" -- x={0}", x);
                           return x*x*x;
                       }, 1, 2))));
 -- x=1
 -- x=2
 -- x=1.5
 -- x=1.25
 -- x=1.75
 -- r=3.75, b=3.9375, n=2
 -- r=3.75, b=3.796875, n=2
3.75

期待どおり圧倒的に速く収束する。

同様にして π/4 = 0.785398… の近似列を二番目の要素まで計算してみよう。

       Console.WriteLine
           (string.Join(", ", improve(integrate(x => {
                           Console.WriteLine(" -- x={0}", x);
                           return 1.0 / (1 + x * x);
                       }, 0, 1)).Take(2)));
 -- x=0
 -- x=1
 -- x=0.5
 -- x=0.25
 -- x=0.75
 -- r=0.783333333333333, b=0.775, n=2
 -- r=0.785392156862745, b=0.782794117647059, n=2
0.783333333333333, 0.785392156862745

「なぜ関数プログラミングは重要か」の記述どおり,二番目の近似値は 1/(1+x*x) を 5 回しか評価せずに小数点以下 5 桁まで正しい値を得ていることが分かる。

新年早々今年1月9日にこの感動的な結果を得て「なぜ関数プログラミングは重要か」を「なぜ LINQ プログラミングは重要か」のように読み替えられると確信したのが本稿をまとめた動機でした。 次章の内容はゲームツリーであり,一見するとさらに LINQ とは無縁に思えますが,勢いで進めてみると応用は案外簡単でした。
(1/3) へ (3/3) へ

Copyright (c) 2015 OKI Software Co., Ltd.