first 10-digit prime found in consecutive digits of e

一道谷歌面试题：

超越数e的小数点后，连续的10位数字，有些能构成一个质数。找到最先出现的这样一个10位质数。

参考了以下链接：
https://blog.csdn.net/huaxi1902/article/details/39155515
https://blog.csdn.net/witnessai1/article/details/68547820
https://zhidao.baidu.com/question/1540327644955312227.html
http://numbers.computation.free.fr/Constants/TinyPrograms/tinycodes.html#tth_sEc2
http://blog.chinaunix.net/uid-21712186-id-1818141.html

一个想法：
假如全部的10位自然数里有4.5%的质数，假设对于e的小数点后任意连续10位，都是一个统计意义上的随机数，那么前100位里出现质数的概率是1-(1-4.5%)^100，非常接近于1了，也就是说第一个质数出现位置的期望其实挺早的。

思路：原问题可以分成两个问题分别求解。

1.高效生成超越数e的n位有效数字
2.高效判断一个10位数是不是质数。

对于问题1，没找到太imba的方法，本质上都是泰勒展开的各种变形。
$e^x=1+x+\frac{x^2}{2!}+\frac{x^3}{3!}+...$
令x=1代入上式，可以知道求e的小数等于求自然数阶乘的倒数和。为了保证精度，直观的一个思路是每计算一位，整体乘10。另外泰勒展开式也可以化成一些其他利于计算的形式。详见链接三。
代码比较短的是第四个链接里写的方法：

main(){
int N=9009,n=N,a[9009],x;
while(--n)a[n]=1+1/n;
for(;N>9;printf("%d",x))
for(n=N--;--n;a[n]=x%n,x=10*a[n-1]+x/n);
}

对于问题2，因为10万乘10万就是11位数了，所以对于一个10位数，只要验证其没有10万以内的约数，就能证明它是质数。比较直接的思路是生成一个10万以内的质数表，然后拿那些10位数挨个除。

接下来介绍一些真正有趣的东西，用来验证一个超大的数是不是质数。
http://blog.codinglabs.org/articles/prime-test.html
这个链接说得比较明白，我复述得很乱。

先介绍一个算法：蒙哥马利快速幂模
$(X^Y) \%Z = [(X* X)\% Z] * [(X* X)\%Z] * ......$
（代码来自链接2）

unsigned Montgomery(unsigned n,unsigned p,unsigned m)  
{ //快速计算(n^p)%m的值  
      unsignedk=1;  
      n%=m;  
     while(p!=1)  
     {  
         if(0!=(p&1))k=(k*n)%m;  
         n=(n*n)%m;  
         p>>=1;  
    }  
    return(n*k)%m;  
}

类似原理，利用二进制特性的快速幂模法：重复平方法。因为对一个数取平方，等于把它化为2的幂指数形式之后，在指数部分结尾加一个0；取平方再乘以自身，相当于再幂指数的结尾加一个1。每次运算时取模，即可得到幂模。

def compute_power(a, p, m):
    result = 1
    p_bin = bin(p)[2:]
    length = len(p_bin)
    for i in range(0, length):
        result = result**2 % m
        if p_bin[i] == '1':
            result = result * a % m
 
    return result

算法复杂度是对数级的，很快。后边的检测法得以实现的前提，就是计算机可以快速用这种算法取模。

费马小定理：
任意正整数N，素数P，N不能被P整除（两数互质），则有：
$N^P\%P=N$
如果N也是素数，则有：
$N^{(P-1)}\%P=1$
现在我们要判断N是不是素数，只要找其他的素数，看看费马测试的结果是不是1即可。利用蒙哥马利快速幂模可以快速计算这个式子的值。比如我们随便找50个质数P，都通过了，那么N大概率是一个质数。要注意，质数一定能通过费马测试，但通过了费马测试的数不一定是质数。有一类合数可以通过所有P值的费马测试，这种合数被称为Carmichael数。卡尔麦克数很稀少，但也有无穷多个。

老师，能不能更给力一些呢？
由费马小定理扩展的更有效的测试法：拉宾米勒检验。
这一算法用于“工业级”检验素数，优点是可控、简单粗暴。每多迭代一次，失误几率就乘以1/4，比如迭代8次那么失误的几率就小于10的负五次方了。控制使用的随机数的个数，即可控制迭代次数。通过了这一测试算法的数字，即使不是质数，在某些情形下也可以当作质数来用。现在展开讲讲这个拉宾米勒测试。

引理：
素数，x是小于p的正整数，且x^2 = 1 mod p，则x要么为1，要么为p-1。

