第100,000,000,000个素数

2020-05-27 07:47:17

这个专栏讲述了一个故事,它有一个寓意。它本身与J没有直接关系,J是这些专栏的表面原因,但它可以被证明是合理的,因为J的直接前身之一是A语言,这是阿瑟·惠特尼(Arthur Whitney)在摩根士丹利(Morgan Stanley)投资银行工作期间开发的。这是惠特尼一天下午在安大略省肯·艾弗森的窑场(Kiln Farm)编写的一页A类口译器的C代码,为罗杰·许(Roger Huge)提供了开始工作所需的方向,J·罗杰在艾弗森软件公司(Iverson Software Inc)出版的J&34;实现附录A中展示了这段代码。惠特尼不再在摩根士丹利工作:他开始是一名自由撰稿人,正在开发一种名为K的语言,这种语言与A和J有一定的亲和力。许仕仁现在摩根士丹利工作,故事讲述的是他冒险寻找一个大质数的故事,而许志永将A作为他瞄准大质数的武器。

故事开始于惠氏在摩根士丹利的老板向他发起挑战,他说:“你认为你很聪明,但你甚至不知道10到10的质数是多少。”惠氏的第一反应是,你是从0开始数还是从1开始数?这位老板大吃一惊,有一分钟左右他都听不懂罗杰的问题。

老板的挑战似乎是一个理论上可以实现但实际上不可能完成的计算任务的例子。这篇文章讲述了罗杰是如何去实现不可能的事情的。我认为它是客户端-服务器技术的胜利。本专栏与其说是一篇关于编程的文章,不如说是一篇关于计算机物流的文章。编程方面虽然很重要,但与Roger如何组织解决方案的故事相比,编程方面是次要的。

惠的解决方案的萌芽是设想一个长度为k的布尔向量p,如果i是素数,那么p的第i个元素是1,否则是0。只需对这个非常长的向量进行求和扫描,并在其中查找1e10的索引。这样的载体应该有多长呢?素数定理说,小于1k的素数的数目大约是1k%(Logk)。在方程式1e10=k%(Log K)中求解k,得到的值约为22.63e11。罗杰出于谨慎,使用了2.7e11的值。*在计算机内存的当前状态下,由2.7E11个元素组成的向量是无法实现的,特别是因为A没有布尔类型:布尔向量需要与整数向量相同的空间。4x2.7e11或111e12字节长的向量根本不在卡片中。即使每个元素只有1位的布尔向量也必须超过3E10个字节的长度,因此显然必须对问题进行划分。

惠的中央程序计算m和n之间的素数,使用筛子方法,消除2,3,5,7,11,13,17等的倍数。这可以在许多小间隔上独立地并行完成,这些小间隔构成了更大的感兴趣的间隔,如果分配到可以与公共中心文件通信的计算机上,将允许在比只有一台计算机处理问题的情况下在更短的时间内解决问题。

分区的方法是由摩根士丹利惠氏部分的150多个工作站提出的。它们都是互连的Unix机器,任何机器都可以在任何其他机器上使用。只需相对较少的工作,就可以建立一个多处理器解决方案,使用这些机器并行运行。这些机器在周末不是很常用,实验是在一个周末进行的。

策略是让机器一次处理10亿个整数。惠氏在摩根士丹利的同事塞思·布莱巴特(Seth Breitbart)建议创建271个文件,命名为0,1,2,.,270,每个文件表示命名间隔为1e9个连续数字,每个文件都是空的。

一旦设置好,机器将做的是查看文件列表,随机挑选一个(命名为1m),擦除它,在间隔t(1e9*m)+i.1e9上工作,一次处理一百万个数字,完成后,编写一个文件,其中包含一个记录,给出该1e9内每个百万号间隔中的质数(有一千个)。完成后,它重复该过程,仅在文件/间隔列表为空时停止。这些机器被设置为一次处理一百万个数字,因为可用的最小的机器有足够的内存一次处理这么多的数字。罗杰指出,如果两台机器意外地选择了相同的文件/间隔,不会造成很大的伤害。在灵活的Unix世界中,任何时候都可以将更多的机器投入使用。

