In [1]:
import numpy as np
import scipy.stats as spstat
import pandas as pd
import bokeh
import bokeh.plotting as bplt
import bokeh.palettes as bpal
# 工具配置
TOOLS = "pan,wheel_zoom,box_zoom,box_select"
# 作图配置
plot_config = dict(plot_width=800, plot_height=500, tools=TOOLS)
bplt.output_notebook()
In [3]:
from IPython.display import Image
Image(url="https://upload.wikimedia.org/wikipedia/commons/b/b9/Sieve_of_Eratosthenes_animation.gif")
Out[3]:
In [2]:
def Sieve(N):
res = np.int64(np.ones(N))
# res的第i个元素为1时，表示其为质数，为0时，表示其不为质数
res[0] = 0
res[1] = 0
# 仅需要检查到sqrt(N)，因为若N非质数，则其必有一个大于1小于等于sqrt(N)的因子
upbound = np.int64(np.sqrt(N))
for i in range(2,upbound):
if res[i]!=0:
# i此时未被其因数标记为0，但我们要把从i^2开始，所有i的倍数都标记为0
res[i**2::i] = 0
else:
continue
return res
In [3]:
for idx,isPrime in enumerate(Sieve(20)):
if isPrime==1:
print idx,
In [4]:
%timeit Sieve(5000)
%timeit Sieve(500000)
%timeit Sieve(50000000)
我们来计算一下筛法的复杂度：若要求出N以内的质数，我们对每一个质数$p$都进行了最多$N/p$次调整
因此$$\sum_{\text{p is prime},p<N} \frac{N}{p} < \sum_{i=2}^N \frac{N}{p} = O(N log N)$$
In [8]:
fiftymilprime = pd.Series(Sieve(50000000).cumsum())[1:]
primeconverge = fiftymilprime*np.log(fiftymilprime.index)/fiftymilprime.index
indexer = [2**i for i in range(1,25)]
p1 = bplt.figure(title="质数数量与N/ln(N)之比",
background_fill="#FFFFFF",x_axis_type="log",**plot_config)
p1.line(primeconverge.index[indexer], primeconverge.values[indexer], line_width=3, color='#3288bd', legend='序列收敛到1')
p1.grid.grid_line_alpha=0.3
p1.xaxis.axis_label = 'N'
p1.yaxis.axis_label = '比值'
p1.ygrid.band_fill_color="steelblue"
p1.ygrid.band_fill_alpha = 0.1
bplt.show(p1)
这样当我们再想使用特定数目个质数的时候，我们就知道该在多大的N内寻找这些质数了。
在知道了质数的数量、密度之后，数学家们进一步寻找质数分布的数学规律就是理所应当的了。
（笔者尚未发现孪生素数的工程应用，请读者指正。）
In [11]:
print 91%13 == 0
print 444%38 == 0
当然，用程序来计算两个数是否互质有更好的办法：辗转相除法，也称欧几里得算法：
首先用$a$除以$b$得到余数$a'$，$a,b$ 互质等价于 $a',b$ 互质。
这是因为，
$$a = br + a'\\ brx+a'x+by = 1 \\ a'x + b(rx+y) = 1$$现在我们很确信的知道$a'<b$，那么再用$b$除以$a'$得到余数$b'$，这样循环往复的操作。
如果这两个数反复相除取余数，在某一步其中一个变成了0，那么另外一个非0的数就是$a,b$的最大公约数。
如果这个数恰好等于1的话，那么$a,b$就是互质的。因为最大的公共约数就是1了，也不可能再有别的公共约数了。
In [84]:
def Euclidean(i1,i2):
a = i1
b = i2
while (a*b!=0):
a = a%b
## 交换一下a和b，确保交换后a比b大。
a,b = b,a
return a+b
In [85]:
print Euclidean(4,9)
print Euclidean(7,14)
In [98]:
def ExtendEuclidean(i1,i2):
a = i1
b = i2
s_2 = 1
t_2 = 0
s_1 = 0
t_1 = 1
while (a*b!=0):
q = a/b
a = a%b
a,b = b,a
s_2 = s_2-q*s_1
t_2 = t_2-q*t_1
s_2,s_1 = s_1,s_2
t_2,t_1 = t_1,t_2
return a+b, s_2,t_2
print ExtendEuclidean(4,9)
第一个1表示$G.C.D(4,9)=1$，而后面的$-2,1$表示$x,y$。
刚才提到，欧拉算法能够得到两个数的最大公约数（Greatest Common Divisor, GCD）。
我们把思路变换一下：如果我突然想求两个数的最小公倍数（Least Common Multiple, LCM），应该如何求呢？
我们先用个粗糙的办法：
因为$ab$是$a,b$的公倍数，但可能不是最小的一个公倍数。遍历$i=1,2,...b$，来测试$ai$是否是$b$的倍数即可：
In [25]:
def LCMBrute(i1,i2):
a = i1
b = i2
for i in range(1,i2+1):
if (a*i)%b==0:
return a*i
print LCMBrute(4,9)
print LCMBrute(7,14)
很显然，当$a,b$很大的时候，上述方案很不好。
我们再次不加证明的给出一个定理，使用本次介绍的内容完全可以证明这一公式：
$$(\forall a,b \in I), ab = G.C.D(a,b)* L.C.M(a,b)$$辗转相除法可以相当快的给出最大公约数，从而可以用如上的公式计算出最小公倍数。
根据这个结果，我们可以轻松的知道，$a,b$互质时，$L.C.M(a,b)=ab$
先简单介绍一个概念：同余
$$p | x_1 - x_2 则 x_1 \equiv x_2 (\text{mod p})$$同余对于加法、乘法而言都是闭的：
$$x_1 \equiv a (\text{mod p}),x_2 \equiv b (\text{mod p}), \text{ then } $$$$x_1+x_2 \equiv a+b (\text{mod p}) , x_1x_2 \equiv ab (\text{mod p})$$结合上面的知识我们进一步了解中国古代一个好玩的数学定理。国外常称之为『中国剩余定理』。我们叫孙子定理就可以了。
曾有一个问题：今有物不知其数,三三数之剩二,五五数之剩三,七七数之剩二,问物几何?
而这个问题的解答可以用一首诗概括：
『三人同行七十稀/五树梅花廿一枝/七子团圆正半月/除百零五便得知』
这就暗含着对此类同余问题的一种解法：
$$ x \equiv 2 (\text{mod 3}) \\ x \equiv 3 (\text{mod 5}) \\ x \equiv 2 (\text{mod 7}) $$那么 $$x \equiv 2*70+3*21+2*15 \equiv 23 (\text{mod 105})$$
因为对于加法、乘法而言，同余性质都能够良好的保持住，仔细观察还可以发现：
这样我们就有： $$ x_1 \equiv 1 (\text{mod 3}) \\ x_1 \equiv 0 (\text{mod 5}) \\ x_1 \equiv 0 (\text{mod 7}), $$
$$ x_2 \equiv 0 (\text{mod 3}) \\ x_2 \equiv 1 (\text{mod 5}) \\ x_2 \equiv 0 (\text{mod 7}), $$$$ x_3 \equiv 0 (\text{mod 3}) \\ x_3 \equiv 0 (\text{mod 5}) \\ x_3 \equiv 1 (\text{mod 7}). $$的一组解：$x_1=70,x_2=21,x_3=15$，这就是诗句中提示的几个数字。
最后，在这几个数字的辅助下，求得：
$$x \equiv 2*x_1+3*x_2+2*x_3 (\text{mod 105})$$我们可以试验一下，105分解成为357之后，$3,5,7$是两两互质的。我们从0遍历到104：
In [27]:
for i in range(105):
print [i,(i%3,i%5,i%7)]
我们先随便试验一下：
In [32]:
print (2**7-2)%7 == 0
print (34**5-34)%5 == 0
print (21**3-21)%3 == 0
换个写法：$$a (a^{p-1} - 1) = a(a-1)(\sum_{i=0}^{p-2} a^i) = (a-1)(\sum_{i=1}^{p-1} a^i) \equiv 0 \text{ (mod p) }$$
当 $a \equiv 0 \text{ or } 1 \text{ ( mod p ) }$是结果是显然的。
不然，$a,2a,3a,...,(p-1)a$ 除以$p$的余数肯定是$1,2,...,p-1$的重排列（可以证明每一个都不可能同余0，而每一对之间都不可能同余）
通过乘法可得 $$a^{p-1} (p-1)! \equiv (p-1)! \text{ (mod p) } $$
此时，$(p-1)!,p$ 二者互质，因此 $$a^{p-1} \equiv 1 \text{ (mod p) } $$
RSA加密原理由以下关键点组成：
假设$N$是两个质数$p,q$的乘积，同时能够找到$e_1$与$(p-1)(q-1)$互质，也就是能找到
$$e_1 e_2 \equiv 1 \text{ ( mod (p-1)(q-1) )}$$此时，在仅仅知道$N,e_1$的情况下，我们可以把一个数字A进行一定的转换：
$$B = A^{e_1} \% N$$只要$A<N$，仅使用$N,e_2$的信息就可以把A给恢复出来（否则，恢复出来的$A'$将是$A$除以$N$的余数）：
$$A = B^{e_2} \% N$$如果你还记得刚才的费马小定理的话：
$$A^{p-1} \equiv 1 \text{ (mod p) } $$$$A^{q-1} \equiv 1 \text{ (mod q) } $$将上述的两个式子分别乘以$q-1,p-1$次方，套用孙子定理：
$$A^{(p-1)(q-1)} \equiv 1 \text{ ( mod N ) }$$于是
$$A^{e_1e_2} = A^{1+k(p-1)(q-1)} \equiv A^1 = A \text{ ( mod N ) }$$我们来实验一下...受限于对大整数的支持，我们只能用小整数玩玩了...
In [151]:
tenmilprimes = pd.Series(Sieve(1000))
primeonly = tenmilprimes[tenmilprimes==1].index
print primeonly[-1],primeonly[-2],primeonly[-3]
In [152]:
p = primeonly[-1]
q = primeonly[-2]
N = p*q
M = (p-1)*(q-1)
e1 = primeonly[-3]
res = ExtendEuclidean(e1,M)
e2 = res[1]
y = res[2]
print e1,e2,(e1*e2)%M
In [154]:
def expmod(a,e,N):
eres = e
modres = 1
curmod = a
while eres>0:
r = eres%2
eres = eres/2
modres = (modres*(curmod**r))%N
curmod = (curmod**2)%N
return modres
A = 101
B = expmod(A,e1,N)
A = expmod(B,e2,N)
print A,B
我们把$(N,e_1)$这组信息成为公钥，也就是大家都能用的钥匙。把$A$称为明文，然后通过公钥加密$A$得到的$B$我们称为密文。
把$A$加密成为$B$很容易，但是除非我们知道$N,e_2$的信息，否则从$B$恢复成为$A$是不容易的。
为了知道$e_2$是什么，黑客们就希望从公开信息$N$中分解出$p$和$q$。
而这项工作在$N$足够大(比如$N$来自于两个超过100位大质数相乘)时就变得很困难。
因此，RSA加密算法的安全性来源于大数的质因子分解十分困难这一现象。
本节内容欲直接了解细节，请参考：
在质数中有一类形式很简洁的质数叫做梅森质数(Mersenne Prime)，它有这种性质：
$$p = 2^k - 1 , k\in I$$人们对它苦苦寻找，因为利用$k$他能够通过某种迭代的方式不断的产生一定长度的数字，而且拥有很多很好的性质：
这种性质很适合产生随机数！
举例说明，C语言最常用的用于产生32bit整数的随机数生成器MT19937参数如下：
而原本是二进制矩阵与向量的乘法操作、矩阵的加法操作或者向量的截取操作，人们可以通过使用计算机最原始的指令集来完成！(通过异或操作、与操作、或操作、向左向右位移操作完成)。
步骤如下：
不加修改的给出代码Wikipedia上的代码如下：
In [156]:
def _int32(x):
# Get the 32 least significant bits.
return int(0xFFFFFFFF & x)
class MT19937:
def __init__(self, seed):
# Initialize the index to 0
self.index = 624
self.mt = [0] * 624
self.mt[0] = seed # Initialize the initial state to the seed
for i in range(1, 624):
self.mt[i] = _int32(
1812433253 * (self.mt[i - 1] ^ self.mt[i - 1] >> 30) + i)
def extract_number(self):
if self.index >= 624:
self.twist()
y = self.mt[self.index]
# Right shift by 11 bits
y = y ^ y >> 11
# Shift y left by 7 and take the bitwise and of 2636928640
y = y ^ y << 7 & 2636928640
# Shift y left by 15 and take the bitwise and of y and 4022730752
y = y ^ y << 15 & 4022730752
# Right shift by 18 bits
y = y ^ y >> 18
self.index = self.index + 1
return _int32(y)
def twist(self):
for i in range(624):
# Get the most significant bit and add it to the less significant
# bits of the next number
y = _int32((self.mt[i] & 0x80000000) +
(self.mt[(i + 1) % 624] & 0x7fffffff))
self.mt[i] = self.mt[(i + 397) % 624] ^ y >> 1
if y % 2 != 0:
self.mt[i] = self.mt[i] ^ 0x9908b0df
self.index = 0
In [159]:
myMT = MT19937(1)
for i in xrange(20):
print myMT.extract_number()
使用这样的随机数生成器，你就可以很均匀的生成从0到$2^{32}-1$的数字了。如果拿这些数直接除以$2^{32}$，那我们就得到了相当好的从0~1的 均匀分布。关于均匀分布为什么很重要，在下一章中你就知道啦！
同时我们把这种被随机数种子唯一决定的随机数序列称为伪随机数，它的特点如下：
所以它更适合用于后者的场合~
我们回到主题上来：随机数生成器帮助我们做了无穷无尽的实验，都要归功于$2^{19937}-1$是一个梅森质数，所以人类一直在追寻越来越大的梅森质数，反正闲着也是闲着！
原本用来全球范围内分布式搜索梅森质数的程序GIMPS开发的PRIME95，在十几年间其实是变成了一个测试PC Server稳定与否的程序，这些CPU超频的极限玩家很可能压根不关心什么数学问题，但有意思的是，最近数学家们验证的梅森质数无一例外是借助PRIME95程序找到的。2016年1月，人类发现并验证了截止笔者写作时最大的梅森质数：
$$2^{74207281} - 1$$这比$2^{19937}-1$可不知道高到哪里去了！
今天，你对奇妙的质数又有了新的认识吗？
In [ ]: