2015年に徳之島にプチ移住をしたとき、ひとときの我が家となった家にはテレビもインターネット接続もなかった。実のところ残念だったのはネット接続がないことのほうで、もしそれが可能だったらネット中継で地元のサッカーチームがJ1でプレーする勇姿を見られたはずだった。
一人で酒を飲むことも控えていたので柄にもなくストイックな生活だったと思う。徳之島方言のコーパスを作ることが目的の移住だったのでコーパス作成の作業に時間を使うようにしていた。作業の内容は形態素に番号を振ることで、そのために多少のプログラミングをしたりすることはあったが、時間の大部分は機械的な作業に費やした。
あまり単調なことばかりしていると頭を使いたくなる。あのころ熱中していたのはちょっと数学的なことをするプログラミングだった。
「タクシー数」というのがある。どうしてこんな変な名前がついているのか、面白いけれど説明が長くなるのでそれは省略するが、
a^2+b^2=c^2+d^2 （^で累乗を表す）
この式を満たす整数の組で、合計が最小のものを言う。式を満たすa,b,c,dをかたはしから探してみるとこれがいくらでも見つかる。もうタクシー数ではなくなっているけれど。
二乗だと簡単すぎるので、a^4+b^4=c^4+d^4 を満たす整数の組を探すことにすると、このような組はいくらでもあるらしいのだが、プログラミングが急に難しくなる。aが４つの整数の中で最大だとするとa>c>d>bとなるのだが、aの値が大きくなるにつれて私がそのころ使っていたVBAでは手に負えなくなってくる。aが100のとき、aの4乗はVBAで浮動小数点表示になる。つまり整数として扱えなくなるのだ。
100以上の数の4乗を整数として扱うためには、8桁以下の数と9桁から16桁までの数の組として数を表すようにする。この整数の組でかけ算をしたり、足し算をしたり、大小比較をしたりする。これを多倍長計算というのだが、これは単純な計算とくらべて計算時間をはるかに使う。
それから、一番単純なアルゴリズムだとループを4重に回すことになる。aの4乗に比例する回数、上の数式が成り立つかどうかチェックすることになる。これは大変なことで、aが100だったら、100000000回よりは少ないがこれの何分の1かの回数だけチェックを行うことになる。aが大きくなるとチェック回数は爆発的に増えていく。
最初はaが100以下の数のときのb,c,dの値がどうなるかを調べていたのだが、あるところから多倍長計算をしなければならなくなったのと計算量が爆発的に増えたせいで、目に見えてプログラムの実行に時間がかかるようになった。
100<a<150 のように区間を区切ってb,c,dの値を調べるのだが、30分ぐらいかけないとプログラムが終了しない。大体30分ぐらいかかるので、プログラムを起動させてから自転車で散歩に出かけて戻ってきて出力結果を見るということを繰り返していた。
一番単純なその代わりに分かりやすいアルゴリズムでは全く歯が立たないので、毎日新しいやりかたを考えて試すことをしていた。タクシー数探索はアルゴリズムの玉手箱のようなもので、本当にいろいろなやり方がある。出発点がaの4乗に比例という最悪なものなので、aの3乗に比例に改善すればかなり劇的な変化になる。
我が家に帰ってから、まずやってみたのはプログラム言語をVBAからVisual Basicに変えることだった。VBAはインタープリターでVisual Basicはコンパイラーなのでこれだけでスピードが100倍ぐらいに上がる。整数表示で扱える桁数が爆増する。配列要素の数がはるかに大きくなる（扱える数が増える）。
アルゴリズムとして一番効果的だと思われたのが、a^4+b^4の全ての組み合わせの合計をソートして、頭から順番に見ていくことだった。n番目とn+1番目が同じだったらa^4+b^4=c^4+d^4が成立していることになる。
n個の要素を最も効率的にソートしたときの計算回数はnlognに比例する。aの最大値がnのとき、対象となる要素の数はn^2/2なので、タクシー数の計算はn^2lognに比例することになる。一番素朴な方法はn^4に比例する計算量だったので、n^2に比例は格段に改良したことになる。念押しになるが、ソートはアルゴリズムの優劣がはっきり存在するので、「最も効率的な」方法を使った場合の話である。
ところが、いかにVisual Basicと言っても、nが20000のときの対象となる4億もの要素を扱うことはできない。うろ覚えで申し訳ないが、nが5000ぐらいでパンクしてしまうはずだ。
そこで登場するのが2分木である。2分木ではn番目の要素は2n番目と2n+1番目の要素より必ず小さい。この2分木を利用する。
分かりやすくするためにお皿を重ねた山をたとえとして用いる。まず、お皿を重ねた山をイメージしてもらいたい。お皿の山は左から右にだんだん大きくなるように並んでいる。一番左の山はお皿1枚だけで、そこには1,1,1と書いてある。「aの値、bの値、合計」の意味である。n番目の山は、n個の皿が積み重ねてあって、一番上の皿にはn,1,n^4+1、一番下の皿にはn,n,n^4+n^4と書いてある。bの値が下に行くほど大きくなっている。
このお皿の山の行列から1枚ずつ皿を取って2分木を構成する。2分木がうまく作れていれば、先頭はその2分木のなかで一番小さい。
分かりやすくするために一番左側の山から10番目の山まで1枚ずつ皿を取り出す。このままで2分木ができているのだが、これからあとは2分木の先頭の皿を取り出していく。皿を取るときに以下の2つの決まりで先頭の皿の補充と2分木の拡張を行う。
1．先頭の皿を取ったときにその皿がもともとあった山から皿を先頭に補充し、2分木を構成する。もし山の皿がなくなっていたら、隣の山から補充する。
2．一番右の山から補充して、もしそれが2枚目の皿だったら、その右の山から1枚目の皿を2分木の一番下に加えて、2分木を構成する。結果として2分木を下に伸ばすことになる。
「2分木を構成」というのは、新しい要素が加わったときにそれが2分木の条件に合致しているかどうか調べ、もしそうでなければ皿を交換するということを条件が満たされるまで行うことを言う。この手間は2分木の長さがnとするとlognとなる（logの底は2）
このようにして、2分木を適宜大きくしながら2分木の先頭の皿を取り出していって、直前に取り出したものと合計が同じかどうかチェックし、もし同じであれば直前の皿とその皿のデータを記録する。これを繰り返す。
たぶん、これが最速のアルゴリズムではないかと思う。あとはできるだけ多倍長計算をサボるぐらいしか、実行時間を短縮する方法がない。
数式を極力使わずにアルゴリズムを記述しようとしたらかえってわかりにくくしてしまったような気がする。本当はすごく単純なことなんです。頭が悪いせいでまわりくどい説明になったのか。
信大のもといた講座の卒論発表会で卒業生のＴ君と再会し、このブログの話をしたら、早速見てくれた形跡があった。Ｔ君にひまつぶしの話もしていた。Ｔ君、あのとき私が言っていたのはこういうことです。