偶然間,看到吳軍老師的《浪潮之巔(第四版)》裡有講到這麼一個故事。
Google曾經在加州的101高速公路上用大廣告牌登了這樣的廣告:
{ 無理數e中前十位連續的素數 }.com
你如果知道這個答案(7427466391.com),就可以通過上述網址進入到Google的招聘網站。而能夠計算出這道題,要很聰明。
“很聰明"?吳軍老師的這句話倒是讓人來興趣了,自己也湊性藉助計算機的力量琢磨琢磨下這個證明題。
思路
大家都知道,e是無理數,換言之,無限不迴圈小數。世間總有些數學愛好人士樂此不疲地用計算機來推算這些無理數的小數位,π如此,e更不例外。所以,
第一步,在一些權威的數學網站上找到e的小數位資料;
這時候,就是谷歌大法好了,來到了e的維基英文百科頁面,留意到了下面e的小數表達形式後面的“A001113”超連結,有蹊蹺,便點進去看下。
接著便是這個OEIS站點的A001113頁面了,原來是數論方面的一個很權威的資料庫網站。
帶著探寶的眼光在這個網頁上下左右掃視,總算發現了目標。LINKS(連結)的第一個條目便是了。
N. J. A. Sloane, Table of 50000 digits of e labeled from 1 to 50000 [based on the ICON Project link below]
作者也是創辦了這個OEIS組織的N.J.A Sloane,看來是個此領域的大人物了,抱歉有眼不識泰山。
點進去這個Table of 50000 digits of e labeled from 1 to 50000,頁面結果讓人很愉悅。直接的純文字資料,左列為小數位數,右列為對應數值。到時候拿到資料時簡單處理下即可展開後續工作,希望這道證明題的答案就在這50000個數中,不然得擴大範圍了。
到此,e無理數的小數位資料有了,便可以展開下一步工作了。
第二步,如何判斷一個數是素數?
上小學的時候開始,數學老師就教導我們,素數的定義是指一個除了1和該數自身外,無法被其他自然數整除的大於1的自然數。所以自然,判斷思路便是對一個大於1的自然數n,依次判斷2 → n能否整除n,如果發現一個數能整除n,那麼n不是素數,否則是。另外,考慮到對稱性,我們也不必一直遞增到n,如對於2*3和3*2,6/2和6/3皆判定6為合數,但是遞增到2時就已經是充分了, 故不必再考慮3。
Python程式碼如下:
def isPrime(n: int) -> bool:
if n <= 1:
return False
i = 2
# Make use of symmetry. For example, 2*3=6, 3*2=6
while i * i < n:
if n % i == 0:
return False
i += 1
return True
在網上查閱素數相關的資料的時候,發現數論裡有個素數分佈規律也可以拿來判斷素數。來源——素數判定演算法
素數分佈規律:
當 n >= 5 時,如果n為素數,那麼 n % 6 = 1 || n % 6 = 5, 即 n 一定出現在 6x(x ≥ 1)兩側。換句話說,任意一個素數都可以被表示為 6x ± 1, x ∈ N 的形式。
證明:
把6x附近的數用以下方式表示:
......(6x - 1), 6x, 6x + 1, 2(3x + 1), 3(2x + 1), 2(3x + 2), 6x + 5, 6(x+1)......
不在6x兩側的數為:2(3x + 1), 3(2x + 1), 2(3x + 2), 它們不是素數,所以素數出現在6x的兩側。
用Python程式碼實現如下,時間複雜度上比前一種相差無幾,不過對於我們的證明來說,是夠用了。
def isPrime(n: int) -> bool:
if n <= 1:
return False
if n <= 3:
return True
# For case of 2(3x + 1), 3(2x + 1), 2(3x + 2)
if n % 2 == 0 or n % 3 == 0:
return False
# For case of the multiplication of prime numbers
i = 5
while i * i < n:
if n % i == 0 or n % (i + 2) == 0:
return False
i += 6
return True
此外,瞭解到密碼學裡面判定素數有很大的用處,比如著名的RSA演算法。在判斷素數演算法方面,並沒有採用上面的時間複雜度較高的簡單取餘演算法,而是Fermat小定理、Miller-Rabin演算法及Solovay-Strassen演算法等,理解起來較為吃力,具體可參考下這篇文章——PyCrypto密碼學庫原始碼解析(一)隨機數和大素數生成, 以及上面那篇。
至此,必要的材料都準備好了,可以進行最後一步了。
第三步,for迴圈e的小數位資料判斷第一個10位長的素數。
開門見山,直接原始碼丟擲來先。
具體思路:先使用requests庫獲取e小數位資料,然後轉存為檔案便於逐行讀取,for迴圈逐行讀取每一小數位資料,進行切片操作,整理成證明所需的10位整數,得到總數量為49991的有序列表,再借助素數判定函式逐個判定這些10位整數,最後得到答案——7427466391。
import requests
import re
response = requests.get('https://oeis.org/A001113/b001113.txt')
# Save sequence to a file for later use
out_file = open('digits_of_e.txt', 'w')
print(response.text, file=out_file)
queue = []
container = ''
counter = 0
in_file = open('digits_of_e.txt', 'r')
list_in_file = list(in_file)
for index, line in enumerate(list_in_file):
segments = list_in_file[index:index+10]
# get lines at a batch of 10 lines
for segment in segments:
matchObj = re.match(r'([\d]*) (\d).*', segment, re.I)
counter += 1
if counter <= 10:
container += matchObj.group(2) if matchObj != None else ''
else:
counter = 1
if len(container) == 10:
queue.append(container)
container = matchObj.group(2) if matchObj != None else ''
in_file.close()
print(len(queue)) # 49991 indeed
def isPrime(n: int) -> bool:
# Prime number definition version:
'''
if n <= 1:
return False
i = 2
# Make use of symmetry. For example, 2*3=6, 3*2=6
while i * i < n:
if n % i == 0:
return False
i += 1
return True
'''
# Distribution pattern of prime number version:
if n <= 1:
return False
if n <= 3:
return True
# For case of 2(3x + 1), 3(2x + 1), 2(3x + 2)
if n % 2 == 0 or n % 3 == 0:
return False
# For case of the multiplication of prime numbers
i = 5
while i * i < n:
if n % i == 0 or n % (i + 2) == 0:
return False
i += 6
return True
result = None
for num in queue:
if isPrime(int(num)):
print(num)
result = num
break
print(result == '7427466391')
print(isPrime(7427466391))
執行結果:
賓果!
結束
證明題解答完畢,我們得走個儀式感,訪問訪問這個網站——7427466391.com。
結果是502錯誤……
行吧,看來這個站點早被人家拋棄了,畢竟這個高速公路廣告也是谷歌在2004年搞的惡作劇。
最後,把原始碼整理了個Kaggle Notebook版本。歡迎查閱!
First10DigitPrimeFoundInConsecutiveDigitsOfE | Kaggle