Python向左,数学向右:乌拉姆的素数研究


1. 引子


       我是一只猫。听起来有点惊悚,但这是千真万确的事实。

       那一年——这样说,并不一定意味着一个历史故事的开始,因为我没有过去和将来的概念。我可以在时间维度上任意移动,就像从东走到西或者从西走到东那样自然。

       那一年,应该是公元1963年。当时我正在美国阿拉莫斯国家试验室的一间会议室的窗台上享受着温暖的午后秋阳,会议室里一群人正在开会。这里处于新墨西哥州杰姆斯山森林的包围之中,空气清新且日照充足,因此我经常从欧洲或者亚洲来这里度假。

       会议主题冗长且毫无新意,有些听众开始打起了哈欠。期间,有个一直埋头写写画画的家伙引起了我的注意。他叫乌拉姆,波兰裔物理学家和数学家,我很早就认识他。

       乌拉姆在这间国家实验室的最大成就,是和氢弹之父泰勒合作提出泰勒-乌拉姆设计(Teller–Ulam design),为核聚变武器研发提供了基础支撑,并最终导致第一颗氢弹的成功试爆。


2. 灵光乍现


       “喵——”

       我悄悄走到乌拉姆身后,打了一个招呼,用胡须轻轻碰了一下乌拉姆右侧的面颊。乌拉姆先用左手拍了拍趴在肩头的我,然后回过头来。

       “嘘——小声点。你知道吗,我可能发现了一个天大的秘密,这绝对比泰勒-乌拉姆设计更有意义。让泰勒那个讨厌的家伙见鬼去吧。”

       乌拉姆面色潮红,两眼放光,看起来有点亢奋。“曼哈顿计划”计划之后,由于内心抗拒继续研制热核武器,乌拉姆在氢弹研制中几乎被人遗忘,心里多少有一点忿忿不平。

       “哦?说来听听。”

       “你看,从1开始,逆时针转圈将连续的整数写在10x10的格子里。”乌拉姆拿起桌上的那张草纸,右手食指指向中间位置,然后开始转圈。纸上是他画的方格、数字和斜线,就像下图这个样子。
在这里插入图片描述

乌拉姆手绘的100以内的素数分布图(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)

       运行结束,打开图片之前,乌拉姆的右手一直放在我的后背上。我能感觉到他的紧张和期待。看到图片的时候,他先是盯着看了三秒钟,然后双手猛然击掌,兴奋地大喊:“果然存在斜线!快,再来几张更大的。”


在这里插入图片描述

1万以内的素数分布图(100x100),右图为放大3倍后的效果

       代码运行效率比想象的要好很多,轻松生成了9万、36万和100万以内的素数分布图。

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

       现在我一点儿都不担心乌拉姆会提出更大的要求,因为这个速度已经远远超出了全世界的计算机算力。而乌拉姆没有想到,他无意中发现了被后人称为“乌拉姆现象”的素数分布特点。

在这里插入图片描述

9万以内的素数分布图(300x300)

在这里插入图片描述

36万以内的素数分布图(600x600)

在这里插入图片描述

1百万以内的素数分布图(1000x1000)


4. 悖论


       短暂的激动之后,乌拉姆盯着这张1百万以内的素数分布图,陷入了沉思。良久,他开口说话了。

       “能否再帮我算一下1千万内和1亿以内的素数的个数?”

       “当然。”

       我把查找素数的函数find_prime_list()稍微改造了一下,将查找范围分成多个区块。这样做主要是为了降低内存消耗,另外也便于并行计算——实际上我并没有启用并行计算,因为逐个区块计算就已经足够快了。

---------------------------------
查找1千万以内的素数耗时0.365秒
共找到664579个素数
---------------------------------
查找1亿以内的素数耗时2.862秒
共找到5761455个素数

       乌拉姆瞅了一眼屏幕,低头在一张草纸上写下了这样一个公式。

