Pollard Rho and AMM算法
Pollard Rho and AMM算法
一.Pollard Rho算法
在O（n ** 1/4）时间复杂度内计算合数n的某个非平凡因子（非平凡因子：除了1和自身之外的因子）
前置知识：
1.快速幂
以求a的b次幂为例，首先将b转换成二进制
如b = 11,即b =1011 =2 ** 0 + 2 ** 1+2**3
求a的b次幂，转换成
完整代码：
# 定义常量 mod（模数）
mod = 10**7
# 快速幂函数 a ^ b % mod
def ksm(a, b, mod):
ans = 1
while b != 0:
if b & 1 != 0:
ans = (ans * a) % mod
a = (a * a) % mod
b >>= 1
return ans
快速乘与快速幂相似，代码如下：
# 快速乘函数 a * b % mod
def ksc(a, b, mod):
ans = 0
while b != 0:
if b & 1 != 0:
ans = (ans + a) % mod
a = (a + a) % mod
b >>= 1
return ans
2.Miller-Rabin 质数检验算法
这个算法可以在O（log n）复杂度中判断一个数是否为质数，运用了费马小定理和二次探测定理
费马小定理判断质数：
a ** (p-1) ≡ 1 mod p
当p为质数时，该定理成立,对于任意数a（a与p互素）,计算a **（p-1）mod p ，结果如果不是1，则一定不是质数。如果结果是1，不一定是质数（例如，卡迈克尔数（Carmichael numbers），是合数，但是满足费马小定理）
二次探测定理：
定义：若p是一个奇素数，且0 < x < p, 则方程x ** 2 ≡ 1 mod p 的解为x =1或x = p-1
Miller-Rabin算法：
设待测试的数为n，假设n为素数且为偶数，那么只有2，很容易判断。如果n为奇数，则n-1为偶数，而偶数可以表示成2 ** s * d的形式，那么n-1可以表示为 2 **s *d的形式，其中d为奇数，s>=0。
随机取一个整数a(1< a <n-1)
如果n为素数，结合费马小定理，则a ** （n-1） ≡ 1 mod n
又因为n- 1 = 2 ** s * d,那么，a ** (2 **s * d) ≡ 1 mod n,
令a * *d = k，k **(2 * *s) ≡ 1 mod n
（k * k）**( 2 * * s -1) ≡ 1 mod n
依次类推
结合Miller-Rabin算法，我们要找满足x ** 2 ≡ 1 mod n 的x
而当(k * k * k...* k)** 2 ≡ 1mod n
x = (k * k * k...* k)
所以我们只需要在2 ** s 变化过程中不断的判断，只要x = 1 或者n-1 则n为素数
但是，如果x一直不等于1或者n-1，不能说明n不是素数，但是，如果尝试几个k值，均不成立，那么n为素数的概率极小，可以将其视为合数
代码实现如下：
import random
def miller_rabin(n, k=5):
if n == 2 or n == 3:
return True
if n <= 1 or n % 2 == 0:
return False
# 将 n - 1 表示为 2^s * d 的形式
s, d = 0, n - 1
while d % 2 == 0:
s += 1
d //= 2
for _ in range(k):
a = random.randint(2, n - 2)
x = pow(a, d, n)
if x == 1 or x == n - 1:
continue
for _ in range(s - 1):
x = pow(x, 2, n)
if x == n - 1:
break
else:
return False
return True
在介绍Pollard Rho算法之前，先介绍一下生日悖论
如果一年只有365天（不计算闰年），且每个人的生日是任意一天的概率均相等，那么只要随机选取23人，就有50%以上的概率有两人同一天生日
解释：第一个人不会和前面的人重生日（因为前面没有人），第二个人不重生日的概率为364/365，第三个人363/365……以此类推，那么只要到第23个人，就有
，说明这时就有50%以上的概率有两人同生日
更多的，当一年有n天时，只要人数到达Θ(sqrt(n))的数量级，有至少两个人同一天生日的概率就可以达到50%以上
利用生日悖论来因数分解，重要的思想就是随机。
已知我们随机地从[2,N-1]中选择出一个数是N的因数的概率是极小的，这也就意味着需要重复随机选择来提高正确率。
那么如果我们不是只选择一个数，而是选择k个数，保证其中两个数的差值是N的因数。
而如果其中两个数x，y的差值为N的因数，就一定会有gcd(|x-y|,N)>1
假设N只有两个因数(除去自己和1)p和q的情况下，那么就意味着此时只有这两个数能整除N
但是如果我们要求的是有多少个数x保证gcd(x,N)>1，此时答案就很多了，有p+q-2个
于是我们就得出了一个简单的策略：
1.在区间[2,N-1]中随机选取k个数，x1~xk
2.判断是否存在gcd(xi-xj,N)>1，若存在，则显然gcd(xi-xj,N)是N的一个因数
同时也出现了一个问题，就是我们需要选取大约（N¼）个数，内存显然是不可能够的，那么又要如何解决这个问题呢？于是Pollard Rho算法就出现了
Pollard Rho算法：
因为数字太多内存不够，所以Pollard Rho只把连续的两个数存在内存中。
也就是说，我们并不会选出k个数，而是一个一个地生成伪随机数，并检查连续的两个数是否符合条件
我们需要一个模 n 的多项式，例如：g(x）=x ** 2 +1 mod n，用于生成“伪随机”数列并且选择一个起始值
代码实现：
def gcd(a, b):
while b:
a, b = b, a%b
return a
def mapx(x):
x=(x*x+1)%n
return x
def pollard_rho (x1,x2):
while True:
x1=mapx(x1)
x2=mapx(mapx(x2))
p=gcd(x1-x2,n)
if (p == n):
print("fail")
return
elif (p!=1):
print("p: "+str(p))
print("q: "+str(n/p))
break
例题：
from Crypto.Util.number import *
flag='XXXXXXXXXX'
e=0x10001
n=3464115689260819392935656139231271022088014497600959975252672820761470484618617542699739764705620767566046150296286140279466041905437740736319886548924058066340624280173573039937440574809212831672936265975209972857747738693288788052694888593629820783280181774594890101819911390468376042288685444703477639361249422054925155575158181408046447773183467474182159096249846479461475003039637685547191529769071424947165685996675043741359728960138725130116665515880652680244470002603320184043266997163009799067135481330853926800023087208636366543210276325733789567957712475079676808820490428973148590780103404729397985019089646026534758042652163974037179260930147000267787012874887703048185189387666199961350969478997330352491403611551096779887936063702892671572230523920277869824082318417366635732798078805070773662305098133013800559413976855639167085996705265872418286939474978058429877355305699191892248323534097124680592084808397263288943661630283254265898940430505371463275516870915261799702701453110383723660477322192658118815861974573497885759332587457883589126213868868281682733732394640188383458580555618181208722172011645063817887462511925363883111475043389367021569982583145912277082528397480756465818166640691294267829613428532584923122279261412957652411789780285394451850006055483598427256129336822790446219260722224551892755568723667049099823707012236618539013908897459508493477299241358343681962043777720375898238069411747593247677719659288861028038609681592049271388766040068274682329498410270047730060645409396168652996740411436504098243588534151795927613112273405281530700440843961363083171583053608594479571127638861391879160062136043924515492933057094406001740725150570943724333872057035272515676895138632004436149126542560186729495087753846670118167613304943153665660470977150177815354160112481310498255965933157902965011317272734556307544837058247794721427834401390455487872367113056168812990496943458920406382650794656884209032284257748450399719067686010156664485086500184078244358783696803659745365860507822269113477617988541658390715640119076036723560146145256418191555730429
m=bytes_to_long(flag)
c=pow(m,e,n)
# print(c)
assert c==1348180992713820685594359831085361247035354818267759694375727771380213764105201303900290589048476671838784929302557807429509301123621642066527762448168417430350203900517746161928291953862336621873868684964993481406952303460424386080800907869577832154132021624995347121475139439850224598911049824567475034545196149320206620924509546294914370633958171840525667753751028800452370972967842349036965813468093121305475860652163573748880159896028384459277781803606279654121781761535513864212749202340439422384362197321739535686506006254445152569135807838677803047115947240251817017395748631234023225248868640103807357800695487118067385969130665294711887439861169737031598274672286881054564295019370600398079829687943670066246766386800967129670975899437309257181151390273864539453238400821713398980792443187238155205199387468000884390778225928491109565539772444103430764411852162036780535068932010631521096301305200494731040799554079996593223419379540630567412342494359961515157700054387356888773321804429041668331430841954949061238002355270616266853430676140296078177898764114420743905073496118755793091902224563802923183897542996809393707165837427396760691616506592877656939455806776333105152244073459704842848097376929745632888315642056267589901527329593161439752243318158600747713335265681764274061086952626629147275106502354039293870614872849851780761537660716667003562253735923584765118319171207551032289722689929226947535724389015974764265109970513715649126340709399701638880824570967389428525179746112494690443526404899052482884526335312491500994658519105204100875596292743155076199283666094556235485722975240358390560161706163569618302627055238860096406627340558428533728688904747040262153892623865784320600845613247625784683300197845026108394831724602330903266980995146046956218083954561639114586880305451066793209385702886944430548686040360906088479432852330611586273784509957389880104225534635731572089138568407411433675077699382999593384336070450058824669351046095205147266720901917218217564113964074957585348227506587531149811564371155522575778721899999227020238022970334707476977210019
题目只给了n,e,c的值，很明显要分解n得到p和q,但是n太大，又6948bit,无法用yafu和在线网站[factordb.com]分解，并且可以猜测n不可能只有p，q两个因数，而是由若干个小素数乘积得到的
这里用Pollard Rho算法来分解n
import gmpy2
def mapx(x):
x =(x * x + 1) % n
return x
def pollard_rho(x1,x2):
while True:
x1 = mapx(x1)
x2 = mapx(mapx(x2))
p =gmpy2.gcd(x1-x2,n)
if p == n:
print("fail")
elif p != 1:
return p
n = 3464115689260819392935656139231271022088014497600959975252672820761470484618617542699739764705620767566046150296286140279466041905437740736319886548924058066340624280173573039937440574809212831672936265975209972857747738693288788052694888593629820783280181774594890101819911390468376042288685444703477639361249422054925155575158181408046447773183467474182159096249846479461475003039637685547191529769071424947165685996675043741359728960138725130116665515880652680244470002603320184043266997163009799067135481330853926800023087208636366543210276325733789567957712475079676808820490428973148590780103404729397985019089646026534758042652163974037179260930147000267787012874887703048185189387666199961350969478997330352491403611551096779887936063702892671572230523920277869824082318417366635732798078805070773662305098133013800559413976855639167085996705265872418286939474978058429877355305699191892248323534097124680592084808397263288943661630283254265898940430505371463275516870915261799702701453110383723660477322192658118815861974573497885759332587457883589126213868868281682733732394640188383458580555618181208722172011645063817887462511925363883111475043389367021569982583145912277082528397480756465818166640691294267829613428532584923122279261412957652411789780285394451850006055483598427256129336822790446219260722224551892755568723667049099823707012236618539013908897459508493477299241358343681962043777720375898238069411747593247677719659288861028038609681592049271388766040068274682329498410270047730060645409396168652996740411436504098243588534151795927613112273405281530700440843961363083171583053608594479571127638861391879160062136043924515492933057094406001740725150570943724333872057035272515676895138632004436149126542560186729495087753846670118167613304943153665660470977150177815354160112481310498255965933157902965011317272734556307544837058247794721427834401390455487872367113056168812990496943458920406382650794656884209032284257748450399719067686010156664485086500184078244358783696803659745365860507822269113477617988541658390715640119076036723560146145256418191555730429
while (n!=1):
p=pollard_rho(2,2)
print(p)
n=n//p
运行一段时间，有如下结果，当然，还没有分解完，这只是一部分结果
观察结果，前6个均为740160703366837，长度为50bit
flag的大小不是很大，就将740160703366837的6次方作为新的n值
from gmpy2 import *
from libnum import *
p = 740160703366837
if is_prime(p):
n = pow(p,6)
e = 0x10001
c = 1348180992713820685594359831085361247035354818267759694375727771380213764105201303900290589048476671838784929302557807429509301123621642066527762448168417430350203900517746161928291953862336621873868684964993481406952303460424386080800907869577832154132021624995347121475139439850224598911049824567475034545196149320206620924509546294914370633958171840525667753751028800452370972967842349036965813468093121305475860652163573748880159896028384459277781803606279654121781761535513864212749202340439422384362197321739535686506006254445152569135807838677803047115947240251817017395748631234023225248868640103807357800695487118067385969130665294711887439861169737031598274672286881054564295019370600398079829687943670066246766386800967129670975899437309257181151390273864539453238400821713398980792443187238155205199387468000884390778225928491109565539772444103430764411852162036780535068932010631521096301305200494731040799554079996593223419379540630567412342494359961515157700054387356888773321804429041668331430841954949061238002355270616266853430676140296078177898764114420743905073496118755793091902224563802923183897542996809393707165837427396760691616506592877656939455806776333105152244073459704842848097376929745632888315642056267589901527329593161439752243318158600747713335265681764274061086952626629147275106502354039293870614872849851780761537660716667003562253735923584765118319171207551032289722689929226947535724389015974764265109970513715649126340709399701638880824570967389428525179746112494690443526404899052482884526335312491500994658519105204100875596292743155076199283666094556235485722975240358390560161706163569618302627055238860096406627340558428533728688904747040262153892623865784320600845613247625784683300197845026108394831724602330903266980995146046956218083954561639114586880305451066793209385702886944430548686040360906088479432852330611586273784509957389880104225534635731572089138568407411433675077699382999593384336070450058824669351046095205147266720901917218217564113964074957585348227506587531149811564371155522575778721899999227020238022970334707476977210019
phi = (p-1)*(p**5)
d = invert(e,phi)
m = pow(c,d,n)
print(n2s(m))
这里将n替换成了新的值，前提条件：m < n, c ≡ m ** e mod n ,并且gcd(e, phi(n)) =1,一般情况不允许篡改n
对于这句（phi = (p-1)*(p**5)）代码，一开始没弄明白，后来做其他的题，遇到类似的问题，尝试询问AI，得到解决
一般情况：
对于n = p ** k
在1到p ** k中，由于p为质数，所以当某个数是p的倍数时
与n不互质，一共有p ** (k-1)个
所以
二.AMM算法
前置知识：
1.二次剩余
x ** 2 ≡ a mod m, gcd(a,m)=1
如果上式x有解，则a叫做模m的二次剩余，否则，非二次剩余
2.欧拉判别条件
设p为奇素数，并且（a,p）=1
(1)a是模p的二次剩余的充分必要条件 ：a **(p-1)/2 ≡ 1 mod p
(2)a是模p的非二次剩余的充分必要条件 ： a **(p-1)/2 ≡ -1 mod p
AMM2（在有限域下开平方根）
x ** 2 ≡ a mod p
p是奇素数，p-1表示成2 ** t * s, s为奇数
根据欧拉判别条件，可以有以下式子：
1.当t-1 = 0 ，即t=1时
2.t-1>=2时
算法代码：
def AMM2(residue,p):
# 计算t,s
t = 0
s = p - 1
while s % 2 == 0:
t += 1
s = s // 2
if t==1:
return pow(residue,(s+1)//2,p)
#生成非二次剩余
noresidue=random.randint(1,p)
while pow(noresidue,(p-1)//2,p)==1:
noresidue = random.randint(1,p)
h=1
a=pow(noresidue,s,p)
b=pow(residue,s,p)
for i in range(1,t):
zhishu=2**(t-1-i)
d=pow(b,zhishu,p)
if d==1:
k=0
else:
k=1
b=b * pow(a, 2*k, p)%p
h=pow(a,k,p)*h%p
a=pow(a,2,p)
return h*pow(residue,(s+1)//2,p)%p
AMM（在有限域下开n次方根）
两种情况：
（1）gcd (r , q-1) = 1
（2）r | q-1
第一种情况，因为逆元存在的条件为gcd(a,b) = 1
所以r存在逆元，
第二种情况：
p-1 = r ** t * s,并且有gcd(r,s)=1
很容易找到一个非负整数α，满足s|r* α - 1，相当于存在r * α -1 = k * s
对于模p的r次剩余δ有δ ** （p-1）/ r ≡ 1 mod p
代码如下：
def AMM(residue,r,p):
# 计算t,s
t = 0
s = p - 1
while s % r == 0:
t += 1
s = s // r
#计算alpha
k=1
while ( k * s + 1 ) % r != 0:
k+=1
alpha=(k*s+1)//r
#当t=1时的情况
if t==1:
return pow(residue,alpha,p)
#t>=2的情况
#计算一个r次非剩余
noresidue = random.randint(1, p)
while pow(noresidue, (p - 1) // r, p) == 1:
noresidue = random.randint(1, p)
a=pow( noresidue , ( r ** ( t-1 )) * s, p )
b=pow( residue , r * alpha - 1, p )
c=pow( noresidue , s , p )
h=1
for i in range(1,t):
d=pow( b , r ** ( t-1-i ) , p )
if d==1:
j=0
else:
j=-math.log(d,a)
b=pow( pow ( c , r , p ) , j , p ) * b % p
h=pow( c , j , p ) * h % p
c=pow( c , r , p )
return pow( residue , alpha , p ) * h % p
浙公网安备 33010602011771号