因为这两种算法都是随机化算法且都与数论问题有关,而且还有许多微妙的联系,因此放在一起整理.

素性检验问题

(主要参考资料:【朝夕的ACM笔记】数论-Miller Rabin素数判定 - 知乎 (zhihu.com))

(不完善的)Fermat素性检验:

由Fermat小定理可知,对于素数$p$,所有$a\in[1,p-1]$,$a^{p-1} =1 \mod p$

其逆命题是不成立的:

对于$p$,若存在$a\in[1,p-1]$,使得$a^{p-1} =1 \mod p$,则$p$是素数.

实际上,$p=341=11 \times 31$时,$1^{340},2^{340},4^{340},8^{340},15^{340},16^{340},23^{340},27^{340},29^{340},30^{340},32^{340},35^{340},39^{340},46^{340},47^{340},54^{340},58^{340},60^{340},61^{340},63^{340},64^{340},70^{340},78^{340},85^{340},89^{340},91^{340},92^{340},94^{340},95^{340},97^{340},101^{340},108^{340},109^{340},116^{340},120^{340},122^{340},123^{340},125^{340},126^{340},128^{340},139^{340},140^{340},147^{340},151^{340},153^{340},156^{340},157^{340},159^{340},163^{340},170^{340},171^{340},178^{340},182^{340},184^{340},185^{340},188^{340},190^{340},194^{340},201^{340},202^{340},213^{340},215^{340},216^{340},218^{340},219^{340},221^{340},225^{340},232^{340},233^{340},240^{340},244^{340},246^{340},247^{340},249^{340},250^{340},252^{340},256^{340},263^{340},271^{340},277^{340},278^{340},280^{340},281^{340},283^{340},287^{340},294^{340},295^{340},302^{340},306^{340},309^{340},311^{340},312^{340},314^{340},318^{340},325^{340},326^{340},333^{340},337^{340},339^{340},340^{340},$取模$341$均等于$1$

这类合数被称为伪Fermat数.

若强化此命题的前提条件,依然不成立:

对于$p$,若除了那些使得$gcd(a,p)\neq 1$以外,所有$a\in[1,p-1]$,均使得$a^{p-1} =1 \mod p$,则$p$是素数.

因为对于$p=561=3\times 11\times 17$时,$1^{560},2^{560},4^{560},5^{560},7^{560},8^{560},10^{560},13^{560},14^{560},16^{560},19^{560},20^{560},23^{560},25^{560},26^{560},28^{560},29^{560},31^{560},32^{560},35^{560},37^{560},38^{560},40^{560},41^{560},43^{560},46^{560},47^{560},49^{560},50^{560},52^{560},53^{560},56^{560},58^{560},59^{560},61^{560},62^{560},64^{560},65^{560},67^{560},70^{560},71^{560},73^{560},74^{560},76^{560},79^{560},80^{560},82^{560},83^{560},86^{560},89^{560},91^{560},92^{560},94^{560},95^{560},97^{560},98^{560},100^{560},101^{560},103^{560},104^{560},106^{560},107^{560},109^{560},112^{560},113^{560},115^{560},116^{560},118^{560},122^{560},124^{560},125^{560},127^{560},128^{560},130^{560},131^{560},133^{560},134^{560},137^{560},139^{560},140^{560},142^{560},145^{560},146^{560},148^{560},149^{560},151^{560},152^{560},155^{560},157^{560},158^{560},160^{560},161^{560},163^{560},164^{560},166^{560},167^{560},169^{560},172^{560},173^{560},175^{560},178^{560},179^{560},181^{560},182^{560},184^{560},185^{560},188^{560},190^{560},191^{560},193^{560},194^{560},196^{560},197^{560},199^{560},200^{560},202^{560},203^{560},205^{560},206^{560},208^{560},211^{560},212^{560},214^{560},215^{560},217^{560},218^{560},223^{560},224^{560},226^{560},227^{560},229^{560},230^{560},232^{560},233^{560},235^{560},236^{560},239^{560},241^{560},244^{560},245^{560},247^{560},248^{560},250^{560},251^{560},254^{560},256^{560},257^{560},259^{560},260^{560},262^{560},263^{560},265^{560},266^{560},268^{560},269^{560},271^{560},274^{560},277^{560},278^{560},280^{560},281^{560},283^{560},284^{560},287^{560},290^{560},292^{560},293^{560},295^{560},296^{560},298^{560},299^{560},301^{560},302^{560},304^{560},305^{560},307^{560},310^{560},311^{560},313^{560},314^{560},316^{560},317^{560},320^{560},322^{560},325^{560},326^{560},328^{560},329^{560},331^{560},332^{560},334^{560},335^{560},337^{560},338^{560},343^{560},344^{560},346^{560},347^{560},349^{560},350^{560},353^{560},355^{560},356^{560},358^{560},359^{560},361^{560},362^{560},364^{560},365^{560},367^{560},368^{560},370^{560},371^{560},373^{560},376^{560},377^{560},379^{560},380^{560},382^{560},383^{560},386^{560},388^{560},389^{560},392^{560},394^{560},395^{560},397^{560},398^{560},400^{560},401^{560},403^{560},404^{560},406^{560},409^{560},410^{560},412^{560},413^{560},415^{560},416^{560},419^{560},421^{560},422^{560},424^{560},427^{560},428^{560},430^{560},431^{560},433^{560},434^{560},436^{560},437^{560},439^{560},443^{560},445^{560},446^{560},448^{560},449^{560},452^{560},454^{560},455^{560},457^{560},458^{560},460^{560},461^{560},463^{560},464^{560},466^{560},467^{560},469^{560},470^{560},472^{560},475^{560},478^{560},479^{560},481^{560},482^{560},485^{560},487^{560},488^{560},490^{560},491^{560},494^{560},496^{560},497^{560},499^{560},500^{560},502^{560},503^{560},505^{560},508^{560},509^{560},511^{560},512^{560},514^{560},515^{560},518^{560},520^{560},521^{560},523^{560},524^{560},526^{560},529^{560},530^{560},532^{560},533^{560},535^{560},536^{560},538^{560},541^{560},542^{560},545^{560},547^{560},548^{560},551^{560},553^{560},554^{560},556^{560},557^{560},559^{560},560^{560},$取模$561$均等于$1$,注意看,和$561$互质的数全在这里了.

