RSA

现代密码

一些基础数学

整除

如果$a=kb(k\geq1)$，称作$a$可被$b$整除，且称$a$是$b$的倍数，$b$是$a$的约数，记作$b\mid a$。

同余

如果$a-b=km$，即$m\mid (a-b)$，则称$a$与$b$同余，记作$a\equiv b(mod\ m)$。

互素

如果$gcd(a,b)=1$，则称$a,b$互素。

欧拉函数

任意整数$n$，$\forall x<n,x\in Z_n^*$，$Z_n$表示模$n$时形成的集合，即从$0$到$n-1$的集合，$Z_n^*$则相对$Z_n$又去掉$0$这个元素。欧拉函数便是在$Z_n^*$这个集合中，与$n$互素的元素个数。记作$\varphi(n)$。

情景1

如果$n=1$，则$\varphi(n)=1$。

情景2

如果$n$是质数，则$\varphi(n)=n-1$。因为质数与小于它的每个数都互素。

情景3

$n=p^{k}$，则$\varphi(n)=p^{k-1}\varphi(p)$。

情景4

$n=pq$，$p,q$互素，则$\varphi(n)=\varphi(p)\varphi(q)$。

标准分解式

设有标准分解式$n=p_1^{\alpha_1}p_2^{\alpha_2}\cdots p_k^{\alpha_k}$,那么

$\varphi(n)=\prod_{i=1}^k(p_i^{\alpha_i}-p_i^{\alpha_i-1})=n\prod_{i=1}^k(1-\frac{1}{p_i})$

拓展情景

复数域下的欧拉函数

$\varphi(n)=(p^2-1)*(q^2-1)$

三阶矩阵群的阶

阶与欧拉函数有相似性，看不懂先跳过即可

$order=p(p+1)(p-1)(p^2+p+1)$

逆元

如果两个数存在$ab\equiv 1(mod\ m)$的关系，就称$a$是$b$在模$m$下的逆元，或者$b$是$a$在模$m$下的逆元。对这个式子进行变换即可得到$ab-1=km$。

欧拉定理

在$a$与$p$互素的情况下，存在这样的关系

$a^{\varphi(p)}\equiv1(mod\ p)$

这个式子便称为欧拉定理。在$p$是素数的情况下，式子即为$a^{p-1}\equiv1(mod\ p)$，这个特殊式子又称费马小定理。关于欧拉定理的证明有些繁琐，就当结论记住吧。

中国剩余定理(crt)

理论了解推荐Bintou，讲得非常清楚明白；crt在库中实际上集成好，只要能调用即可，from sympy.ntheory.modular import crt这是sympy库中的crt，sage中自带有crt，但两者使用有点区别，sympy库中crt模数在前，形如crt([p, q], [i, j])，而sage模数在后，形如crt([i, j], [p, q])

RSA

RSA是一种模数下的幂运算，具体运算规则如下

$c\equiv m^e\mod p\\ ed\equiv1\mod \varphi(p)\\ m\equiv c^d\mod p$

关于ed在指数上为什么是模欧拉函数而不是模数本身的证明

$m^e\equiv c\mod n,\quad m^{\varphi(n)}\equiv 1\mod n\\ log_me=c,\quad log_m\varphi(n)=1,\quad log_m(e+\varphi(n))=c\\ m^{e+\varphi(n)}=c$

快速幂，快速幂算法应该熟悉吧，在大整数运算时不可能先算乘方再取模，在acm里面就会用到快速幂，而python的pow函数有模形式，看这种算法表达形式要记住，在ECC里面会出现。

ll qpow(ll a,ll b)
{
    ll res=1;
    while(b)
    {
        if(b&1) res=res*a%mod;
        a=a*a%mod;
        b>>=1;
    }
    return res;
}

最简单的加密程序

from Crypto.Util.number import *
m = b'hello world'
m = bytes_to_long(m)
p = getPrime(256)
e = 65537
c = pow(m, e, p)
print(f"p={p}, c={c}")
# p=104572425491036779976425338195598909418460679281450078773956190724671461860829, c=15413370856595181721630593414126186156343530569374585851966215761918789311844

m是明文，c是密文，p是模数，e是指数(e与$\varphi(p)$要互素，上面提到求逆元时需要互素的条件，不然无法算出逆元，若不互素恢复时就不能常规算出d，需要其他方法入手)。

更加官话点的RSA如下：任意选取两个素数$p,q$，计算出$n=pq$，选取一个$e$，满足$gcd(e,\varphi(n))=1$，$e$作为公钥，同时生成私钥d，满足$ed\equiv1\mod\varphi(n)$，$e,d$是模$\varphi(n)$下的逆元组。

在指数运算中存在欧拉定理满足

$m^{\varphi(n)}\equiv1\mod n,\quad ed=k*\varphi(n)+1\\ m^{k*\varphi(n)+1}=m^{k*\varphi(n)}*m=(m^{\varphi(n)})^{k}*m\equiv 1^k*m\mod n\\$

在模n条件下，$(m^{\varphi(n)})^{k}$必为1，即得到m

加密$c\equiv m^{e}(mod\ n)$

解密$m\equiv c^{d}(mod\ n)$

from Crypto.Util.number import long_to_bytes, inverse
p = 
q = 
e = 
c = 
n = p * q
phi = (p - 1) * (q - 1)

d = inverse(e, phi)
print(long_to_bytes(pow(c, d, n)))

关于$RSA$的背景说明，个人比较推荐阮一峰的博客；攻击脚本总结，还得是Lazzaro；各类攻击理论证明和脚本运用，独奏的文章更好。

python中一些常用的库

pip install pycryptodome gmpy2 libnum sympy pwn

# Crypto.Util.number
long_to_bytes(a)	# 将整数a转换成bytes
bytes_to_long(str.encode())	# 将bytes类型转换成整数
#	不推荐用long2str和str2long,是上面两函数早期版本,虽然也能实现功能,但会警告说已经被这两函数取代
inverse(e, n)	# 求e在模n下的逆元
isPrime(m)		# 判断m是否为质数
getPrime(512)	# 生成一个二进制(bits)长度为512的质数，一般审题会出现

# gmpy2
gcd()
lcm()
invert(e, n)	# 与inverse()功能一样
iroot(m, e)		# 对m开e次方，形成一个列表，列表中第二项表示开e次方是否整数
next_prime(m)	# 求m的下一个质数
gcd_ext			# 拓展欧几里得算法

# libnum(Crypto安装不成功可以作为替代)
s2n(m)	# 字符m转成十进制,与bytes_to_long()一样
s2b(m)	# 字符m转成二进制
n2s(m)	# 整数m转成字符,与long_to_bytes()一样
b2s(m)	# 二进制转成字符
invmod(e, n)	# e在模n下的逆元
xgcd()	# 拓展欧几里得算法

# sympy
nextprime(x)		# x的下一个质数
prevprime(x)		# x的上一个质数

