素性检验问题和模平方根问题

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

素性检验问题

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

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

由Fermat小定理可知,对于素数p,所有a[1,p1],ap1=1modp

其逆命题是不成立的:

对于p,若存在a[1,p1],使得ap1=1modp,则p是素数.

实际上,p=341=11×31时,1340,2340,4340,8340,15340,16340,23340,27340,29340,30340,32340,35340,39340,46340,47340,54340,58340,60340,61340,63340,64340,70340,78340,85340,89340,91340,92340,94340,95340,97340,101340,108340,109340,116340,120340,122340,123340,125340,126340,128340,139340,140340,147340,151340,153340,156340,157340,159340,163340,170340,171340,178340,182340,184340,185340,188340,190340,194340,201340,202340,213340,215340,216340,218340,219340,221340,225340,232340,233340,240340,244340,246340,247340,249340,250340,252340,256340,263340,271340,277340,278340,280340,281340,283340,287340,294340,295340,302340,306340,309340,311340,312340,314340,318340,325340,326340,333340,337340,339340,340340,取模341均等于1

这类合数被称为伪Fermat数.

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

对于p,若除了那些使得gcd(a,p)1以外,所有a[1,p1],均使得ap1=1modp,则p是素数.

因为对于p=561=3×11×17时,1560,2560,4560,5560,7560,8560,10560,13560,14560,16560,19560,20560,23560,25560,26560,28560,29560,31560,32560,35560,37560,38560,40560,41560,43560,46560,47560,49560,50560,52560,53560,56560,58560,59560,61560,62560,64560,65560,67560,70560,71560,73560,74560,76560,79560,80560,82560,83560,86560,89560,91560,92560,94560,95560,97560,98560,100560,101560,103560,104560,106560,107560,109560,112560,113560,115560,116560,118560,122560,124560,125560,127560,128560,130560,131560,133560,134560,137560,139560,140560,142560,145560,146560,148560,149560,151560,152560,155560,157560,158560,160560,161560,163560,164560,166560,167560,169560,172560,173560,175560,178560,179560,181560,182560,184560,185560,188560,190560,191560,193560,194560,196560,197560,199560,200560,202560,203560,205560,206560,208560,211560,212560,214560,215560,217560,218560,223560,224560,226560,227560,229560,230560,232560,233560,235560,236560,239560,241560,244560,245560,247560,248560,250560,251560,254560,256560,257560,259560,260560,262560,263560,265560,266560,268560,269560,271560,274560,277560,278560,280560,281560,283560,284560,287560,290560,292560,293560,295560,296560,298560,299560,301560,302560,304560,305560,307560,310560,311560,313560,314560,316560,317560,320560,322560,325560,326560,328560,329560,331560,332560,334560,335560,337560,338560,343560,344560,346560,347560,349560,350560,353560,355560,356560,358560,359560,361560,362560,364560,365560,367560,368560,370560,371560,373560,376560,377560,379560,380560,382560,383560,386560,388560,389560,392560,394560,395560,397560,398560,400560,401560,403560,404560,406560,409560,410560,412560,413560,415560,416560,419560,421560,422560,424560,427560,428560,430560,431560,433560,434560,436560,437560,439560,443560,445560,446560,448560,449560,452560,454560,455560,457560,458560,460560,461560,463560,464560,466560,467560,469560,470560,472560,475560,478560,479560,481560,482560,485560,487560,488560,490560,491560,494560,496560,497560,499560,500560,502560,503560,505560,508560,509560,511560,512560,514560,515560,518560,520560,521560,523560,524560,526560,529560,530560,532560,533560,535560,536560,538560,541560,542560,545560,547560,548560,551560,553560,554560,556560,557560,559560,560560,取模561均等于1,注意看,和561互质的数全在这里了.

这类合数被称为Carmichael数.

再次强化命题的条件:

对于p,若所有a[1,p1],使得ap1=1modp,则p是素数.

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

因此需要考虑别的思路.