这类合数被称为Carmichael数.

再次强化命题的条件:

对于$p$,若所有$a\in[1,p-1]$,使得$a^{p-1} =1 \mod p$,则$p$是素数.

这次倒是成立了,但将这个命题实际应用,相比暴力枚举判断质数有任何优势吗......

因此需要考虑别的思路.

二次探测定理:

对于奇素数$p$,同余方程$x^2=1 (\mod p)$在$[0,p-1]$上的解有且仅有两个:$x_1=1,x_2=p-1$

证明:

$$\begin{aligned}
x^2-1=&0 (\mod p) \\
(x+1)(x-1)=&0 (\mod p)
\end{aligned}$$

因此$p|(x+1)(x-1)$,又因为$p$是素数,所以要么$(x+1)(x-1)=0$,要么$(x+1)(x-1)$是$p$的倍数,因此在在$[0,p-1]$上的解有且仅有两个:$x_1=1,x_2=p-1$

这个定理的用处在于其逆否命题:

若同余方程$x^2=1 (\mod p)$在$[0,p-1]$上存在着$1,p-1$以外的解,则$p$为合数.

按如下方法改进Fermat素性检验:

改进的Fermat素性检验:

输入:奇数$p$

1. 任取$a$,计算$a^{p-1}(\mod p)$,
   1. 若结果不为$1$,则根据Fermat小定理,$p$必然是合数
   2. 否则进行下一步
2. 由于$p$是奇数,那么$p-1$必然是偶数,记$t=\frac{p-1}{2}$
3. 因为$(a^{t})^2 =1 (\mod p)$,计算$a^t (\mod p)$,
   1. 若$a^t (\mod p)$不为$1$或$p-1$,则同余方程$$x^2=1 (\mod p)$$存在着非$1,p-1$的解$x=a^t$,根据二次探测定理,$p$必然为合数.
   2. 若$a^t=1$,
      1. 当$t$为偶数时,令$t'=\frac{t}{2}$,由于$(a^{t'})^2=a^t=1 (\mod p)$,可回到步骤3继续判断
      2. 当$t$为奇数时,无法判断$p$的素性,结束程序
   3. 若$a^t=p-1$,此时无法判断$p$的素性,但可以不直接结束程序,而是继续计算$a^{t/2},a^{t/4},a^{t/8},\cdots$直到$t/2^{k}$除不开
      1. 最终没能使得$a^{t/2^k}=1$,则结束程序
      2. 找到了$a^{t/2^k}=1$,则令$t'=t/2^k$,回到步骤3.

由此可以得到一个高效的素性检测算法:

Miller-Rabin素性检测算法:

输入待检测的数$p\geq 2$

