無限数列、無限級数、円周率

SICP ネタです。

SICP ３章のストリームの節に、ライプニッツ級数というものが紹介されています。

pi/4 = 1 - 1/3 + 1/5 - 1/7 + ...

という式です。項毎に符号が入れ替わる奇数の逆数の無限列、その総和が、パイの４分の１に収束する、というものです。

一見理解しがたい話ですが、このことの理論的背景については Wikipedia の解説を参照していただきたいと思います。

理論はさて置き、本を読んでいて「なんか簡単に JavaScript でも出来そう」と思ったので、やってみました (本文最後にデモページへのリンクがあります)。

以下、コードは全て Scheme からの翻訳です。

Stream in JavaScript

まずストリームの実装を考えましょう。

ストリームとは、遅延リストとも言い、コンピュータの有限リソースの中で無限列の表現を可能にするものです (有限列でも構いませんが)。

Haskell 等、遅延評価をベースとする言語においてはリストはネイティブに遅延リストなんですが、Scheme, JavaScript その多の言語では、プログラマが遅延評価のための仕組みを用意する必要があります。

マクロがあれば割合簡単に実現できるんですが、そうでない場合は少々煩わしいです。

まず、コンストラクタ及びセレクタをこのように定義します:

function cons_stream(x, p) {
  return [x, p];
}

function stream_car(s) {
  return s[0];
}

function stream_cdr(s) {
  return s[1]();
}

つまり、cons_stream の第２引数として、遅延したい計算を引数ゼロの関数として渡す (プログラマが明示的に遅延させる必要がある)、ということです。

さらにマッピング、反復その他、必要なユーティリティ関数を用意します:

var the_empty_stream = [];

function stream_ref(s, n) {
  if (n == 0) {
    return stream_car(s);
  } else {
    return stream_ref(stream_cdr(s), n-1);
  }
}

function stream_nullp(s) {
  return s.length == 0;
}

function stream_map(f) {
  var argstreams = Array.prototype.slice.call(arguments,1);
  if (stream_nullp(car(argstreams))) {
    return the_empty_stream;
  } else {
    return cons_stream(apply(f, map(stream_car, argstreams)),
                       function() {
                         return apply(stream_map,
                                      cons(f, map(stream_cdr,
                                                  argstreams)));
                       });
  }
}

function stream_for_each(f, s) {
  if (stream_nullp(s)) {
    return "done";
  } else {
    f(stream_car(s));
    setTimeout(function() {
      stream_for_each(f, stream_cdr(s));
    }, 100);
  }
}

stream_for_each において、setTimeout を使って再帰していますが、これはストリームの表示が出来るようにするためです。map, apply 等の定義は省略しています。

例として、自然数列は次のように書けます:

function integers_starting_from(n) {
  return cons_stream(n,
                     function() {
                       return integers_starting_from(n+1);
                     });
}

var integers = integers_starting_from(1);

stream_for_each(function(x){ window.status = x },
                integers);

数列 (sequence) と級数 (series)

では、ライプニッツ級数を実装していきます。

function pi_summands(n) {
  return cons_stream(1.0/n,
                     function() {
                       return stream_map(function(x){ return -x },
                                         pi_summands(n+2));
                     });
}
pi_summands(1); // ==> 1, -1/3, 1/5, -1/7, ...

この関数で生成されるのはライプニッツ級数、ではなくて、数列です (数列を足し合わせたものを級数と言います)。そこで、一般的な数列に対してその和を求める関数を次に作ります。

function add_streams(s, t) {
  return stream_map(function(x,y){ return x+y }, s, t);
}

function partial_sums(s) {
  return cons_stream(stream_car(s),
                     function() {
                       return add_streams(stream_cdr(s),
                                          partial_sums(s))
                     });
}

partial_sums が級数の部分和を生成する関数です。たとえば、1,2,3,4,... というストリームを与えられると 1,3,6,10,... というストリームを返します。ライプニッツ級数においては、項を限りなく大きく取った時の部分和が円周率の４分の１に収束する、とされているわけです。

ということで、部分和のストリームを４倍するための関数も必要です:

function scale_stream(s, factor) {
  return stream_map(function(x){ return x*factor }, s);
}

これにより、円周率に収束するストリームは

var pi_stream = scale_stream(partial_sums(pi_summands(1)), 4);

として表すことができます。

さて、実際にこれを表示してみると、3.14 近辺で値が上下に振動する様子が分かります。が、その触れ幅がなかなか小さくならず、もどかしいです。収束がとんでもなく遅い級数なんですね。

そこで SICP ではさらに、このような収束の遅い級数を加速する (しかも劇的に) 工夫が紹介されています。

収束の加速装置 (Euler Transform)

次に示すのが、項毎に正負が入れ替わるような級数の収束を加速するという、オイラーによる変換関数です:

function square(n) {
  return n * n;
}

function euler_transform(s) {
  var s0 = stream_ref(s, 0);
  var s1 = stream_ref(s, 1);
  var s2 = stream_ref(s, 2);
  return cons_stream(s2 - (square(s2 - s1) / (s0 + (-2 * s1) + s2)),
                     function() {
                       return euler_transform(stream_cdr(s));
                     });
}

euler_transform に pi_stream を渡して表示してみると、確かに元のストリームよりも相当速く収束しそうな気配が感じられます。

これだけでも凄いんですが、もっと劇的に加速する方法があるんです:

function make_tableu(transform, s) {
  return cons_stream(s,
                     function() {
                       return make_tableu(transform, transform(s))
                     });
}

function accelerated_sequence(transform, s) {
  return stream_map(stream_car, make_tableu(transform, s));
}

accelerated_sequence(euler_transform, pi_stream);

euler_transform によって加速されたストリームにさらに euler_transform を適用し、それを無限に再帰していく、というものです。

ストリームのストリーム（無限列の無限列）が出来上がるわけなんですが、それぞれの１項目を取ることで、数値の無限列を得ることができます。

以上の結果はこちらをご覧下さい: pi-stream.htm (倍精度での演算では小数点以下 14 桁まで合ったところで NaN と表示されると思います)