1. 引子

我是一只貓，聽起來有點驚悚，但這是千真萬確的事實，

那一年——這樣說，并不一定意味著一個歷史故事的開始，因為我沒有過去和將來的概念，我可以在時間維度上任意移動，就像從東走到西或者從西走到東那樣自然，

那一年，應該是公元1963年，當時我正在美國阿拉莫斯國家試驗室的一間會議室的窗臺上享受著溫暖的午后秋陽，會議室里一群人正在開會，這里處于新墨西哥州杰姆斯山森林的包圍之中，空氣清新且日照充足，因此我經常從歐洲或者亞洲來這里度假，

會議主題冗長且毫無新意，有些聽眾開始打起了哈欠，期間，有個一直埋頭寫寫畫畫的家伙引起了我的注意，他叫烏拉姆，波蘭裔物理學家和數學家，我很早就認識他，

烏拉姆在這間國家實驗室的最大成就，是和氫彈之父泰勒合作提出泰勒-烏拉姆設計（Teller–Ulam design），為核聚變武器研發提供了基礎支撐，并最終導致第一顆氫彈的成功試爆，

2. 靈光乍現

“喵——”

我悄悄走到烏拉姆身后，打了一個招呼，用胡須輕輕碰了一下烏拉姆右側的面頰，烏拉姆先用左手拍了拍趴在肩頭的我，然后回過頭來，

“噓——小聲點，你知道嗎，我可能發現了一個天大的秘密，這絕對比泰勒-烏拉姆設計更有意義，讓泰勒那個討厭的家伙見鬼去吧，”

烏拉姆面色潮紅，兩眼放光，看起來有點亢奮，“曼哈頓計劃”計劃之后，由于內心抗拒繼續研制熱核武器，烏拉姆在氫彈研制中幾乎被人遺忘，心里多少有一點忿忿不平，

“哦？說來聽聽，”

“你看，從1開始，逆時針轉圈將連續的整數寫在10x10的格子里，”烏拉姆拿起桌上的那張草紙，右手食指指向中間位置，然后開始轉圈，紙上是他畫的方格、數字和斜線，就像下圖這個樣子，

“把其中的素數用紅筆標上顏色，我突然發現，這些素數居然都出現在斜線上！”

這時我才明白，烏拉姆為什么會如此激動了，千百年來，從古希臘的歐幾里得、古埃及的埃拉托色尼，到近代德國的高斯和黎曼、法國的勒讓德，人類一直在嘗試解決素數的多少以及如何分布的問題，并試圖找到一個可以計算出所有素數的簡單公式，盡管這些數學天才們窮盡了可能的手段，但是并沒有得到期望的結果，烏拉姆也許覺得自己找到了素數的分布規律，

“這只是100以內的素數，說明不了什么吧？”

“是的，這張草紙太小了，要是給我一張足夠大的紙，也許可以從中發現什么，”

“是個好主意，不過，一個數字一個數字的填寫，效率太低了，為什么不試試Python呢？”

在那一瞬間，忽然想起來我在北大混日子的時光，那時候，人們喜歡叫我哲學貓，其實，除了喜歡旁聽哲學和藝術，我偶爾也去聽聽計算機課，在那里，我多少了解了一點Python的知識，

“Python？蟒蛇也研究數學嗎？”

看著烏拉姆一臉懵逼，我意識到自己犯了一個錯誤：Python直到1989年才誕生，那時烏拉姆已經去世五年了，想到這一點，我心里微微有一些傷感，

“這是很多年后一個名叫吉多·范羅蘇姆的荷蘭人搞出來的編程語言，可以協助處理一些重復的、繁瑣的事務性作業，這個語言最初無人問津，沒想到幾十年后卻風靡全球，火得燙手，沒關系，我可以提前拿來用一下，幫你搞搞這個，”我指了一下他畫的草圖，

對于我能從未來臨時借用Pyhton，烏拉姆并未表現出任何的驚奇或不解，這倒不是因為他是頂尖的物理學家，可以理解四維生物的生存狀態，而是早在我協助薛定諤做量子實驗的時候，他就已經認識我了，薛定諤是奧地利人，烏拉姆生于波蘭，從格拉茨開車到華沙只有800公里，這段距離，擱在中國就算是同鄉了，加上又是同行，盡管年齡相差二十幾歲，兩人卻是忘年之交，

