β関数 (ゲーデルのβ関数のほう) を SymPy で計算してみます。
その際に、できるだけ結果のゲーデル数が小さくなるように工夫してみます。
なお、本記事は著者の NumPy 配列の練習を兼ねているので、無理やりおしゃれベクトル計算しているところがあります。
参考文献
- 林晋, 八杉満利子 訳・解説「ゲーデル 不完全性定理」岩波文庫 (青 944-1), 2006.
- (原典) K. Gödel. プリンキピア・マテマティカおよび関連した体系の形式的に決定不能な命題について I, Monatshefte für Mathematik und Physik (38), pp. 173-198, 1931.
- 前原昭二「数学基礎論入門」基礎数学シリーズ26, 朝倉書店, 2000.
- Wikipedia "Gödel numbering for sequences" https://en.wikipedia.org/wiki/G%C3%B6del_numbering_for_sequences
- Wikipedia「中国の剰余定理」 https://ja.wikipedia.org/wiki/%E4%B8%AD%E5%9B%BD%E3%81%AE%E5%89%B0%E4%BD%99%E5%AE%9A%E7%90%86
1. ゲーデルのβ関数とは
自然数の列を自然数に符号化 (エンコード) する際に、復号 (デコード) のために使われた関数のことです。ゲーデルの不完全性定理 (1931) の証明に出てきたことから、こう呼ばれます。なお、素数の冪を使う、いわゆるゲーデル数とは別の部分に出てきたものです。
以下を仮定します。
- 自然数とは 0 を含むとします。
- 数列の添字は 0 から始まるものとします。ゲーデルの論文と同じです。
自然数 $l,m,i$ ($m>0,0\le i\le n-1$) に対して関数 $\beta$ を
\beta(l,m,i) = \text{rem}(l,(i+1)m + 1)
で定義します (rem は通常の剰余演算。典型的には mod とか % が使われる)1。
このとき、具体的に与えられた、長さ $n\ge 1$ の自然数列 $a=\langle a_0, \dots, a_{n-1}\rangle$ に対して、ある自然数 $l,m$ を計算することができて、$0\le i\le n-1$ に対して $\beta(l,m,i) = a_i$ となるようにできます2。すなわち、$(l,m)$ がその数列の自然数表現 (つまりゲーデル数) であるとみなせます。自然数のペアを単一の自然数で符号化する方法もある (補足参照) ので、結果として、自然数列を自然数として符号化できることになります。なお、以下の議論ではペアの符号化/復号は明示せずに進めます。
$l$ の計算には中華剰余定理 (Chinese Remainder Theorem; CRT) が必要なので、関数 crt を持つ Python のライブラリである SymPy を使用します。
具体例が必要なので、長さ n=5 の配列 a = [31, 5, 51, 0, 2] を符号化する例を考えます。
>>> import numpy as np
>>> from sympy.ntheory.modular import crt
>>> a = [31, 5, 51, 0, 2]
>>> a
[31, 5, 51, 0, 2]
>>> n = len(a)
>>> n
5
何かと便利なので、NumPy 配列 a_ に変換しておきましょう。あと、a_ の添字の配列も作っておきます。
>>> a_ = np.array(a)
>>> a_
array([31, 5, 51, 0, 2])
>>> a_idx = np.arange(n)
>>> a_idx
array([0, 1, 2, 3, 4])
2. β関数による符号化の計算
符号化のための計算について、方針を述べた後、実際に Python で計算してみます。
上記のように数列が与えられたとき、符号化の方針は以下のようになります。
- 各 $m_i$ が互いに素になるような (＋ある条件を満たすような) $m$ を何らかの方法で計算する。
- 各 $m_i$ が互いに素であることを用いて、中華剰余定理で $l$ を計算する。
一般的には階乗を使って $m$ を計算します。しかし、この方法では $m$ が非常に大きくなります。そうすると $l$ も大きくなりがちです。本記事ではできるだけ $m$、ひいては $l$ が小さくなる方法で計算してみます。ただし、これで最小になっているかどうかは私にはわかっていません。
以下、$m_i = (i+1)m + 1$ と表します。つまり、β関数の定義は $\beta(l,m,i) = \text{rem}(l, m_i)$ となります。
m を求める
Wikipediaの "Gödel numbering for sequences" の記事はちょっと読みにくいのですが、$m$ が満たすべきできるだけ弱い条件について説明されており、参考になります (最弱か否かは不明)。それによると $m$ が以下の条件を満たすならば、正しい符号化のために採用できます。
i\mid m \quad (1\le i\le n-1) \\
a_i < m_i \quad (0\le i \le n-1)
第1式は $m$ が $1,\dots,n-1$ の公倍数であることを意味しています。よって、$m$ の候補として $1,\dots,n-1$ の最小公倍数 (LCM) $m_0$ の倍数を列挙し、そのうちで第2式を満たす最小のものを見つけるとよいでしょう。ではやってみます。
m0 を求める
NumPy には2引数関数 lcm があります。これはいわゆる ufunc なので、reduce メソッドを使うことで配列に対して LCM を計算できます。やってみます。さきほど添字の配列 $[0,1,\dots,n-1]$ を作ったのでそのスライスとしての $[1,\dots,n-1]$ を使用します。
>>> m0 = np.lcm.reduce(a_idx[1:])
>>> m0
12
$1,2,3,4$ の LCM を手計算すると確かに $12$ になります。この $m_0$ の倍数 $m$ は必ず第1式を満たします。
$m_i$ を生成する
次の節で、いろいろな $m$ に対して $\{m_i\}$ を計算します。ここではその方法について先に説明します。
各 $m$ に対して、$\{m_i\}$ を求めて $\{a_i\}$ と比較する必要があります。
この前者は $\{ m_i \mid 0\le i \le n-1 \} = \{ (i+1)m +1 \mid 0\le i \le n-1 \}$ です。$i$ の範囲が先に定義しておいた a_idx と一致するので、配列計算で一気に求めることができます。求める関数を m_i_seq として定義します。
>>> def m_i_seq(m): return (a_idx + 1) * m + 1
ためしに m=12 に対して計算するとこうなります。
>>> mis = m_i_seq(12)
>>> mis
array([13, 25, 37, 49, 61])
(mis は $m_i$s のイメージ)
条件を満たす m を探す
では、条件の第2式を満たす $m$ ($m_0$ の倍数) を探してみましょう。まず $m=m_0=12$ の場合についての $\{m_i\}$ はさきほど計算したとおり mis = [13 25 37 49 61] です。これと $\{a_i\}$ の比較はこれも NumPy 配列のベクトル計算機能で一気に計算できます。
>>> a_ < mis
array([False, True, False, True, True])
これは2つの $i$ について成り立ってないことが確認できます。
次に、では $m_0$ を2倍にして $m=24$ としてみます。このとき、
>>> mis = m_i_seq(24)
>>> mis
array([ 25, 49, 73, 97, 121])
>>> a_ < mis
array([False, True, True, True, True])
まだ成り立っていません。$m=36$ だとこうなります。
>>> mis
array([ 37, 73, 109, 145, 181])
>>> a_ < mis
array([ True, True, True, True, True])
すべて True になったので OK です。この $m=36$ を採用します。
(オプション) $m_i$ の coprimality を確認
前述のように $m$ を決めると、$\{m_i\}$ の相異なる2要素は互いに素になることが証明できます。ここでは省きますが、さきほどの条件のうち第1式を使います。
本記事では、これを計算で確認してみます。NumPy の gcd は ufunc なので、直積集合に関数を適用する outer メソッドを持っています。これで最大公約数が1に等しいか否かを確認できます。
>>> np.gcd.outer(mis, mis)
array([[ 37, 1, 1, 1, 1],
[ 1, 73, 1, 1, 1],
[ 1, 1, 109, 1, 1],
[ 1, 1, 1, 145, 1],
[ 1, 1, 1, 1, 181]])
たしかに対角成分を除いて1になっていることが確認できました (対角成分は gcd(x,x) なので結果は x になる)3。
$l$ を計算
あとは中華剰余定理をふつうに適用できます。各 $m_i$ が互いに素であることから、合同式系
l\equiv a_i \pmod {m_i}
は $m_0\dots m_{i-1}$ を法として一意な解 $l$ を持ちます。crt を使って解きます。crt の引数には NumPy 配列を指定できないようなので、List に戻した/変換したものを指定します。
>>> l, M = crt(mis.tolist(), a)
>>> l, M
(mpz(2496662055), 7726764205)
l が得られました。これは gmpy2 による多倍長整数になっています。
なお、M は $m_0\dots m_{i-1}$ の値です。いちおう確認しておきます。
>>> np.multiply.reduce(mis)
7726764205
最後に、求めた $l$ は $l\equiv a_i \pmod {m_i}$ を満たしますが、$a_i < m_i$ であるように設定してあることから、$a_i = \text{rem}(l, m_i)$ が成り立ちます。
以上で題意を満たす $(l, m) = (2496662055, 36)$ が得られました。(なお、$\{m_i\}$ は $[36\cdot 1+1, 36\cdot 2+1, 36\cdot 3+1, 36\cdot 4+1, 36\cdot 5+1] = [ 37, 73, 109, 145, 181]$ です。)
3. β関数による復号 (概要)
逆に、$l$ を各 $m_i$ で割った余りが $a_i$ になることを確認します。
np.mod は ufunc なので、配列やスカラーを指定すればよしなに計算してくれます。
>>> decoded = np.mod(l, mis)
>>> decoded
array([mpz(31), mpz(5), mpz(51), mpz(0), mpz(2)], dtype=object)
>>> a_
array([31, 5, 51, 0, 2])
これで、復号結果がもとの a と一致することが確認できました。
なお、復号結果は gmpy2 のクラス mpz でラップされていることがわかります。これは特に気にせず、普通の int 値として、後続の計算で使うことができるはずです。ラッパクラスを外すこともできると思います (未確認)。
>>> decoded[1] * 6
mpz(30)
パッケージ化
Python の特殊メソッドを実装することにより、実際に配列であるかのようにパッケージ化できます。
(コードはまたこんど)
まとめ
というわけでβ関数による数列の表現方法について紹介しました。Python の多倍長整数と数論的関数のサポートにより、現実的な (もっと大きい) 長さの数列に関しても、符号化を実際に計算することができるはずです。また NumPy のベクトル演算により、効率の良い計算ができます。
実際に計算することは、抽象的な理論のより実用的な理解に役立つはずです。ではまた。
補足
補足 1: β関数に関してもう少し
エンコードとしてもっとシンプルな方法は存在します。例えばゲーデル論文の前半に出てくるように
\langle a_0, \dots, a_{n-1}\rangle \mapsto 2^{a_0}\cdot 3^{a_1}\cdot\dots \cdot p_{n-1}^{a_{n-1}}
とする方法もあります。β関数による方法の特徴は、復号が剰余演算だけで行えることです。これによりゲーデルは、論文の前半で提示した決定不能な論理式が、算術的論理式 ($+,\cdot, =$ のみによる1階述語論理式) としても表現できることを示しました。
補足 2: 数列の長さ
β関数を用いる場合、符号化された値を見ても、元の数列の長さはわかりません。ゲーデルの証明においては数列の長さを復号できる必要はありませんでしたが、もし長さが必要ならば
- 符号化のとき: 数列 $a$ の先頭に $a$ の長さ $n$ を挿入した列 $a' = \langle n, a_0, \dots, a_{n-1}\rangle$ を符号化する。
- 復号のとき: 数列の先頭要素を復号して $n$ を得る。$1\le i\le n$ について $a'[i]$ を復号し、これを $a[i-1]$ とする。
とすることにより正しく処理できます。
補足 3: 自然数のペアの符号化
多様な方法があると思います。
- 高校数学でよく出てくる、(階段状に数字を並べる) 全単射を使う方法: $(x,y) \mapsto (x+y)(x+y+1)/2+y$
- 中華剰余定理において、合同式が2個のバージョンを使う方法
もありますが、ここでは、復号が比較的容易なものとして以下のような方法を紹介します。参考: https://slideplayer.com/slide/16991262/
(x,y) \mapsto 2^x(2y+1)-1
$2^x$ は偶数、$2y+1$ は奇数です。よって、符号化後の値を2で割れるだけ割ることにより復号ができます。