代码急转弯——Tupper（塔珀自指公式）

1、塔珀自指公式

塔珀自指公式（Tupper’s Self-Referential Formula）来自于Jeff Tupper在2001年发表的一篇论文，涉及到一个函数 $f(x,y)$ ,一个参数 $k$，和一个图片 $a$。
$f(x,y)$ 对 $a$ 中的每一个点 $(x,y)$ 进行判断，决定该点的颜色是0还是1，因此决定了 $a$ 的图案。如果给定不同的 $k$，就会得到不同的图案。令人感到惊奇的是，有一个特别的 $k$，它会使得所得 $a$ 竟然是 $f$ 用来做判断的公式本身。下面就是塔珀公式，
$$\frac{1}{2}<\lfloor mod(\lfloor \frac{y}{17}\rfloor 2^{-17\lfloor x\rfloor -mod(\lfloor y\rfloor ,17)},2)\rfloor$$
其中，$\lfloor x\rfloor$ 是向下取整函数，$mod(a,b)$ 计算 $a$ 除以 $b$ 的余数，$0\le x\le 105$ ，$k\le y\le k+16$ 。$k$ 会是一个什么样的数？
k = 960939379918958884971672962127852754715004339660129306651505519271702802395266424689642842174350718121267153782770623355993237280874144307891325963941337723487857735749823926629715517173716995165232890538221612403238855866184013235585136048828693337902491454229288667081096184496091705183454067827731551705405381627380967602565625016981482083418783163849115590225610003652351370343874461848378737238198224849863465033159410054974700593138339226497249461751545728366702369745461014655997933798537483143786841806593422227898388722980000748404719

print(len(str(k))
543

好家伙，足足有543位。我们来看看图片 $a$，体会其中的“自指”的含义。

2、如何绘制图片 $a$ ？

import numpy as np
import matplotlib.pyplot as plt

def Tupper_self_referential_formula(k): 
    aa = np.zeros((17,106))
    def f(x, y):
        y += k
        a1 = 2**-(-17*x - y%17)
        a2 = (y // 17) // a1
        return 1 if a2 % 2 > 0.5 else 0
 
    for y in range(17):
        for x in range(106):
            aa[y, x] = f(x, y)
    
    return aa[:,::-1]

k = 960939379918958884971672962127852754715004339660129306651505519271702802395266424689642842174350718121267153782770623355993237280874144307891325963941337723487857735749823926629715517173716995165232890538221612403238855866184013235585136048828693337902491454229288667081096184496091705183454067827731551705405381627380967602565625016981482083418783163849115590225610003652351370343874461848378737238198224849863465033159410054974700593138339226497249461751545728366702369745461014655997933798537483143786841806593422227898388722980000748404719
aa = Tupper_self_referential_formula(k)
plt.figure(figsize=(15,10))
plt.imshow(aa,origin='lower')

运行结果就是上面的图片。如果删除 $k$ 前面的第一位数字9，那么将得到另一个完全不同的图片。

注意，不能直接按照公式编码计算，

a1 = 2**(-17*x - y%17)
a2 = (y // 17) * a1

这样做会溢出。

3、如何计算参数 $k$ ？

在给定图片 $a$ 后，如何找到 $k$ ？
只要对 $a$ 中的各元素进行排列，得到一个2进制数（因为 $a$ 的元素都是1或者0），转换为10进制，再乘以17就得到了 $k$。排列顺序为，第0列自下而上，然后是第1列，也是自下而上，把各元素排列到一起，成为一个长长的二进制数字的字符串，如下例所示，

array([['a', 'b', 'c', 'd'],
       ['e', 'f', 'g', 'h'],
       ['i', 'j', 'k', 'l']], dtype='<U1')

排序的结果为 ieajfbkgclhd，从中可观察出排列的规则。下面是计算 $k$ 代码，

b1 = aa[::-1] # 上下颠倒
b2 = b1.T.reshape(-1) # 降维成一维数组
b3 = [str(int(k)) for k in b2] # 各元素转变为0/1字母
b4 = ''.join(b3) # 拼接成二进制串
b5 = int(b4,2)*17 # 二进制串变为十进制数，再乘以17

''' 验证一下 '''
print('k - b5 =', k - b5)

现在用计算所得的 $b5$ 与 $k$ 相减，结果为零。

k - b5 = 0

【下一个谜题】life游戏——最简元胞自动机