MillerRabin 素性检测

$ 2^{256} - 2^{32} - 2^9 - 2^8 - 2^7 - 2^6 - 2^4 -1=$

$FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFEFFFFFC2F$

115792089237316195423570985008687907853269984665640564039457584007908834671663 是个质数吗？

$2^{255}-19$
57896044618658097711785492504343953926634992332820282019728792003956564819949 是个质数吗？

### 测试

测试100是否为质数
只需要测试 `[2, sqrt(100) ]` 区间的数能被100整除否

100/10 = 10
100/5 = 20
100/20 = 5
这算是某种 “对称性”吧

更快速的方法是：求得[2, sqrt(100)]区间的所有质数，能否被100整除否，既然取质数，偶数除了2以外，都不是质数

`rane(start, end, step)`， 所有step可以为2

记得有一种筛选法寻找质数，
剔除2的倍数，3的倍数，5的倍数，7的倍数 ...... 依次类推

### ruby 有特殊的技巧

```ruby
require 'prime'

Pcurve = 2**256 - 2**32 - 2**9 - 2**8 - 2**7 - 2**6 - 2**4 -1

Prime.each(10).to_a  #10以内的质数表
# [2, 3, 5, 7]

Prime.each(10).to_a.size #10以内的质数数量
# 4

Prime.each(100000000).to_a.size() #100000000以内的质数数量
# 5761455

Prime.prime?(7) # 7是质数么
# true

Prime.prime?(Math.sqrt(Pcurve).to_i) 
# 340282366920938463463374607431768211456 是质数么
# false, 速度很快

Prime.prime?(Pcurve)
# 10分钟后还在计算
```

**Pcurve 难以想象的大**

### python code

按照以上想法做个实现

```python
#coding:utf-8

import math

Pcurve = 2**256 - 2**32 - 2**9 - 2**8 - 2**7 - 2**6 - 2**4 -1
 
def _sqrt(x):
    r = int(math.sqrt(x))
    return r

def is_prime():
    for i in range(2, _sqrt(Pcurve), 2):
        if (Pcurve % i) == 0:
            print(Pcurve," is not a prime number")
            print(i," times",Pcurve//i,"is",Pcurve)
            break
        else:
            print(Pcurve,"is a prime number")


if __name__ == '__main__':
    is_prime()
```

代码不能运行，提示 `OverflowError: range() result has too many items`，直接溢出了。

这就是感觉上很简单，验证却比较困难的问题。

**Pcurve大的超出了想象**

### 下一步的方案

我们可以找到  $\sqrt{Pcurve} $ 内的质数，
欲找到  $\sqrt{Pcurve} $内的质数，先找到$ \sqrt{ \sqrt{Pcurve} }$内的质数， 
 欲找到$ \sqrt{ \sqrt{Pcurve} }$内的质数，先找到$ \sqrt{\sqrt{ \sqrt{Pcurve}}}$内的质数， 

如果想更快，就需要了解数论方面的研究了

😶😶😶😶😶😶

### MillerRabin 素性检测

不满足 $2^{(n-1)} \mod n = 1$ 的n一定不是素数；如果满足的话则多半是素数。这样，一个比试除法效率更高的素性判断方法出现了：制作一张伪素数表，记录某个范围内的所有伪素数，那么所有满足2^(n-1) mod n = 1且不在伪素数表中的n就是素数。之所以这种方法更快，是因为我们可以使用二分法快速计算2^(n-1) mod n 的值，这在计算机的帮助下变得非常容易；在计算机中也可以用二分查找有序数列、Hash表开散列、构建Trie树等方法使得查找伪素数表效率更高。
    有人自然会关心这样一个问题：伪素数的个数到底有多少？换句话说，如果我只计算2^(n-1) mod n的值，事先不准备伪素数表，那么素性判断出错的概率有多少？研究这个问题是很有价值的，毕竟我们是OIer，不可能背一个长度上千的常量数组带上考场。统计表明，在前10亿个自然数中共有50847534个素数，而满足2^(n-1) mod n = 1的合数n有5597个。这样算下来，算法出错的可能性约为0.00011。这个概率太高了，如果想免去建立伪素数表的工作，我们需要改进素性判断的算法。