3. 烏拉姆現象

“這是全美目前最快的電腦，每秒鐘可以完成3000次計算，”

烏拉姆將我帶進機房，指著一臺美國資料設備公司（DEC）生產的POP-5小型機說，

“你這里，也許只有電源是可用的，”

我苦笑一下，打開了自己的電腦，邊敲鍵盤，邊和正在背身操作咖啡機的烏拉姆討論如何繪制素數分布圖，

“我打算這樣做：先搞一個可以計算n以內素數的函式，再搞一個將1到n逆時針轉圈排列成生成方陣的函式，呼叫這兩個函式，分別得到素數表和數字方陣，然后將方陣中的素數置為255，合數置為0……”

“哇嗚，這主意聽起來相當美妙！最后得到一張灰度圖，白色像素點表示素數，黑色像素點表示合數，素數的分布規律——如果有的話，就會一目了然，”

這家伙腦子的確好使，沒上過計算機課，也能搞得門兒清，我抬頭斜著眼瞅了烏拉姆一下，他知趣地用咖啡杯堵住了嘴，

既然他都明白了，就不用繼續說我的計劃了，直接寫代碼吧，先把需要的模塊匯入，萬一要是沒有安裝所需模塊，還得回到幾十年后才能找到可用的模塊，

import numpy as np
from PIL import Image

幸好我的電腦上早已安裝了NumPy和pillow這兩個模塊，開局非常順利，再定義一個函式，找出不大于n的所有素數，

def find_prime_list(n):
    """回傳不大于n的素陣列成的numpy陣列"""
    
    nums = np.arange(n+1) # 生成0到n的陣列
    nums[1] = 0 # 陣列第1和元素置0，從2開始，都是非0的
    k, m = 0, pow(n+1, 0.5)
    
    while True: # 回圈
        primes = nums[nums>0] # 找出所有非0的元素
        if primes.any(): # 如果存在不為0的元素
            p = primes[k] # 則第一個元素為素數
            k += 1
            nums[2*p::p] = 0 # 這個素數的所有倍數，都置為0
            if p > m: # 如果找到的素數大于上限的平方根
                break # 則退出回圈
        else:
            break # 全部0，也退出回圈
    
    return nums[nums>0]

將1到n逆時針轉圈排列成生成方陣，稍微有點一點難度，看我想得出神，烏拉姆建議說：“為了簡化設計，我們可以限定n必須是一個平方數，這樣最終得到的方陣的行列數都是n的算術平方根，你寫程式也容易一點，”

我接受了烏拉姆的建議，有了這個約束條件，思路一下就清晰了：假定n的算術平方根為side，不考慮效率的話，可以先生成一個side*side的二維串列，其元素全部為None；如果side為奇數，則n位于二維串列的右下角，n-1位于n的左側，否則，n位于二維串列的左上角，n-1位于n的右側；然后從n開始逆序將數字按照順時針方向寫入串列，

def get_square(n):
    """將從1開始的n個連續的整數排成一個串列"""
    
    side = int(pow(n, 1/2)) # 方陣邊長
    
    if side%2:
        row, col = side-1, side-1
        direct = 'left'
    else:
        row, col = 0, 0
        direct = 'right'
        
    result = [[None for j in range(side)] for i in range(side)]
    
    for i in range(n, 0, -1):
        result[row][col] = i
        
        if direct == 'right':
            if col+1 == side or result[row][col+1]: # 如果不能繼續向右，則向下:
                row += 1
                direct = 'down'
            else: # 否則繼續向右
                col += 1
        elif direct == 'down':
            if row+1 == side or result[row+1][col]: # 如果不能繼續向下，則向左
                col -= 1
                direct = 'left'
            else: # 否則繼續向下
                row += 1
        elif direct == 'left':
            if col-1 < 0 or result[row][col-1]: # 如果不能繼續向左，則向上
                row -= 1
                direct = 'up'
            else: # 否則繼續向左
                col -= 1
        elif direct == 'up':
            if row-1 < 0 or result[row-1][col]: # 如果不能繼續向上，則向右
                col += 1
                direct = 'right'
            else: # 否則繼續向上
                 row -= 1
    
    return np.array(result)

有了素數表，有了方陣，只剩下替換操作，就簡單多了，

