0
1

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 1 year has passed since last update.

MillerRabin 素性检测

Last updated at Posted at 2022-03-08

$ 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 有特殊的技巧

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

按照以上想法做个实现

#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表示成d2^r(其中d是一个奇数)。那么我们需要计算的东西就变成了a的d2^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表示成d2^r,如果n是一个素数,那么或者a^d mod n=1,或者存在某个i使得a^(d2^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或高精度计算。

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;
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

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 算法

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

0
1
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
0
1

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?