最简单的想法就是，我们刚才只考虑了a=2的情况。对于式子 a^(n-1) mod n ，取不同的a可能导致不同的结果。一个合数可能在a=2时通过了测试，但a=3时的计算结果却排除了素数的可能。于是，人们扩展了伪素数的定义，称满足a^(n-1) mod n = 1的合数n叫做以a为底的伪素数(pseudoprime to base a)。前10亿个自然数中同时以2和3为底的伪素数只有1272个，这个数目不到刚才的1/4。这告诉我们如果同时验证a=2和a=3两种情况，算法出错的概率降到了0.000025。容易想到，选择用来测试的a越多，算法越准确。通常我们的做法是，随机选择若干个小于待测数的正整数作为底数a进行若干次测试，只要有一次没有通过测试就立即把
这个数扔回合数的世界。这就是Fermat素性测试。
人们自然会想，如果考虑了所有小于n的底数a，出错的概率是否就可以降到0呢？没想到的是，居然就有这样的合数，它可以通过所有a的测试（这个说法不准确，详见我在地核楼层的回复）。Carmichael第一个发现这样极端的伪素数，他把它们称作Carmichael数。你一定会以为这样的数一定很大。错。第一个Carmichael数小得惊人，仅仅是一个三位数，561。前10亿个自然数中Carmichael数也有600个之多。Carmichael数的存在说明，我们还需要继续加强素性判断的算法。

Miller和Rabin两个人的工作让Fermat素性测试迈出了革命性的一步，建立了传说中的Miller-Rabin素性测试算法。新的测试基于下面的定理：如果p是素数，x是小于p的正整数，且x^2 mod p = 1，那么要么x=1，要么x=p-1。这是显然的，因为x^2 mod p = 1相当于p能整除x^2-1，也即p能整除(x+1)(x-1)。由于p是素数，那么只可能是x-1能被p整除(此时x=1)或x+1能被p整除(此时x=p-1)。
我们下面来演示一下上面的定理如何应用在Fermat素性测试上。前面说过341可以通过以2为底的Fermat测试，因为2^340 mod 341=1。如果341真是素数的话，那么2^170 mod 341只可能是1或340；当算得2^170 mod 341确实等于1时，我们可以继续查看2^85除以341的结果。我们发现，2^85 mod 341=32，这一结果摘掉了341头上的素数皇冠，面具后面真实的嘴脸显现了出来，想假扮素数和我的素MM交往的企图暴露了出来。

**这就是Miller-Rabin素性测试的方法**。不断地提取指数n-1中的因子2，把n-1表示成d*2^r（其中d是一个奇数）。那么我们需要计算的东西就变成了a的d*2^r次方除以n的余数。于是，a^(d * 2^(r-1))要么等于1，要么等于n-1。如果a^(d * 2^(r-1))等于1，定理继续适用于a^(d * 2^(r-2))，这样不断开方开下去，直到对于某个i满足a^(d * 2^i) mod n = n-1或者最后指数中的2用完了得到的a^d mod n=1或n-1。
这样，Fermat小定理加强为如下形式：尽可能提取因子2，把n-1表示成d*2^r，如果n是一个素数，那么或者a^d mod n=1，或者存在某个i使得a^(d*2^i) mod n=n-1 ( 0<=i<r ) （注意i可以等于0，这就把a^d mod n=n-1的情况统一到后面去了）
**Miller-Rabin素性测试同样是不确定算法**，我们把可以通过以a为底的Miller-Rabin测试的合数称作以a为底的强伪素数(strong pseudoprime)。第一个以2为底的强伪素数为2047。第一个以2和3为底的强伪素数则大到1 373 653。

Miller-Rabin算法的代码也非常简单：计算d和r的值（可以用位运算加速），然后二分计算a^d mod n的值，最后把它平方r次。

### Pascal code

函数IsPrime返回对于特定的底数a，n是否是能通过测试。如果函数返回False，那说明n不是素数；如果函数返回True，那么n极有可能是素数。注意这个代码的数据范围限制在longint，你很可能需要把它们改成int64或高精度计算。