def plot_prime(side):
    """繪制不大于side*side的質數分布圖"""
    
    n = side * side
    square = get_square(n)
    primes = find_prime_list(n)
    
    im_arr = np.isin(square, primes).astype(np.uint8) * 255
    im = Image.fromarray(im_arr)
    im.save('%dx%d.jpg'%(side, side))

寫到這里的時候，我抬起頭，注意到烏拉姆一臉驚訝的表情，

“這個方陣里的素數是怎么被替換成255的？我無法理解np.isin()這個神奇的東西，感覺創造Python的那個荷蘭人，智商比愛因斯坦還高！”

烏拉姆平生最服氣愛因斯坦，一旦遇到牛人，總喜歡和愛因斯坦作比較，

“哦，你說的是NumPy，那的確是一個很牛叉的東西，不過你不用因此而自卑，NumPy不是吉多·范羅蘇姆寫的，而是由一個偉大的團隊研發的，”

“我就說嘛，這怎么可能是一個人的力量可以做到的呢，其實大家都明白，氫彈也不是泰勒那家伙一個人搞出來的，”

看來烏拉姆對于泰勒貪功一事一直無法釋懷，我沒有接烏勒姆的話茬，而是轉向筆記本電腦，開始了第一次的運行測驗，先從100x100的方陣開始，

plot_prime(100)

運行結束，打開圖片之前，烏拉姆的右手一直放在我的后背上，我能感覺到他的緊張和期待，看到圖片的時候，他先是盯著看了三秒鐘，然后雙手猛然擊掌，興奮地大喊：“果然存在斜線！快，再來幾張更大的，”

代碼運行效率比想象的要好很多，輕松生成了9萬、36萬和100萬以內的素數分布圖，

plot_prime(300)
plot_prime(600)
plot_prime(1000)

現在我一點兒都不擔心烏拉姆會提出更大的要求，因為這個速度已經遠遠超出了全世界的計算機算力，而烏拉姆沒有想到，他無意中發現了被后人稱為“烏拉姆現象”的素數分布特點，

4. 悖論

短暫的激動之后，烏拉姆盯著這張1百萬以內的素數分布圖，陷入了沉思，良久，他開口說話了，

“能否再幫我算一下1千萬內和1億以內的素數的個數？”

“當然，”

我把查找素數的函式find_prime_list()稍微改造了一下，將查找范圍分成多個區塊，這樣做主要是為了降低記憶體消耗，另外也便于并行計算——實際上我并沒有啟用并行計算，因為逐個區塊計算就已經足夠快了，

---------------------------------
查找10000000以內的素數耗時0.365秒
共找到664579個素數
---------------------------------
查找1億以內的素數耗時2.862秒
共找到5761455個素數

烏拉姆瞅了一眼螢屏，低頭在一張草紙上寫下了這樣一個公式，

lim ? x → + ∞ π ( x ) x = 0 \lim\limits_{x\to+\infty}{\frac{\pi(x)}{x}} = 0 x→+∞lim?xπ(x)?=0

“這里的 π \pi π可不是圓周率，而是表示不大于自然數 x x x的素數的個數，”他左手端起早已涼透的半杯咖啡，用右手中指的關節輕扣草紙，

“歐拉和勒讓德早已證明了當 x x x趨于無窮大時， π ( x ) \pi(x) π(x)與 x x x之比等于0——這意味，即便素數有無窮多個，也不應該出現在自然數中，因為素數在自然數中出現的概率為0，”

“這可真是一個有趣的悖論，”

我不知道烏拉姆要表達什么，只好敷衍地回應了一句，

沉默有頃，烏拉姆幽幽地說："Python在幾秒鐘內就可以計算出，當 x x x趨近1億時，素數出現的概率降低到5.76%，這是一個強大的工具，可以瞬間完成單憑人力無法完成的計算任務，但是，卻剝奪了人類享受在黑暗中探索和思考的樂趣，”

“伸手不見五指，沒有方向，沒有終點，那不是痛苦嗎？”

“不，那是一種享受，天堂向左，數學向右，”

“好吧，我想我能理解你的意思，時間的不可逆，和對時光逆轉的渴望，令人類創造了一些有趣的概念，比如天堂，還有Python，為什么不是Python向左，數學向右呢，烏拉姆博士？”