全宇宙の素粒子の数を超えて…C++で巨大数に挑戦!

 早速ですが、本ブログ開設1周年を記念して、C++で、巨大数の計算を行なってみたいと思います。
 といっても今回は、多倍長整数、いわゆるBigIntとかそういう話ではなく、巨大数を表現するための表記……具体的には「コーンウェイのチェーン表記」の話です。
 ちなみに、ちょっとした式を計算しようと思っただけでもBigIntはおろか、全宇宙の素粒子を使っても到底表現できないほどの値になりますので、実用性は皆無。

コーンウェイのチェーン表記とは

 まず、コーンウェイのチェーン表記とは何か、Wikipediaを見てみましょう。

タワー表記の拡張による。まず、3つ組の場合を定義する。0でない自然数 a, b, c について
a \rightarrow \ b \rightarrow \ c = a \uparrow^{c} b = a \underbrace{\uparrow \cdots \uparrow}_{c} b
次に、2つ組の場合を定義する。2つ組の場合は次のように単純な累乗と定めるものとする。
a \rightarrow b = a \uparrow b = a^b
4つ以上の0でない自然数を連ねて書いたときは、以下の規則によって変形できるものとする。

  1. チェーンの数の最後及び最後から2番目が共に 1 でない場合
    • タワー表記における定義をコンウェイのチェーン表記で表すと次のようになる
      • a→b→n = a→(a→(b-1)→n)→(n-1)
    • よって、4つ以上の場合もこれに習って以下のように、最後の数は変えずに最後から2番目の数を 1 だけ減らしたチェーンと、最後の数を 1 だけ減らしたチェーンで表すものとする
        • a→b→...→x→y→z = a→b→...→x→(a→b→...→x→(y-1)→z)→(z-1)
  2. チェーンの数の最後が 1 の場合
    • タワー表記における定義をコンウェイのチェーン表記で表すと次のようになる
      • a→b→1 = a→b
    • よって、4つ以上の場合もこれに習って以下のように最後の 1 を消去できるものとする
        • a→b→...→x→y→1 = a→b→...→x→y
  3. チェーンの数の最後から2番目が 1 の場合
    • タワー表記における定義によると、 a→1→n は次のようにnの値に関係なく a となる
      • a \rightarrow 1 \rightarrow n = a\uparrow^n 1 = a\uparrow^{n-1} (a\uparrow^n 0) = = a\uparrow^{n-1} 1 = ... = a\uparrow 1 = a \rightarrow 1 = a
http://ja.wikipedia.org/wiki/%E3%82%B3%E3%83%B3%E3%82%A6%E3%82%A7%E3%82%A4%E3%81%AE%E3%83%81%E3%82%A7%E3%83%BC%E3%83%B3%E8%A1%A8%E8%A8%98

 よくわかりませんが、兎に角そういう事みたいです。
 とりあえずは、このルールに従ってライブラリを作ってみましょうか。

ハイパー演算子の実装

 まずは、最初に出てきた「タワー表記」とやらを実装してみるとしましょう。こちらも、Wikipediaによると次のようになっています。

ここでクヌースは、二重矢印をテトレーション(指数計算の反復)を表す演算子として定義した。
a\uparrow\uparrow b= \underbrace{a \uparrow a \uparrow \dots \uparrow a}_ {b\mathrm{ \ copies \ of \ }a } = \underbrace{a_{}^{a^{{}^{.\,^{.\,^{.\,^a}}}}}}_ {b\mathrm{ \ copies \ of \ }a }
(中略)
だがクヌースはこれに飽き足らず、「2重矢印」による演算を反復する演算子として、「3重矢印」を定義した。
a\uparrow\uparrow\uparrow b= \underbrace{a \uparrow\uparrow a\uparrow\uparrow \dots \uparrow\uparrow a}_ {b\mathrm{ \ copies \ of \ }a }
同様に、「4重矢印」演算子も定義できる。
a\uparrow\uparrow\uparrow\uparrow b= \underbrace{a \uparrow\uparrow\uparrow a \uparrow\uparrow\uparrow \dots \uparrow\uparrow\uparrow a}_ {b\mathrm{ \ copies \ of \ }a }
これを一般的に述べると、n 重の矢印演算子は、(n − 1) 重の矢印演算子の反復として表すことができる。
a \underbrace{\uparrow\uparrow\cdots\uparrow}_n b = \underbrace{a \underbrace{\uparrow\cdots\uparrow}_{n-1} a \underbrace{\uparrow\cdots\uparrow}_{n-1} a \dots a \underbrace{\uparrow\cdots\uparrow}_{n-1} a}_ {b\mathrm{ \ copies \ of \ }a }
なお、矢印を使った指数の記法 a \uparrow b = a ^ b も、クヌースの矢印記号の特殊例(一重矢印)として再解釈される。

