背景

两个月前的DASCTF2022大爆零，RSA的题目一个下午了只做出来了拆解二平方和一个步骤。恰好前几天UNCTF的群里有人问类似的题目，故在此稍作归纳。

以DASCTF2022-easyrsa为例，观察前几行代码：

1
2
3
4
5
6
7
8
from gmpy2 import *
from Crypto.Util.number import *

p1 = getPrime(128)
q = getPrime(128)
n = p1**2+q**2
print('n = ',n)
# n = 86073852484226203700520112718689325205597071202320413471730820840719099334770

某乎上真是什么问题都有：一个整数可以拆成两个整数的平方和，5201314可以拆成哪两个数的平方和？

定理1和定理2非常重要，它们把这个题目简化为了以下两个步骤：

  1. 分解质因数
  2. 把每个质因数拆解成二平方和
  3. 合并

接下来介绍每一步的具体实现：

分解质因数

用yafu或者factordb。

4k+1型质因数的二平方和拆解

重复定理1：

对素数 p>2p>2，pp 能表示成两个整数的平方和当且仅当 p≡1 (mod 4)p \equiv 1 \space (mod \space 4).

因此，第一步分解出的所有质因数都必须是4k+1型的方可求解。

这里给出一个4k+1型质因数的二平方和拆解方案。

脚本如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
import gmpy2

def lim(a,m): # 限定a在[-m/2, m/2]范围内
    if a>m//2:
        a -= m
    elif a<-(m//2):
        a+=m
    return a

def prime_solver(p):
    if p == 2: # 对p==2的情况进行特判
        return 1,1
    a = 2
    f = 0
    while not f: # 二次剩余求解
        m = (p-1)//2
        if pow(a,m,p)>1:
            f = 1
            x=pow(a,(m//2),p)
        else:
            a = int(gmpy2.next_prime(a))
    a = x
    b = 1
    while 1: # 求a,b
        m1 = (a**2+b**2)//p
        if m1 == 1:
            break
        u = a % m1
        v = b % m1
        v = lim(v,m1)
        u = lim(u,m1)
        a1 = (u*a+v*b)//m1
        b1 = (u*b-v*a)//m1
        a = a1
        b = b1
    assert a*a+b*b == p
    return abs(a),abs(b)

合并

根据定理2：

(a2+b2)(c2+d2)=(ac+bd)2+(ad−bc)2(a^2+b^2)(c^2+d^2)=(ac+bd)^2+(ad-bc)^2

可以将第二步分解的结果依次合并。然而，合并时可能出现多解，我们需要用深搜和约束条件找出我们需要的那一组解。

脚本如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
import gmpy2

def judge(n):# 最终分解出的p,q的约束条件
    return gmpy2.is_prime(n) or gmpy2.is_prime(n)

def search(p,q,cnt): # 深搜求解
    if cnt == len(factors):
        if judge(p) and judge(q):
            print(p,q)
        return
    search(p*roots[cnt][0]+q*roots[cnt][1],abs(p*roots[cnt][1]-q*roots[cnt][0]),cnt+1)
    search(p*roots[cnt][1]+q*roots[cnt][0],abs(p*roots[cnt][0]-q*roots[cnt][1]),cnt+1)

factors = []
roots = []

for f in factors:
    if f%4 == 1 or f == 2:
        roots.append(prime_solver(f))
    else:
        raise Exception("4k+3!")
print(roots)
search(roots[0][0],roots[0][1],1)
search(roots[0][1],roots[0][0],1)

总结

汇总上面的程序得出的最终脚本如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
import gmpy2

def lim(a,m): # 限定a在[-m/2, m/2]范围内
    if a>m//2:
        a -= m
    elif a<-(m//2):
        a+=m
    return a

def prime_solver(p):
    if p == 2:
        return 1,1
    a = 2
    f = 0
    while not f: # 二次剩余求解
        m = (p-1)//2
        if pow(a,m,p)>1:
            f = 1
            x=pow(a,(m//2),p)
        else:
            a = int(gmpy2.next_prime(a))
    a = x
    b = 1
    while 1: # 求a,b
        m1 = (a**2+b**2)//p
        if m1 == 1:
            break
        u = a % m1
        v = b % m1
        v = lim(v,m1)
        u = lim(u,m1)
        a1 = (u*a+v*b)//m1
        b1 = (u*b-v*a)//m1
        a = a1
        b = b1
    assert a*a+b*b == p
    return abs(a),abs(b)

def judge(n):# 最终分解出的p,q的约束条件
    return gmpy2.is_prime(n) or gmpy2.is_prime(n)

def search(p,q,cnt): # 深搜求解
    if cnt == len(factors):
        if judge(p) and judge(q):
            print(p,q)
        return
    search(p*roots[cnt][0]+q*roots[cnt][1],abs(p*roots[cnt][1]-q*roots[cnt][0]),cnt+1)
    search(p*roots[cnt][1]+q*roots[cnt][0],abs(p*roots[cnt][0]-q*roots[cnt][1]),cnt+1)
    
# factordb.com
factors = [2,5,2341,191429240956028917,18808725088317359977,1021179744308298062426500817712456533]
roots = []

for f in factors:
    if f%4 == 1 or f == 2:
        roots.append(prime_solver(f))
    else:
        raise Exception("4k+3!")
print(roots)
search(roots[0][0],roots[0][1],1)
search(roots[0][1],roots[0][0],1)

运行后得出p1和q的解。