1. 若$p=2$,直接输出"质数",若为其他偶数,直接输出"合数"
2. 在$p=u\times 2^{t}$的形式下找到最大的$t$
3. 选取多个底数$a_1,a_2,\cdots,a_n$,对于每个底数$a_i$
   1. 计算$b_1=a_i^u,b_2=a_i^{u\times 2},b_3=a_i^{u\times 2^2},b_4=a_i^{u\times 2^3},\cdots,b_t=a_i^{u \times 2^{t}}$
   2. 若$b_t \neq 1$,输出"合数"
   3. 若存在$j$,使得$b_{j+1}=1,b_{j}\neq 1,b_j\neq p-1$,则输出"合数".
4. 全部步骤完成后,输出"质数"

注意,该算法存在着微乎其微的错把合数当作质数的可能性

(据说)检测范围为int时,取底数为$\{2,7,61\}$可保证不出错,检测范围为long long 时,取底数为$\{2,325,9375,28178,450775,9780504,1795265022\}$可保证不出错.

代码:

import time
def divi(a): #分解质因数
    i=2
    ans=[]
    while(i*i<a):
        while(a%i==0):
            ans.append(i)
            a=a//i
        i=i+1
    if(a>1):
        ans.append(a)
    return ans

def gcd(a,b): #辗转相除法
    if(b == 0):
        return a
    else:
        return gcd(b,a%b)
    
def exgcd(a,b): #扩展欧几里得算法
    if(b == 0):
        x = 1
        y = 0
        return a,x,y
    ans,t,x = exgcd(b,a%b)
    y = t-a//b*x
    return ans,x,y

def qpow(base,exp,mod): #快速幂
    ans = 1
    while(exp>0):
        if(exp%2 == 1):
            ans = ans*base%mod
        exp = exp//2
        base = base*base%mod
    return ans

def rev(a,mod,type ='not_prime'):
    if(type =='prime'): #费马小定理求逆元
        return qpow(a,mod-2,mod)
    else: #辗转相除法求逆元
        g,x,y = exgcd(a,mod)
        if(g!= 1):
            return -1 #逆元不存在
        else:
            return (x+mod)%mod



arr=[2,325,9375,28178,450775,9780504,1795265022]
arr.append(int(time.time()*10000000%1000000007))# 添加随机的底数
p=int(input())
for i in range(len(arr)):
    arr[i]=arr[i]%p
print("base:",arr)
if(p==2):
    print("Prime")
    exit()
if(p % 2==0):
    print("Not prime")
    exit()

t=0
t_pow_2=1
u=p-1
while((p-1)%(t_pow_2*2)==0):
    t=t+1
    t_pow_2=t_pow_2*2
    u=u//2
print(p-1,'=',u,'*','2 ^ (',t,')')
for a in arr:
    last=qpow(a,u*qpow(2,t,p),p)
    print(a,'^ (',u,'*','2 ^ (',t,')',')','=',a,'^','(',u*qpow(2,t,p),')','=',last)
    if(last!=1):
        print("Not prime")
    for j in range(t-1,-1,-1):
        b=qpow(a,u*qpow(2,j,p),p)
        print(a,'^ (',u,'*','2 ^ (',j,')',')','=',a,'^','(',u*qpow(2,j,p),')','=',b)
        if(last==1 and b!=1 and b!=p-1):
            print("Not prime")
            exit()
        last=b
print("Prime")

求模意义下的平方根问题:

二次剩余与模平方根:

(主要参考资料:算法学习笔记(41): 二次剩余 - 知乎 (zhihu.com),二次剩余 - OI Wiki (oi-wiki.org))

回顾密码协议学习笔记(1):密码协议引论与密码学基础 - Isakovsky - 博客园 (cnblogs.com)中阶,原根的概念.

注意,这一部分只讨论模数是奇素数的情况

已知$a\in[0,p-1]$,在$[0,p-1]$上寻找$x$,使得$a=x^2 \mod p$

需要注意的是,并非所有模平方根问题都有解,如果上述方程有非零解,则称$a$是$p$的二次剩余,如果无解,则称$a$是$p$的二次非剩余.

为方便讨论,引入Legendre符号:

