素性检验问题和模平方根问题
因为这两种算法都是随机化算法且都与数论问题有关,而且还有许多微妙的联系,因此放在一起整理.
素性检验问题
(主要参考资料:【朝夕的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$
- 任取$a$,计算$a^{p-1}(\mod p)$,
- 若结果不为$1$,则根据Fermat小定理,$p$必然是合数
- 否则进行下一步
- 由于$p$是奇数,那么$p-1$必然是偶数,记$t=\frac{p-1}{2}$
- 因为$(a^{t})^2 =1 (\mod p)$,计算$a^t (\mod p)$,
- 若$a^t (\mod p)$不为$1$或$p-1$,则同余方程$$x^2=1 (\mod p)$$存在着非$1,p-1$的解$x=a^t$,根据二次探测定理,$p$必然为合数.
- 若$a^t=1$,
- 当$t$为偶数时,令$t'=\frac{t}{2}$,由于$(a^{t'})^2=a^t=1 (\mod p)$,可回到步骤3继续判断
- 当$t$为奇数时,无法判断$p$的素性,结束程序
- 若$a^t=p-1$,此时无法判断$p$的素性,但可以不直接结束程序,而是继续计算$a^{t/2},a^{t/4},a^{t/8},\cdots$直到$t/2^{k}$除不开
- 最终没能使得$a^{t/2^k}=1$,则结束程序
- 找到了$a^{t/2^k}=1$,则令$t'=t/2^k$,回到步骤3.
由此可以得到一个高效的素性检测算法:
Miller-Rabin素性检测算法:
输入待检测的数$p\geq 2$
- 若$p=2$,直接输出"质数",若为其他偶数,直接输出"合数"
- 在$p=u\times 2^{t}$的形式下找到最大的$t$
- (随机)选取多个底数$a_1,a_2,\cdots,a_n$,对于每个底数$a_i$
- 计算$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}}$
- 若$b_t \neq 1$,输出"合数"
- 若存在$j$,使得$b_{j+1}=1,b_{j}\neq 1,b_j\neq p-1$,则输出"合数".
- 全部步骤完成后,输出"质数"
注意,该算法存在着微乎其微的错把合数当作质数的可能性
(据说)检测范围为int时,取底数为$\{2,7,61\}$可保证不出错,检测范围为long long 时,取底数为$\{2,325,9375,28178,450775,9780504,1795265022\}$可保证不出错.
代码:
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
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")
exit()
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]$上一半的数都是二次非剩余,也可以使用随机化方法枚举.
代码:
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 qpow_for_complex(base,unit,exp,mod): #对复数的快速幂,unit=i*i
ans = [1,0]
while(exp>0):
if(exp%2 == 1):
ans = mult_for_complex(ans,base,unit,mod)
exp = exp//2
base = mult_for_complex(base,base,unit,mod)
return ans
def mult_for_complex(a,b,unit,mod):
c=[(a[0]*b[0]+a[1]*b[1]*unit)%mod,(a[0]*b[1]+a[1]*b[0])%mod]
return c
a=int(input())
p=int(input())
if(a==0):
print(0)
elif(qpow(a,(p-1)//2,p)==p-1):
print(-1)
else:
for b in range(1,p):
if(qpow((b*b-a)%p,(p-1)//2,p)==p-1):
print('b=',b)
print('i*i=',(b*b-a)%p)
x=qpow_for_complex([b,1],(b*b-a)%p,(p+1)//2,p)
print(x[0],p-x[0])
break
