Close

每一個圓周率裏都書寫了無數卷漢姆萊特——但你需要一些耐心

之前看到一則 math.stackexchange 上面的問題:圓周率 \pi 是否窮盡了所有的數字的排列組合,這個答案是否定的,所有數字的排列組合是 \aleph_1 的,這個證明可以用對角線方法來實現:記 \pi 的第n位是P_n,可以構造這樣的一個序列 S_n 滿足 S_n\neq P_{2n},於是我們可以得知 S_n 必然不是 P_n 的一個子列。這是由於自第k位開始的 P_n 的子列的第k位必然和 S_k 不相等。

但是換一種說法,\pi 是否窮盡了所有有限長度的數字的排列組合,答案是肯定的。事實上,所有有限長度的數字序列的個數是可以枚舉的(秩數 \aleph_0 的)。

又有一種說法,如果一隻猴子可以有無窮多的時間區敲打打字機,那麼也會有一天它恰好打出了漢姆萊特這樣的鉅著。理由和上面一樣,因爲漢姆萊特是一本有限長度的書(這種說法不得不讓人想到博爾赫斯的《沙之書》不是嘛)。

如果把上述兩個結論綜合起來的話,就是:

在圓周率的深處沉眠了無數卷的漢姆萊特……

當然,前提是我們能夠把數字轉義成拉丁字母。

第三種說法是,問題的核心往往不是有沒有,而是在哪裏,「每年有一半的廣告開銷是有實際價值的,而問題是我不知道那是哪一半。」亨利·福特是這樣表達的。同樣的,圓周率中的漢姆萊特究竟在哪裏,趁着夜色迷離,讓我們做一次不太有意義的嘗試。

胡思亂想時間

計算 \pi 一直是各種算法書當中津津樂道的話題,而一個真正的數學愛好者選擇的方法,當然是到 PiDay 網站上直接下載圓周率前 100 萬位啦。用能幫助你準時下班的 Python,我們這麼做

import requests
from bs4 import BeautifulSoup
page = requests.get('https://www.piday.org/million/').content
soup = BeautifulSoup(page).find('div', {'id': 'million_pi'})
PIstr = soup.getText().replace('.', '')

考慮一下如何將數字轉譯成英文字母的話,ASCII 碼可能是一個比較合理簡單的方式,排除掉所有無法打印的字符之後,我們可以把所有字符壓縮到 0~94 當中,爲了方便處理,我們也把所有超過 94 的數字轉化爲空白,那麼:

def digits_to_char(d: str) -> str:
	num: int = int(d)
	if num > 94:
		return ' '
	return chr(num + 32)

另一方面我們可以把前述得到的 PIstr 打斷成兩位一組的序列

def chunk_string(s: str) -> List[str]:
	return list(s[i: i + 2] for i in range(0, len(s), 2))

然後試着讀一下可能包含着漢姆萊特的大作:

pirnt(''.join(digits_to_char(i) for i in chunk_string(PIstr)))

# ?I[:UZ }7t^`As; "xI 0} ES%r)j~M|>n0H^v4y ^pBrUJ+fc...

???跟標題講得好像不太一樣???事實上,在轉譯出來的前 50 萬位字符中別說 Hamlet 了,連火腿(ham)也沒有,當然也沒有 spam,雖然明明我們用了 Python(用梗過度)。

假裝理性時間

難以卒讀的原因主要有以下兩個問題:

  • 由於 \pi 中每個數字出現的比例大致相同,導致在轉譯出的文字中每個字符出現的比率非常接近,這與自然英語文本中的字符頻率完全不一致。
  • 在自然語言中,每個字符與其後繼的字符並非相互獨立的事件,例如,在字符q後面,99.95% 的概率是字符u(我有一個叫 Ashfaq 的朋友來着的……),這點在轉譯序列中也沒有得到體現。

但是莫要緊,統計漲落會爲我們處理這些問題,只要有耐心(可能還是要劃掉),也是會出現某個局部乖離均值的區域,而我們的漢姆萊特就在那裏等待我們。

數學的 TL; DR

那麼,考慮一本由 a 個字母組成的字典,其中存在一個長度爲 k 的單詞,現在在一本字數爲 n 的書中找到這個單詞的概率是多少呢?不妨說,讀完 n 個字才找到這個單詞的概率是多少呢?單詞在第 n 位首次出現的概率記爲 P_{(n, k)},單詞在 n 位時已經出現過的概率記爲 S_{(n,k)},顯然我們有:

P_{(n, k)}=(1-S_{(n -k, k)})(\frac1a)^k\\ S_{(n, k)}=\sum\limits_{i=0}^nP_{(i,k)}\\ P_{(n, k)} = 0\iff n<k

或者 S_n-S_{n-1}=(\frac1a)^k(1-S_{n-k}),這個算式用代碼表示的話:

def s(n: int, k: int) -> float:
	if n < k:
		return 0
	return s(n - 1, k) + 0.01 ** k * (1 - s(n - k, k))

在我們的例子裏,字典的大小是 100 個字符,單詞的長度大約是130,000 字符(Hamlet 的長度),爲了確保 95% 的置信度找到漢姆萊特,可能真的需要很久。由於上面的代碼並沒有Tail call optimized,無法幫助我們實際計算出需要讀取的圓周率的位數,只能說非常遺憾了。

因此,今天晚上也是既沒有讀書,也沒有學習的晚上,乾杯!

Leave a Reply

Your email address will not be published. Required fields are marked *

This site uses Akismet to reduce spam. Learn how your comment data is processed.

沪ICP备19039127号-1