偶然间,看到吴军老师的《浪潮之巅(第四版)》里有讲到这么一个故事。

“很聪明"?吴军老师的这句话倒是让人来兴趣了,自己也凑性借助计算机的力量琢磨琢磨下这个证明题。


思路

大家都知道,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

在网上查阅素数相关的资料的时候,发现数论里有个素数分布规律也可以拿来判断素数。来源——素数判定算法

用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

参考资料汇总

  1. 吴军老师的《浪潮之巅(第四版)》P44
  2. 无理数e的维基英文百科
  3. OEIS站点的A001113页面
  4. Table of 50000 digits of e labeled from 1 to 50000
  5. 素数判定算法
  6. PyCrypto密码学库源码解析(一)随机数和大素数生成
  7. What Google’s Genius Billboard From 2004 Can Teach Us About Solving Problems
03-05 15:19