m.to_bytes((m.bit_length() + 7) // 8, 'big')	# 大端序
m.to_bytes((m.bit_length() + 7) // 8, 'little')	# 小端序
int.from_bytes(m, 'big')
int.from_bytes(m, 'little')

RSA各类变量及其方法

$d\equiv e^{-1}\mod (p-1)*(q-1)\\ dp\equiv d\mod p-1\\ dq\equiv d\mod q-1$

低指数爆破

当$e$非常小时，$m^e$也会非常小，而如果$n$较大，那么$m^e-c=kn$这个恒等式中$k$也不会很大。甚至很多题目中$m^e$直接比$n$小。

$m^e=kn+c\\ m=\sqrt[e]{kn+c}$

from gmpy2.gmpy2 import iroot
from Crypto.Util.number import *
n = 
e = 
c = 
for i in range(100000):
	t = iroot(i * n + c, e)		# 可开e次方时为1,否则为0
	if t[1] == 1:
		print(long_to_bytes(t[0]))
		break

例题

n = 10888751337932558679268839254528888070769213269691871364279830513893837690735136476085167796992556016532860022833558342573454036339582519895539110327482234861870963870144864609120375793020750736090740376786289878349313047032806974605398302398698622431086259032473375162446051603492310000290666366063094482985737032132318650015539702912720882013509099961316767073167848437729826084449943115059234376990825162006299979071912964494228966947974497569783878833130690399504361180345909411669130822346252539746722020515514544334793717997364522192699435604525968953070151642912274210943050922313389271251805397541777241902027
e = 3
c = 2449457955338174702664398437699732241330055959255401949300755756893329242892325068765174475595370736008843435168081093064803408113260941928784442707977000585466461075146434876354981528996602615111767938231799146073229307631775810351487333

共模攻击

$c_1\equiv m^{e_1}(mod\ n),\quad c_2\equiv m^{e_2}(mod\ n)$

如果$gcd(e_1,e_2)=1$，拓展欧几里得公式可得$x\cdot e_1+y\cdot e_2=1$

$c_1^{x}\cdot c_2^{y}\equiv m^{x\cdot e_1}\cdot m^{y\cdot e_2}\equiv m^{x\cdot e_1+y\cdot e_2}\equiv m(mod\ n)$

import gmpy2
n = 
c1 = 
c2 = 
e1 = 
e2 = 
_, x, y = gmpy2.gcdext(e1, e2)
m = pow(c1, x, n) * pow(c2, y, n) % n
print(m.to_bytes((m.bit_length() + 7) // 8, 'big'))

e不互素+低指数

$RSA$解密本质便是构造逆元，使得$m^{ed}\equiv1(mod\ n)$，而当$gcd(e,phi)\neq1$，便无法直接构造出这样的逆元。如若可讲$e$分解为素数的乘积$e=ab$形式(方便表示便只取两个元素)，如$gcd(a,phi)=1$，可构造$a$逆元，使原式化简$m^{ab\cdot a^{-1}}=m^{b}$，此时$gcd(b,phi)\neq1$，不能再次化简。如果$b$较小，且$m^b$长度能小于$n$，便能直接用开方形式算出$m$确切值。如$b$仍相对较大，可以考虑用其他方法进行不断化简。

$c\equiv m^{e}(mod\ n),c\equiv m^{ab}(mod\ n)$

设$a$关于$phi$逆元为$t$

$c^{t}\equiv m^{b}(mod\ n)$

此时$c^t$便是$m^b$在模$n$下确值，如$m^b$位数小于$n$​,直接开方便可

$m=\sqrt[t]{c^t}$

exp

from Crypto.Util.number import *
from gmpy2 import iroot
def rsa_not_gcd(p, q, e, c):
    n = p*q
    phi = (p-1)*(q-1)
    g = GCD(e, phi)
    c = pow(c, inverse(e//g, phi), n)
    m = iroot(c, g)[0]
    return m

例题

p = 134261118796789547851478407090640074022214132682000430136383795981942884853000826171189906102866323044078348933419038543719361923320694974970600426450755845839235949167391987970330836004768360774676424958554946699767582105556239177450470656065560178592346659948800891455240736405480828554486592172443394370831
q = 147847444534152128997546931602292266094740889347154192420554904651813340915744328104100065373294346723964356736436709934871741161328286944150242733445542228293036404657556168844723521815836689387184856871091025434896710605688594847400051686361372872763001355411405782508020591933546964183881743133374126947753
n = 19850163314401552502654477751795889962324360064924594948231168092741951675262933573691070993863763290962945190372400262526595224437463969238332927564085237271719298626877917792595603744433881409963046292095205686879015029586659384866719514948181682427744555313382838805740723664050846950001916332631397606277703888492927635867870538709596993987439225247816137975156657119509372023083507772730332482775258444611462771095896380644997011341265021719189098262072756342069189262188127428079017418048118345180074280858160934483114966968365184788420091050939327341754449300121493187658865378182447547202838325648863844192743
c = 13913396366755010607043477552577268277928241319101215381662331498046080625902831202486646020767568921881185124894960242867254162927605416228460108399087406989258037017639619195506711090012877454131383568832750606102901110782045529267940504471322847364808094790662696785470594892244716137203781890284216874035486302506042263453255580475380742959201314003788553692977914357996982118328587119124144181290753389394149235381045389696841471483947310663329993873046123134587149661347999774958105091103806375702387084149309542351541021140111048408248121408401601979108510758891595550054699719801708646232427198902271953673874
e = 28

d泄露分解

若给定$d$，我们有$k=e∗d−1$，通过定义我们可以知道，$k$是$\varphi(n)$的倍数，我们知道$\varphi(n)$是偶数，同时有$k=2^t∗r$，（$r$是奇数）。对于$g\in Z_n^*$，有$g^k≡1\mod n$，并且$g^\frac{k}{2}$是模$n$的平方根，由中国剩余定理我们知道，对于$n=p*q$，这样的数有4个。其中，有两个是正负一，还有两个是$\pm x$，其中$x≡1\mod p$,$x≡−1\mod q$。使用后两者，我们就可以通过计算$gcd(x-1, n)$得到$n$的分解。有简单的论证表明，如果$g$是在$Z^∗_n$中随机取，那么序列$g^\frac{k}{2},g^\frac{k}{4},\cdots,g^\frac{k}{2^t}\mod n$的其中一个是平方根且能分解$n$的概率至少为$1/2$。序列中的所有元素都可以在$O(log_2n)^3)$的时间内有效计算。

将$e,d$写成如下形式

$ed-1=k\varphi(n)=2^st\\$

只要一半的$a\in Z_n$，存在一个$i\in[1,s]$满足

$r=a^{2^{i-1}t}\ne1\mod n\\ a^{2^it}\equiv 1\mod n$

那么此时，$p=gcd(r-1,n)$是$n$的一个因子

dp leak

dp与$\varphi(p)$位数接近，e与k位数接近，在e比较小时，可以用爆破法求解

$dp\equiv d\mod\varphi(p)\\ e*dp-1=k\varphi(p)\\ \varphi(p)=\frac{e*dp-1}{k}\\ gcd(\frac{e*dp-1}{k}+1,n)\ne1$

from Crypto.Util.number import *
def dp_leak(e, dp, n, c):
    for k in range(1, e):
        phi = (e*dp - 1)//k
        p = GCD(phi + 1, n)
        if p != 1:
            q = n//p
            phi = (p-1)*(q-1)
            d = inverse(e, phi)
            m = pow(c, d, n)
            return m

coppersmith下的dp泄露

看Tover博客时看到的一种方法，在e大时无法如上述进行爆破

$ed_p-1=k(p-1)\\ ed_p-1+k=kp$

令$A=ed_p-1$，令$x=k$，就是解同余方程：

$A+x\equiv0\mod p$

求出小量$k$，$k$无论如何都小于$e$，如果$e$过大那么$k$就不是一个相对小量，那么就无法适用，求出$k$后自然能恢复$p$，进而得到所有的量，可求解范围大概是$e<n^{1/4}$

#!/usr/bin/env sage

e  = 
n  = 
dp = 
c  = 

PR.<x> = PolynomialRing(Zmod(n))
f = e * dp + x - 1
f = f.monic()
roots = f.small_roots(X=e, beta=0.47)
assert roots

k = Integer(roots[0])
assert (e * dp + k - 1) % k == 0

p = (e * dp + k - 1) // k
assert n % p == 0
q = n // p

phi = (p-1) * (q-1)
d = Integer(e).inverse_mod(phi)
m = Integer(pow(c, d, n))
flag = bytes.fromhex(m.hex())
print(flag.decode())

例题

from secret import flag
import libnum

bits = 1024
p, q = [random_prime(2^bits) for _ in range(2)]
n = p * q
e = 2*10^76-3
d = e.inverse_mod((p-1) * (q-1))
dp = d % (p-1)

m = libnum.s2n(flag)
c = pow(m, e, n)

print('e  = %d' % e)
print('n  = %d' % n)
print('dp = %d' % dp)

print('c  = %d' % c)

'''
e  = 19999999999999999999999999999999999999999999999999999999999999999999999999997
n  = 7195506839435218889565105541674965483194164483027741709706696451513641438345177472634371310250998546706062462270851552911697354605048972081656931006641878545036542923897114647393564522132057589249800431430995780074871171268958056358251827104531889348948541240686274977093185746573748206617663459128090693743840574459752890533065398493485714768878646999590143805843490432318539260302521682823958290340460403361801534822098048095280034600065200137857346827560676300256938953222718633375808719441534702981763523406056651752914141143665893462943582116716812913462656214604870428310720751101481210148746546806273965485289
dp = 34961801811050613471700883525108632060492526395401334090302835931304663757529660746363964830407055340550990256271716811099606849841913560556222756478612800702209651907866303152581107449312861896692310607989826809665245295483724533775337076019316812377921373194504440845718347150919782506437242366281376701299
c  = 3014636373048664939954772778404195986026862165799593915685719641505606570670923436003664110094703916031096486273947905494103538805486521321522443488182065845367347589071783679908494724693530639371358965655992560909299314626568439587755874253430614726720724608456333450258184012429367293386944954388615812902809362326474915645899324083994448117282677622943580354006160302366855350193039875335543211982510928721395526768129547143054319585071252781483346116972611571317425047748862917945459911485505200762492537496489429730213393936533514665994680707861503489288913062785427211743828345144957201996243444547648085230048
'''

dp问题通解

无论e范围多大，任意取一个数t，求出$t^{e*dp}-t$，再将结果与$n$取公因数，代入上题确实可以解决

$dp\equiv d\mod p-1\\ e*dp-1=k(p-1)\\ 2^{e*dp}\equiv 2^{k(p-1)+1}\equiv2^1\mod p\\ 2^{e*dp}-2^{1}=kp$

from Crypto.Util.number import *
def common_dp(e, dp, n, c):
    t = pow(2, e*dp, n)
    p = GCD(t - 2, n)
    q = n//p
    phi = (p-1)*(q-1)
    m = pow(c, inverse(e, phi), n)
    return m

e=2(Rabin算法)

算法的攻击条件：$e=2$

$m^2\equiv c(mod\ n)\Rightarrow m^2\equiv c(mod\ p),m^2\equiv c(mod\ q)$

即$c$是二次剩余，二次剩余具有一些特别的性质，当$c$是二次非剩余时同余为$-1$：

$c^{\frac{p-1}{2}}\equiv1(mod\ p)$

因为$p\equiv3(mod\ 4)$和$q\equiv3(mod\ 4)$，

$m_p^2\equiv c(mod\ p)=c*c^{\frac{p-1}{2}}(mod\ p)=c^{\frac{p+1}{2}}(mod\ p)\\ m_p\equiv\sqrt{c}(mod\ p)=c^\frac{p+1}{4}(mod\ p)\\ m_q=\sqrt{c}(mod\ q)=c^{\frac{q+1}{4}}(mod\ q)$

利用拓展欧几里得算出$y_p和y_q$：

$y_pp+y_qq=1\\$

存在方程

$m\equiv\begin{cases} m_q(mod\ q) \\ m_p(mod\ p) \end{cases}$

中国剩余定理求和公式：

$m=\sum_{i=1}^j n_{i}^{*} (n_i^{*})^{-1} a_{i}\\ m_p^2\equiv c^{\frac{p+1}{2}}(mod\ p)\Rightarrow m_p\equiv\pm c^\frac{p+1}{4}(mod\ p)\\ m_q\equiv\pm c^{\frac{q+1}{4}}(mod\ q)$

根据中国剩余定理，解出四个明文：

$a=(y_ppm_q+y_qqm_p)(mod\ n)\\b=n-a\\c=(y_ppm_q-y_qqm_p)(mod\ n)\\d=n-c$

import gmpy2
def rabin(cipher, p, q):
    n = p * q
    mp = pow(cipher, (p + 1) // 4, p)
    mq = pow(cipher, (q + 1) // 4, q)
    yp = gmpy2.invert(p, q)
    yq = gmpy2.invert(q, p)
    a = (yp * p * mq + yq * q * mp) % n
    b = n - a
    c = (yp * p * mq - yq * q * mp) % n
    d = n - c
    return (a, b, c, d)

p-1光滑

光滑数($smooth\ number$)，是一个可以因数分解为小质数乘积的正整数。

如果整数的所有素因子都不大于$B$，我们称这个整数是$B-Smooth$数。

$p-1=q_1*q_2*q_3*\cdots*q_m,且q_1,q_2,q_3,\cdots,q_m\le B$

这时候可以通过从小素数开始依次枚举，直到将其完全分解，每一个$q_n$项可以表示质数的任意幂次方。

关于算法原理可以参考这里

跟$p-1$光滑相关还有$p+1$光滑，在原理上$p+1$理解难度更大，两者推荐直接使用脚本。python中有直接的库可以调用，看到python -m知道在cmd中运行。

python -m primefac -vs -m=p-1 xxxxxx
python -m primefac -vs -m=p+1 xxxxxx

p-1光滑分解脚本如下，p+1脚本找到的不太能用，就不给出了

from gmpy2 import gcd
a = 2
n = 2
N = 
while True:
    a = pow(a, n, N)
    res = gcd(a - 1, N)
    if res != 1 and res != N:
        q = N // res
        print("p =", res)
        print("q =", q)
        break
    n += 1

例题:

from random import choice
from Crypto.Util.number import isPrime, sieve_base as primes
from flag import flag


def myPrime(bits):
    while True:
        n = 2
        while n.bit_length() < bits:
            n *= choice(primes)
        if isPrime(n + 1):
            return n + 1

e = 0x10001
m = int.from_bytes(flag.encode(), 'big')
p = myPrime(2048)
q = getPrime(2048)
n = p * q
c = pow(m, e, n)

# n = 1224542620373200232525165378018470774334801515191193204875465445916504222883935188890019876845076388385357911730689547124547346252951158814249284724565588433721828377715469374541007756509231399263095022024229078845538543233785364809917406108015271780070196140628158427541531563472532516237632553353655535922926443707617182025475004547531104052989085765070550150028833424395972178427807901747932614235844448614858629761183210600428438018388051958214596857405813088470933109693499438012040822262549119751099008671892966082341548512112435591881692782766559736840448702039918465573051130405935280702181505538733234675792472428666968900055706926735800561218167237812066851519973807203332801575980055838563085817664973968944323258406789203078387708964307931318918136664885818917720073433998810127482159223895026085726623747340692196977140382318293090736558135980651252533606603312148824142669800602887109353065489282386215179238458743567166284295855288783740314247952124965482197632971993708775190564519250754150756867653033527903848903210074426177258586450311109023467944412194124015505951966140443860862968311560843608415723549525497729679097936310538451467530605937684408079363677707513923579164067808729408365886209340192468399685190639
# c = 145742860621666495489510776707734134231023214235535481878205099324276369445463746101469487674333600296204530932386373415987357363515200117271393133347844479863240936801112306080456942844796779477817786176831015954410967693647534326733641573842953783193563678040093734579772976410574013857063137696465850300484753282472377882118892522844694078667622111244886303620349388556315704648609353412177123230438077637042880490566244740468503369707900343076369151796123461132932226563486870411965536062339169788331659119981901553536009275158600580698576110294775989992794065611215170351808698605911258789407992833170968332058255364527244293283228694886707241979238145252395651417561576433516407782575454294499521347378058366557950770592472271985004818847838711060048422015207674862177145761946560579360220239667890707135827136815780729363013864130107808776517514214310689477005999830284272130148939734935547341627208913181919190392205389452185597444280635342938046191904062547803917870268485346888653569349729643793041018550170090471310374856687407102762116819004790791936814214507908374380597027347007448114684844276041116955473180015221164545212550832233007714133699817366745648092776901013502840540012912660742166994968977400188176557657864

综合练习题

from Crypto.Util.number import isPrime, bytes_to_long, sieve_base
from random import choice
from secret import flag
 
m=bytes_to_long(flag)
def uniPrime(bits):
    while True:
        n = 2
        while n.bit_length() < bits:
            n *= choice(sieve_base)
        if isPrime(n + 1):
            return n + 1
 
 
p=uniPrime(512)
q=uniPrime(512)
n=p*q
e= 196608
c=pow(m,e,n)
 
print("n=",n)
print("c=",c)
 
'''
n= 3326716005321175474866311915397401254111950808705576293932345690533263108414883877530294339294274914837424580618375346509555627578734883357652996005817766370804842161603027636393776079113035745495508839749006773483720698066943577445977551268093247748313691392265332970992500440422951173889419377779135952537088733
c= 2709336316075650177079376244796188132561250459751152184677022745551914544884517324887652368450635995644019212878543745475885906864265559139379903049221765159852922264140740839538366147411533242116915892792672736321879694956051586399594206293685750573633107354109784921229088063124404073840557026747056910514218246
'''

Schmidt-Samoa

$N=p^2\cdot q$，$d=invert(N, \varphi(pq))$,公钥对$(N,d)$
加密

$C=m^N(mod\ N)$

解密

$m=C^d(mod\ pq)$

$d$和$N$在模$\varphi(pq)$下互为逆元，故$Nd\equiv1(mod\ \varphi(pq))$
任取一个数字$a$

$a^{Nd}\equiv a\mod pq\Rightarrow\quad a^{Nd}-a=k\cdot pq\\ pq=gcd(a^{nd}-a, N)$

exp

from Crypto.Util.number import *
from gmpy2 import iroot

def solve(N, d, c):
    a = pow(2, N*d, N)
    n = GCD(a - 2, N)
    m = pow(c, d, n)
    return m

例题:

from secret import flag
from Crypto.Util.number import *

p = getPrime(1024)
q = getPrime(1024)

N = p*p*q

d= inverse(N, (p-1)*(q-1)//GCD(p-1, q-1))

m = bytes_to_long(flag)

c = pow(m, N, N)

print('c =', c)
print('N =', N)
print('d =', d)

# c = 1653396627113549535760516503668455111392369905404419847336187180051939350514408518095369852411718553340156505246372037811032919080426885042549723125598742783778413642221563616358386699697645814225855089454045984443096447166740882693228043505960011332616740785976743150624114653594631779427044055729185392854961786323215146318588164139423925400772680226861699990332420246447180631417523181196631188540323779487858453719444807515638025771586275969579201806909799448813112034867089866513864971414742370516244653259347267231436131850871346106316007958256749016599758599549180907260093080500469394473142003147643172770078092713912200110043214435078277125844112816260967490086038358669788006182833272351526796228536135638071670829206746835346784997437044707950580087067666459222916040902038574157577881880027391425763503693184264104932693985833980182986816664377018507487697769866530103927375926578569947076633923873193100147751463
# N = 1768427447158131856514034889456397424027937796617829756303525705316152314769129050888899742667986532346611229157207778487065194513722005516611969754197481310330149721054855689646133721600838194741123290410384315980339516947257172981002480414254023253269098539962527834174781356657779988761754582343096332391763560921491414520707112852896782970123018263505426447126195645371941116395659369152654368118569516482251442513192892626222576419747048343942947570016045016127917578272819812760632788343321742583353340158009324794626006731057267603803701663256706597904789047060978427573361035171008822467120148227698893238773305320215769410594974360573727150122036666987718934166622785421464647946084162895084248352643721808444370307254417501852264572985908550839933862563001186477021313236113690793843893640190378131373214104044465633483953616402680853776480712599669132572907096151664916118185486737463253559093537311036517461749439
# d = 20650646933118544225095544552373007455928574480175801658168105227037950105642248948645762488881219576174131624593293487325329703919313156659700002234392400636474610143032745113473842675857323774566945229148664969659797779146488402588937762391470971617163496433008501858907585683428652637958844902909796849080799141999490231877378863244093900363251415972834146031490928923962271054053278056347181254936750536280638321211545167520935870220829786490686826062142415755063724639110568511969041175019898031990455911525941036727091961083201123910761290998968240338217895275414072475701909497518616112236380389851984377079

Wilson定理

设$p$是素数，$r_1,\cdots,r_{p-1}$是模$p$的既约剩余系，我们有：

$\prod_{r\ mod\ p}r\equiv r_{1}…r_{p-1} \equiv-1(mod\ p)$

在模质数$p$下，$gcd(r_i,p)=1$，每个元素都有逆元，$r_1$和$r_{p-1}$的逆元都是其本身，其他元素逆元两两相消，特别地，有：

$(p-1)!\equiv-1(mod\ p)\\ -1(mod\ p)=(p-1)(mod\ p)\\ (p-2)!\equiv1(mod\ p)$

例题(RoarCTF 2019)

import sympy
import random

def myGetPrime():
    A= getPrime(513)
    print(A)
    B=A-random.randint(1e3,1e5)
    print(B)
    return sympy.nextPrime((B!)%A)
p=myGetPrime()
#A1=21856963452461630437348278434191434000066076750419027493852463513469865262064340836613831066602300959772632397773487317560339056658299954464169264467234407
#B1=21856963452461630437348278434191434000066076750419027493852463513469865262064340836613831066602300959772632397773487317560339056658299954464169264467140596

q=myGetPrime()
#A2=16466113115839228119767887899308820025749260933863446888224167169857612178664139545726340867406790754560227516013796269941438076818194617030304851858418927
#B2=16466113115839228119767887899308820025749260933863446888224167169857612178664139545726340867406790754560227516013796269941438076818194617030304851858351026

r=myGetPrime()

n=p*q*r
#n=85492663786275292159831603391083876175149354309327673008716627650718160585639723100793347534649628330416631255660901307533909900431413447524262332232659153047067908693481947121069070451562822417357656432171870951184673132554213690123308042697361969986360375060954702920656364144154145812838558365334172935931441424096270206140691814662318562696925767991937369782627908408239087358033165410020690152067715711112732252038588432896758405898709010342467882264362733
c=pow(flag,e,n)
#e=0x1001
#c=75700883021669577739329316795450706204502635802310731477156998834710820770245219468703245302009998932067080383977560299708060476222089630209972629755965140317526034680452483360917378812244365884527186056341888615564335560765053550155758362271622330017433403027261127561225585912484777829588501213961110690451987625502701331485141639684356427316905122995759825241133872734362716041819819948645662803292418802204430874521342108413623635150475963121220095236776428
#so,what is the flag?

$A$是$513bits$的大整数，而$B$与$A$相差很小，要求出$B!$，这么大的整数，直接算明显不合理，威尔逊定理刚好可以解决这个问题。首先$(A-2)!\equiv1(mod\ A)$成立，那么从这开始反推

$(A-2)!\equiv1(mod\ A)\\ (A-3)!\equiv(A-2)^{-1}(mod\ A)\\ (A-4)!\equiv(A-3)^{-1}(A-2)^{-1}(mod\ A)\\ \cdots$

直接推到$B!$为止，那么再使用$nextprime$便是一个因子。

def willion(a, b):
    o = 1
    for i in range(b + 1, a - 1):
        o = o * inverse(i, a) % a
    return sympy.nextprime(o)

不给e,但d非常大(Wiener’s attack)

连分数

一个有理数可以化成如下形式

$\large x_0+\frac{1}{x_1+\frac{1}{x_2+\frac{1}{\cdots+\frac{1}{x_n}}}}$

该形式即为连分数的表达形式，通常简记为$<x_0,x_1,x_2,\cdots,x_n>$或者$[x_0,x_1,x_2,\cdots,x_n]$因为上式数位有限，所以是有限简单连分数。形如下方的式子则称为无限简单连分数

$\large x_0+\frac{1}{x_1+\frac{1}{x_2+\frac{1}{\cdots}}}$

无限简单连分数通常简记为$<x_0,x_1,x_2,\cdots>$或者$[x_0,x_1,x_2,\cdots]$。

定理1

任何一个有理数可以写成有限连分数形式

有理数必然可以用分数表示，且分子分母都是固定的数字，根据辗转相除法，必然只执行有限次除法取余数。例如取有理数$\frac{b_0}{r_0}$

\begin{align} &b_0=a_0r_0+r_1\\ &r_0=a_1r_1+r_2\\ &r_1=a_2r_2+r_3\\ &\cdots\\ \end{align}

\begin{align} \frac{b_0}{r_0}&=a_0+\frac{r_1}{r_0}=a_0+\frac{1}{\frac{r_0}{r_1}}\\ &=a_0+\frac{1}{a_1+\frac{r_1}{r_2}}\\ &=a_0+\frac{1}{a_1+\frac{1}{a_2+\frac{r_2}{r_3}}}\\ &\cdots\cdots\\ &=[a_0,a_1,a_2,\cdots,a_n] \end{align}

渐进分数

对于连分数的简写形式，我们也可以将其还原成分数。将其代入分数形式的定义式中，便能求出其值。

对于$[a_0,a_1,a_2,\cdots,a_n]$这个简写形式。设$0\leq k\leq n$，我们把有限分数$[a_0,a_1,\cdots,a_k]$称为有限连分数的第$k$个渐进分数，把$a_k$称为它的第$k$个部分商。取的$k$越大，还原出的分数越逼近原值。

$\large |\frac{b_0}{r_0}-a_0|>|\frac{b_0}{r_0}-(a_0+\frac{1}{a_1})|\\ \large |\frac{b_0}{r_0}-(a_0+\frac{1}{a_1})|>|\frac{b_0}{r_0}-(a_0+\frac{1}{a_1+\frac{1}{a_2}})|$

对于有限连分数，当数目取到极值时，便能还原成正确分数，而对于无限连分数，便是一个不断接近真理的过程。

对于无限连分数$[a_0,a_1,a_2,\cdots]$一定是收敛的，也就是说，设$r^{(n)}=[a_0,\cdots,a_n]$是它的第$n$个渐进分数，那么一定存在极限

$\lim_{n\rightarrow\infin}r^{(n)}=\theta\\$

此外还有

$r^{(0)}<r^{(2)}<\cdots<r^{(2t)}<\theta<\cdots<r^{(2s-1)}<\cdots<r^{(3)}<r^{(1)}$

以及$\theta$一定是无理数。

我们拿$\sqrt{11}$来举例

\begin{align} 3.3166247903554\approx\sqrt{11}&=3+\sqrt{11}-3 &&=3<\sqrt{11}\\ &=3+\frac{1}{3+\frac{\sqrt{11}-3}{2}} &&\approx3.3\mathbf{3}33333333333333>\sqrt{11}\\ &=3+\frac{1}{3+\frac{1}{6+(\sqrt{11}-3)}} &&\approx3.31\mathbf{5}7894736842106<\sqrt{11}\\ &=3+\frac{1}{3+\frac{1}{6+\frac{1}{3+\frac{\sqrt{11}-3}{2}}}} &&\approx3.3166\mathbf{6}6666666667>\sqrt{11}\\ &\cdots &&\approx3.31662\mathbf{2}691292876<\sqrt{11}\\ &\cdots &&\approx3.316624\mathbf{8}95572264>\sqrt{11} \end{align}

定理2前言

$[x_0]=\frac{x_0}{1},\quad[x_0,x_1]=x_0+\frac{1}{x_1}=\frac{x_0x_1+1}{x_1}\\ [x_0,x_1,x_2]=[x_0,x_1+\frac{1}{x_2}]=\frac{x_0(x_1+\frac{1}{x_2})+1}{x_1+\frac{1}{x_2}}=\frac{(x_0x_1+1)x_2+x_0}{x_1x_2+1}\\ [x_0,x_1,x_2,x_3]=[x_0,x_1,x_2+\frac{1}{x_3}]=\frac{((x_0x_1+1)x_2+x_0)x_3+(x_0x_1+1)}{(x_1x_2+1)x_3+x_1}$

如果将$x_0,x_1,\cdots$看做实变数，那么

$[x_0,\cdots,x_{n-1},x_n]=\frac{P_n}{Q_n},\quad n\geq0$

定理2

设$x_0,x_1,x_2,\cdots$是无穷实数列，$x_j>0,j\geq1$；再设

$P_{-2}=0,\quad P_{-1}=1,\quad Q_{-2}=1,\quad Q_{-1}=0$

那么可以得出
$\begin{cases} P_n=x_nP_{n-1}+P_{n-2}\\ Q_n=x_nQ_{n-1}+Q_{n-2} \end{cases}$(形式在上面展开时可以通过规律发现)

根据$P_n,Q_n$的递推关系可以得出(使用数学归纳法)

$[x_0,\cdots,x_n]=\frac{P_n}{Q_n},\quad n\geq0\\ P_nQ_{n-1}-P_{n-1}Q_n=(-1)^{n+1},\quad n\geq-1\\ P_nQ_{n-2}-P_{n-2}Q_n=(-1)^{n}x_n,\quad n\ge0\\ [x_0,\cdots,x_{n-1},x_n]-[x_0,\cdots,x_{n-1}]=(-1)^{n+1}(Q_nQ_{n-1})^{-1},\quad n\geq1\\ [x_0,\cdots,x_{n-2},x_{n-1},x_n]-[x_0,\cdots,x_{n-2}]=(-1)^nx_n(Q_nQ_{n-2})^{-1},\quad n\ge2$

证明(数学归纳法)

$[x_0,\cdots,x_{k-1},x_k,x_{k+1}]=[x_0,\cdots,x_{k-1},x_k+\frac{1}{x_{k+1}}]\\$

\begin{align} [x_0,\cdots,x_{k-1},x_k,x_{k+1}]&=\frac{(x_k+\frac{1}{x_{k+1}})P_{k-1}+P_{k-2}}{(x_k+\frac{1}{x_{k+1}})Q_{k-1}+Q_{k-2}}\\ &=\frac{x_{k+1}(x_kP_{k-1}+P_{k-2})+P_{k-1}}{x_{k+1}(x_kQ_{k-1}+Q_{k-2})+Q_{k-1}}\\ &=\frac{x_{k+1}P_k+P_{k-1}}{x_{k+1}Q_k+Q_{k-1}}=\frac{P_{k+1}}{Q_{k+1}} \end{align}

\begin{align} P_kQ_{k-1}-P_{k-1}Q_k&=(x_kP_{k-1}+P_{k-2})Q_{k-1}-P_{k-1}(x_kQ_{k-1}+Q_{k-2})\\ &=-(P_{k-1}Q_{k-2}-P_{k-2}Q_{k-1})\\ &=-(-1)^{k}=(-1)^{k+1}\\ P_kQ_{k-2}-P_{k-2}Q_k&=(x_kP_{k-1}+P_{k-2})Q_{k-2}-P_{k-2}(x_kQ_{k-1}+Q_{k-2})\\ &=x_k(P_{k-1}Q_{k-2}-P_{k-2}Q_{k-1})\\ &=x_k(-1)^{k} \end{align}

$P_nQ_{n-1}-P_{n-1}Q_n=(-1)^{n+1}\Rightarrow \frac{P_n}{Q_n}-\frac{P_{n-1}}{Q_{n-1}}=(-1)^{n+1}(Q_nQ_{n-1})^{-1}\\ P_nQ_{n-2}-P_{n-2}Q_n=(-1)^nx_n\Rightarrow\frac{P_n}{Q_{n}}-\frac{P_{n-2}}{Q_{n-2}}=(-1)^nx_n(Q_nQ_{n-2})^{-1}$

定理3

设$[a_0,a_1,a_2,\cdots]$是无限简单连分数，记极限为

$\xi_n=[a_n,a_{n+1},\cdots],\quad n\geq0$

那么有

$a_n=[\xi_n]$(向下取整)，$\quad\xi_{n+1}=\{\xi_n\}^{-1}$(特定运算)
$[a_0,a_1,a_2,\cdots]=[a_0,\cdots,a_n,\xi_{n+1}],\quad n\geq0$
推导特定运算：$\xi_n=a_n+\frac{1}{\xi_{n+1}}\Rightarrow\xi_{n+1}=(\xi_n-a_n)^{-1}=\{\xi_n\}^{-1}$

定理4

$\frac{1}{Q_n(Q_n+Q_{n+1})}<|\xi_0-r^{(n)}|=|\xi_0-\frac{P_n}{Q_n}|<\frac{1}{Q_nQ_{n+1}}$

证明

$\xi_0=[a_0,\cdots,a_n,\xi_{n+1}]=\frac{P_n\xi_{n+1}+P_{n-1}}{Q_n\xi_{n+1}+Q_{n-1}}$

\begin{align} \xi_0-r^{(n)}&=\frac{Q_nP_{n-1}-Q_{n-1}P_n}{Q_n(Q_n\xi_{n+1}+Q_{n-1})}\\ &=\frac{(-1)^{n}}{Q_n(Q_n\xi_{n+1}+Q_{n-1})} \end{align}

由于$a_{n+1}<\xi_{n+1}<a_{n+1}+1$

$\begin{aligned} Q_{n+1}&=a_{n+1}Q_n+Q_{n-1}<Q_n\xi_{n+1}+Q_{n-1}\\ &<a_{n+1}Q_n+Q_{n-1}+Q_n=Q_{n+1}+Q_n \end{aligned}$

$\frac{1}{Q_n(Q_n+Q_{n+1})}<\frac{1}{Q_n(Q_n\xi_{n+1}+Q_{n-1})}<\frac{1}{Q_nQ_{n+1}}$，式子成立

定理5

设$\xi_0$是无理数，若有有理分数$\frac{a}{b},b\ge1$，使得

$|\xi_0-\frac{a}{b}|<\frac{1}{2b^2}$

那么$\frac{a}{b}$一定是$\xi_0$的某个渐进分数。

attack algorithm

$N=pq$，其中$q<p<2q$，若$d<\frac{1}{3}N^{\frac{1}{4}}$时，给定公钥$(N,e)$，且$ed\equiv 1(mod\ \varphi(N))$，那么可以有效得到$d$。

$ed\equiv1(mod\ \varphi(N)),N-\varphi(N)=p+q-1\\ ed-k\varphi(N)=1,p+q<3\sqrt{N}\\ e<\varphi(N)\Rightarrow k<d<\frac{1}{3}N^{\frac{1}{4}}$

\begin{align} |\frac{e}{N}-\frac{k}{d}|&=|\frac{ed-kN-k\varphi(N)+k\varphi(N)}{Nd}|\\ &=|\frac{1-k(N-\varphi(N))}{Nd}|<|\frac{3k\sqrt{N}}{Nd}|=\frac{3k}{d\sqrt{N}}\\ &<\frac{3d}{d(3d)^2}=\frac{1}{3d^2}<\frac{1}{2d^2} \end{align}

$|\frac{e}{N}-\frac{k}{d}|<\frac{1}{2d^2}$

根据定理6，$\frac{k}{d}$一定是$\frac{e}{N}$的某个渐进分数

fra = continued_fraction(e/N)
for i in range(1, len(fra)):
    k = fra.numerator(i)
    d = fra.denominator(i)
    if (e*d-1)%k == 0:
            print(k, d)

对于连分数展开成列表形式，使用$Sagemath$中的集成函数$continued\_fraction()$。既然$\frac{k}{d}$一定是$\frac{e}{N}$的某个渐进分数，那么在连分数展开长度内，对每一个渐进分数赋值给$\frac{k}{d}$，直到出现$ed-1=kn$这样形式出现时，便是渐进出正确的$k,d$。将式子改写也便是$(ed-1)\%k = 0$。

例题

N = 0xb8af3d3afb893a602de4afe2a29d7615075d1e570f8bad8ebbe9b5b9076594cf06b6e7b30905b6420e950043380ea746f0a14dae34469aa723e946e484a58bcd92d1039105871ffd63ffe64534b7d7f8d84b4a569723f7a833e6daf5e182d658655f739a4e37bd9f4a44aff6ca0255cda5313c3048f56eed5b21dc8d88bf5a8f8379eac83d8523e484fa6ae8dbcb239e65d3777829a6903d779cd2498b255fcf275e5f49471f35992435ee7cade98c8e82a8beb5ce1749349caa16759afc4e799edb12d299374d748a9e3c82e1cc983cdf9daec0a2739dadcc0982c1e7e492139cbff18c5d44529407edfd8e75743d2f51ce2b58573fea6fbd4fe25154b9964d
e = 0x9ab58dbc8049b574c361573955f08ea69f97ecf37400f9626d8f5ac55ca087165ce5e1f459ef6fa5f158cc8e75cb400a7473e89dd38922ead221b33bc33d6d716fb0e4e127b0fc18a197daf856a7062b49fba7a86e3a138956af04f481b7a7d481994aeebc2672e500f3f6d8c581268c2cfad4845158f79c2ef28f242f4fa8f6e573b8723a752d96169c9d885ada59cdeb6dbe932de86a019a7e8fc8aeb07748cfb272bd36d94fe83351252187c2e0bc58bb7a0a0af154b63397e6c68af4314601e29b07caed301b6831cf34caa579eb42a8c8bf69898d04b495174b5d7de0f20cf2b8fc55ed35c6ad157d3e7009f16d6b61786ee40583850e67af13e9d25be3
c = 0x3f984ff5244f1836ed69361f29905ca1ae6b3dcf249133c398d7762f5e277919174694293989144c9d25e940d2f66058b2289c75d1b8d0729f9a7c4564404a5fd4313675f85f31b47156068878e236c5635156b0fa21e24346c2041ae42423078577a1413f41375a4d49296ab17910ae214b45155c4570f95ca874ccae9fa80433a1ab453cbb28d780c2f1f4dc7071c93aff3924d76c5b4068a0371dff82531313f281a8acadaa2bd5078d3ddcefcb981f37ff9b8b14c7d9bf1accffe7857160982a2c7d9ee01d3e82265eec9c7401ecc7f02581fd0d912684f42d1b71df87a1ca51515aab4e58fab4da96e154ea6cdfb573a71d81b2ea4a080a1066e1bc3474

Coppersmith攻击

Coppersmith proved that any root $x_0$ with $|x_0|<N^{1/\delta}$ can be found in polynomial time, where $\delta=deg\ p$. The technique consists in building a lattice that contains the solutions of the modular polynomial equation; all small solutions are shown be belong to an hyperplane of the lattice; an equation of this hyperplane is obtained by considering that last vector of an $LLL-reduced$ basis; this gives a polynomial $h(x)$ such that $h(x_0)=0$ over the integers, from which one can recover $x_0$. The method can be extended to handle multivariate modular polynomial equations, but the extension is heuristic only.

参考

有一个e阶多项式$f$，那么可以

• 在模$n$意义下，快速求出$n^{\frac{1}{e}}$以内的根
• 给定$\beta$，快速求出模某个$b$意义下较小根，其中$b\ge n^{\beta}$，是$n$的因数

p高位

coppersmith可以解决多项式在模$n$的某个因数下的根

$p_h+x=0\quad(mod\ sth\ divides\ n)$

# Sage
p>>128<<128类型
t = p >> 128 << 128
kbits = 128
R.<x> = PolynomialRing(Zmod(n))
f = x + t
x0 = f.small_roots(X=2^kbits, beta=0.4)[0]
p = p + int(x0)

$f(x)=x+t,\quad f(x_0)=x_0+t\equiv0(mod\ n)\\ x_0+t=k*n$

用$X=2^{kbits}$将求解范围上限界定，对coppersmith攻击的范围稍微测试了一下，在X=2^kbits, beta=0.4条件下kbits为227时，即已知高位的285，大概率是能恢复出p，测试了很多组只有一次没有恢复出来，而如果高位仅知道284位，大概率是恢复不出来，测试了十组没有一组成功。那么如果看到p >> 227 << 227这样的形式，只要里面的数字小于等于227，大概率便能直接使用此方法恢复p，下面是测试时写的代码。

def check(kbits):
    p = getPrime(512)
    q = getPrime(512)
    n = p * q
    leak = p >> kbits
	
    temp = leak << kbits
    R.<x> = PolynomialRing(Zmod(n))
    f = x + temp
    res = f.small_roots(X=2^kbits, beta=0.4)
    print(res, end='    ')
    print(res[0]+temp==p)

small_roots默认epsilon为0.05，将epsilon调小可以在已知更少位数时算出答案，但所需时间也越长，经个人测试在kbits为248时大概率算出，即已知高位264位，kbits到249时大概率算不出，起码测试时没有成功算出的例子，这样通过修改epsilon参数使得coppersmith攻击范围进一步扩大了21位。

def check(kbits):
    p = getPrime(512)
    q = getPrime(512)
    n = p * q
    leak = p >> kbits
	
    temp = leak << kbits
    R.<x> = PolynomialRing(Zmod(n))
    f = x + temp
    res = f.small_roots(X=2^kbits, beta=0.4, epsilon=0.01)
    print(res, end='    ')
    print(res[0]+temp==p)

[2024-SICTF-Round3]铜匠

from Crypto.Util.number import *
from enc import flag

def Decimal_conversion(num):
    if num == 0:
        return '0'
    digits = []
    while num:
        digits.append(str(num % 5))
        num //= 5
    return ''.join(reversed(digits))

m = bytes_to_long(flag)
p = getPrime(512)
q = getPrime(512)
e = 65537
n = p*q
c = pow(m,e,n)
print(f"leak = {Decimal_conversion(p)[:112]}")
print(f"n = {n}")
print(f"e = {e}")
print(f"c = {c}")
'''
leak = 2011133132443111302000224204142244403203442000141102312242343143241244243020003333022112141220422134444214010012
n = 85988668134257353631742597258304937106964673395852009846703777410474172989069717247424903079500594820235304351355706519069516847244761609583338251489134035212061654870087550317540291994559481862615812258493738064606592165529948648774081655902831715928483206013332330998262897765489820121129058926463847702821
e = 65537
c = 64708526479058278743788046708923650158905888858865427385501446781738669889375403360886995849554813207230509920789341593771929287415439407977283018525484281064769128358863513387658744063469874845446480637925790150835186431234289848506337341595817156444941964510251032210939739594241869190746437858135599624562
'''

首先对题目中代码进行分析，leak是转换成五进制的p的高位，用getPrime(512)生成的五进制数看刚好都是221位，知道五进制的高位，后面补零，这样不就能大致恢复出高位了，但因为补零存在误差，需要先自己生成几组数字观察，多组测试发现每一次最后两位大概率不同，这时是248位，是调整epsilon后的最大范围，那么还是不能直接求出。

from Crypto.Util.number import *

def Decimal_conversion(num):
    if num == 0:
        return '0'
    digits = []
    while num:
        digits.append(str(num % 5))
        num //= 5
    return ''.join(reversed(digits))

p = getPrime(512)
t = Decimal_conversion(p)
temp = t[:112]
print(p>>248)
print(int(temp + '0' * (221 - 112), 5)>>248)

"""
19695928862172181276359271675291516196118182510955101279750210725275456117222121
19695928862172181276359271675291516196118182510955101279750210725275456117222100
"""

可以对后两位进行爆破处理，总有一个答案是正确的高位数据，爆破也就只是100种情况。这样还是能算出答案的，因为有了高264位，可以恢复出p，之后便是常见的解密。

leak = '2011133132443111302000224204142244403203442000141102312242343143241244243020003333022112141220422134444214010012'
n = 85988668134257353631742597258304937106964673395852009846703777410474172989069717247424903079500594820235304351355706519069516847244761609583338251489134035212061654870087550317540291994559481862615812258493738064606592165529948648774081655902831715928483206013332330998262897765489820121129058926463847702821
e = 65537
c = 64708526479058278743788046708923650158905888858865427385501446781738669889375403360886995849554813207230509920789341593771929287415439407977283018525484281064769128358863513387658744063469874845446480637925790150835186431234289848506337341595817156444941964510251032210939739594241869190746437858135599624562

# print(int(leak + '0' * (221 - 112), 5)>>248)
# 26907901789126261702200127621264495244443424222236698342075135000543694885363317

def Decimal_conversion(num):
    if num == 0:
        return '0'
    digits = []
    while num:
        digits.append(str(num % 5))
        num //= 5
    return ''.join(reversed(digits))

def check(leak, de):
    if de < int(2):
        for i in range(10):
            check(leak + str(i), de + 1)
    elif de == 2:
        temp = Integer(int(leak))
        kbits = 512 - temp.nbits()

        temp = temp << kbits
        R.<x>= PolynomialRing(Zmod(n))
        f = x + temp
        res = f.small_roots(X=2 ^ kbits, beta=0.4, epsilon=0.01)
        if res:
            print(res[0] + temp)
# 12170789707638469557363249767228204966074269830454332436369564884472290413677820876733675261757160662428920138919536054003455470370040518528641001253487453
# SICTF{1d213ffc-d5dc-42d0-b206-e9a0c1a3cb69}

m高位

$m^e-c=(m_h+x)^e-c=0$

不止可以恢复$p$，甚至可以用来已知$m$高位恢复$m$，构造的式子$f$则需要变化。

from Crypto.Util.number import getPrime,bytes_to_long,long_to_bytes
from random import randint
from secret import flag

p = getPrime(1024)
q = getPrime(1024)
n = p*q
print(n)

m = bytes_to_long(long_to_bytes(randint(0,30))*208+flag)
assert(m.bit_length()==2044)
print((m>>315)<<315)
c = pow(m,3,n)
print(c)

exp

# Sage
from Crypto.Util.number import *
n = 14113948189208713011909396304970377626324044633561155020366406284451614054260708934598840781397326960921718892801653205159753091559901114082556464576477585198060530094478860626532455065960136263963965819002575418616768412539016154873800614138683106056209070597212668250136909436974469812231498651367459717175769611385545792201291192023843434476550550829737236225181770896867698281325858412643953550465132756142888893550007041167700300621499970661661288422834479368072744930285128061160879720771910458653611076539210357701565156322144818787619821653007453741709031635862923191561438148729294430924288173571196757351837
t = 1520800285708753284739523608878585974609134243280728660335545667177630830064371336150456537012842986526527904043383436211487979254140749228004148347597566264500276581990635110200009305900689510908049771218073767918907869112593870878204145615928290375086195098919355531430003571366638390993296583488184959318678321571278510231561645872308920917404996519309473979203661442792048291421574603018835698487725981963573816645574675640357569465990665689618997534740389987351864738104038598104713275375385003471306823348792559733332094774873827383320058176803218213042061965933143968710199376164960850951030741280074168795136
c = 6635663565033382363211849843446648120305449056573116171933923595209656581213410699649926913276685818674688954045817246263487415328838542489103709103428412175252447323358040041217431171817865818374522191881448865227314554997131690963910348820833080760482835650538394814181656599175839964284713498394589419605748581347163389157651739759144560719049281761889094518791244702056048080280278984031050608249265997808217512349309696532160108250480622956599732443714546043439089844571655280770141647694859907985919056009576606333143546094941635324929407538860140272562570973340199814409134962729885962133342668270226853146819
e = 3
R.<x> = PolynomialRing(Zmod(n))
f = (t+x) ^ e - c
m = f.small_roots(X=2^315, beta=0.4)[0]
print(long_to_bytes(t+m))

d低位

$ed\equiv 1\mod\varphi(n)\\ ed_l\equiv 1+k(n-p-q+1)\mod 2^{dbits}\\ ed_l*p\equiv p+k(np-p^2-n+p)\mod 2^{dbits}$

遍历k小于e，解这个一元方程可以得到p低位，在用低位恢复p

线性攻击

$\begin{cases} x^e-c_1=0\\(x+1)^e-c_2=0 \end{cases}$

$m$是两个方程的共解，那么$(x-m)$是这两多项式的公因式，用gcd即可求

二维Coppersmith

[2021 Zer0pts CTF]easy pseudo random

from Crypto.Util.number import*
from flag import flag

nbits = 256
p = random_prime(1 << nbits)
Fp = Zmod(p)
P.<v> = PolynomialRing(Fp)

b = randrange(p)
d = 2
F = v^2 + b

v0 = randrange(p)
v1 = F(v0)

k = ceil(nbits * (d / (d + 1)))
w0 = (v0 >> (nbits - k))
w1 = (v1 >> (nbits - k))

# encrypt
m = bytes_to_long(flag)
v = v1
for i in range(5):
    v = F(v)
    m ^^= int(v)

print(f"p = {p}")
print(f"b = {b}")
print(f"m = {m}")
print(f"w0 = {w0}")
print(f"w1 = {w1}")

已知$w_0,w_1$的高位，二元$Coppersmith$攻击恢复$w_0,w_1$。$Sage$中的$small\_roots$方法仅支持一元恢复，二元可以借用$github$上的拓展。关于拓展问题，需要将其放在sage可引用的目录下，输入__import__('os').system('pwd')可以看到当前目录，如果是wsl或者windows本身，可以直接将文件放进去

$v_0+b=v_1\\ (w_0<<85+x)^2+b=(w_1<<85+y)\\ f(x,y)=(w_0<<85+x)^2+b-(w_1<<85+y)$

$small\_roots$方法求出规定范围内使$f(x,y)=0$的合适解

# Sage
from Crypto.Util.number import *
p = 86160765871200393116432211865381287556448879131923154695356172713106176601077
b = 71198163834256441900788553646474983932569411761091772746766420811695841423780
m = 88219145192729480056743197897921789558305761774733086829638493717397473234815
w0 = 401052873479535541023317092941219339820731562526505
w1 = 994046339364774179650447057905749575131331863844814

load("coppersmith.sage") # amazing https://github.com/defund/coppersmith
x, y = PolynomialRing(GF(p), names='x, y').gens()

t = 85
pol = (w0 * 2**t + x)**2 + b - (w1 * 2**t + y)
rs = small_roots(pol, [2**t, 2**t])
x, y = rs[0]
v0 = int(w0*2**t + x)
v1 = int(w1*2**t + y)
assert (v0**2 + b) % p == v1

v = v1
for i in range(5):
    v = (v**2+b) % p
    m ^^= int(v)
long_to_bytes(m)

coppersmith.sage源文件

import itertools

def small_roots(f, bounds, m=1, d=None):
	if not d:
		d = f.degree()

	if isinstance(f, Polynomial):
		x, = polygens(f.base_ring(), f.variable_name(), 1)
		f = f(x)

	R = f.base_ring()
	N = R.cardinality()
	
	f /= f.coefficients().pop(0)
	f = f.change_ring(ZZ)

	G = Sequence([], f.parent())
	for i in range(m+1):
		base = N^(m-i) * f^i
		for shifts in itertools.product(range(d), repeat=f.nvariables()):
			g = base * prod(map(power, f.variables(), shifts))
			G.append(g)

	B, monomials = G.coefficient_matrix()
	monomials = vector(monomials)

	factors = [monomial(*bounds) for monomial in monomials]
	for i, factor in enumerate(factors):
		B.rescale_col(i, factor)

	B = B.dense_matrix().LLL()

	B = B.change_ring(QQ)
	for i, factor in enumerate(factors):
		B.rescale_col(i, 1/factor)

	H = Sequence([], f.parent().change_ring(QQ))
	for h in filter(None, B*monomials):
		H.append(h)
		I = H.ideal()
		if I.dimension() == -1:
			H.pop()
		elif I.dimension() == 0:
			roots = []
			for root in I.variety(ring=ZZ):
				root = tuple(R(root[var]) for var in f.variables())
				roots.append(root)
			return roots

	return []

Boneh and Durfee attack

boneh attack攻击范围比wiener更广，如果wiener解不出来，可以尝试使用这个攻击，当$d$较小时，满足$d<N^{0.292}$。接下来便是原理分析：

$ed\equiv1(mod\ \varphi(N)/2)\\ ed+k\varphi(N)/2=1\\ k\varphi(N)/2\equiv1(mod\ e)\\ \varphi(N)=(p-1)(q-1)=pq-p-1+1=N-p-q+1\\ k(N-p-q+1)/2\equiv1(mod\ e)$

假设$A=\frac{N+1}{2}，y=\frac{-p-1}{2}$，原式可化为

$f(k, y)=k(A+y)\equiv1(mod\ e)$

其中

$|k|<\frac{2ed}{\varphi(N)}<\frac{3ed}{N}=3*\frac{e}{N}*d<3*\frac{e}{N}*N^{delta}\\ |y|<2*N^{0.5}$

$y$的估计用到了$p,q$比较均匀的假设。这里$delta$为预估的小于$0.292$的值。如果我们求得了该二元方程的根，那么我们自然也就可以解一元二次方程$N=pq，p+q=-2y$来得到$p$与$q$。虽然原理也和没说一样，github上有基于sagemath的拓展脚本，虽然原理看起来并不难懂，算法实现确实异常困难，亦是基于格的攻击。以下是项目文件的源码，同时将其改成可调用版本，使用方法便是example(e, n)，实际上后面还有两个参数，第三个参数是用来确定上界，第四个参数是用来确定格的维度。example(e, n, delta=.18)，delta默认是0.18，传入时不能大于0.292，example(e, n, delta=.18, m=4)，格的维度越大，解的速度越慢，实际上控制好delta的大小便能迅速出答案，m调大虽然也能拿到答案但速度就慢，当delta如何调都跑不出来，而且理论分析的界满足攻击情况时，可以适当将m调大看看，m默认是$4$。

import time

############################################
# Config
##########################################

"""
Setting debug to true will display more informations
about the lattice, the bounds, the vectors...
"""
debug = True

"""
Setting strict to true will stop the algorithm (and
return (-1, -1)) if we don't have a correct
upperbound on the determinant. Note that this
doesn't necesseraly mean that no solutions
will be found since the theoretical upperbound is
usualy far away from actual results. That is why
you should probably use `strict = False`
"""
strict = False

"""
This is experimental, but has provided remarkable results
so far. It tries to reduce the lattice as much as it can
while keeping its efficiency. I see no reason not to use
this option, but if things don't work, you should try
disabling it
"""
helpful_only = True
dimension_min = 7 # stop removing if lattice reaches that dimension

############################################
# Functions
##########################################

# display stats on helpful vectors
def helpful_vectors(BB, modulus):
    nothelpful = 0
    for ii in range(BB.dimensions()[0]):
        if BB[ii,ii] >= modulus:
            nothelpful += 1

    print(nothelpful, "/", BB.dimensions()[0], " vectors are not helpful")

# display matrix picture with 0 and X
def matrix_overview(BB, bound):
    for ii in range(BB.dimensions()[0]):
        a = ('%02d ' % ii)
        for jj in range(BB.dimensions()[1]):
            a += '0' if BB[ii,jj] == 0 else 'X'
            if BB.dimensions()[0] < 60:
                a += ' '
        if BB[ii, ii] >= bound:
            a += '~'
        print(a)

# tries to remove unhelpful vectors
# we start at current = n-1 (last vector)
def remove_unhelpful(BB, monomials, bound, current):
    # end of our recursive function
    if current == -1 or BB.dimensions()[0] <= dimension_min:
        return BB

    # we start by checking from the end
    for ii in range(current, -1, -1):
        # if it is unhelpful:
        if BB[ii, ii] >= bound:
            affected_vectors = 0
            affected_vector_index = 0
            # let's check if it affects other vectors
            for jj in range(ii + 1, BB.dimensions()[0]):
                # if another vector is affected:
                # we increase the count
                if BB[jj, ii] != 0:
                    affected_vectors += 1
                    affected_vector_index = jj

            # level:0
            # if no other vectors end up affected
            # we remove it
            if affected_vectors == 0:
                print("* removing unhelpful vector", ii)
                BB = BB.delete_columns([ii])
                BB = BB.delete_rows([ii])
                monomials.pop(ii)
                BB = remove_unhelpful(BB, monomials, bound, ii-1)
                return BB

            # level:1
            # if just one was affected we check
            # if it is affecting someone else
            elif affected_vectors == 1:
                affected_deeper = True
                for kk in range(affected_vector_index + 1, BB.dimensions()[0]):
                    # if it is affecting even one vector
                    # we give up on this one
                    if BB[kk, affected_vector_index] != 0:
                        affected_deeper = False
                # remove both it if no other vector was affected and
                # this helpful vector is not helpful enough
                # compared to our unhelpful one
                if affected_deeper and abs(bound - BB[affected_vector_index, affected_vector_index]) < abs(bound - BB[ii, ii]):
                    print("* removing unhelpful vectors", ii, "and", affected_vector_index)
                    BB = BB.delete_columns([affected_vector_index, ii])
                    BB = BB.delete_rows([affected_vector_index, ii])
                    monomials.pop(affected_vector_index)
                    monomials.pop(ii)
                    BB = remove_unhelpful(BB, monomials, bound, ii-1)
                    return BB
    # nothing happened
    return BB

""" 
Returns:
* 0,0   if it fails
* -1,-1 if `strict=true`, and determinant doesn't bound
* x0,y0 the solutions of `pol`
"""
def boneh_durfee(pol, modulus, mm, tt, XX, YY):
    """
    Boneh and Durfee revisited by Herrmann and May
    
    finds a solution if:
    * d < N^delta
    * |x| < e^delta
    * |y| < e^0.5
    whenever delta < 1 - sqrt(2)/2 ~ 0.292
    """

    # substitution (Herrman and May)
    PR.<u, x, y> = PolynomialRing(ZZ)
    Q = PR.quotient(x*y + 1 - u) # u = xy + 1
    polZ = Q(pol).lift()

    UU = XX*YY + 1

    # x-shifts
    gg = []
    for kk in range(mm + 1):
        for ii in range(mm - kk + 1):
            xshift = x^ii * modulus^(mm - kk) * polZ(u, x, y)^kk
            gg.append(xshift)
    gg.sort()

    # x-shifts list of monomials
    monomials = []
    for polynomial in gg:
        for monomial in polynomial.monomials():
            if monomial not in monomials:
                monomials.append(monomial)
    monomials.sort()
    
    # y-shifts (selected by Herrman and May)
    for jj in range(1, tt + 1):
        for kk in range(floor(mm/tt) * jj, mm + 1):
            yshift = y^jj * polZ(u, x, y)^kk * modulus^(mm - kk)
            yshift = Q(yshift).lift()
            gg.append(yshift) # substitution
    
    # y-shifts list of monomials
    for jj in range(1, tt + 1):
        for kk in range(floor(mm/tt) * jj, mm + 1):
            monomials.append(u^kk * y^jj)

    # construct lattice B
    nn = len(monomials)
    BB = Matrix(ZZ, nn)
    for ii in range(nn):
        BB[ii, 0] = gg[ii](0, 0, 0)
        for jj in range(1, ii + 1):
            if monomials[jj] in gg[ii].monomials():
                BB[ii, jj] = gg[ii].monomial_coefficient(monomials[jj]) * monomials[jj](UU,XX,YY)

    # Prototype to reduce the lattice
    if helpful_only:
        # automatically remove
        BB = remove_unhelpful(BB, monomials, modulus^mm, nn-1)
        # reset dimension
        nn = BB.dimensions()[0]
        if nn == 0:
            print("failure")
            return 0,0

    # check if vectors are helpful
    if debug:
        helpful_vectors(BB, modulus^mm)
    
    # check if determinant is correctly bounded
    det = BB.det()
    bound = modulus^(mm*nn)
    if det >= bound:
        print("We do not have det < bound. Solutions might not be found.")
        print("Try with highers m and t.")
        if debug:
            diff = (log(det) - log(bound)) / log(2)
            print("size det(L) - size e^(m*n) = ", floor(diff))
        if strict:
            return -1, -1
    else:
        print("det(L) < e^(m*n) (good! If a solution exists < N^delta, it will be found)")

    # display the lattice basis
    if debug:
        matrix_overview(BB, modulus^mm)

    # LLL
    if debug:
        print("optimizing basis of the lattice via LLL, this can take a long time")

    BB = BB.LLL()

    if debug:
        print("LLL is done!")

    # transform vector i & j -> polynomials 1 & 2
    if debug:
        print("looking for independent vectors in the lattice")
    found_polynomials = False
    
    for pol1_idx in range(nn - 1):
        for pol2_idx in range(pol1_idx + 1, nn):
            # for i and j, create the two polynomials
            PR.<w,z> = PolynomialRing(ZZ)
            pol1 = pol2 = 0
            for jj in range(nn):
                pol1 += monomials[jj](w*z+1,w,z) * BB[pol1_idx, jj] / monomials[jj](UU,XX,YY)
                pol2 += monomials[jj](w*z+1,w,z) * BB[pol2_idx, jj] / monomials[jj](UU,XX,YY)

            # resultant
            PR.<q> = PolynomialRing(ZZ)
            rr = pol1.resultant(pol2)

            # are these good polynomials?
            if rr.is_zero() or rr.monomials() == [1]:
                continue
            else:
                print("found them, using vectors", pol1_idx, "and", pol2_idx)
                found_polynomials = True
                break
        if found_polynomials:
            break

    if not found_polynomials:
        print("no independant vectors could be found. This should very rarely happen...")
        return 0, 0
    
    rr = rr(q, q)

    # solutions
    soly = rr.roots()

    if len(soly) == 0:
        print("Your prediction (delta) is too small")
        return 0, 0

    soly = soly[0][0]
    ss = pol1(q, soly)
    solx = ss.roots()[0][0]

    #
    return solx, soly

def example(e, N, delta = .18, m = 4):
    ############################################
    # How To Use This Script
    ##########################################

    #
    # The problem to solve (edit the following values)
    #

    # the modulus
    # N = 0xc2fd2913bae61f845ac94e4ee1bb10d8531dda830d31bb221dac5f179a8f883f15046d7aa179aff848db2734b8f88cc73d09f35c445c74ee35b01a96eb7b0a6ad9cb9ccd6c02c3f8c55ecabb55501bb2c318a38cac2db69d510e152756054aaed064ac2a454e46d9b3b755b67b46906fbff8dd9aeca6755909333f5f81bf74db
    # the public exponent
    # e = 0x19441f679c9609f2484eb9b2658d7138252b847b2ed8ad182be7976ed57a3e441af14897ce041f3e07916445b88181c22f510150584eee4b0f776a5a487a4472a99f2ddc95efdd2b380ab4480533808b8c92e63ace57fb42bac8315fa487d03bec86d854314bc2ec4f99b192bb98710be151599d60f224114f6b33f47e357517

    # the hypothesis on the private exponent (the theoretical maximum is 0.292)
    # delta = .18 # this means that d < N^delta

    #
    # Lattice (tweak those values)
    #

    # you should tweak this (after a first run), (e.g. increment it until a solution is found)
    # m = 4 # size of the lattice (bigger the better/slower)

    # you need to be a lattice master to tweak these
    t = int((1-2*delta) * m)  # optimization from Herrmann and May
    X = 2*floor(N^delta)  # this _might_ be too much
    Y = floor(N^(1/2))    # correct if p, q are ~ same size

    #
    # Don't touch anything below
    #

    # Problem put in equation
    P.<x,y> = PolynomialRing(ZZ)
    A = int((N+1)/2)
    pol = 1 + x * (A + y)

    #
    # Find the solutions!
    #

    # Checking bounds
    if debug:
        print("=== checking values ===")
        print("* delta:", delta)
        print("* delta < 0.292", delta < 0.292)
        print("* size of e:", int(log(e)/log(2)))
        print("* size of N:", int(log(N)/log(2)))
        print("* m:", m, ", t:", t)

    # boneh_durfee
    if debug:
        print("=== running algorithm ===")
        start_time = time.time()

    solx, soly = boneh_durfee(pol, e, m, t, X, Y)

    # found a solution?
    if solx > 0:
        print("=== solution found ===")
        if False:
            print("x:", solx)
            print("y:", soly)

        d = int(pol(solx, soly) / e)
        print("private key found:", d)
    else:
        print("=== no solution was found ===")

    if debug:
        print(("=== %s seconds ===" % (time.time() - start_time)))

[2023 基地"楚慧杯"]so-large-e

from Crypto.Util.number import *
from Crypto.PublicKey import RSA
from flag import flag
import random

m = bytes_to_long(flag)

p = getPrime(512)
q = getPrime(512)
n = p*q
e = random.getrandbits(1024)
assert size(e)==1024
phi = (p-1)*(q-1)
assert GCD(e,phi)==1
d = inverse(e,phi)
assert size(d)==269

pub = (n, e)
PublicKey = RSA.construct(pub)
with open('pub.pem', 'wb') as f :
    f.write(PublicKey.exportKey())

c = pow(m,e,n)
print('c =',c)
print(long_to_bytes(pow(c,d,n)))

# c = 6838759631922176040297411386959306230064807618456930982742841698524622016849807235726065272136043603027166249075560058232683230155346614429566511309977857815138004298815137913729662337535371277019856193898546849896085411001528569293727010020290576888205244471943227253000727727343731590226737192613447347860
# e = 113449247876071397911206070019495939088171696712182747502133063172021565345788627261740950665891922659340020397229619329204520999096535909867327960323598168596664323692312516466648588320607291284630435682282630745947689431909998401389566081966753438869725583665294310689820290368901166811028660086977458571233
# n = 116518679305515263290840706715579691213922169271634579327519562902613543582623449606741546472920401997930041388553141909069487589461948798111698856100819163407893673249162209631978914843896272256274862501461321020961958367098759183487116417487922645782638510876609728886007680825340200888068103951956139343723

可以尝试wiener attack，是跑不出答案的，d和n的位数都已知，算一下269/1024可以得到0.26的答案，这很明显超出wiener的范围，而这时就得用到boneh attack。

c = 6838759631922176040297411386959306230064807618456930982742841698524622016849807235726065272136043603027166249075560058232683230155346614429566511309977857815138004298815137913729662337535371277019856193898546849896085411001528569293727010020290576888205244471943227253000727727343731590226737192613447347860
e = 113449247876071397911206070019495939088171696712182747502133063172021565345788627261740950665891922659340020397229619329204520999096535909867327960323598168596664323692312516466648588320607291284630435682282630745947689431909998401389566081966753438869725583665294310689820290368901166811028660086977458571233
n = 116518679305515263290840706715579691213922169271634579327519562902613543582623449606741546472920401997930041388553141909069487589461948798111698856100819163407893673249162209631978914843896272256274862501461321020961958367098759183487116417487922645782638510876609728886007680825340200888068103951956139343723
load('boneh_durffe.sage')
example(e, n)			# 可以看到是无d解出
example(e, n, 0.26)		# 将delta调到0.26,就能格约束出来

AMM算法

有限域上的高次开根AMM算法可以解决ctf部分rsa中的e和phi不互素的问题。

在$e=2$情况下：

$m^2\equiv c(mod\ q),q-1=2^t*s\\ x^{q-1}\equiv1(mod\ q),\quad x^{\frac{q-1}{2}}=\begin{cases}x^{2^{t-1}*s}\equiv1(mod\ q) \\ x^{2^{t-1}*s}\equiv-1(mod\ q) \end{cases}$

$t=1$：对式子两边乘上$x$并进行开方。

$x^s\equiv1(mod\ q)\Rightarrow x^{s+1}\equiv x(mod\ q)\Rightarrow x^{\frac{s+1}{2}}\equiv x^{\frac{1}{2}}(mod\ q)$

将$c$代入$x^{\frac{s+1}{2}}$：

$c^{\frac{s+1}{2}}\equiv c^{\frac{1}{2}}(mod\ q)=m(mod\ q)$

$t>=2$：

$\begin{cases}x^{2^{t-1}*s}\equiv1(mod\ q)\\x^{2^{t-1}*s}\equiv-1(mod\ q)\end{cases}$

对其开根会有剩余与非剩余情况，为方便计算，对非剩余的情况再次引入一个非剩余量相乘，使得右侧恒为$1$。

$x^{2^{t-2}*s}*y^{2^{t-1}*s}\equiv1(mod\ q)$

将两种情况合并，引入变量$k$，控制是否需要乘上这个非二次剩余。$k=0$，为剩余类，$k=1$，为非剩余类，此后不断对$x$进行开根重复操作：

$x^{2^{t-2}*s}*y^{2^{t-1}*s*k}\equiv1(mod\ q)\\ x^s*y^{s*(2^1k_1+2^2k_2+2^3k_3+…+2^{t-1}k_{t-1})}\equiv1(mod\ q)$

之后将式子乘上$x$再开方代入$c$：

$x^{\frac{s+1}{2}}*y^{s*(k_1+2k_2+2^2k_3+…+2^{t-2}k_{t-1})}\equiv1(mod\ q)\\ c^{\frac{s+1}{2}}*y^{s*(k_1+2k_2+2^2k_3+…+2^{t-2}k_{t-1})}\equiv1(mod\ q)$

在$e>2$情况下：

$x^{e}\equiv\delta(mod\ q)，q-1=e^t*s$

构造以下式子：

$\delta^{\frac{q-1}{e}}\equiv(x^e)^{\frac{q-1}{e}}\equiv x^{q-1}\equiv1(mod\ q)$

$gcd(e,s)=1$，则可以构造$e\alpha-s\beta=1$，拓展欧几原式应为$+$，方便式子推去则取相反数。推出$s\beta=e\alpha-1$,$s|(e\alpha-1)$。

$(\delta^{\frac{q-1}{e}})^\beta\equiv(\delta^{\frac{e^t*s}{e}})^\beta\equiv\delta^{e^{t-1}*s\beta}\equiv\delta^{e^{t-1}*(e\alpha-1)}\equiv1(mod\ q)$

​ $t=1$：

$\delta^{e\alpha}\equiv\delta(mod\ q),(\delta^{\alpha})^e\equiv\delta(mod\ q)$

当$t>2$时：

取$r$次非剩余$\rho$:

$\rho^{\frac{q-1}{e}}\neq1(mod\ q)\Rightarrow不存在x 使得(x^e)^{\frac{q-1}{e}}=1(mod\ q)$

构造集合

$K_i=\rho^{i*\frac{q-1}{e}}=\rho^{i*s*e^{t-1}},i的范围[0,e-1]\\ K_i*K_{e-i}=\rho^{i*\frac{q-1}{e}}*\rho^{(e-i)*{\frac{q-1}{e}}}=\rho^{e*\frac{q-1}{e}}=\rho^{q-1}\equiv1$

即两者互为逆元，$\rho$是之前取的一个$e$次非剩余，我们可以定义一个生成元为$\rho$的模$q$的循环群$G$。

$G=\{1,\rho^2,\rho^3,…,\rho^{q-2},\rho^{q-1} \}$

在该循环群$G$中，我们再定义$\rho_1=\rho^{\frac{q-1}{e}}$，所以可以得到一个生成元为$\rho_1$的子群$H$为：

$H = \{1,\rho^{\frac{q-1}{e}},\rho^{\frac{q-1}{e}*2},\rho^{\frac{q-1}{e}*3},…,\rho^{\frac{q-1}{e}*(e-1)} \}$

在这个子群里，所有的元素都是由$\rho_1=\rho^{\frac{q-1}{e}}$生成，所以所有元素的$e$次方都为$1$，因为子群$H$里正好有$e$个元素，所以$H$是模$q$开$e$次方根的结果的集合。同理，$K$中每个元素对应一个开根结果。

$\delta^{\frac{q-1}{e}}\equiv\delta^{e^{t-1}*(e\alpha-1)}\equiv(\delta^{e^{t-2}*(e\alpha-1)})^e$

设开$e$次方下对应集合中的第$j$个数：

$\delta^{e^{t-2}*(e\alpha-1)}\equiv K_j\\ \delta^{e^{t-2}*(e\alpha-1)}*K_{e-j}\equiv1$

之后对式子不断开$e$次方，乘上逆元，直到把$t$消掉，最终得到式子：

$\delta^{e\alpha-1}*(\rho^s)^{j_1e+j_2e^2+…+j_{t-1}e^{t-1}}\equiv1\\ (\delta^\alpha*\rho^{s*(j_1+j_2e+…+j_{t-1}e^{t-2})})^e\equiv\delta$

这样算出来的式子是其中一个解，根据有生成元的循环群性质，则可以推出剩余的解：

$(x^\frac{q-1}{e})^e\equiv x^{q-1}\equiv1\\ m_1=(m_0*x^{\frac{q-1}{e}})^e\equiv c(mod\ q)$

import gmpy2
from libnum import *
import random
import math
import time

def onemod(e, q):
    p = random.randint(1, q - 1)
    while (pow(p, (q - 1) // e, q) == 1):  # (r,s)=1
        p = random.randint(1, q)
    return p

def AMM_rth(o, r, q):  # r|(q-1)
    print('start to calculate primitive root...')
    start = time.time()
    assert ((q - 1) % r == 0)
    p = onemod(r, q)
    print('p:%d' % p)
    t = 0
    s = q - 1
    while (s % r == 0):
        s = s // r
        t += 1
    k = 1
    while ((s * k + 1) % r != 0):
        k += 1
    alp = (s * k + 1) // r
    print('s:%d,t:%d r:%d alp:%d' % (s, t, r, alp))
    a = pow(p, r ** (t - 1) * s, q)
    b = pow(o, r * a - 1, q)
    c = pow(p, s, q)
    h = 1
    for i in range(1, t - 1):
        d = pow(int(b), r ** (t - 1 - i), q)
        print('d:%d' % d)
        if d == 1:
            j = 0
        else:
            j = (-math.log(d, a)) % r
        b = (b * (c ** (r * j))) % q
        h = (h * c ** j) % q
        c = (c * r) % q
    result = (pow(o, alp, q) * h)
    end = time.time()

    print('result:%d' % result)
    print('Finish in {} seconds'.format(end - start))
    return result

# 2020.4.3 finish
def ALL_Solution(m, q, rt, cq, e):
    print('start to calculate all root...')
    start = time.time()
    mp = []
    for pr in rt:
        r = (pr * m) % q
        assert (pow(r, e, q) == cq)
        mp.append(r)
    end = time.time()
    print('Finish in {} seconds'.format(end - start))
    return mp

def calc(mp, mq, e, p, q):
    print('satrting crt... ')
    i = 1
    j = 1
    start = time.time()
    t1 = gmpy2.invert(q, p)
    t2 = gmpy2.invert(p, q)
    for mp1 in mp:
        if (j % 1000 == 0):
            print(j)
        j += 1
        for mq1 in mq:
            ans = (mp1 * t1 * q + mq1 * t2 * p) % (p * q)
            if check(ans):
                print('Finish in {} seconds'.format(time.time() - start))
                return
    return


def check(m):
    try:
        a = n2s(m)
        if a.startswith('NCTF'):
            print(a)
            return True
        else:
            return False
    except:
        return False


def ALL_ROOT2(r, q):  # use function set() and .add() ensure that the generated elements are not repeated
    print('start to find all primitive root...')
    start = time.time()
    li = set()
    while (len(li) < r):
        p = pow(random.randint(1, q - 1), (q - 1) // r, q)
        li.add(p)
    end = time.time()
    print('Finish in {} seconds'.format(end - start))
    return li


if __name__ == '__main__':
    c = 10562302690541901187975815594605242014385201583329309191736952454310803387032252007244962585846519762051885640856082157060593829013572592812958261432327975138581784360302599265408134332094134880789013207382277849503344042487389850373487656200657856862096900860792273206447552132458430989534820256156021128891296387414689693952047302604774923411425863612316726417214819110981605912408620996068520823370069362751149060142640529571400977787330956486849449005402750224992048562898004309319577192693315658275912449198365737965570035264841782399978307388920681068646219895287752359564029778568376881425070363592696751183359
    p = 199138677823743837339927520157607820029746574557746549094921488292877226509198315016018919385259781238148402833316033634968163276198999279327827901879426429664674358844084491830543271625147280950273934405879341438429171453002453838897458102128836690385604150324972907981960626767679153125735677417397078196059
    q = 112213695905472142415221444515326532320352429478341683352811183503269676555434601229013679319423878238944956830244386653674413411658696751173844443394608246716053086226910581400528167848306119179879115809778793093611381764939789057524575349501163689452810148280625226541609383166347879832134495444706697124741
    e = 0x1337
    cp = c % p
    cq = c % q

    mp = AMM_rth(cp, e, p)
    mq = AMM_rth(cq, e, q)

    rt1 = ALL_ROOT2(e, p)
    rt2 = ALL_ROOT2(e, q)

    amp = ALL_Solution(mp, p, rt1, cp, e)
    amq = ALL_Solution(mq, q, rt2, cq, e)

    calc(amp, amq, e, p, q)

多项式环上的RSA

任意取两个多项式$g(p),g(q)$

$g(n)=g(p)*g(q)\\ \varphi(g(n))=phi\\ g(c)\equiv g(m)^{e}(mod\ g(n))\\ g(m)\equiv g(c)^{d}(mod\ g(n))$

对于素数,$\varphi(x)=x-1$，但是对于不可约多项式$p(x)$，$\varphi(p(x))=x-1$是不成立的
不可约多项式$p(x)$，除了$0$，长度为$n$每一个多项式都与$p(x)$互素，因此$\varphi(p(x))=2^{n}-1$

$P(x)=a_nx^n+a_{n-1}x^{n-1}+\cdots+a_1x+a_0\\ Q(x)=b_mx^m+b_{m-1}x_{m-1}+\cdots+b_1x+b_0\\ N(x)=P(x)Q(x)\\$

度小于$m+n$的多项式$K(x)=\sum_{i=0}^{m+n-1}k_ix^i,k_i\in\Z_p$。系数有$m+n$个($k_0,k_1,\cdots,k_{m+n-1}$)，那么组合情况共有$p^{m+n}$种。

而对于多项式$P(x)$，找到$P(x)$的倍数的话，即$P(x)A(x)$存在，求$A(x)$共有多少种情况。$A(x)=\sum_{i=0}^{m-1}k_ix^i,k_i\in\Z_p$，系数有$m$个($k_0,k_1,\cdots,k_{m-1}$)，组合情况有$p^m$种。

对于多项式$Q(x)$也类似，组合情况有$p^n$种。

$p^{m+n}$是所有度小于$m+n$的多项式总数，$p^m$则是$P(x)$多项式倍数个数，$p^n$则是$Q(x)$多项式倍数个数。除去$P(x),Q(x)$的倍数，剩下都是与$N(x)$互素的多项式，那么这些多项式个数便是多项式的欧拉函数。

$s=(p^m-1)(p^n-1)$

例题

flag = bytearray(raw_input())
flag = list(flag)
length = len(flag)
bits = 16

## Prime for Finite Field.
p = random_prime(2^bits-1, False, 2^(bits-1))

file_out = open("downloads/polynomial_rsa.txt", "w")
file_out.write("Prime: " + str(p) + "\n")

## Univariate Polynomial Ring in y over Finite Field of size p
R.<y> = PolynomialRing(GF(p))

## Analogous to the primes in Z
def gen_irreducable_poly(deg):
    while True:
        out = R.random_element(degree=deg)
        if out.is_irreducible():
            return out


## Polynomial "primes"
P = gen_irreducable_poly(ZZ.random_element(length, 2*length))
Q = gen_irreducable_poly(ZZ.random_element(length, 2*length))

## Public exponent key
e = 65537

## Modulus
N = P*Q
file_out.write("Modulus: " + str(N) + "\n")

## Univariate Quotient Polynomial Ring in x over Finite Field of size 659 with modulus N(x)
S.<x> = R.quotient(N)

## Encrypt
m = S(flag)
c = m^e

file_out.write("Ciphertext: " + str(c))
file_out.close()

"""
Prime: 43753
Modulus: 34036*y^177 + 23068*y^176 + 13147*y^175 + 36344*y^174 + 10045*y^173 + 41049*y^172 + 17786*y^171 + 16601*y^170 + 7929*y^169 + 37570*y^168 + 990*y^167 + 9622*y^166 + 39273*y^165 + 35284*y^164 + 15632*y^163 + 18850*y^162 + 8800*y^161 + 33148*y^160 + 12147*y^159 + 40487*y^158 + 6407*y^157 + 34111*y^156 + 8446*y^155 + 21908*y^154 + 16812*y^153 + 40624*y^152 + 43506*y^151 + 39116*y^150 + 33011*y^149 + 23914*y^148 + 2210*y^147 + 23196*y^146 + 43359*y^145 + 34455*y^144 + 17684*y^143 + 25262*y^142 + 982*y^141 + 24015*y^140 + 27968*y^139 + 37463*y^138 + 10667*y^137 + 39519*y^136 + 31176*y^135 + 27520*y^134 + 32118*y^133 + 8333*y^132 + 38945*y^131 + 34713*y^130 + 1107*y^129 + 43604*y^128 + 4433*y^127 + 18110*y^126 + 17658*y^125 + 32354*y^124 + 3219*y^123 + 40238*y^122 + 10439*y^121 + 3669*y^120 + 8713*y^119 + 21027*y^118 + 29480*y^117 + 5477*y^116 + 24332*y^115 + 43480*y^114 + 33406*y^113 + 43121*y^112 + 1114*y^111 + 17198*y^110 + 22829*y^109 + 24424*y^108 + 16523*y^107 + 20424*y^106 + 36206*y^105 + 41849*y^104 + 3584*y^103 + 26500*y^102 + 31897*y^101 + 34640*y^100 + 27449*y^99 + 30962*y^98 + 41434*y^97 + 22125*y^96 + 24314*y^95 + 3944*y^94 + 18400*y^93 + 38476*y^92 + 28904*y^91 + 27936*y^90 + 41867*y^89 + 25573*y^88 + 25659*y^87 + 33443*y^86 + 18435*y^85 + 5934*y^84 + 38030*y^83 + 17563*y^82 + 24086*y^81 + 36782*y^80 + 20922*y^79 + 38933*y^78 + 23448*y^77 + 10599*y^76 + 7156*y^75 + 29044*y^74 + 23605*y^73 + 7657*y^72 + 28200*y^71 + 2431*y^70 + 3860*y^69 + 23259*y^68 + 14590*y^67 + 33631*y^66 + 15673*y^65 + 36049*y^64 + 29728*y^63 + 22413*y^62 + 18602*y^61 + 18557*y^60 + 23505*y^59 + 17642*y^58 + 12595*y^57 + 17255*y^56 + 15316*y^55 + 8948*y^54 + 38*y^53 + 40329*y^52 + 9823*y^51 + 5798*y^50 + 6379*y^49 + 8662*y^48 + 34640*y^47 + 38321*y^46 + 18760*y^45 + 13135*y^44 + 15926*y^43 + 34952*y^42 + 28940*y^41 + 13558*y^40 + 42579*y^39 + 38015*y^38 + 33788*y^37 + 12381*y^36 + 195*y^35 + 13709*y^34 + 31500*y^33 + 32994*y^32 + 30486*y^31 + 40414*y^30 + 2578*y^29 + 30525*y^28 + 43067*y^27 + 6195*y^26 + 36288*y^25 + 23236*y^24 + 21493*y^23 + 15808*y^22 + 34500*y^21 + 6390*y^20 + 42994*y^19 + 42151*y^18 + 19248*y^17 + 19291*y^16 + 8124*y^15 + 40161*y^14 + 24726*y^13 + 31874*y^12 + 30272*y^11 + 30761*y^10 + 2296*y^9 + 11017*y^8 + 16559*y^7 + 28949*y^6 + 40499*y^5 + 22377*y^4 + 33628*y^3 + 30598*y^2 + 4386*y + 23814
Ciphertext: 5209*x^176 + 10881*x^175 + 31096*x^174 + 23354*x^173 + 28337*x^172 + 15982*x^171 + 13515*x^170 + 21641*x^169 + 10254*x^168 + 34588*x^167 + 27434*x^166 + 29552*x^165 + 7105*x^164 + 22604*x^163 + 41253*x^162 + 42675*x^161 + 21153*x^160 + 32838*x^159 + 34391*x^158 + 832*x^157 + 720*x^156 + 22883*x^155 + 19236*x^154 + 33772*x^153 + 5020*x^152 + 17943*x^151 + 26967*x^150 + 30847*x^149 + 10306*x^148 + 33966*x^147 + 43255*x^146 + 20342*x^145 + 4474*x^144 + 3490*x^143 + 38033*x^142 + 11224*x^141 + 30565*x^140 + 31967*x^139 + 32382*x^138 + 9759*x^137 + 1030*x^136 + 32122*x^135 + 42614*x^134 + 14280*x^133 + 16533*x^132 + 32676*x^131 + 43070*x^130 + 36009*x^129 + 28497*x^128 + 2940*x^127 + 9747*x^126 + 22758*x^125 + 16615*x^124 + 14086*x^123 + 13038*x^122 + 39603*x^121 + 36260*x^120 + 32502*x^119 + 17619*x^118 + 17700*x^117 + 15083*x^116 + 11311*x^115 + 36496*x^114 + 1300*x^113 + 13601*x^112 + 43425*x^111 + 10376*x^110 + 11551*x^109 + 13684*x^108 + 14955*x^107 + 6661*x^106 + 12674*x^105 + 21534*x^104 + 32132*x^103 + 34135*x^102 + 43684*x^101 + 837*x^100 + 29311*x^99 + 4849*x^98 + 26632*x^97 + 26662*x^96 + 10159*x^95 + 32657*x^94 + 12149*x^93 + 17858*x^92 + 35805*x^91 + 19391*x^90 + 30884*x^89 + 42039*x^88 + 17292*x^87 + 4694*x^86 + 1497*x^85 + 1744*x^84 + 31071*x^83 + 26246*x^82 + 24402*x^81 + 22068*x^80 + 39263*x^79 + 23703*x^78 + 21484*x^77 + 12241*x^76 + 28821*x^75 + 32886*x^74 + 43075*x^73 + 35741*x^72 + 19936*x^71 + 37219*x^70 + 33411*x^69 + 8301*x^68 + 12949*x^67 + 28611*x^66 + 42654*x^65 + 6910*x^64 + 18523*x^63 + 31144*x^62 + 21398*x^61 + 36298*x^60 + 27158*x^59 + 918*x^58 + 38601*x^57 + 4269*x^56 + 5699*x^55 + 36444*x^54 + 34791*x^53 + 37978*x^52 + 32481*x^51 + 8039*x^50 + 11012*x^49 + 11454*x^48 + 30450*x^47 + 1381*x^46 + 32403*x^45 + 8202*x^44 + 8404*x^43 + 37648*x^42 + 43696*x^41 + 34237*x^40 + 36490*x^39 + 41423*x^38 + 35792*x^37 + 36950*x^36 + 31086*x^35 + 38970*x^34 + 12439*x^33 + 7963*x^32 + 16150*x^31 + 11382*x^30 + 3038*x^29 + 20157*x^28 + 23531*x^27 + 32866*x^26 + 5428*x^25 + 21132*x^24 + 13443*x^23 + 28909*x^22 + 42716*x^21 + 6567*x^20 + 24744*x^19 + 8727*x^18 + 14895*x^17 + 28172*x^16 + 30903*x^15 + 26608*x^14 + 27314*x^13 + 42224*x^12 + 42551*x^11 + 37726*x^10 + 11203*x^9 + 36816*x^8 + 5537*x^7 + 20301*x^6 + 17591*x^5 + 41279*x^4 + 7999*x^3 + 33753*x^2 + 34551*x + 9659
"""

p = 43753
R.<y> = PolynomialRing(GF(p))
N = R("34036*y^177 + 23068*y^176 + 13147*y^175 + 36344*y^174 + 10045*y^173 + 41049*y^172 + 17786*y^171 + 16601*y^170 + 7929*y^169 + 37570*y^168 + 990*y^167 + 9622*y^166 + 39273*y^165 + 35284*y^164 + 15632*y^163 + 18850*y^162 + 8800*y^161 + 33148*y^160 + 12147*y^159 + 40487*y^158 + 6407*y^157 + 34111*y^156 + 8446*y^155 + 21908*y^154 + 16812*y^153 + 40624*y^152 + 43506*y^151 + 39116*y^150 + 33011*y^149 + 23914*y^148 + 2210*y^147 + 23196*y^146 + 43359*y^145 + 34455*y^144 + 17684*y^143 + 25262*y^142 + 982*y^141 + 24015*y^140 + 27968*y^139 + 37463*y^138 + 10667*y^137 + 39519*y^136 + 31176*y^135 + 27520*y^134 + 32118*y^133 + 8333*y^132 + 38945*y^131 + 34713*y^130 + 1107*y^129 + 43604*y^128 + 4433*y^127 + 18110*y^126 + 17658*y^125 + 32354*y^124 + 3219*y^123 + 40238*y^122 + 10439*y^121 + 3669*y^120 + 8713*y^119 + 21027*y^118 + 29480*y^117 + 5477*y^116 + 24332*y^115 + 43480*y^114 + 33406*y^113 + 43121*y^112 + 1114*y^111 + 17198*y^110 + 22829*y^109 + 24424*y^108 + 16523*y^107 + 20424*y^106 + 36206*y^105 + 41849*y^104 + 3584*y^103 + 26500*y^102 + 31897*y^101 + 34640*y^100 + 27449*y^99 + 30962*y^98 + 41434*y^97 + 22125*y^96 + 24314*y^95 + 3944*y^94 + 18400*y^93 + 38476*y^92 + 28904*y^91 + 27936*y^90 + 41867*y^89 + 25573*y^88 + 25659*y^87 + 33443*y^86 + 18435*y^85 + 5934*y^84 + 38030*y^83 + 17563*y^82 + 24086*y^81 + 36782*y^80 + 20922*y^79 + 38933*y^78 + 23448*y^77 + 10599*y^76 + 7156*y^75 + 29044*y^74 + 23605*y^73 + 7657*y^72 + 28200*y^71 + 2431*y^70 + 3860*y^69 + 23259*y^68 + 14590*y^67 + 33631*y^66 + 15673*y^65 + 36049*y^64 + 29728*y^63 + 22413*y^62 + 18602*y^61 + 18557*y^60 + 23505*y^59 + 17642*y^58 + 12595*y^57 + 17255*y^56 + 15316*y^55 + 8948*y^54 + 38*y^53 + 40329*y^52 + 9823*y^51 + 5798*y^50 + 6379*y^49 + 8662*y^48 + 34640*y^47 + 38321*y^46 + 18760*y^45 + 13135*y^44 + 15926*y^43 + 34952*y^42 + 28940*y^41 + 13558*y^40 + 42579*y^39 + 38015*y^38 + 33788*y^37 + 12381*y^36 + 195*y^35 + 13709*y^34 + 31500*y^33 + 32994*y^32 + 30486*y^31 + 40414*y^30 + 2578*y^29 + 30525*y^28 + 43067*y^27 + 6195*y^26 + 36288*y^25 + 23236*y^24 + 21493*y^23 + 15808*y^22 + 34500*y^21 + 6390*y^20 + 42994*y^19 + 42151*y^18 + 19248*y^17 + 19291*y^16 + 8124*y^15 + 40161*y^14 + 24726*y^13 + 31874*y^12 + 30272*y^11 + 30761*y^10 + 2296*y^9 + 11017*y^8 + 16559*y^7 + 28949*y^6 + 40499*y^5 + 22377*y^4 + 33628*y^3 + 30598*y^2 + 4386*y + 23814")
C = R("5209*y^176 + 10881*y^175 + 31096*y^174 + 23354*y^173 + 28337*y^172 + 15982*y^171 + 13515*y^170 + 21641*y^169 + 10254*y^168 + 34588*y^167 + 27434*y^166 + 29552*y^165 + 7105*y^164 + 22604*y^163 + 41253*y^162 + 42675*y^161 + 21153*y^160 + 32838*y^159 + 34391*y^158 + 832*y^157 + 720*y^156 + 22883*y^155 + 19236*y^154 + 33772*y^153 + 5020*y^152 + 17943*y^151 + 26967*y^150 + 30847*y^149 + 10306*y^148 + 33966*y^147 + 43255*y^146 + 20342*y^145 + 4474*y^144 + 3490*y^143 + 38033*y^142 + 11224*y^141 + 30565*y^140 + 31967*y^139 + 32382*y^138 + 9759*y^137 + 1030*y^136 + 32122*y^135 + 42614*y^134 + 14280*y^133 + 16533*y^132 + 32676*y^131 + 43070*y^130 + 36009*y^129 + 28497*y^128 + 2940*y^127 + 9747*y^126 + 22758*y^125 + 16615*y^124 + 14086*y^123 + 13038*y^122 + 39603*y^121 + 36260*y^120 + 32502*y^119 + 17619*y^118 + 17700*y^117 + 15083*y^116 + 11311*y^115 + 36496*y^114 + 1300*y^113 + 13601*y^112 + 43425*y^111 + 10376*y^110 + 11551*y^109 + 13684*y^108 + 14955*y^107 + 6661*y^106 + 12674*y^105 + 21534*y^104 + 32132*y^103 + 34135*y^102 + 43684*y^101 + 837*y^100 + 29311*y^99 + 4849*y^98 + 26632*y^97 + 26662*y^96 + 10159*y^95 + 32657*y^94 + 12149*y^93 + 17858*y^92 + 35805*y^91 + 19391*y^90 + 30884*y^89 + 42039*y^88 + 17292*y^87 + 4694*y^86 + 1497*y^85 + 1744*y^84 + 31071*y^83 + 26246*y^82 + 24402*y^81 + 22068*y^80 + 39263*y^79 + 23703*y^78 + 21484*y^77 + 12241*y^76 + 28821*y^75 + 32886*y^74 + 43075*y^73 + 35741*y^72 + 19936*y^71 + 37219*y^70 + 33411*y^69 + 8301*y^68 + 12949*y^67 + 28611*y^66 + 42654*y^65 + 6910*y^64 + 18523*y^63 + 31144*y^62 + 21398*y^61 + 36298*y^60 + 27158*y^59 + 918*y^58 + 38601*y^57 + 4269*y^56 + 5699*y^55 + 36444*y^54 + 34791*y^53 + 37978*y^52 + 32481*y^51 + 8039*y^50 + 11012*y^49 + 11454*y^48 + 30450*y^47 + 1381*y^46 + 32403*y^45 + 8202*y^44 + 8404*y^43 + 37648*y^42 + 43696*y^41 + 34237*y^40 + 36490*y^39 + 41423*y^38 + 35792*y^37 + 36950*y^36 + 31086*y^35 + 38970*y^34 + 12439*y^33 + 7963*y^32 + 16150*y^31 + 11382*y^30 + 3038*y^29 + 20157*y^28 + 23531*y^27 + 32866*y^26 + 5428*y^25 + 21132*y^24 + 13443*y^23 + 28909*y^22 + 42716*y^21 + 6567*y^20 + 24744*y^19 + 8727*y^18 + 14895*y^17 + 28172*y^16 + 30903*y^15 + 26608*y^14 + 27314*y^13 + 42224*y^12 + 42551*y^11 + 37726*y^10 + 11203*y^9 + 36816*y^8 + 5537*y^7 + 20301*y^6 + 17591*y^5 + 41279*y^4 + 7999*y^3 + 33753*y^2 + 34551*y + 9659")
e = 65537
P, Q = factor(N)
P, Q = P[0], Q[0]
phi = (p**P.degree()-1)*(p**Q.degree()-1)
d = inverse_mod(e, phi)
m = pow(C, d, N)
bytes(m.coefficients())

复数域下RSA

[Crypto CTF 2019]Complex RSA

from gmpy2 import invert,lcm,is_prime
import sys
sys.setrecursionlimit(2047)

f = (378122348642214690905411683807377396279362526734068034297186493255873563264253248095197659138069664908850853277799239471404715546747080714581633876343291058815268009173602365587463522797812239900814474973941995698621021680129530453282192603316731832323767320307650941745085796583822798379896337325, 205549984221850341303682190742446959375043769671555741781145106776498798455293849755553794941345135675348633855787880355726112506553703853175271830126908219803257875131387292937296039523359975622001865593227631119497003599166091642640101746534401068507972834895023584843991641629874709898672559632)
e = 59107
p = 228517792080140341
q = 1675909164550923263854591345270445396052847869117231939809062226222204253885693425526434134321712288675268468398852452684029376569327518089966506865838909486699078280423099271324646863671350838232140981094611254627738568184261530942469845202934677427234062382272736609418198352717

n = p*q

def cadd(a,b,n):
    return (a[0]+b[0]%n,a[1]+b[1]%n)

def cmul(a,b,n):
    return ((a[0]*b[0]-a[1]*b[1])%n,(a[0]*b[1]+a[1]*b[0])%n)

def cpow(a,k,n):
    if(k==0):
        return (1,0)
    if(k==1):
        return a
    if(k%2==0):
        a=cmul(a,a,n)
        return cpow(a,k/2,n)
    else:
        return cmul(a,cpow(cmul(a,a,n),(k-1)/2,n),n)

o=lcm((p*p-1),(q*q-1))
fm=cpow(f,invert(e,o),n)
print(fm)

在复数域下欧拉函数发生了变化

$\varphi(p)=p^2-1\quad \varphi(q)=q^2-1\\ \varphi(n)=(p^2-1)(q^2-1)$

在$C/pC*$复数域下，实部和虚部都在0到p-1范围内，除了0+0i无法构成群内元素外，其他都与p互素，那么共有$p^2-1$种情况

结式

给定任意两个多项式f(x)和g(x)，如果两个多项式的结式为0，那么两个多项式有公因式。

$f(x)=2x^2+x+1,\quad g(x)=3x^3+4x^2+5x+6\\ R(f,g)=\left|\begin{matrix}2&1&1&&\\&2&1&1&\\&&2&1&1\\3&4&5&6&\\&3&4&5&6\end{matrix}\right|=158$

将f(x)的系数按g(x)的度进行3次重复，g(x)的系数按f(x)的度进行2次重复，列出的行列式自然阶数是f(x)的度加上g(x)的度。结式不为0，自然两个多项式没有公因式，通过sage进行gcd(f, g)验证也可得出结果是1。

结式可以求两个一元多项式的曲线方程

$\begin{cases}x=cos\theta\\y=sin\theta\end{cases} \Rightarrow t=tan\frac{\theta}{2}\quad \begin{cases}x=\frac{1-t^2}{1+t^2}\\y=\frac{2t}{1+t^2}\end{cases} \Rightarrow\begin{cases}(1+t^2)x-(1-t^2)=0\\(1+t^2)y-2t=0\end{cases} \Rightarrow\begin{cases}(x+1)t^2+(x-1)=0\\yt^2-2t+y=0\end{cases}$

$R=\left|\begin{matrix}x+1&0&x-1&\\&x+1&0&x-1\\y&-2&y\\&y&-2&y\end{matrix}\right|=4x^2+4y^2-4=0\\$

那么可以得到圆方程为$x^2+y^2=1$。

from Crypto.Util.number import *
from secret import flag

l = len(flag)
assert l == 56
x = bytes_to_long(flag[:l//2])
y = bytes_to_long(flag[l//2:])

p = getPrime(1024)
e = 65537

x = pow(x, e, p)
y = pow(y, e, p)

a = (7 * x + x * y + 77 * y ** 7) % p
b = (x ** 7 + 777 * y) % p
"""
print(f'p = {p}')
print(f'a = {a}')
print(f'b = {b}')
"""

p = 160676801612994301361202519503059426958636739446670462398261976532859847492256822690640058297338763725128097587993428329580105931247817467950370089691908132361316857330836120708767594061772979871315614755470773991633234068651435625372887767258609941208307491359777513843529144444836847722372845148836203335627
a = 30318995909014771647618268716833486449002423009996671727903532973647046764624121316716790986592523978549131384964872198795285872746623966910764159262479160147876027157581577141632378119375701270068263640642243000011932466519579133761464923463402462812787531220639360431295348786697861069940729757964584951972
b = 51036630170491152581994259808984114372634216659979376101433163181132141957563047348137651942358538069256102718534893846618166559129391336639526588292370462975735415885732360576961407017238385374280336346614960555565504032093702784952402038043052556719843691506943605133036720410419999467125928578673380637828

from sage.matrix.matrix2 import Matrix
from Crypto.Util.number import long_to_bytes

def resulant(f1, f2, var):
    return Matrix.determinant(f1.sylvester_matrix(f2, var))

p = 160676801612994301361202519503059426958636739446670462398261976532859847492256822690640058297338763725128097587993428329580105931247817467950370089691908132361316857330836120708767594061772979871315614755470773991633234068651435625372887767258609941208307491359777513843529144444836847722372845148836203335627
a = 30318995909014771647618268716833486449002423009996671727903532973647046764624121316716790986592523978549131384964872198795285872746623966910764159262479160147876027157581577141632378119375701270068263640642243000011932466519579133761464923463402462812787531220639360431295348786697861069940729757964584951972
b = 51036630170491152581994259808984114372634216659979376101433163181132141957563047348137651942358538069256102718534893846618166559129391336639526588292370462975735415885732360576961407017238385374280336346614960555565504032093702784952402038043052556719843691506943605133036720410419999467125928578673380637828
e = 65537

P.<x, y> = PolynomialRing(Zmod(p))
f1 = 7*x + x*y + 77*y**7 - a
f2 = x**7 + 777*y - b

hx = resulant(f1, f2, y)
rx = hx.univariate_polynomial().roots()
x, _ = zip(*rx)
y = [((b - i^7) * inverse_mod(777, p)) % p for i in x]

d = inverse_mod(e, p - 1)

for i in range(len(x)):
    m1 = pow(x[i], d, p)
    m2 = pow(y[i], d, p)
    try:
        print(bytes.fromhex(hex(m1)[2:]).decode(), end='')
        print(bytes.fromhex(hex(m2)[2:]).decode())
    except:
        pass

在sage中，直接使用long_to_bytes这类函数，会因为类型原因报错，因为通过pow运算出来的类型是<class 'sage.rings.finite_rings.integer_mod.IntegerMod_gmp'>，而函数所需的类型是int，这就需要在pow之外套上层int，或者就像题解中所使用的bytes.fromhex，hex进行层转换，好像也没多大差别。但decode这个函数的使用非常精妙，因为非flag中有着不可显字符的错在，会爆'utf-8' codec can't decode byte 0xb1 in position 0: invalid start byte错误，这样通过try能排掉错误答案，不再需要去人眼辨识。

[2024XYCTF]Complex_rsa

from Crypto.Util.number import *
from secrets import flag


class Complex:
    def __init__(self, re, im):
        self.re = re
        self.im = im

    def __mul__(self, c):
        re_ = self.re * c.re - self.im * c.im
        im_ = self.re * c.im + self.im * c.re
        return Complex(re_, im_)

    def __str__(self):
        if self.im == 0:
            return str(self.re)
        elif self.re == 0:
            if abs(self.im) == 1:
                return f"{'-' if self.im < 0 else ''}i"
            else:
                return f"{self.im}i"
        else:
            return f"{self.re} {'+' if self.im > 0 else '-'} {abs(self.im)}i"


def complex_pow(c, exp, n):
    result = Complex(1, 0)
    while exp > 0:
        if exp & 1:
            result = result * c
            result.re = result.re % n
            result.im = result.im % n
        c = c * c
        c.re = c.re % n
        c.im = c.im % n
        exp >>= 1
    return result


m = bytes_to_long(flag)
key = getRandomNBitInteger(m.bit_length())
c = m ^ key
com = Complex(key, c)
p = getPrime(512)
q = getPrime(512)
e = 9
enc = complex_pow(com, e, p * q)
print(enc)
print(Complex(p, q) * Complex(p, q))
# 66350931528185981323649477263900844564494528747802437244889229343520648562164950914433406604580038018765783183569276743239888668912948977370163046257917321742455772852779551569446155827368453262479370103326286297164105599131090881306108546341785251895116423206455175290083968296281375908109039893280371271943 + 65266730684129269806656018828265187384002656633231286337130178390517924611751697965395744944541329793503617856896706439809241745206839328124348693486741656130890593895436661857688522977866438805549144904296596887218275440542852887148071837153436265722614658566275517205322945316112048487893204059562830581004i
# -28814875173103880290298835537218644402443395484370652510062722255203946330565951328874411874019897676900075613671629765922970689802650462822117767927082712245512492082864958877932682404829188622269636302484189627580600076246836248427780151681898051243380477561480415972565859837597822263289141887833338111120 + 235362412848885579543400940934854106052672292040465052424316433330114813432317923674803623227280862945857543620663672974955235166884830751834386990766053503640556408758413592161122945636548462064584183165189050320898315823173824074873376450569212651128285746330837777597290934043912373820690250920839961482862i

$(p+qi)^2=p^2-q^2+2pqi$

那么可以通过sympy库解方程得到p,q

p = 10205509456040823453782883291824829705816298450992336115902525676510986341532486274067943978039013674207011185602314114359146043975207543018267242312660911
q = 11531144714660489617244423924607904114688972598683439999377362467380544567879231460623526800789918614728790840508257630983753525432337178000220918002499321

$enc\equiv com^9\mod p$

在指数方法，因为gcd(e,p**2-1)=3，不能求逆元，而且e与其不存在互素部分，不能进行降次，题解中提到构造如下式子

$9t\equiv 3\mod p^2-1$

证明其正确性，因为p**2-1是三的倍数，模数必然是整数

$9t-3=k(p^2-1)\\ 3t-1=\frac{k(p^2-1)}{3}\\ t\equiv 3^{-1}\mod \frac{p^2-1}{3}$

那么式子可以化成

$com^3\equiv enc^{t}\mod p$

求三次方比求九次方方便多了

$(k+mi)^2=(k^2-m^2)+2kmi\\ (k+mi)^3=(k^3-3km^2)+(3k^2m-m^3)i$

接下来也是结式的运算

from Crypto.Util.number import *
from gmpy2 import *

e = 9
pq2 = 235362412848885579543400940934854106052672292040465052424316433330114813432317923674803623227280862945857543620663672974955235166884830751834386990766053503640556408758413592161122945636548462064584183165189050320898315823173824074873376450569212651128285746330837777597290934043912373820690250920839961482862
p2_q2 = -28814875173103880290298835537218644402443395484370652510062722255203946330565951328874411874019897676900075613671629765922970689802650462822117767927082712245512492082864958877932682404829188622269636302484189627580600076246836248427780151681898051243380477561480415972565859837597822263289141887833338111120

class Complex:
    def __init__(self, re, im):
        self.re = re
        self.im = im

    def __mul__(self, c):
        re_ = self.re * c.re - self.im * c.im
        im_ = self.re * c.im + self.im * c.re
        return Complex(re_, im_)

    def __str__(self):
        if self.im == 0:
            return str(self.re)
        elif self.re == 0:
            if abs(self.im) == 1:
                return f"{'-' if self.im < 0 else ''}i"
            else:
                return f"{self.im}i"
        else:
            return f"{self.re} {'+' if self.im > 0 else '-'} {abs(self.im)}i"


def complex_pow(c, exp, n):
    result = Complex(1, 0)
    while exp > 0:
        if exp & 1:
            result = result * c
            result.re = result.re % n
            result.im = result.im % n
        c = c * c
        c.re = c.re % n
        c.im = c.im % n
        exp >>= 1
    return result

p = 10205509456040823453782883291824829705816298450992336115902525676510986341532486274067943978039013674207011185602314114359146043975207543018267242312660911
q = pq2 // 2 // p

com = Complex(66350931528185981323649477263900844564494528747802437244889229343520648562164950914433406604580038018765783183569276743239888668912948977370163046257917321742455772852779551569446155827368453262479370103326286297164105599131090881306108546341785251895116423206455175290083968296281375908109039893280371271943 , 65266730684129269806656018828265187384002656633231286337130178390517924611751697965395744944541329793503617856896706439809241745206839328124348693486741656130890593895436661857688522977866438805549144904296596887218275440542852887148071837153436265722614658566275517205322945316112048487893204059562830581004)
re = int(complex_pow(com, inverse(3,(p^2-1)//3), p).re)
im = int(complex_pow(com, inverse(3,(p^2-1)//3), p).im)


PR.<a,b> = PolynomialRing(Zmod(p))
f1 = a^3 - 3*a*b^2 - re
f2 = 3*a^2*b - b^3 - im
h = f1.sylvester_matrix(f2, a).det()
res1 = h.univariate_polynomial().monic().roots()
print(res1)


re = int(complex_pow(com, inverse(3,(q^2-1)//3), q).re)
im = int(complex_pow(com, inverse(3,(q^2-1)//3), q).im)
PR.<a,b> = PolynomialRing(Zmod(q))
f1 = a^3 - 3*a*b^2 - re
f2 = 3*a^2*b - b^3 - im
h = f1.sylvester_matrix(f2, a).det()
res1 = h.univariate_polynomial().monic().roots()
print(res1)

[GeekChallenge-2023]EzComplex

from Crypto.Util.number import *
flag = b'FAKE{Do_You_konw_Complex_numbers}'
p = random_prime(1 << 384)
q = random_prime(1 << 384)
n = p * q
e = 0x10001
N = pow(p, 2) + pow(q, 2)
m = bytes_to_long(flag)
c = pow(m, e, n)

print(c)
print(N)
'''
122977267154486898127643454001467185956864368276013342450998567212966113302012584153291519651365278888605594000436279106907163024162771486315220072170917153855370362692990814276908399943293854077912175867886513964032241638851526276
973990451943921675425625260267293227445098713194663380695161260771362036776671793195525239267004528550439258233703798932349677698127549891815995206853756301593324349871567926792912475619794804691721625860861059975526781239293017498
'''

$N=p^2+q^2=(p+qi)*(p-qi)$

在复数域上可以看作两个复数相乘，那么可以用高斯整数分解，试过$two\_squares$跑不出正确解，跑出来的两个数并不是质数，只能用高斯整数分解方法

c = 122977267154486898127643454001467185956864368276013342450998567212966113302012584153291519651365278888605594000436279106907163024162771486315220072170917153855370362692990814276908399943293854077912175867886513964032241638851526276
N = 973990451943921675425625260267293227445098713194663380695161260771362036776671793195525239267004528550439258233703798932349677698127549891815995206853756301593324349871567926792912475619794804691721625860861059975526781239293017498
GI = GaussianIntegers()
d = divisors(GI(N))
for k in d:
    if 'I' in str(k):
        p = int(abs(k.imag()))
        q = int(abs(k.real()))
        if p.bit_length() in range(382,385) and is_prime(p) and q.bit_length() in range(382,385) and is_prime(q):
            print((p,q))

p,q = (8732781022306464325787401448517171026218291389436971731700810979177651389459896422549428444142746055523338740248707, 29962125885196559918101088622575501736433575381042696980660846307183241725227137854663856022170515177120773072848343)
e = 0x10001
c = 122977267154486898127643454001467185956864368276013342450998567212966113302012584153291519651365278888605594000436279106907163024162771486315220072170917153855370362692990814276908399943293854077912175867886513964032241638851526276
f = (p-1)*(q-1)
d = inverse_mod(e,f)
m = pow(c,d,p*q)
print(bytes.fromhex(hex(m)[2:]))