プログラミングの備忘録

プログラムをつくる過程を残すもの

タッパーの自己言及式の謎を解く

こんにちは。
今回は「タッパーの自己言及式 (Tupper's self-referential formula)」を取りあげようと思います。
(記事タイトルを AI が生成してくれるようになったので、さっそく使ってみました。)


例のごとく、少し前にこんなツイートを見かけました。

(Fermat's Library は理系的な雑学のツイートが多く、知見が広がるのでよく見ています。)

どこかで見たことがあるという方もいるかもしれません。
数式を図にすると数式そのものが浮かび上がる。
これは掘り下げない訳にはいきません。


目次


まずは実装してみる

さっそくですが、その真偽を確かめるべく実装してみます。


先程のツイートには載っていませんでしたが、$k$ という値が必要になります。

詳しくはまた後で触れますが、$k$ の値は…

$$ k = 960939379918958884971672962127852754715004339660129306651505519271702802395266424689642842174350718121267153782770623355993237280874144307891325963941337723487857735749823926629715517173716995165232890538221612403238855866184013235585136048828693337902491454229288667081096184496091705183454067827731551705405381627380967602565625016981482083418783163849115590225610003652351370343874461848378737238198224849863465033159410054974700593138339226497249461751545728366702369745461014655997933798537483143786841806593422227898388722980000748404719 $$

543桁にもなる大きな数です。


図は、ある座標 $(x, y)$ について右辺を計算し、その計算結果が $\displaystyle \frac{1}{2}$ と比べて大きいか小さいかで色分けをすることによって描くことができます。
(そのまま図示すると上下左右が逆になった図になるため、少し調整が必要になります。)


また、$\lfloor a \rfloor$ は $a$ を超えない最大の整数、$\mathrm{mod} (a, b)$ は $a$ を $b$ で割った余りを指します。


ここでは、Processing と Python で実装してみようと思います。

Processing

巨大数を扱うにあたって

実は、素直に実装しようとしても、いつも通りやっていてはうまくいきません。

普段は全く気にしませんが、扱える数値には下限上限があります。
int型での上限は $2147483647$ とされており、到底543桁には及びません。
(int型だ、と宣言するとある程度のデータ量をその数値に割くようにしてくれて、その「ある程度」が関係するとか)

long型という、もっと大きな (小さな) 数を扱えるものもありますが、それでも上限は $9223372036854775807$ と足りません。


そこで使うのが、「BigDecimal型 (BigDecimal (Java Platform SE 8))」です。
上限は $2147483647$ くらいだとされているようで、とんでもなく大きな数まで扱えることがわかると思います。

そんなBigDecimal型は Java のライブラリで、

import java.math.BigDecimal;

で読み込むことができます。

実装

あとは (関数は BigDecimal のものを使わなければいけませんが) 普通に計算していけば実装できます。

import java.math.BigDecimal;

String strk = "960939379918958884971672962127852754715004339660129306651505519271702802395266424689642842174350718121267153782770623355993237280874144307891325963941337723487857735749823926629715517173716995165232890538221612403238855866184013235585136048828693337902491454229288667081096184496091705183454067827731551705405381627380967602565625016981482083418783163849115590225610003652351370343874461848378737238198224849863465033159410054974700593138339226497249461751545728366702369745461014655997933798537483143786841806593422227898388722980000748404719";
BigDecimal k = new BigDecimal(strk);
int a = 17; // 縦
int b = 106; // 横

void setup(){
  size(800, 200);
  background(255);
  
  float l = width / b; // セル1辺の長さ
  
  for(int i = 0; i < b; i ++){
  for(int j = 0; j < a; j ++){
    BigDecimal x = new BigDecimal(i);
    BigDecimal y = new BigDecimal(j).add(k);
    
    BigDecimal v1 = new BigDecimal(a).multiply(x);
    BigDecimal v2 = y.remainder(new BigDecimal(a));
    BigDecimal v3 = v1.add(v2);
    BigDecimal v4 = new BigDecimal(2).pow(v3.intValue());
    BigDecimal v5 = new BigDecimal(1).divide(v4);
    BigDecimal v6 = y.divideToIntegralValue(new BigDecimal(a));
    BigDecimal v7 = v6.multiply(v5);
    BigDecimal v8 = v7.remainder(new BigDecimal(2));
    
    if(v8.intValue() > 0.5){
      fill(0);
    }
    else{
      fill(255);
    }
    rect(width-i*l-2*l, j*l+2*l, l, l);
  }
  }
  print("end"); // 計算が終わったら「end」と出力
}


変数に v シリーズがいくつも出てきてややこしかったり、描画位置が微妙な位置になっていて気持ち悪かったりしますが、実行すると一応描画できました。
(終了までに20秒程度かかりました。)

Python

Python では、int型の上限はなんとなし
しかし、float型は上限 $308$ 桁くらいと足りません。
(それでも十分大きいですが…)

結局、「decimal (decimal --- 十進固定及び浮動小数点数の算術演算 — Python 3.12.1 ドキュメント)」をインポートして巨大数を扱うことになります。

from decimal import Decimal, getcontext
import matplotlib.pyplot as plt

def tupper(k, a, i, j): # タッパーの自己言及式
    x = i
    y = j + k

    t = int(((Decimal(y) // 17) * (2 ** Decimal(- a * x - y % a))) % 2)

    if t > 0.5:
        return 1
    else:
        return 0

getcontext().prec = 1000 # 精度

k = 960939379918958884971672962127852754715004339660129306651505519271702802395266424689642842174350718121267153782770623355993237280874144307891325963941337723487857735749823926629715517173716995165232890538221612403238855866184013235585136048828693337902491454229288667081096184496091705183454067827731551705405381627380967602565625016981482083418783163849115590225610003652351370343874461848378737238198224849863465033159410054974700593138339226497249461751545728366702369745461014655997933798537483143786841806593422227898388722980000748404719
a = 17 # 縦
b = 106 # 横

print('a = {}, b = {}'.format(a, b))

xx = []
yy = []
for x in range(0, b):
    for y in range(0, a):
        if tupper(k, a, x, y) == 1:
            xx.append(-x)
            yy.append(-y)

# 散布図の描画
fig = plt.figure()
ax = fig.add_subplot()
ax.scatter(xx, yy, s=25, c="black", marker=",")
plt.show()


ネタばらし

改めて、タッパーの自己言及式を見てみます。

$$ \frac {1}{2} < \left\lfloor \mathrm{mod} \left( \left\lfloor {\frac {y}{17}} \right\rfloor 2^{-17 \lfloor x \rfloor - \mathrm{mod} \left( \lfloor y \rfloor ,17 \right)}, 2 \right) \right\rfloor $$

この数式では何が起こっているのでしょうか。
(いろいろありまして、先程つくったコードと以下の説明で使っている文字が一部対応していません。悪しからず。)


ざっくり言えば、この数式は $k$ をグラフ上に図示する機械で、$k$ は図示するとその数式そのものの形になるように調整されたもの、です。
つまり数式が先にあって、$k$ をうまいことつくって自己言及的になるようにした、という感じです。


さて、ある2進数 $a$ を $2^{0}$ の位から座標平面に $(0, 0), (0, 1), (0, 2), \dots$ と $h$ 個並べ、$(1, 0)$ に戻ってきてまた縦に並べていく、という操作を繰り返し、縦 $h$、横 $w$ の四角形をつくることを考えます。
($a$ をうまく調整して、この四角形 (ビットマップ) が何かしらの絵 (上下左右が逆) になるようにしておきます。)

$a$ の $n$ 桁目 ($2^{0}$ の位を $n = 0$ とする) の数は

$$ \left\lfloor \mathrm{mod} \left( a \times 2^{-n}, 2 \right) \right\rfloor $$

で求められ、これが $0$ か $1$ かを判定したいので $\displaystyle \frac{1}{2}$ を基準として、

$$ \frac{1}{2} < \left\lfloor \mathrm{mod} \left( a \times 2^{-n}, 2 \right) \right\rfloor $$

ここで、四角形における位置 $(i, j)$ が $a$ の何桁目にあたるかを求めるには、$n = h \times i + j$ を計算すれば良いので、

$$ \frac{1}{2} < \left\lfloor \mathrm{mod} \left( a \times 2^{- h \times i - j}, 2 \right) \right\rfloor $$


座標からその位置でのビットを計算したいので、数式から $a$ を消去します。
これは数学的なトリック (?) なのですが、$(i, j)$ から $(x, y) = (i, j + k)$ ($k = h \times a$) に平行移動させます。
すると、$\displaystyle a = \left\lfloor \frac{y}{h} \right\rfloor$、$j = \mathrm{mod} (y, h)$ と表せるので、

$$ \frac{1}{2} < \left\lfloor \mathrm{mod} \left( \left\lfloor \frac{y}{h} \right\rfloor \times 2^{- h \times x - \mathrm{mod} (y, h)}, 2 \right) \right\rfloor $$


$x$、$y$ は実数なので、四角形に描画するために整数に変換して、

$$ \frac{1}{2} < \left\lfloor \mathrm{mod} \left( \left\lfloor \frac{y}{h} \right\rfloor \times 2^{- h \times \lfloor x \rfloor - \mathrm{mod} (\lfloor y \rfloor, h)}, 2 \right) \right\rfloor $$


以上で、一般化した式が求められました。
タッパーの自己言及式では縦 $h$ を $17$ にしているため、初めに挙げたような式となります。
横 $w$ は式に含まれず、2進数 $a$ の桁数を $d$ として、$\displaystyle \left\lceil \frac{d}{h} \right\rceil$ で求めることができます。


ということは…?

$$ \frac {1}{2} < \left\lfloor \mathrm{mod} \left( \left\lfloor {\frac {y}{17}} \right\rfloor 2^{-17 \lfloor x \rfloor - \mathrm{mod} \left( \lfloor y \rfloor ,17 \right)}, 2 \right) \right\rfloor $$

この数式は $k$ をグラフ上に図示する機械で、$k$ は図示するとその数式そのものの形になるように調整されたものだった、ということでした。

つまり、$k$ を変えることで、好きなように絵を描くことができるということになります。


これに関して、便利なサイトがありました。

Tupper's Formula Tools

ドット絵を $k$ に変えてくれます。(逆も可)


例えば、k = 37067602316440724456404475121574976829800827292700312151932062052691617906507927381280234914856650067959341707704037256202764341388719745773546243699883854050921017111451676703126563717998971858347996166958438528749534028200461401835215861451434273715380484400478716115151572056421758253407955986722944213998525692367560978044283475116945590335229170642234077735071680796317197511645338572277873548247607374220806113516110376740106816215430532077738073620942940121330994950836541296197937646843827484372199351437787467973352862443913609216 とすると、


こうしてしまうと、自己言及でも何でもない、ただのお絵かきマシーンになってしまいますが。


まとめ

以上、「タッパーの自己言及式」を Processing と Python で実装し、そのカラクリの解明をしてみました。


参考