SICP 3.5節 にある無限ストリームを C++ で実装した。
それを使ってフィボナッチ数列や素数、Aitken 加速などを実装できた。
ラムダ式とか auto とかをいい感じに使えた、と思う。
あとこの記事を描くのに noweb を使った文芸的プログラミングを初めて使った。
実装
SICP では無限ストリームは car の値と cdr 部分のストリームを計算する関数の組として実装されている。
たとえばnから始まる整数を並べた無限ストリームは
;; n から始まる整数ストリーム
(define (integers-starting-from n)
(cons-stream n (integers-starting-from (+ n 1))))となる。
これと同じように C++ でストリームを実装しようとするとすごく面倒臭くなる。
C++ にはラムダ式と auto を使って関数を定義するときには再帰が使えないという制限があるので、上の例と似た関数を
// コンパイルできない
auto integers_starting_from = [](auto n) {return cons(n, integers_starting_from(n + 1)); };
と定義することができず、型(いまの例なら std::functional<Stream<int>(int)> )を毎回書かなければならない。
流石にこれは面倒臭すぎるので、ストリームを配列っぽく実装した。
T 型のストリームは第n項を計算する関数によって初期化され、
一度計算した値は std::unordered_map<size_t, T> に保存される。
さらに、第n項を計算するときに他の項の値を参照できるように、ストリームを初期化するときに与える
関数には自分自身 (*this) が引数として与えられる。
#include <functional>
#include <iostream>
#include <unordered_map>
template <typename T>
struct Stream {
private:
// genrator(self, n) はストリームの第n項を返す
// 第一引数には自分自身 (*this) が与えられるので、第n項を計算するときに他の項を参照することができる
std::function<T(Stream &, size_t)> generator;
// 計算された値を保存するハッシュ
std::unordered_map<size_t, T> memo;
public:
Stream() {}
explicit Stream(std::function<T(Stream &, size_t)> generator_)
: generator(generator_) {}
// ストリームの第n項を返す
T &operator[](size_t n) {
if (memo.count(n) == 1) {
// 計算済みなら、その値を返す
return memo[n];
} else {
// 計算してないなら、計算を行い、結果を保存してから、結果を返す
memo[n] = generator(*this, n);
return memo[n];
}
}
};
以下のプログラムではフィボナッチ数列を計算している。
<<Streamの定義>>
// フィボナッチ数列
static Stream<size_t> fib([](auto &self, auto n) -> size_t {
if (n < 2) {
return 1;
} else {
// 他の項を参照する。メモ化は適切に行われる。
return self[n - 1] + self[n - 2];
}
});
int main() {
for (size_t i = 0; i < 10; i++) {
std::cout << "fib[" << i << "] = " << fib[i] << "\n";
}
std::cout << "fib[" << 100 << "] = " << fib[100] << "\n";
return 0;
}
fib[0] = 1
fib[1] = 1
fib[2] = 2
fib[3] = 3
fib[4] = 5
fib[5] = 8
fib[6] = 13
fib[7] = 21
fib[8] = 34
fib[9] = 55
fib[100] = 1298777728820984005その他の例:
#include <cassert>
<<Streamの定義>>
// 整数
static Stream<size_t> integers([](auto &, auto n) { return n; });
// 素数
static Stream<size_t> primes([](auto &self, auto n) {
static int i = 2;
for (;; i++) {
auto flag = true;
// n番目の素数は、n-1番目までの素数全てで割り切れない最小の自然数
for (size_t j = 0; j < n; j++) {
flag &= (i % self[j]) != 0;
}
if (flag) {
return i;
}
}
});
// 双子素数
static Stream<std::pair<size_t, size_t>> twin_primes(
[](auto &, auto n) -> std::pair<size_t, size_t> {
// primes[i] と primes[i+1] が双子素数であるようなiを返す
static Stream<size_t> twin_primes_idx([](auto &self, auto n) -> size_t {
if (n == 0) {
return 1;
}
auto i = self[n - 1] + 1;
for (;; i++) {
if (primes[i] + 2 == primes[i + 1]) {
return i;
}
}
});
auto younger = primes[twin_primes_idx[n]];
return {younger, younger + 2};
});
// コラッツ予想の操作をしたときに何回で1になるかを返す
static Stream<size_t> Collatz([](auto &self, auto n) -> size_t {
if (n == 1) {
return 0;
} else if ((n % 2) == 0) {
return self[n / 2] + 1;
} else {
return self[3 * n + 1] + 1;
}
});
int main() {
assert(integers[0] == 0);
assert(integers[1] == 1);
assert(integers[2] == 2);
assert(integers[3] == 3);
assert(integers[100] == 100);
assert(primes[0] == 2);
assert(primes[1] == 3);
assert(primes[2] == 5);
assert(primes[3] == 7);
assert(primes[4] == 11);
assert(primes[99] == 541);
assert(twin_primes[0].first == 3);
assert(twin_primes[1].first == 5);
assert(twin_primes[2].first == 11);
assert(twin_primes[3].first == 17);
assert(twin_primes[99].first == 3821);
assert(Collatz[1] == 0);
assert(Collatz[11] == 14);
assert(Collatz[27] == 111);
assert(Collatz[871] == 178);
std::cout << "OK\n";
return 0;
}
OK数列の加速
SICP にある Aitken 加速の例も実装できたので解説する。
まず \(\frac{\pi}{4}\) に収束する数列
\begin{align*} a_{n} = \sum_{i=0}^{n} \frac{(-1)^{n}}{2n+1} \end{align*}
を計算するストリームを作った。
この数列は
\begin{align*}
a_{n} =
\begin{cases}
1 & (n = 0) \\
a_{n-1} + \frac{(-1)^{n}}{2n+1} & (\text{otherwise})
\end{cases}
\end{align*}
と書き直せるので、この数列の第 n 項を計算するストリームを以下のように定義した。
#include <cmath>
// pi / 4 に収束する
static Stream<double> A([](auto &self, auto n) -> double {
if (n == 0) {
return 1;
} else {
return self[n - 1] + std::pow(-1, n) / (2 * n + 1);
}
});
#include <iostream>
<<Streamの定義>>
<<数列Aの定義>>
int main() {
for (size_t i = 0; i < 10; i++) {
std::cout << "4*A[" << i << "] = " << 4*A[i] << "\n";
}
std::cout << "4*A[" << 100 << "] = " << 4*A[100] << "\n";
std::cout << "4*A[" << 1000 << "] = " << 4*A[1000] << "\n";
std::cout << "4*A[" << 10000 << "] = " << 4*A[10000] << "\n";
return 0;
}
4*A[0] = 4
4*A[1] = 2.66667
4*A[2] = 3.46667
4*A[3] = 2.89524
4*A[4] = 3.33968
4*A[5] = 2.97605
4*A[6] = 3.28374
4*A[7] = 3.01707
4*A[8] = 3.25237
4*A[9] = 3.04184
4*A[100] = 3.15149
4*A[1000] = 3.14259
4*A[10000] = 3.14169かなりゆっくりではあるが A が \(\frac{\pi}{4}\) に収束しているのが分かる。
この数列 A に対して Aitken 加速を適用することができる。 数列 \(x_{n}\) に対して Aitken 加速を適用した数列 \(x_{n}^{\prime}\) は
\begin{align*}
x^{\prime}_{n} &= x_{n} - \frac{(\Delta x_{n})^{2}}{\Delta^{2}x_{n}} \\
\Delta x_{n} &= x_{n+1} - x_{n} \\
\Delta^{2} x_{n} &= \Delta x_{n+1} - \Delta x_{n} \\
\end{align*}
と定義され、\(x_{n}^{\prime}\) は \(x_{n}\) に比べて収束がとても速くなる(ことがある)。
Aitken 加速をストリームを受け取ってストリームを返す関数として実装した。
auto Aitken = [](auto &x) {
return Stream<double>([&x](auto &, auto n) {
auto ddx = (x[n + 2] - x[n + 1]) - (x[n + 1] - x[n]);
auto dx = x[n + 1] - x[n];
if (ddx == 0) {
// 分母が0になるときの例外処理
return x[n];
}
return x[n] - dx * dx / ddx;
});
};
Aitken 加速によって収束が速くなることを確認した:
#include <iostream>
<<Streamの定義>>
<<数列Aの定義>>
<<Aitken加速の定義>>
int main() {
auto AA = Aitken(A);
for (size_t i = 0; i < 10; i++) {
std::cout << "4*AA[" << i << "] = " << 4*AA[i] << "\n";
}
return 0;
}
4*AA[0] = 3.16667
4*AA[1] = 3.13333
4*AA[2] = 3.14524
4*AA[3] = 3.13968
4*AA[4] = 3.14271
4*AA[5] = 3.14088
4*AA[6] = 3.14207
4*AA[7] = 3.14125
4*AA[8] = 3.14184
4*AA[9] = 3.14141第10項で小数以下5桁まで正しく計算できているので、収束が加速されたことが分かる。
Aitken 加速を適用した数列に更に Aitken 加速を適用することができ、さらにその数列を加速することもできる。 この再帰的な Aitken 加速をストリームのストリームを使って実装した。
// 再帰的に Aitken 加速されたストリームのストリームを返す
auto AitkenTableau = [](auto &x) {
return Stream<Stream<double>>([&x](auto &self, auto n) {
if (n == 0) {
return Aitken(x);
} else {
return Aitken(self[n - 1]);
}
});
};
auto AitkenRecursive = [](auto &x) {
auto tableau = AitkenTableau(x);
// 計算した tableau をクロージャの中に入れる
// tableau は計算の中で更新されるので、mutable が必須
return Stream<double>(
[tableau](auto &, auto n) mutable { return tableau[n][0]; });
};
#include <iostream>
#include <iomanip>
<<Streamの定義>>
<<数列Aの定義>>
<<Aitken加速の定義>>
<<再帰的Aitken加速の定義>>
int main() {
auto AAA = AitkenRecursive(A);
std::cout << std::setprecision(15);
for (size_t i = 0; i < 10; i++) {
std::cout << "4*AAA[" << i << "] = " << 4*AAA[i] << "\n";
}
return 0;
}
4*AAA[0] = 3.16666666666667
4*AAA[1] = 3.14210526315789
4*AAA[2] = 3.141599357319
4*AAA[3] = 3.14159271403378
4*AAA[4] = 3.14159265397529
4*AAA[5] = 3.14159265359118
4*AAA[6] = 3.14159265358978
4*AAA[7] = 3.1415926535898
4*AAA[8] = 3.14159265358979
4*AAA[9] = 3.14159265358979第9項の値は小数以下14桁まで正しい。