$$(\frac{a}{p}):=\left\{\begin{matrix}
1  &  \text{$a$是模$p$的二次剩余}\\
0  & a=0 \\
-1  & \text{$a$是模$p$的非二次剩余}
\end{matrix}\right.$$

为什么要这样定义呢?这是因为如下定理的存在

Euler准则:

对于奇质数$p$,$a\in[0,p-1]$,总有

$$(\frac{a}{p})=a^{\frac{p-1}{2}} (\mod p)$$

(神奇吧!)

证明:

首先证明$a^\frac{p-1}{2}$为什么只有$1,0,-1$三个值,

当$a=0$时,$a^\frac{p-1}{2}=0$是显然的.

当$a\neq 0$时根据Fermat小定理,$a^{p-1}=1 \mod p$,而$a^\frac{p-1}{2}$的值实际上是$x^2=1 (\mod p)$的一个解,根据二次探测定理,$a^\frac{p-1}{2}$要么是$1$,要么是$-1$.

下证,$a$是模$p$的二次剩余$\Leftrightarrow$ $a^\frac{p-1}{2}=1 (\mod p)$

$\Rightarrow$: 既然$a$是模$p$的二次剩余,那么设$a=x^2 (\mod p)$,则$a^\frac{p-1}{2}=x^{p-1}=1 (\mod p)$

$\Leftarrow$:

设$g$是模$p$的一个原根(奇素数必然存在原根),记$a=g^k$(原根能生成$[1,p-1]$中所有元素),$a^\frac{p-1}{2}=g^{\frac{k}{2}(p-1)}=1 (\mod p)$,因此,$\frac{k}{2}(p-1)$是$\varphi(p)=p-1$的倍数.因此,$k$是偶数,此时令$x=g^\frac{k}{2}$,可令$a=x^2 (\mod p)$,通过构造证明了$a$是模$p$的二次剩余.

那么,既然

$a^\frac{p-1}{2}$只有$1,0,-1$三个值,

$a$是模$p$的二次剩余$\Leftrightarrow$ $a^\frac{p-1}{2}=1 (\mod p)$,

而$a=0$时,$a^\frac{p-1}{2}=0$,

那剩下一种情况,$a$是模$p$的非二次剩余,自然有$a^\frac{p-1}{2}=-1 (\mod p)$

证毕.

Euler定理的推论(奇变偶不变):

当$\frac{p-1}{2}$是奇数时,

$a$是模$p$的二次剩余,则$p-a$是$p$的非二次剩余.

$a$是模$p$的非二次剩余,则$p-a$是$p$的二次剩余.

当$\frac{p-1}{2}$是偶数时,

$a$是模$p$的二次剩余,则$p-a$是$p$的二次剩余.

$a$是模$p$的非二次剩余,则$p-a$是$p$的非二次剩余.

二次剩余方程:

如果$a$是模$p$的二次剩余,该如何计算$a=x^2 (\mod p)$的解呢?

显然地,该方程有两个互为相反数(在模$p$意义下即相加得$p$)的解(实际上,模$[1,p-1]$中共有$\frac{p-1}{2}$个二次剩余,它们的模平方根没有重复,正好是$[1,p-1]$中的所有整数,证明略)

Cipolla算法:

首先找到$b$使得$(b^2-a)^\frac{p-1}{2}=-1$,那么根据Euler准则,$b^2-a$就是二次非剩余

虽然它是二次非剩余,不存在模平方根,但假设

$$i^2=b^2-a (\mod p)$$

(可以理解为在模意义下定义了虚数,$i$在计算的最后肯定是会被消掉的)

然后计算

$$(b+i)^\frac{p+1}{2}$$

实际上:

$$\begin{aligned}
 &(b+i)^\frac{p+1}{2}\\
=&((b+i)^{p+1})^\frac{1}{2}\\
=&((b+i)^p\cdot(b+i))^\frac{1}{2}\\
=&((b^p+i^p)\cdot(b+i))^\frac{1}{2} \qquad (\text{二项式展开后,其他项因为含$p$都被模掉了})\\
=&((b-i)\cdot(b+i))^\frac{1}{2} \qquad (\text{根据Fermat小定理,$b^p=b(\mod p),i^{p-1}=(b^2-b)^\frac{p-1}{2}=-1$})\\
=&((b-i)(b+i))^\frac{1}{2} \\
=&(b^2-i^2)^\frac{1}{2} \\
=&(b^2-b^2+a)^\frac{1}{2} \\
=&a^\frac{1}{2}
\end{aligned}$$

因为$a$是二次剩余,所以$(b+i)^\frac{p+1}{2}=a^\frac{1}{2}$必然在$[1,p-1]$上存在,这便是模平方根的一个解,另一个解是$p-(b+i)^\frac{p+1}{2}$

那么,接下来的问题是,怎么找到$b$?暴力枚举肯定可以,毕竟$[1,p-1]$上一半的数都是非二次剩余.