```Pascal
function pow( a, d, n:longint ):longint;
begin
   if d=0 then exit(1)
   else if d=1 then exit(a)
   else if d and 1=0 then exit( pow( a*a mod n, d div 2, n) mod n)
   else exit( (pow( a*a mod n, d div 2, n) * a) mod n);
end;
```

```Pascal
function IsPrime( a,n:longint ):boolean;
var
   d,t:longint;
begin
   if n=2 then exit(true);
   if (n=1) or (n and 1=0) then exit(false);
   d:=n-1;
   while d and 1=0 do d:=d shr 1;
   t:=pow( a, d, n );
   while ( d<>n-1 ) and ( t<>1 ) and ( t<>n-1 ) do
   begin
      t:=(t * t)mod n;
      d:=d shl 1;
   end;
   exit( (t=n-1) or (d and 1=1) );
end;
```

对于大数的素性判断，目前Miller-Rabin算法应用最广泛。一般底数仍然是随机选取，但当待测数不太大时，选择测试底数就有一些技巧了。比如，如果被测数小于4 759 123 141，那么只需要测试三个底数2, 7和61就足够了。当然，你测试的越多，正确的范围肯定也越大。如果你每次都用前7个素数(2, 3, 5, 7, 11, 13和17)进行测试，所有不超过341 550 071 728 320的数都是正确的。如果选用2, 3, 7, 61和24251作为底数，那么10^16内唯一的强伪素数为46 856 248 255 981。这样的一些结论使得Miller-Rabin算法在OI中非常实用。通常认为，Miller-Rabin素性测试的正确率可以令人接受，随机选取k个底数进行测试算法的失误率大概为4^(-k)。

Miller-Rabin算法是一个RP算法。RP是时间复杂度的一种，主要针对判定性问题。一个算法是RP算法表明它可以在多项式的时间里完成，对于答案为否定的情形能够准确做出判断，但同时它也有可能把对的判成错的（错误概率不能超过1/2）。RP算法是基于随机化的，因此多次运行该算法可以降低错误率。还有其它的素性测试算法也是概率型的，比如Solovay-Strassen算法。另外一些素性测试算法则需要预先知道一些辅助信息（比如n-1的质因子），或者需要待测数满足一些条件（比如待测数必须是2^n-1的形式）。前几年AKS算法轰动世界，它是第一个多项式的、确定的、无需其它条件的素性判断算法。当时一篇论文发表出来，题目就叫PRIMES is in P，然后整个世界都疯了，我们班有几个MM那天还来了初潮。算法主要基于下面的事实：n是一个素数当且仅当(x-a)^n≡(x^n-a) (mod n)。注意这个x是多项式中的未知数，等式两边各是一个多项式。举个例子来说，当a=1时命题等价于如下结论：当n是素数时，杨辉三角的第n+1行除两头的1以外其它的数都能被n整除。

### Python code

```python
import random

def miller_rabin(n, k):
    if n == 2:
        return True

    if n % 2 == 0:
        return False

    r, s = 0, n - 1
    while s % 2 == 0:
        r += 1
        s //= 2
    for _ in range(k):
        a = random.randrange(2, n - 1)
        x = pow(a, s, n)
        if x == 1 or x == n - 1:
            continue
        for _ in range(r - 1):
            x = pow(x, 2, n)
            if x == n - 1:
                break
        else:
            return False
    return True


if __name__ == "__main__":

  secp256k1 = 2**256 - 2**32 - 2**9 - 2**8 - 2**7 - 2**6 - 2**4 - 1
  ed25519 = 2*255 - 19

  miller_rabin(secp256k1, 61)
  miller_rabin(ed25519, 24251)
  miller_rabin(27742317777372353535851937790883648493, 24251)
```

### 使用 openssl 检测质数

openssl 内置了质数检查器，使用的就是 Miller-Rabin 算法

```python
openssl prime 115792089237316195423570985008687907853269984665640564039457584007908834671663

# FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFEFFFFFC2F is prime

openssl prime 57896044618658097711785492504343953926634992332820282019728792003956564819949

# 7FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFED is prime
```

---

参考：

https://github.com/ruby/prime
http://www.matrix67.com/blog/archives/234
https://security.stackexchange.com/a/176396