二次探测定理:

对于奇素数p,同余方程x2=1(modp)[0,p1]上的解有且仅有两个:x1=1,x2=p1

证明:

x21=0(modp)(x+1)(x1)=0(modp)

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

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

若同余方程x2=1(modp)[0,p1]上存在着1,p1以外的解,则p为合数.

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

改进的Fermat素性检验:

输入:奇数p

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

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

Miller-Rabin素性检测算法:

输入待检测的数p2

  1. p=2,直接输出"质数",若为其他偶数,直接输出"合数"
  2. p=u×2t的形式下找到最大的t
  3. (随机)选取多个底数a1,a2,,an,对于每个底数ai
    1. 计算b1=aiu,b2=aiu×2,b3=aiu×22,b4=aiu×23,,bt=aiu×2t
    2. bt1,输出"合数"
    3. 若存在j,使得bj+1=1,bj1,bjp1,则输出"合数".
  4. 全部步骤完成后,输出"质数"

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

(据说)检测范围为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[0,p1],在[0,p1]上寻找x,使得a=x2modp

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

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

(ap):={1a是模p的二次剩余0a=01a是模p的二次非剩余

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

Euler准则:

对于奇质数p和整数a[0,p1],总有

(ap)=ap12(modp)

(神奇吧!)

证明:

首先证明ap12为什么只有1,0,1三个值,

a=0时,ap12=0是显然的.

a0时根据Fermat小定理,ap1=1modp,而ap12的值实际上是x2=1(modp)的一个解,根据二次探测定理,ap12要么是1,要么是1.

下证,a是模p的二次剩余 ap12=1(modp)

: 既然a是模p的二次剩余,那么设a=x2(modp),则ap12=xp1=1(modp)

:

g是模p的一个原根(奇素数必然存在原根),记a=gk(原根能生成[1,p1]中所有元素),ap12=gk2(p1)=1(modp),因此,k2(p1)φ(p)=p1的倍数.因此,k是偶数,此时令x=gk2,可令a=x2(modp),通过构造证明了a是模p的二次剩余.

那么,既然

ap12只有1,0,1三个值,

a是模p的二次剩余 ap12=1(modp),

a=0时,ap12=0,

那剩下一种情况,a是模p的二次非剩余,自然当且仅当ap12=1(modp)

证毕.

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

p12是奇数时,

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

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

p12是偶数时,

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

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

二次剩余方程:

如果a是模p的二次剩余,该如何计算a=x2(modp)的解呢?

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

Cipolla算法:

首先找到b使得(b2a)p12=1,那么根据Euler准则,b2a就是二次非剩余

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

i2=b2a(modp)

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

然后参考复数的运算,计算

(b+i)p+12

实际上:

(b+i)p+12=((b+i)p+1)12=((b+i)p(b+i))12=((bp+ip)(b+i))12(二项式展开后,其他项因为含p都被模掉了)=((bi)(b+i))12(根据Fermat小定理,bp=b(modp),ip1=(b2b)p12=1)=((bi)(b+i))12=(b2i2)12=(b2b2+a)12=a12

因为a是二次剩余,所以(b+i)p+12=a12必然在[1,p1]上存在,这便是模平方根的一个解,另一个解是p(b+i)p+12

那么,接下来的问题是,怎么找到b?暴力枚举肯定可以,毕竟[1,p1]上一半的数都是二次非剩余,也可以使用随机化方法枚举.

代码:

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

 

posted @   Isakovsky  阅读(82)  评论(0编辑  收藏  举报
相关博文:
阅读排行:
· TypeScript + Deepseek 打造卜卦网站:技术与玄学的结合
· Manus的开源复刻OpenManus初探
· AI 智能体引爆开源社区「GitHub 热点速览」
· 从HTTP原因短语缺失研究HTTP/2和HTTP/3的设计差异
· 三行代码完成国际化适配,妙~啊~
点击右上角即可分享
微信分享提示