偶然间,看到吴军老师的《浪潮之巅(第四版)》里有讲到这么一个故事。
Google曾经在加州的101高速公路上用大广告牌登了这样的广告:
{ 无理数e中前十位连续的素数 }.com
你如果知道这个答案(http://7427466391.com),就可以通过上述网址进入到Google的招聘网站。而能够计算出这道题,要很聪明。
“很聪明"?吴军老师的这句话倒是让人来兴趣了,自己也凑性借助计算机的力量琢磨琢磨下这个证明题。
大家都知道,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密码学库源码解析(一)随机数和大素数生成, 以及上面那篇。
至此,必要的材料都准备好了,可以进行最后一步了。
开门见山,直接源码抛出来先。
具体思路:先使用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))
运行结果:
宾果!
证明题解答完毕,咱得走个仪式感,访问访问这个网站——http://7427466391.com。
结果是502错误……
行吧,看来这个站点早被人家抛弃了,毕竟这个高速公路广告也是谷歌在2004年搞的恶作剧。
最后,把源码整理了个Kaggle Notebook版本。欢迎查阅!
First10DigitPrimeFoundInConsecutiveDigitsOfE | Kaggle