非素数模下的二次剩余

非素数模下的二次剩余

$Instruction$

二次剩余定义为，一个数$a$，如果不是$p$的倍数且模$p$同余于某个数的平方，则称 $a$为模$p$的二次剩余
$$x^2\equiv a(\mathrm{mod}p)$$
对于普通的二次剩余来说，要应用正常求解二次剩余的方法有很多种，但前提都是模数为素数，当模数为素数时，可应用求解的方法有Atkin 算法，Cipolla 算法，Legendre 算法，Tonelli-Shanks 算法

算法的具体内容可见：二次剩余 - OI Wiki (oi-wiki.org)

另外关于Tonelli-Shanks算法的Python实现：(10条消息) Tonelli-Shanks算法_python_M3ng@L的博客-CSDN博客

---

但是这里我们并不重点看模数为素数的情况，而是看模数为非素数的情况，条件是在整数域的情况下

若模数不是素数，那么必为合数；将模数按质因数分解为素数后，把原来的同余式变换成同余式组；求解一元一次同余式组当然是应用孙子定理，但是这里并不是一元一次同余式组$\{\begin{array}{rl}& x\equiv a_1(\mathrm{mod}{p}_1)\\ & x\equiv a_2(\mathrm{mod}{p}_2)\\ & \cdots \\ & x\equiv a_n(\mathrm{mod}{p}_n)\end{array}$，而是一元二次同余式组$\{\begin{array}{rl}& x^2\equiv a_1(\mathrm{mod}{p}_1)\\ & x^2\equiv a_2(\mathrm{mod}{p}_2)\\ & \cdots \\ & x^2\equiv a_n(\mathrm{mod}{p}_n)\end{array}$

所以要先求解同余式组中每个同余式的二次剩余，得到一元二次同余式组中的每个同余式的解$x_i$，再代入写作新的一元一次同余式组$\{\begin{array}{rl}& x\equiv x_1(\mathrm{mod}{p}_1)\\ & x\equiv x_2(\mathrm{mod}{p}_2)\\ & \cdots \\ & x\equiv x_n(\mathrm{mod}{p}_n)\end{array}$，应用孙子定理求解出满足整个同余式组的$x$；至此，对于非素数模下的二次剩余就求解完成了

另外有两个点需要注意，一是将初始$p$分解质因数之后很可能是以$p_1^{q_1}\cdot p_2^{q_2}\cdot \cdots \cdot p_n^{q_n}$的形式存在，但是按$p_i^{q_i}$分隔开来，也满足孙子定理的条件：模数之间两两互素。二是求解一元二次同余式组中的各个二次剩余，单个二次剩余可能会有多个解满足该二次同余式，那么就将所有解依次作为该同余式的唯一解构造新的一元一次同余式组求解

$Example$

$Problem$

来源于Crypto CTF 2022 starter_ecc

from Crypto.Util.number import *
from secret import n, a, b, x, flag

y = bytes_to_long(flag.encode('utf-8'))

assert y < n
E = EllipticCurve(Zmod(n), [a, b])

try:
    G = E(x, y)
    print(f'x = {x}')
    print(f'a = {a}')
    print(f'b = {b}')
    print(f'n = {n}')
    print('Find the flag :P')
except:
    print('Ooops, ERROR :-(')
"""
x = 10715086071862673209484250490600018105614048117055336074437503883703510511249361224931983788156958581275946729175531468251871452856923140435984577574698574803934567774824230985421074605062371141877954182153046477020617917601884853827611232355455223966039590143622792803800879186033924150173912925208583
a = 31337
b = 66826418568487077181425396984743905464189470072466833884636947306507380342362386488703702812673327367379386970252278963682939080502468506452884260534949120967338532068983307061363686987539408216644249718950365322078643067666802845720939111758309026343239779555536517718292754561631504560989926785152983649035
n = 117224988229627436482659673624324558461989737163733991529810987781450160688540001366778824245275287757373389887319739241684244545745583212512813949172078079042775825145312900017512660931667853567060810331541927568102860039898116182248597291899498790518105909390331098630690977858767670061026931938152924839936
Find the flag :P
"""

$Analysis$

检验n发现是非素数，而flag作为y位于椭圆曲线上，其中已知椭圆曲线的参数以及y对应的横坐标x；根据函数EllipticCurve()传入的参数可知构造的椭圆曲线方程为
$$y^2\equiv x^3+ax+b(\mathrm{mod}n)$$
那么相当于同余式右侧已知，求左侧的$y$值大小；也就是一个非素数模下的二次剩余$x^2\equiv a(\mathrm{mod}n)$，按照之前分析非素数模下的二次剩余求解即可

先将n分解为质因数乘积的形式

factors_of_n = ecm.factor(n)
# 也可以用factor()
factors_of_n = factor(n)
"""
2^63 * 690712633549859897233^6 * 651132262883189171676209466993073^5
"""

可以使用ecm.factor()的原因是在特定条件下（其中一个因子大小大约为80bits 或者说是 25个十进制数字），ta的质因数分解速度是最快的

接着将原来的二次同余式拆分成一元二次同余式组$\{\begin{array}{rl}& x^2\equiv a_1(\mathrm{mod}{p}_1)\\ & x^2\equiv a_2(\mathrm{mod}{p}_2)\\ & \cdots \\ & x^2\equiv a_n(\mathrm{mod}{p}_n)\end{array}$，单独对每个同余式进行二次剩余求解

from sympy.ntheory.residue_ntheory import _sqrt_mod_prime_power, _sqrt_mod1
x = 10715086071862673209484250490600018105614048117055336074437503883703510511249361224931983788156958581275946729175531468251871452856923140435984577574698574803934567774824230985421074605062371141877954182153046477020617917601884853827611232355455223966039590143622792803800879186033924150173912925208583
a = 31337
b = 66826418568487077181425396984743905464189470072466833884636947306507380342362386488703702812673327367379386970252278963682939080502468506452884260534949120967338532068983307061363686987539408216644249718950365322078643067666802845720939111758309026343239779555536517718292754561631504560989926785152983649035
n = 117224988229627436482659673624324558461989737163733991529810987781450160688540001366778824245275287757373389887319739241684244545745583212512813949172078079042775825145312900017512660931667853567060810331541927568102860039898116182248597291899498790518105909390331098630690977858767670061026931938152924839936

factors_of_n = {2:63, 690712633549859897233:6, 651132262883189171676209466993073:5}
A = x**3 + a*x + b 

# 之后的x是一元二次同余式组中不同二次剩余的解的结合x
x = []  
p = []  

for p_i, e_i in factors_of_n.items():  
    if A % p_i == 0:  
        x_i = _sqrt_mod1(A, p_i, e_i)  # find solution to ``x**2 == a mod p**n`` when ``a % p == 0``
    else:  
        x_i = _sqrt_mod_prime_power(A, p_i, e_i)  # find the solutions to ``x**2 = a mod p**k`` when ``a % p != 0``
    x.append(x_i)  
    p.append(p_i**e_i)

求得每个二次剩余的解$x_i$，将各个$x_i$重新构造出一个新的同余式组$\{\begin{array}{rl}& x\equiv x_1(\mathrm{mod}{p}_1)\\ & x\equiv x_2(\mathrm{mod}{p}_2)\\ & \cdots \\ & x\equiv x_n(\mathrm{mod}{p}_n)\end{array}$；但是这里一个二次剩余会产生多个解，也就是$x_{i_j}$，那么按照之前提到的解决办法，将每个二次剩余的解都先看作是唯一解构造新的同余式组，直到将每个二次剩余的解按照不同组合遍历完即可，求得的结果设置条件判断输出

from sympy.polys.galoistools import gf_crt
from sympy.polys.domains import ZZ
from Crypto.Util.number import *

...

for combin_x in _product(*x):
    r = gf_crt(combin_x, p, ZZ)
    flag = long_to_bytes(int(r))
    if b'CTF' in flag:
        print(flag)
        break

$Solution$

# type: ignore
from sympy.ntheory.residue_ntheory import _sqrt_mod_prime_power, _product
from sympy.polys.galoistools import gf_crt1, gf_crt2, gf_crt
from sympy.polys.domains import ZZ
from Crypto.Util.number import *

x = 10715086071862673209484250490600018105614048117055336074437503883703510511249361224931983788156958581275946729175531468251871452856923140435984577574698574803934567774824230985421074605062371141877954182153046477020617917601884853827611232355455223966039590143622792803800879186033924150173912925208583
a = 31337
b = 66826418568487077181425396984743905464189470072466833884636947306507380342362386488703702812673327367379386970252278963682939080502468506452884260534949120967338532068983307061363686987539408216644249718950365322078643067666802845720939111758309026343239779555536517718292754561631504560989926785152983649035
n = 117224988229627436482659673624324558461989737163733991529810987781450160688540001366778824245275287757373389887319739241684244545745583212512813949172078079042775825145312900017512660931667853567060810331541927568102860039898116182248597291899498790518105909390331098630690977858767670061026931938152924839936


factors_of_n = {2:63, 690712633549859897233:6, 651132262883189171676209466993073:5}
A = x**3 + a*x + b 

# 之后的x是一元二次同余式组中不同二次剩余的解的结合x
x = []  
p = []  

for p_i, e_i in factors_of_n.items():  
    if A % p_i == 0:  
        x_i = _sqrt_mod1(A, p_i, e_i)  # find solution to ``x**2 == a mod p**n`` when ``a % p == 0``
    else:  
        x_i = _sqrt_mod_prime_power(A, p_i, e_i)  # find the solutions to ``x**2 = a mod p**k`` when ``a % p != 0``
    x.append(x_i)  
    p.append(p_i**e_i)

for combin_x in _product(*x):
    r = gf_crt(combin_x, p, ZZ)
    flag = long_to_bytes(int(r))
    if b'CTF' in flag:
        print(flag)
        break
# # 下面这种方法也可以求解，将CRT分成两部分解决，具体使用方法可见http://lidavidm.github.io/sympy/_modules/sympy/polys/galoistools.html
# mm, e, s = gf_crt1(pv ,ZZ)
# for vx in _product(*v):
#     r = gf_crt2(vx, pv, mm, e, s, ZZ)
#     flag = long_to_bytes(int(r))
#     if b'CTF' in flag:
#         print(flag)
#         break

$Reference$

CTFtime.org / Crypto CTF 2022 / starter_ecc / Writeup