说明对于任意一个小于p的正整数x，如果发现1（模p）的非平凡平方根存在，则说明p是合数。

具体算法：
对于待检测的一个很大的奇数p，有偶数p-1，并且必然有一组u和t，其中u是奇数，t是正整数，使得 $p-1=u 2^t$ ，于是有 $a^{p-1}=a^{u 2^t}=(a^u)^{2^t}$
就是说先算出a^u，然后平方t次，在任意一次平方时发现了非平凡平方根，可以判定p是合数。
u和t的值在二进制下非常容易找到。

举个例子：
检测561是不是质数：
560二进制为1000110000，去掉末尾的4个0，剩下100011，十进制为35，所以560=35 * 2^4，u=35，t=4，那么：
2^35 mod 561 = 263
263^2 mod 561 = 166
166^2 mod 561 = 67
67^2 mod 561 = 1
找到了非平凡平方根67，所以561是合数。
如果最后一步不等于1，那么561很大概率就是质数了。

计算机实现：对于待测大数p：

1.选择随机数 $A<p$
2.找到R和M，使得 $p=(2^R)*M+1$ 成立。
快速得到R和M的方式：p用二进制数B来表示，令C=B-1。因为p为奇数（素数都是奇数），所以C的最低位为0，从C的最低位的0开始向高位统计，一直到遇到第一个1。这时0的个数即为R，M为B右移R位的值。
3.如果 $A^M\%N=1$ ，则通过A对于p的测试，然后进行下一个A的测试
4.如果 $A^M\%N\neq 1$ ，那么令i由0迭代至R，进行下面的测试
5.如果 $A^ ((2^i)*M)\%N=N-1$ 则通过A对于N的测试，否则进行下一个i的测试
6.如果i=r，且尚未通过测试，则此A对于N的测试失败，说明p为合数。
7.进行下一个A对p的测试，直到测试完指定个数的A
代码来自第二个链接：

bool RabbinMillerTest( unsigned n )  
{  
    if (n<2)  
    {   
         // 小于2的数即不是合数也不是素数  
        throw 0;  
    }  
    const unsigned nPrimeListSize=sizeof(g_aPrimeList)/sizeof(unsigned);//求素数表元素个数  
    for(int i=0;i<nPrimeListSize;++i)  
    {  
        // 按照素数表中的数对当前素数进行判断  
        if (n/2+1<=g_aPrimeList[i])  
        {  
            // 如果已经小于当前素数表的数，则一定是素数  
            return true;  
        }  
        if (0==n%g_aPrimeList[i])  
        {  
            // 余数为0则说明一定不是素数  
            return false;  
        }  
    }  
    // 找到r和m，使得n = 2^r * m + 1;  
    int r = 0, m = n - 1; // ( n - 1 ) 一定是合数  
    while ( 0 == ( m & 1 ) )  
    {  
        m >>= 1; // 右移一位  
        r++; // 统计右移的次数  
    }  
    const unsigned nTestCnt = 8; // 表示进行测试的次数  
    for ( unsigned i = 0; i < nTestCnt; ++i )  
    {   
        // 利用随机数进行测试，  
        int a = g_aPrimeList[ rand() % nPrimeListSize ];  
        if ( 1 != Montgomery( a, m, n ) )  
        {  
            int j = 0;  
            int e = 1;  
            for ( ; j < r; ++j )  
            {  
                if ( n - 1 == Montgomery( a, m * e, n ) )   
                {  
                    break;  
                }  
                e <<= 1;  
            }  
            if (j == r)  
            {  
                return false;  
            }  
        }  
    }  
    return true;  
}

补充：
1费马小定理由费马于1636-1640年发现并提出，不是费马证明的（不知道是不是因为信纸上的地方不够），莱布尼茨1683年最早证明了它，但手稿没发表。最早发表证明的是欧拉，1736年。
2关于费马小定理的一个浅显而奇妙的证明：Proof by counting bracelets，by Golomb
已知：p是素数，a是整数且1≤ a ≤ p
求证：a^p − a 能被p整除
证明：假设一串由p个珠子组成的项链，项链中珠子的种数为最多为a，那么所有可能的不重复的排列为a^p ，其中种数为1的情况有a种，那么去除这些情况后一共就有a^p− a种。现在把这a^p− a种项链首尾相接，那么就有可能出现重复的情况了。因为p是素数，所以一个包含着p个珠子的环旋转后会出现p种不同的情况，亦即，这a^p− a种情况可以分为一些类，每类有p中情况，亦即，a^p − a 能被p整除。