http://ja.wikipedia.org/wiki/%E3%82%AF%E3%83%8C%E3%83%BC%E3%82%B9%E3%81%AE%E7%9F%A2%E5%8D%B0%E8%A1%A8%E8%A8%98

 ……とのことですが、どうやらより広い範囲で使えるらしい「ハイパー演算子」なるものもあるようです(引用は省略)。ですのでまず、こちらを実装してみます。

template <class T>
T hyper(const T& a, const T& n, const T& b)
{
  if (0 == n) return (T)(b + 1) ;
  if (1 == n) return (T)(a + b) ;
  if (2 == n) return (T)(a * b) ;

  if (0 == b) return (T)1 ;

  T result = a ;
  T n_1 = (T)(n - 1) ;
  for (T i = b - 1 ; i > 0 ; --i) {
    result = hyper(a, n_1, result) ;
  }
  return result ;
}

 はい、簡単ですね。一応、テンプレート化しておいて、どんな整数型に対しても適用できるようにしておきます。

チェーン計算処理の追加

 まず、任意のコンテナをチェーンとして扱う関数を作るとしましょう。

template <class Container>
typename Container::value_type conway_chain(const Container& list)
{
  typedef typename Container::value_type T ;
  const T ONE = (T)1 ;

  auto begin = list.begin() ;
  size_t n = std::find(begin, list.end(), ONE) - begin ; // 「1」が出てくる以降の無視
  if (2 <= n && 2 == *begin && 2 == *(begin + 1)) return (T)4 ;
  switch (n) {
  case 0 : return ONE ;
  case 1 : return *begin ;
  case 2 : return hyper(*begin, (T)3, *(begin + 1)) ;
  case 3 : return hyper(*begin, *(begin + 2) + 2, *(begin + 1)) ;
  default :
    // TODO
  }
}

 まず、チェーンの長さにかかわらず「1」以降は式の値に影響を与えないため、「1」以降は無視します。
 この結果として長さが0になったチェーンの値は1とします*1。長さ1のチェーンは最初の値そのもの、長さ2の場合はただのべき乗なのでhyper(a, 3, b)、長さ3の場合もハイパー演算子でそのまま表現できます。
 ついでに、最初2つの項目がともに2の場合は必ず4になりますので、計算を省略しておきましょう。

 それでは最後に、n≧4の場合の処理を追加しましょう。

std::vector<T> v(begin, begin + n) ;

auto rit0 = v.rbegin() ;
auto rit1 = rit0 + 1 ;
while (true) {
  --*rit1 ;
  *rit1 = conway_chain(v) ;
  if (ONE == --*rit0) {
    v.pop_back() ;
    if (3 == v.size()) return hyper(v[0], v[2] + 2, v[1]) ;

    // イテレータを求め直す
    rit0 = v.rbegin() ;
    rit1 = rit0 + 1 ;
  }
}

 まず、チェーンは処理の途中で変形されるため、チェーンのコピーを作ります。本当は不要な部分は元の値を参照した方が(特にBigIntを使う場合に)都合がよいのでしょうが、今回は割愛します。
 次に、変形規則1に応じた変形を行ないます。その結果最後の値が1になったら削除しますが、削除の結果として要素数が3になれば、あとはハイパー演算子を使って計算できるようになります。そうでなければ後から2番目の値を減らす作業に戻る。
 これを延々と繰り返して行けば、理論上は唯一の値が求まるはずです。

 それでは、やってみましょう。

int main()
{
  std::cout << conway_chain(std::vector<unsigned long long>({3ull, 3ull, 2ull})) << std::endl ;
}
7625597484987

 よし、Wikipedia3\uparrow\uparrow3 の値と一致*2

*1:長さ0=1で始まるチェーン=1のべき乗なので

*2:これ以上長くすると、1や2の特殊処理だけで十分な場合を除いて酷いことになります……