lim ⁡ x → + ∞ π ( x ) x = 0 \lim\limits_{x\to+\infty}{\frac{\pi(x)}{x}} = 0 x+limxπ(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向左,数学向右呢,乌拉姆博士?”

相关推荐
简介 笔者当初为了学习JAVA收集了很多经典源码源码难易程度分为初级、中级、高级等详情看源码列表需要可以直接下载! 这些源码反映了那时那景笔者对未来盲目对代码热情、执着对IT憧憬、往!此时此景笔者只专注Android、Iphone等移动平台开发看着这些源码心中有万分感慨写此文章纪念那时那景! Java 源码包 Applet钢琴模拟程序java源码 2个目标文件提供基本音乐编辑功能。编辑音乐软件朋友这款实例会对你有所帮助。 Calendar万年历 1个目标文件 EJB 模拟银行ATM流程及操作源代码 6个目标文件EJB来模拟银行ATM机流程及操作获取系统属性初始化JNDI取得Home对象引用创建EJB对象并将当前计数器初始化调用每一个EJB对象count()方法保证Bean正常被激活和钝化EJB对象是用完毕从内存中清除从账户中取出amt如果amt>账户余额抛出异常一个实体Bean可以表示不同数据实例我们应该通过主键来判断删除哪个数据实例…… ejbCreate函数用于初始化一个EJB实例 5个目标文件演示Address EJB实现 创建一个EJB测试客户端得到名字上下文查询jndi名通过强制转型得到Home接口getInitialContext()函数返回一个经过初始化上下文用clientgetHome()函数调用Home接口函数得到远程接口引用用远程接口引用访问EJB。 EJB中JNDI使用源码例子 1个目标文件JNDI使用例子有源代码可以下载参考JNDI使用初始化Context,它是连接JNDI树起始点查找你要对象打印找到对象关闭Context…… ftp文件传输 2个目标文件FTP目标是(1)提高文件共享性(计算机程序和/或数据)(2)鼓励间接地(通过程序)使用远程计算机(3)保护用户因主机之间文件存储系统导致变化(4)为了可靠和高效地传输虽然用户可以在终端上直接地使用它但是它主要作用是供程序使用。本规范尝试满足大型主机、微型主机、个人工作站、和TACs 不同需求。例如容易实现协议设计。 Java EJB中有、无状态SessionBean两个例子 两个例子无状态SessionBean可会话Bean必须实现SessionBean获取系统属性初始化JNDI取得Home对象引用创建EJB对象计算利息等;在有状态SessionBean中用累加器以对话状态存储起来创建EJB对象并将当前计数器初始化调用每一个EJB对象count()方法保证Bean正常被激活和钝化EJB对象是用完毕从内存中清除…… Java Socket 聊天通信演示代码 2个目标文件一个服务器一个客户端。 Java Telnet客户端实例源码 一个目标文件演示Socket使用。 Java 组播组中发送和接受数据实例 3个目标文件。 Java读写文本文件示例代码 1个目标文件。 java俄罗斯方块 一个目标文件。 Java非对称加密源码实例 1个目标文件 摘要:Java源码,算法相关,非对称加密   Java非对称加密源程序代码实例本例中使用RSA加密技术定义加密算法可用 DES,DESede,Blowfish等。   设定字符串为“张三你好我是李四”   产生张三密钥对(keyPairZhang)   张三生成公钥(publicKeyZhang)并发送给李四,这里发送是公钥数组字节   通过网络或磁盘等方式,把公钥编码传送给李四李四接收到张三编码后公钥,将其解码李四用张三公钥加密信息并发送给李四张三用自己私钥解密从李四处收到信息…… Java利用DES私钥对称加密代码实例 同上 java聊天室 2个目标文件简单。 java模拟掷骰子2个 1个目标文件输出演示。 java凭图游戏 一个目标文件简单。 java求一个整数因子 如题。 Java生成密钥实例 1个目标文件 摘要:Java源码,算法相关,密钥   Java生成密钥、保存密钥实例源码通过本源码可以了解到Java如何产生单钥加密密钥(myKey)、产生双钥密钥对(keyPair)、如何保存公钥字节数组、保存私钥到文件privateKey.dat、如何用Java对象序列化保存私钥,通常应对私钥加密后再保存、如何从文件中得到公钥编码字节数组、如何从字节数组解码公钥。 Java数据压缩与传输实例 1个目标文件 摘要:Java源码,文件操作,数据压缩,文件传输   Java数据压缩与传输实例可以学习一下实例化套按字、得到文件输入流、压缩输入流、文件输出流、实例化缓冲
<p> <b><span style="background-color:#FFE500;">【超实用课程内容】</span></b> </p> <p> <br /> </p> <p> <br /> </p> <p> 本课程内容包含讲解<span>解读Nginx基础知识</span><span>解读Nginx核心知识、带领学员进行</span>高并发环境下Nginx性能优化实战让学生能够快速将所学融合到企业应用中。 </p> <p> <br /> </p> <p style="font-family:Helvetica;color:#3A4151;font-size:14px;background-color:#FFFFFF;"> <b><br /> </b> </p> <p style="font-family:Helvetica;color:#3A4151;font-size:14px;background-color:#FFFFFF;"> <b><span style="background-color:#FFE500;">【课程如何观看?】</span></b> </p> <p style="font-family:Helvetica;color:#3A4151;font-size:14px;background-color:#FFFFFF;"> PC端<a href="https://edu.csdn.net/course/detail/26277"><span id="__kindeditor_bookmark_start_21__"></span></a><a href="https://edu.csdn.net/course/detail/27216">https://edu.csdn.net/course/detail/27216</a> </p> <p style="font-family:Helvetica;color:#3A4151;font-size:14px;background-color:#FFFFFF;"> 移动端CSDN 学院APP(注意不是CSDN APP哦) </p> <p style="font-family:Helvetica;color:#3A4151;font-size:14px;background-color:#FFFFFF;"> 本课程为录播课课程永久有效观看时长大家可以抓紧时间学习后一起讨论哦~ </p> <p style="font-family:"color:#3A4151;font-size:14px;background-color:#FFFFFF;"> <br /> </p> <p class="ql-long-24357476" style="font-family:"color:#3A4151;font-size:14px;background-color:#FFFFFF;"> <strong><span style="background-color:#FFE500;">【学员专享增值服务】</span></strong> </p> <p class="ql-long-24357476" style="font-family:"color:#3A4151;font-size:14px;background-color:#FFFFFF;"> <b>源码开放</b> </p> <p class="ql-long-24357476" style="font-family:"color:#3A4151;font-size:14px;background-color:#FFFFFF;"> 课件、课程案例代码完全开放给你你可以根据所学知识自行修改、优化 </p> <p class="ql-long-24357476" style="font-family:"color:#3A4151;font-size:14px;background-color:#FFFFFF;"> 下载方式电脑登录<a href="https://edu.csdn.net/course/detail/26277"></a><a href="https://edu.csdn.net/course/detail/27216">https://edu.csdn.net/course/detail/27216</a>播放页面侧点击课件进行资料打包下载 </p> <p> <br /> </p> <p> <br /> </p> <p> <br /> </p>
©️2020 CSDN 皮肤主题: 代码科技 设计师:Amelia_0503 返回首页
实付 15.20元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、C币套餐、付费专栏及课程。

余额充值