如果一个人是Unix系统的超级用户,你可以随意使用这些机器,比如找到所有其他机器的ID,但许志永(我认为这是明智的)不想被这样的功能玷污,所以为了得到机器的名字,他在地板上四处阅读贴在每台机器上的纸条上的机器名称,然后坐在他的机器前,询问他名单上每台机器的状态。如果它是空闲的,他就让它去解决问题。

当他这样做的时候,有一个很好的两难境地要解决。在启动更多机器之前,他的时间应该花在改进算法上,还是应该花时间寻找更多的机器?出于新颖性的考虑,他倾向于后一种方法,最终使用了大约15个IBMRS/6000和60个Sun Sparc 2和Sparc 10。

20个小时后,他有271个文件,每个文件有1000条记录。从这些数据中,他得出了区间1e6*i.271000中素数的271000个元素向量。因此,通过对此进行总和扫描,他知道了包含第10^10次素数的间隔。他的函数psieve返回一个布尔列表,选择m和n之间的素数。将其应用于魔术区间得到实际的第10^10个素数。

关于他的程序的一些细节:如果一个数q不能被小于或等于1sqrt(Q)的任何数整除,则q是素数。因此,要测试一个小于2.7e11的数字的素性,只需使用小于1sqrt(2.7e11)或大致小于6e5的试除数即可。在实践中,惠递归地预计算出了所有小于1e6的78,498个素数的列表,从2 3 5 7开始自举。(这只需要几秒钟。)。当时,确定任何小于1e12的数字是否为素数是例行公事:只要看看它相对于每个小于一百万的素数的残数是否是非零;如果是,那它就是素数。

对于好奇的人来说,这里是前10,000,000,002个素数的简明列表,以及它们的0起源序号。

0 2 1 3 2 5 3 7 4 11.1e10-1252,097,800,623 1e10 252,097,800,629 1e10+1 252,097,800,637。

这样做之后,许发现威廉·贾德森·列维克(William Judson Leveque)在“数论基础”第1.1节中有一张表,上面列出了少于10^3+I.8的素数数量。惠氏的表格与Leveque的表格10^3+I.7一致。然而,对于10^10,Leveque说是455,052,512,许说是455,052,511。事实证明莱维克是错的,惠曾在滑铁卢大学的李·迪基的帮助下检查了他的结果。Dickey告诉Hue,他的同事们推测Leveque可能是从D.N.Lehmer编制的列表中获得了他的数字,其中包括1作为质数,而Leveque可能没有从那个特定的数字中减去1。(1不是素数,因为它不满足素数的定义:正好有两个不同因子1和n的正整数n。)。

现在来看看这个故事的寓意:许志永告诉我,从那以后,他也找到了一些工作,可以更容易地发现任何第n个素数。德国天文学家梅塞尔在19世纪70年代发现了一种计算pi(X)的个体值的方法,即素数t<;:x的计数函数。但他的方法是基于部分筛选函数的递归,他用它来计算pi(1e9),其中,pi是一个计算小于或等于其自变量的素数的函数。D.H.Lehmer简化和扩展了Meissel的方法。最近,包括一些新的筛选技术的Meissel-Lehmer方法的进一步改进已经由Lagaria,Miller和Odlyzko在";Computing pi(X):Meissel-Lehmer method";中报告,Maxims of Computation,卷44,编号170,1985年4,第537到560页。在本文中,作者给出了所得算法的渐近运行时间分析,表明对于每个e>;0,它最多使用O(x^(2r3)+e)算术运算,并且最多使用计算机上的O(x^(1r3)+e)个存储位置,使用长度为1+<;2^.x比特的字。使用并行处理器可以进一步加速算法。它们表明,有一种算法,当给定M个并行处理机时,在每个并行处理机上最多使用O(x^(1r3)+e)个存储位置来计算时间至多为O((%M)*x^(2r3)+e)的pi(X),提供了M<;:x^1r3。*实现了该算法的一个变体,并将其用于计算pi(4e16)。

他们报告称,在IBM 3081K上,pi(4e16)花了他们1730分钟;npi(2e12)花了3分钟;npi(3e12)花了4分钟。如果他知道这种方法,在相当慷慨地使用素数定理缩小搜索区间后,辉可能会使用二进制搜索技术,使用OO(2logn)技术来找到第110个素数。惠计算的值是1pi(2.52e11),他在60-70个工作站上花了20个小时。大脑再次战胜肌肉:一个设计良好、数学知识渊博的算法击败了蛮力!