I participated in zer0pts CTF in Super Guesser, and we reached 1st place.

Since there are many challenges to cover, this writeup will be very concise.

The solution codes may be a bit messy, since I wrote this in a hurry :P

Also, r4j solved signme, so I'll leave that writeup to r4j :D

 

war(sa)mup

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
from Crypto.Util.number import getStrongPrime, GCD
from random import randint
from flag import flag
import os
 
def pad(m: int, n: int):
  # PKCS#1 v1.5 maybe
  ms = m.to_bytes((m.bit_length() + 7// 8"big")
  ns = n.to_bytes((n.bit_length() + 7// 8"big")
  assert len(ms) <= len(ns) - 11
 
  ps = b""
  while len(ps) < len(ns) - len(ms) - 3:
    p = os.urandom(1)
    if p != b"\x00":
      ps += p
  return int.from_bytes(b"\x00\x02" + ps + b"\x00" + ms, "big")
 
 
while True:
  p = getStrongPrime(512)
  q = getStrongPrime(512)
  n = p * q
  phi = (p-1)*(q-1)
  e = 1337
  if GCD(phi, e) == 1:
    break
 
= pad(int.from_bytes(flag, "big"), n)
c1 = pow(m, e, n)
c2 = pow(m // 2, e, n)
 
print("n =", n)
print("e =", e)
print("c1=", c1)
print("c2=", c2)
 
cs

 

Solution

Writing $\lfloor m/2 \rfloor = x$, we see that $m$ is either $2x$ or $2x+1$. If $m = 2x$, we can perform the related message attack with polynomials $$p_1(x) = (2x)^e - c_1, \quad p_2(x) = x^e - c_2$$ Of course, in this case we would have $c_1 \equiv 2^e c_2 \pmod{n}$ and polynomial GCD would be useless. If $m = 2x+1$, we can do it with $$p_1(x) = (2x+1)^e - c_1, \quad p_2(x) = x^e - c_2$$ which works and gives our desired value of $x$. $m$ can be recovered as $2x+1 \pmod{n}$.

 

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
# sage
= 113135121314210337963205879392132245927891839184264376753001919135175107917692925687745642532400388405294058068119159052072165971868084999879938794441059047830758789602416617241611903275905693635535414333219575299357763227902178212895661490423647330568988131820052060534245914478223222846644042189866538583089
= 1337
c1 = 89077537464844217317838714274752275745737299140754457809311043026310485657525465380612019060271624958745477080123105341040804682893638929826256518881725504468857309066477953222053834586118046524148078925441309323863670353080908506037906892365564379678072687516738199061826782744188465569562164042809701387515
c2 = 18316499600532548540200088385321489533551929653850367414045951501351666430044325649693237350325761799191454032916563398349042002392547617043109953849020374952672554986583214658990393359680155263435896743098100256476711085394564818470798155739552647869415576747325109152123993105242982918456613831667423815762
 
P.<x> = PolynomialRing(Zmod(n))
 
def GCD(f, g, n):
    g = g % f
    if g == 0:
        return f
    t = g.lc()
    if gcd(t, n) != 1:
        print(t)
        exit()
    tt = inverse_mod(Integer(t), n)
    g = g * tt
    return GCD(g, f, n)
 
= (2 * x + 1)^e - c1
= x^e - c2
 
print(GCD(f, g, n))
 
# back to python
= -113131683468284119607964196541575421728552328779560800910707074969146190674728631210678514366174908879152073378298687396133007061413453069629241334690374397179449596372765780570923850759894857943624531609494119853246060991087070835093327088258517606688574536921269027653158032526646833175962482267255713310835
= 113135121314210337963205879392132245927891839184264376753001919135175107917692925687745642532400388405294058068119159052072165971868084999879938794441059047830758789602416617241611903275905693635535414333219575299357763227902178212895661490423647330568988131820052060534245914478223222846644042189866538583089
= 1337
c1= 89077537464844217317838714274752275745737299140754457809311043026310485657525465380612019060271624958745477080123105341040804682893638929826256518881725504468857309066477953222053834586118046524148078925441309323863670353080908506037906892365564379678072687516738199061826782744188465569562164042809701387515
c2= 18316499600532548540200088385321489533551929653850367414045951501351666430044325649693237350325761799191454032916563398349042002392547617043109953849020374952672554986583214658990393359680155263435896743098100256476711085394564818470798155739552647869415576747325109152123993105242982918456613831667423815762
 
print(n+m)
 
fin = (2 * (n+m) + 1) % n
print(long_to_bytes(fin))
cs

 

OT or NOT OT

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
import os
import signal
import random
from base64 import b64encode
from Crypto.Util.number import getStrongPrime, bytes_to_long
from Crypto.Util.Padding import pad
from Crypto.Cipher import AES
from flag import flag
 
= getStrongPrime(1024)
 
key = os.urandom(32)
iv = os.urandom(AES.block_size)
aes = AES.new(key=key, mode=AES.MODE_CBC, iv=iv)
= aes.encrypt(pad(flag, AES.block_size))
 
key = bytes_to_long(key)
print("Encrypted flag: {}".format(b64encode(iv + c).decode()))
print("p = {}".format(p))
print("key.bit_length() = {}".format(key.bit_length()))
 
signal.alarm(600)
while key > 0:
    r = random.randint(2, p-1)
    s = random.randint(2, p-1)
    t = random.randint(2, p-1)
    print("t = {}".format(t))
 
    a = int(input("a = ")) % p
    b = int(input("b = ")) % p
    c = int(input("c = ")) % p
    d = int(input("d = ")) % p
    assert all([a > 1 , b > 1 , c > 1 , d > 1])
    assert len(set([a,b,c,d])) == 4
 
    u = pow(a, r, p) * pow(c, s, p) % p
    v = pow(b, r, p) * pow(c, s, p) % p
    x = u ^ (key & 1)
    y = v ^ ((key >> 1& 1)
    z = pow(d, r, p) * pow(t, s, p) % p
 
    key = key >> 2
 
    print("x = {}".format(x))
    print("y = {}".format(y))
    print("z = {}".format(z))
 
cs

 

Solution

The key idea is to make a small subgroup which the values of $x, y \in \mathbb{F}_p$ must belong to.

To do this, we wait until $p \equiv 1 \pmod{4}$, then take any $t \in \mathbb{F}_p$ with order $4$.

This can be done by taking some random $g \in \mathbb{F}_p$ and choosing $t = g^{(p-1)/4}$. 

Checking if order of $t$ is $4$ can now be done simply by verifying $t^2 \neq 1$ in $\mathbb{F}_p$.

 

Now we simply send $a = t$, $b = t^2$, $c = t^3$. Since $a^4 = b^4 = c^4 = 1$, we also have $u^4 = v^4 = 1$.

Therefore, we can see if $x = u$ and if $y= v$ by simply checking if $x^4=1$ and if $y^4=1$.

 

Of course, it should be noted that getStrongPrime does not generate a prime of the form $p=2q+1$ with $q$ prime. :wink:

 

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
= remote('crypto.ctf.zer0pts.com''10130')
 
def recvint():
    s = r.recvline()[:-1]
    s = s.split()[-1]
    return int(s)
 
= r.recvline()[:-1]
= S.split()[-1]
= base64.b64decode(S)
iv = V[:16]
ctxt = V[16:]
 
= recvint()
keylen = recvint()
 
if p % 4 != 1:
    print("BAD PRIME")
    exit()
 
= 0
for g in range(22000):
    t = pow(g, (p-1)//4, p)
    if (t * t) % p != 1:
        break
 
print("Begin Key Finding")
keyv = 0
add = 1
for i in tqdm(range(0128)):
    t = recvint()
    r.sendline(str(t))
    r.sendline(str((t * t) % p))
    r.sendline(str((t * t * t) % p))
    r.sendline(str(5))
    x = recvint()
    y = recvint()
    z = recvint()
    if pow(x, 4, p) != 1:
        keyv += add
    add = add * 2
    if pow(y, 4, p) != 1:
        keyv += add
    add = add * 2
    if i >= 126:
        keyv = long_to_bytes(keyv)
 
        aes = AES.new(key=keyv, mode=AES.MODE_CBC, iv=iv)
        ptxt = aes.decrypt(ctxt)
        print(ptxt)
cs

 

janken vs yoshiking

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
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
import random
import signal
from flag import flag
from Crypto.Util.number import getStrongPrime, inverse
 
HANDNAMES = {
    1"Rock",
    2"Scissors",
    3"Paper"
}
 
def commit(m, key):
    (g, p), (x, _) = key
    r = random.randint(2, p-1)
    c1 = pow(g, r, p)
    c2 = m * pow(g, r*x, p) % p
    return (c1, c2)
 
 
def decrypt(c, key):
    c1, c2 = c
    _, (x, p)= key
 
    m = c2 * inverse(pow(c1, x, p), p) % p
    return m
 
 
def keygen(size):
    p = getStrongPrime(size)
    g = random.randint(2, p-1)
    x = random.randint(2, p-1)
 
    return (g, p), (x, p)
 
 
signal.alarm(3600)
key = keygen(1024)
(g, p), _ = key
print("[yoshiking]: Hello! Let's play Janken(RPS)")
print("[yoshiking]: Here is g: {}, and p: {}".format(g, p))
 
round = 0
wins = 0
while True:
    round += 1
    print("[system]: ROUND {}".format(round))
 
    yoshiking_hand = random.randint(13)
    c = commit(yoshiking_hand, key)
    print("[yoshiking]: my commitment is={}".format(c))
 
    hand = input("[system]: your hand(1-3): ")
    print("")
    try:
        hand = int(hand)
        if not (1 <= hand <= 3):
            raise ValueError()
    except ValueError:
        print("[yoshiking]: Ohhhhhhhhhhhhhhhh no! :(")
        exit()
 
    yoshiking_hand = decrypt(c, key)
    print("[yoshiking]: My hand is ... {}".format(HANDNAMES[yoshiking_hand]))
    print("[yoshiking]: Your hand is ... {}".format(HANDNAMES[hand]))
    result = (yoshiking_hand - hand + 3) % 3
    if result == 0:
        print("[yoshiking]: Draw, draw, draw!!!")
    elif result == 1:
        print("[yoshiking]: Yo! You win!!! Ho!")
        wins += 1
        print("[system]: wins: {}".format(wins))
 
        if wins >= 100:
            break
    elif result == 2:
        print("[yoshiking]: Ahahahaha! I'm the winnnnnnner!!!!")
        print("[yoshiking]: You, good loser!")
        print("[system]: you can check that yoshiking doesn't cheat")
        print("[system]: here's the private key: {}".format(key[1][0]))
        exit()
 
print("[yoshiking]: Wow! You are the king of roshambo!")
print("[yoshiking]: suge- flag ageru")
print(flag)
 
cs

 

Solution

The first observation is that if $g$ is a quadratic residue (QR), the commitment leaks the legendre symbol of $m$.

Indeed, if $g$ is QR, then the legendre symbol of $m$ is simply the legendre symbol of $c_2$.

 

We generate instances until we find a prime $p$ that not all $1, 2, 3$ are QR modulo $p$. 

In this case, by leaking the legendre symbol of $m$, we can reduce the number of candidates for yoshiking's hand from 3 to 2.

Therefore, we can at least force a tie, and occasionally win. This is enough to solve the challenge.

 

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
63
64
65
66
67
= remote('crypto.ctf.zer0pts.com''10463')
r.recvline()
 
= r.recvline()
= s.split()
= int(T[4][:-1])
= int(T[-1].rstrip())
 
 
QR = []
NQR = []
 
if pow(g, (p-1)//2 , p) != 1:
    print("BAD PRIME 1")
    exit()
QR.append(1)
if pow(2, (p-1)//2, p-1== 1:
    QR.append(2)
else:
    NQR.append(2)
if pow(3, (p-1)//2, p-1== 1:
    QR.append(3)
else:
    NQR.append(3)
 
if len(QR) == 3:
    print("lol")
    exit()
 
def get_commit():
    s = r.recvline()
    t = s.split(b'(')[-1]
    c1 = int(t.split(b',')[0])
    c2 = int(t.split(b',')[1][1:-2])
    return c1, c2
 
for i in range(0500):
    print(r.recvline())
    c1, c2 = get_commit()
    if pow(c2, (p-1// 2, p) == 1:
        if len(QR) == 1:
            r.sendline(b"3")
        else:
            if 2 in QR:
                r.sendline(b"1")
            else:
                r.sendline(b"3")
    else:
        if len(NQR) == 1:
            if 2 in NQR:
                r.sendline(b"1")
            if 3 in NQR:
                r.sendline(b"2")
        else:
            r.sendline(b"2")
    r.recvline()
    r.recvline()
    r.recvline()
    s = r.recvline()
    if s.split()[-1].rstrip() == b"draw!!!":
        continue
    else:
        s = r.recvline()
        print(s)
        if s.split()[-1].rstrip() == b"100":
            for t in range(04):
                print(r.recvline())
cs

 

easy pseudorandom

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
from Crypto.Util.number import*
from flag import flag
 
nbits = 256
= random_prime(1 << nbits)
Fp = Zmod(p)
P.<v> = PolynomialRing(Fp)
 
= randrange(p)
= 2
= v^2 + b
 
v0 = randrange(p)
v1 = F(v0)
 
= ceil(nbits * (d / (d + 1)))
w0 = (v0 >> (nbits - k))
w1 = (v1 >> (nbits - k))
 
# encrypt
= bytes_to_long(flag)
= 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}")
 
cs

 

Solution

It's clear that we just need to find the value of $v_0$. Denote $l = 256$. Using the given info, we can write $$\begin{equation*} \begin{aligned} v_0 &= w_0 \cdot 2^{l - k} + x \\ v_1 &= w_1 \cdot 2^{l-k} + y \\ v_1 & \equiv v_0^2 + b \pmod{n} \end{aligned} \end{equation*}$$ Combining, this gives us the key relation $$ w_1 \cdot 2^{l-k} + y \equiv w_0^2 \cdot 2^{2l-2k} + b + w_0 \cdot 2^{l-k+1} \cdot x + x^2 \pmod{n}$$ The key idea is to remove the quadratic term $x^2$. Denote $$\begin{equation*} \begin{aligned} A &= w_0 \cdot 2^{l-k+1} \\ B &= w_1 \cdot 2^{l-k} - w_0^2 \cdot 2^{2l-2k} - b \end{aligned} \end{equation*}$$ Note that $A, B$ are constants that we know. Our relation is now $$ Ax \equiv B + y - x^2 \pmod{n}$$ Since $x, y$ are all below $2^{l-k}$, we can see that $$ B - 2^{2l-2k} \le Ax \pmod{n} \le B + 2^{l-k}$$ which is a type of relation that can be solved. Check out my writeups for PBCTF Special Gift and SECCON crypto01.

 

Note that there are $2^{l-k}$ possible values of $x$, and there are approximately $2^{2l-2k}$ possible values of $Ax \pmod{n}$. 

Since $l-k + 2l-2k = 3l-3k \approx l$, we can decide that there must be fairly small number of $x$ that satisfy the desired inequality.

This type of analysis was done in both of writeups above, and recently highlighted by Mystiz in AeroCTF writeup.

 

Now that there are three challenges I solved with this method, I think I'll add this to my CVP repository as well...

 

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
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
# notation a bit different from solution, but same idea
 
# simplify polynomial in sage
= 86160765871200393116432211865381287556448879131923154695356172713106176601077
= 71198163834256441900788553646474983932569411761091772746766420811695841423780
= 88219145192729480056743197897921789558305761774733086829638493717397473234815
w0 = 401052873479535541023317092941219339820731562526505
w1 = 994046339364774179650447057905749575131331863844814
 
nbits = 256
P.<x, y> = PolynomialRing(GF(p))
= (nbits * 2 + 2// 3
 
= w0 * 2^(nbits - k) + x
= w1 * 2^(nbits - k) + y
 
= g - f^2 - b
print(h)
 
# python
 
def ceil(n, m):
    return (n + m - 1// m
 
def optf(A, M, L, R):
    if L == 0:
        return 0
    if 2 * A > M:
        L, R = R, L
        A = M - A
        L = M - L
        R = M - R
    cc_1 = ceil(L, A)
    if A * cc_1 <= R:
        return cc_1
    cc_2 = optf(A - M % A, A, L % A, R % A)
    return ceil(L + M * cc_2, A)
 
= 18867904637006146022735447
= 19342813113834066795298816
= 86160765871200393116432211865381287556448879131923154695356172713106176601077
= 71198163834256441900788553646474983932569411761091772746766420811695841423780
= 88219145192729480056743197897921789558305761774733086829638493717397473234815
w0 = 401052873479535541023317092941219339820731562526505
w1 = 994046339364774179650447057905749575131331863844814
 
C1 = 55130802749277213576496911760053178817655787149958046010477129311148596128757
C2 = 78083221913223461198494116323396529665894773452683783127339675579334647310194
 
nbits = 256
= (nbits * 2 + 2// 3
 
delt = nbits - k
# C1 x + y + C2 - x^2 == 0 mod p
# C1 x == x^2 - y - C2 mod p
 
= (0 - (1 << (delt)) - C2 + p) % p
= ((1 << (2 * delt)) - C2 + p) % p
 
lst = 0
while lst <= (1 << delt):
    NL = (L - C1 * (lst + 1)) % p
    NR = (R - C1 * (lst + 1)) % p
    if NL > NR:
        lst = lst + 1
    else:
        cc = optf(C1, p, NL, NR)
        lst = lst + 1 + cc
    mm = m
    x = lst
    v0 = w0 * (1 << (nbits - k)) + x
    v1 = (v0 * v0 + b) % p
    v = v1
    # try out!
    for i in range(5):
        v = (v * v + b) % p
        mm ^= v
    print(long_to_bytes(mm))
cs

 

3-AES

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
63
64
65
66
from Crypto.Cipher import AES
from Crypto.Random import get_random_bytes
from binascii import hexlify, unhexlify
from hashlib import md5
import os
import signal
from flag import flag
 
keys = [md5(os.urandom(3)).digest() for _ in range(3)]
 
 
def get_ciphers(iv1, iv2):
    return [
        AES.new(keys[0], mode=AES.MODE_ECB),
        AES.new(keys[1], mode=AES.MODE_CBC, iv=iv1),
        AES.new(keys[2], mode=AES.MODE_CFB, iv=iv2, segment_size=8*16),
    ]
 
def encrypt(m: bytes, iv1: bytes, iv2: bytes) -> bytes:
    assert len(m) % 16 == 0
    ciphers = get_ciphers(iv1, iv2)
    c = m
    for cipher in ciphers:
        c = cipher.encrypt(c)
    return c
 
def decrypt(c: bytes, iv1: bytes, iv2: bytes) -> bytes:
    assert len(c) % 16 == 0
    ciphers = get_ciphers(iv1, iv2)
    m = c
    for cipher in ciphers[::-1]:
        m = cipher.decrypt(m)
    return m
 
signal.alarm(3600)
while True:
    print("==== MENU ====")
    print("1. Encrypt your plaintext")
    print("2. Decrypt your ciphertext")
    print("3. Get encrypted flag")
    choice = int(input("> "))
 
    if choice == 1:
        plaintext = unhexlify(input("your plaintext(hex): "))
        iv1, iv2 = get_random_bytes(16), get_random_bytes(16)
        ciphertext = encrypt(plaintext, iv1, iv2)
        ciphertext = b":".join([hexlify(x) for x in [iv1, iv2, ciphertext]]).decode()
        print("here's the ciphertext: {}".format(ciphertext))
 
    elif choice == 2:
        ciphertext = input("your ciphertext: ")
        iv1, iv2, ciphertext = [unhexlify(x) for x in ciphertext.strip().split(":")]
        plaintext = decrypt(ciphertext, iv1, iv2)
        print("here's the plaintext(hex): {}".format(hexlify(plaintext).decode()))
 
    elif choice == 3:
        plaintext = flag
        iv1, iv2 = get_random_bytes(16), get_random_bytes(16)
        ciphertext = encrypt(plaintext, iv1, iv2)
        ciphertext = b":".join([hexlify(x) for x in [iv1, iv2, ciphertext]]).decode()
        print("here's the encrypted flag: {}".format(ciphertext))
        exit()
 
    else:
        exit()
 
cs

 

Solution

Let's consider plaintexts/ciphertexts of one block first. Denote $E_k$, $D_k$ as ECB encryption/decryption with key $k$.

The encryption process can be written as $$ P \rightarrow E_{k_1}(P) \rightarrow E_{k_2}(E_{k_1}(P) \oplus IV_1) \rightarrow E_{k_2}(E_{k_1}(P) \oplus IV_1) \oplus E_{k_3}(IV_2) = C$$ The key relation is $$E_{k_1}(P) \oplus IV_1 = D_{k_2}(C \oplus E_{k_3}(IV_2))$$ For a fixed value of $C, IV_2$, the right hand side is fixed. Therefore, the left hand side is also fixed.

 

We start by using the decryption oracle with $C, IV_1, IV_2$ being all 16 bytes of \x00. Denote $P_0$ as the resulting plaintext. 

Then, we use the decryption oracle with $C, IV_1, IV_2$ all \x00, but $IV_1$ has the last byte \x01. Denote $P_1$ as the resulting plaintext.

 

Our observations lead to $$\text{bytes_to_long}(E_{k_1}(P_0) \oplus E_{k_1}(P_1)) = 1$$ This can be tested for all $2^{24}$ possible keys, and with that we find $k_1$. 

 

Now that we know $k_1$, we also know the value of $D_{k_2} \circ E_{k_3}$ applied to 16 bytes of \x00.

This is a perfect scenario for a meet-in-the-middle attack, and with this we can find $k_2$, $k_3$. 

 

retrieve data from server

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
# get data from server
= remote('crypto.ctf.zer0pts.com''10929')
 
def get_enc(ptxt):
    r.recvuntil(b">")
    r.sendline(b"1")
    ptxt = bytes.hex(ptxt)
    r.sendline(ptxt)
    s = r.recvline().split()[-1].rstrip()
    iv1, iv2, ctxt = s.split(b":")
    iv1 = bytes.fromhex(iv1.decode())
    iv2 = bytes.fromhex(iv2.decode())
    ctxt = bytes.fromhex(ctxt.decode())
    return iv1, iv2, ctxt
 
def get_dec(ctxt, iv1, iv2):
    r.recvuntil(b">")
    r.sendline(b"2")
    ctxt = bytes.hex(ctxt)
    iv1 = bytes.hex(iv1)
    iv2 = bytes.hex(iv2)
    goal = iv1 + ":" + iv2 + ":" + ctxt
    r.sendline(goal)
    s = r.recvline()
    print(s)
    s = s.split()[-1].rstrip()
    ptxt = bytes.fromhex(s.decode())
    return ptxt
 
def get_flag():
    r.recvuntil(b">")
    r.sendline(b"3")
    s = r.recvline().split()[-1].rstrip()
    iv1, iv2, ctxt = s.split(b":")
    iv1 = bytes.fromhex(iv1.decode())
    iv2 = bytes.fromhex(iv2.decode())
    ctxt = bytes.fromhex(ctxt.decode())
    assert len(iv1) == 16
    assert len(iv2) == 16
    assert len(ctxt) == 48
    return iv1, iv2, ctxt
 
 
= open("lol.txt""w")
ptxt_0 = get_dec(b"\x00" * 16, b"\x00" * 16, b"\x00" * 16)
ptxt_1 = get_dec(b"\x00" * 16, b"\x00" * 15 + b"\x01", b"\x00" * 16)
 
 
f.write(str(bytes_to_long(ptxt_0)) + "\n")
f.write(str(bytes_to_long(ptxt_1)) + "\n")
 
iv1, iv2, ctxt = get_flag()
 
f.write(str(bytes_to_long(iv1)) + "\n")
f.write(str(bytes_to_long(iv2)) + "\n")
f.write(str(bytes_to_long(ctxt)) + "\n")
 
cs

 

find the first key

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
def bytexor(a, b):
    assert len(a) == len(b)
    return bytes(x ^ y for x, y in zip(a, b))
 
for i in tqdm(range(01 << 24)):
    cc = long_to_bytes(i)
    if len(cc) < 3:
        cc = b"\x00" * (3 - len(cc)) + cc
    k1 = hashlib.md5(cc).digest()
    cipher = AES.new(k1, mode = AES.MODE_ECB)
    val_1 = cipher.encrypt(ptxt_0)
    val_2 = cipher.encrypt(ptxt_1)
    det = bytexor(val_1, val_2)
    if bytes_to_long(det) == 1:
        print(i)
cs

 

meet in the middle preparation

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
# we now know the first key
cipher = AES.new(key = key_1, mode = AES.MODE_ECB)
vv = cipher.encrypt(ptxt_0)
 
= open("Data1.txt""w")
for i in tqdm(range(01 << 24)):
    cc = long_to_bytes(i)
    if len(cc) < 3:
        cc = b"\x00" * (3 - len(cc)) + cc
        assert bytes_to_long(cc) == i
    k2 = hashlib.md5(cc).digest()
    cipher = AES.new(key = k2, mode = AES.MODE_ECB)
    res = cipher.encrypt(vv)
    f.write(str(bytes_to_long(res)) + "\n")
 
for i in tqdm(range(01 << 24)):
    cc = long_to_bytes(i)
    if len(cc) < 3:
        cc = b"\x00" * (3 - len(cc)) + cc
    k3 = hashlib.md5(cc).digest()
    cipher = AES.new(key = k3, mode = AES.MODE_ECB)
    res = cipher.encrypt(b"\x00" * 16)
    f.write(str(bytes_to_long(res)) + "\n")
f.close()
cs

 

meet in the middle in C++

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
map<stringint> M;
 
int main(void)
{
    fio; int i, j;
    freopen("Data1.txt""r", stdin);
    for(i=0 ; i < (1<<24) ;  i++)
    {
        if(i % 1000000 == 0cout << i / 1000000 << endl;
        string s; cin >> s;
        M[s] = i;
    }
    for(i=0 ; i<(1<<24) ; i++)
    {
        if(i % 1000000 == 0cout << i / 1000000 << endl;
        string s; cin >> s;
        if(M[s] != 0)
        {
            cout << M[s] << " " << i << endl;
            return 0;
        }
    }
    return 0;
}
cs

 

NOT Mordell primes

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
from Crypto.Util.number import bytes_to_long
from secrets import k, FLAG
 
 
= 13046889097521646369087469608188552207167764240347195472002158820809408567610092324592843361428437763328630003678802379234688335664907752858268976392979073
= 10043619664651911066883029686766120169131919507076163314397915307085965058341170072938120477911396027902856306859830431800181085603701181775623189478719241
= 12964455266041997431902182249246681423017590093048617091076729201020090112909200442573801636087298080179764338147888667898243288442212586190171993932442177
 
= EllipticCurve(GF(p),[a,b])
 
= E(1128360620302355288075151618990689693489224136092325178068938705418318741031525951872324247759313197901044260703591395247778139170748768869166170361843998012748862750577419812619234165922125135009793011470953429653398381275403229335519006908182956425430354120606424111151410237675942385465833703061487938776991)
= k*P
= (k+1)*P
 
= int(Q[0])
= int(R[0])
 
assert is_prime(p)
assert is_prime(q)
 
= 0x10001
= p*q
= bytes_to_long(FLAG)
= pow(m,e,N)
 
print(f'N = {N}')
print(f'c = {c}')
 
cs

 

Solution

Denote $P = (u, v)$, $Q = (x, y)$. We can utilize the sum of points formula to get the x-coordinate of $R$ as $$ \left( \frac{y-v}{x-u} \right)^2 - (x + u)$$ Therefore, we have the relation $$ g(x, y) = x(y-v)^2 - (x+u)(x-u)^2 - N(x-u)^2 \equiv 0 \pmod{p}$$ We now have to combine this relation with $$y^2 \equiv x^3 + ax + b \pmod{p}$$ Initially I tried to calculate the resultant, but it bugged out (maybe due to incorrect parameters given) so I changed my approach.

 

If we print out the monomials of $g(x, y)$, we see that $y$ appears in two monomials, $xy^2$ and $xy$.

We can change $y^2$ as $x^3 + ax + b$ in the first one. Now we can drag the $xy$ term into the right hand side, square the equation, and change $y^2$ to $x^3 + ax + b$ to obtain a equation with only one variable $x$. For more details, please refer to the solution code.

 

This polynomial equation can now be solved, which gives us the candidates for the prime divisor of $N$. 

 

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
pytp = 13046889097521646369087469608188552207167764240347195472002158820809408567610092324592843361428437763328630003678802379234688335664907752858268976392979073
= 10043619664651911066883029686766120169131919507076163314397915307085965058341170072938120477911396027902856306859830431800181085603701181775623189478719241
= 12964455266041997431902182249246681423017590093048617091076729201020090112909200442573801636087298080179764338147888667898243288442212586190171993932442177
 
= EllipticCurve(GF(p),[a,b])
 
= E(1128360620302355288075151618990689693489224136092325178068938705418318741031525951872324247759313197901044260703591395247778139170748768869166170361843998012748862750577419812619234165922125135009793011470953429653398381275403229335519006908182956425430354120606424111151410237675942385465833703061487938776991)
= G[0]
= G[1]
 
= 22607234899418506929126001268361871457071114354768385952661316782742548112938224795906631400222949082488044126564531809419277303594848211922000498018284382244900831520857366772119155202621331079644609558409672584261968029536525583401488106146231216232578818115404806474812984250682928141729397248414221861387
= 15850849981973267982600456876579257471708532525108633915715902825196241000151529259632177065183069032967782114646012018721535909022877307131272587379284451827627191021621449090672315265556221217089055578013603281682705976215360078119427612168005716370941190233189775697324558168779779919848728188151630185987
 
P.<x, y> = PolynomialRing(GF(p))
= y ** 2 - x ** 3 - a * x - b
= x * ((y - v) * (y - v)  - (x + u) * (x-u) * (x-u)) - N * (x-u) * (x-u)
 
print(g.monomials())
# [x^4, x^3, xy^2, x^2, xy, x, 1]
 
= g.coefficients()
P.<x> = PolynomialRing(Zmod(p))
= 0
= 0
+= T[0* x^4
+= T[1* x^3 
+= T[2* x * (x^3 + a* x + b)
+= T[3* x^2 
+= T[5* x
+= T[6]
 
+= T[4* T[4* x^2 * (x^3 + a*+ b)
 
 
= f * f - g 
 
print(h.roots())
 
= 22607234899418506929126001268361871457071114354768385952661316782742548112938224795906631400222949082488044126564531809419277303594848211922000498018284382244900831520857366772119155202621331079644609558409672584261968029536525583401488106146231216232578818115404806474812984250682928141729397248414221861387
= 15850849981973267982600456876579257471708532525108633915715902825196241000151529259632177065183069032967782114646012018721535909022877307131272587379284451827627191021621449090672315265556221217089055578013603281682705976215360078119427612168005716370941190233189775697324558168779779919848728188151630185987
 
# cand = h.roots()
cand = [(52666479036523526653095613318351861523276271632713318115554199785641910004700605665354284976751168870025415689045359043450374250110154575852620226048974511), (42925282481368613878909113199174559468414118724732506754095097356205723116364073618588815566773856095001786294300257105174112147027045971030053962344407371), (112836062030235528807515161899068969348922413609232517806893870541831874103152595187232424775931319790104426070359139524777813917074876886916617036184399802)]
 
for p, t in cand:
    q = N // p
    phi = (p-1* (q-1)
    print(p * q - N)
    d = inverse(65537, phi)
    print(long_to_bytes(pow(c, d, N)))
cs

 

pure division

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
from Crypto.Util.number import *
from flag import flag
 
assert len(flag) == 40
 
= 74894047922780452080480621188147614680859459381887703650502711169525598419741
a1 = 22457563127094032648529052905270083323161530718333104214029365341184039143821
a2 = 82792468191695528560800352263039950790995753333968972067250646020461455719312
Fp = Zmod(p^3)
= EllipticCurve(Fp, [a1, a2])
 
= bytes_to_long(flag)
= E(201395103510950985196528886887600944697931024970644444173327129750000389064102542826357168547230875812115987973230106228243893553395960867041978131850021580112077013996963515239128729448812815223970675917812499157323530103467271226217465854493032911836659600850860977113580889059985393999460199722148747745817726547235063418161407320876958474804964632767671151534736727858801825385939645586103320316229199221863893919847277366752070948157424716070737997662741835)
= m * S
 
with open("output.txt""w"as f:
    f.write("T: {}".format(T))
 
 
cs

 

Solution

The required paper for this challenge is https://arxiv.org/pdf/2010.15543.pdf.

Basically, (like a few challenges this year) we need a homomorphism to a group we can efficiently calculate discrete logarithm on.

It turns out that the given elliptic curve over $\mathbb{Z}_{p^3}$ has a homomorphism to $\mathbb{Z}_{p^2}$.

UPD : It should also be noted that the binary operator in $\mathbb{Z}_{p^2}$ is sum, not product.

 

To calculate this, we need to know how to perform point multiplication in a projective curve/division free setting.

The formula for this can be found online, for example, here. I wonder if this computation can be done in sage...

 

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
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
 
= 74894047922780452080480621188147614680853399736985793708596679454987247185378
# N is the order of the same elliptic curve, except the underlying field is GF(p) instead of ring Z/p^3Z
= 74894047922780452080480621188147614680859459381887703650502711169525598419741
= 22457563127094032648529052905270083323161530718333104214029365341184039143821
= 82792468191695528560800352263039950790995753333968972067250646020461455719312
mod = p * p * p
 
def point_double(x, y, z):
    t = (3 * x * x + a * z * z) % mod
    u = (2 * y * z) % mod
    v = (2 * u * x * y) % mod
    w = (t * t - 2 * v + mod * 3) % mod
    x2 = (u * w) % mod 
    y2 = (t * (v - w) - 2 * (u * u * y * y)) % mod
    y2 = (y2 + mod) % mod
    z2 = (u * u * u) % mod
    return x2, y2, z2
 
def point_add(x0, y0, z0, x1, y1, z1):
    if x0 == 0 and y0 == 0:
        return x1, y1, z1
    if x1 == 0 and y1 == 0:
        return x0, y0, z0
    t0 = (y0 * z1) % mod
    t1 = (y1 * z0) % mod
    t = (t0 - t1 + mod) % mod
    u0 = (x0 * z1) % mod
    u1 = (x1 * z0) % mod
    u = (u0 - u1 + mod) % mod
    u2 = (u * u) % mod
    v = (z0 * z1) % mod
    w = (t * t * v - u2 * (u0 + u1)) % mod 
    w = (w + mod) % mod
    u3 = (u * u2) % mod
    x2 = (u * w) % mod
    y2 = (t * (u0 * u2 - w) - t0 * u3) % mod
    y2 = (y2 + mod) % mod
    z2 = (u3 * v) % mod
    return x2, y2, z2
 
def point_mul(x, y, z, n):
    xr, yr, zr = 001
    while n:
        if n % 2 == 1:
            xr, yr, zr = point_add(xr, yr, zr, x, y, z)
            n -= 1
        n = n // 2
        x, y, z = point_double(x, y, z)
    return xr, yr, zr
 
x1 = 201395103510950985196528886887600944697931024970644444173327129750000389064102542826357168547230875812115987973230106228243893553395960867041978131850021580112077013996963515239128729448812815223970675917812499157323530103467271226
y1 = 217465854493032911836659600850860977113580889059985393999460199722148747745817726547235063418161407320876958474804964632767671151534736727858801825385939645586103320316229199221863893919847277366752070948157424716070737997662741835
 
x2 = 49376632602749543055345783411902198690599351794957124343389298933965298693663616388441379424236401744560279599744281594405742089477317921152802669021421009909184865835968088427615238677007575776072993333868804852765473010336459028 
y2 = 342987792080103175522504176026047089398654876852013925736156942540831035248585067987750805770826115548602899760190394686399806864247192767745458016875262322097116857255158318478943025083293316585095725393024663165264177646858125759
 
print((y1 ** 2 - x1 ** 3 - a * x1 - b) % mod)
print((y2 ** 2 - x2 ** 3 - a * x2 - b) % mod)
 
xf, yf, zf = point_mul(x1, y1, 1, N)
xs, ys, zs = point_mul(x2, y2, 1, N)
 
= xs * yf
= ys * xf
= GCD(A, B)
//= g
//= g
 
gg = (A * inverse(B, p * p)) % p
 
 
for i in range(02):
    cc = bytes_to_long(b"zer0pts{")
    val = (gg - cc * (1 << 256)) % p + i * p
    flag = b"zer0pts{" + long_to_bytes(val)
    print(flag)
cs

 

Tokyo Network

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
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
import base64
import random
from Crypto.Cipher import AES
from Crypto.Random import get_random_bytes
from Crypto.Util.Padding import pad
from qulacs import QuantumState, QuantumCircuit
from qulacs.gate import *
from secret import flag
 
GATES = {
    'CNOT': (CNOT, 2),
    'H': (H, 1), 'X': (X, 1), 'Y': (Y, 1), 'Z': (Z, 1),
    'S': (S, 1), 'SDAG': (Sdag, 1), 'T': (T, 1), 'TDAG': (Tdag, 1)
}
def parse_circuit(asm, qbits=9):
    """ Convert assembly into quantum circuit
    i.e.  q0 ---+--[Z]--
                |        <= CNOT 0,1; Z 0; H 1;
          q1 --[X]-[H]--
    """
    def apply(gate, args):
        return gate(*args)
 
    circuit = QuantumCircuit(qbits)
    cnt = 0
    for instruction in asm.replace('\n''').split(';'):
        t = instruction.strip().split()
        if t == []:
            continue
 
        if len(t) < 2:
            print("[-] Invalid instruction")
            exit(0)
 
        opecode, operand = t[0].upper(), t[1:]
        if opecode not in GATES:
            print("[-] Invalid gate")
            exit(0)
 
        operand = list(map(lambda x: int(x), ''.join(t[1:]).split(',')))
        if not all(map(lambda x: 0 <= x < qbits, operand)):
            print("[-] Invalid quantum bit specified")
            exit(0)
 
        if GATES[opecode][1!= len(operand):
            print("[-] Invalid number of operands")
            exit(0)
 
        gate = apply(GATES[opecode][0], operand)
        circuit.add_gate(gate)
 
        cnt += 1
        if cnt > 100:
            print("[-] Too large circuit")
            exit(0)
 
    return circuit
 
def transfer_and_measure(q):
    # Oops, noise might happen :(
    noise = QuantumCircuit(9)
    idx = random.randrange(09)
    noise.add_gate(DephasingNoise(idx, 0.31337))
    noise.add_gate(DepolarizingNoise(idx, 0.31337))
    noise.add_gate(BitFlipNoise(idx, 0.31337))
    noise.update_quantum_state(q)
 
    # Quantum arrived! You can update (or keep) the state :)
    circuit = parse_circuit(input('[?] Your Circuit: '))
    circuit.update_quantum_state(q)
 
    # Now you measure the quantum state :P
    return random.choices(range(2**9),
                          map(lambda x: abs(x), q.get_vector()))[0& 1
 
if __name__ == '__main__':
    N = 128
    xi, xip = 0.980.98
    p = (xi * (1 + xi))**0.5 - xi
    Np = int(N * (1 + 2*xi + 2*(xi*(1+xi))**0.5 + xip))
    print("[+] Np = " + str(Np))
 
    encoder = parse_circuit("""
    CNOT 0,3; CNOT 0,6;
    H 0; H 3; H 6;
    CNOT 0,1; CNOT 0,2;
    CNOT 3,4; CNOT 3,5;
    CNOT 6,7; CNOT 6,8;
    """)
 
    measured = 0
    ra, ba = 00
    for i in range(Np):
        ra = (ra << 1| random.choice([01])
        ba = (ba << 1| random.choices([10], [p, 1-p])[0]
        q = QuantumState(9)
        q.set_zero_state()
        if ra & 1:
            X(0).update_quantum_state(q)
        if ba & 1:
            H(0).update_quantum_state(q)
        encoder.update_quantum_state(q)
 
        measured = (measured << 1| transfer_and_measure(q)
        del q
 
    print("[+] Measured state: " + bin(measured))
    bb = int(input('[?] bb = '), 2)
    print("[+] ba = " + bin(ba))
 
    Nx, Nz = 00
    for i in range(Np):
        if (ba >> i) & 1 == (bb >> i) & 1 == 1:
            Nx += 1
        elif (ba >> i) & 1 == (bb >> i) & 1 == 0:
            Nz += 1
    if Nx < N * xi or Nz - N < N * xi:
        print("[-] Something went wrong :(")
        exit(0)
 
    xa = 0
    for i in range(Np):
        if (ba >> i) & 1 == (bb >> i) & 1 == 1:
            xa = (xa << 1| ((ra >> i) & 1)
    print("[+] xa = " + bin(xa))
 
    l = []
    for i in range(Np):
        if (ba >> i) & 1 == (bb >> i) & 1 == 0:
            l.append(i)
    chosen = random.sample(l, Nz - N)
    m = 0
    for i in chosen:
        m |= 1 << i
    print("[+] m = " + bin(m))
 
    k = 0
    for i in sorted(list(set(l) - set(chosen))):
        k = (k << 1| ((ra >> i) & 1)
 
    key = int.to_bytes(k, N // 8'big')
    iv = get_random_bytes(AES.block_size)
    cipher = AES.new(key, AES.MODE_CBC, iv)
    ct = cipher.encrypt(pad(flag, AES.block_size))
    print("Result: " + base64.b64encode(iv + ct).decode())
 
cs

 

Solution

The context of this challenge is quite nice - it's quantum error correction with Shor and QKD with BB84.

I realized the BB84 part of the program after I read the flag, but if you look at it knowing it's BB84, it's quite clear.

 

Working backwards, we observe the following and eventually find the solution

  • our final goal is to find the AES key $k$
  • we know $ba$ and $m$, and we are free to choose $bb$, so we know the array $l$ and $chosen$
  • therefore, our goal is to find the bits of $ra$ in the indexes of $l \setminus chosen$
  • note that in those indexes, we have $ba$ and $bb$ bits all $0$
  • we see that this implies we don't need to worry about Hadamard gates
  • we only need to check whether or not X gate was applied at the 0th qubit
  • if there was no encoding and noise, we could just measure the 0th qubit to see whether 0 or 1 pops out
  • however, we have encoding and noise : how do we deal with this?
  • it turns out that the encoding part is simply the first part of Shor's quantum error correction code
  • therefore, we can just send the remaining part of the Shor's quantum error correction code and be done with it
  • note that our noise is single qubit, so we have no worries whatsoever
  • however, we need Toffoli gates to make the circuit : we don't have that in our arsenal of gates
  • this can be resolved by creating Toffoli gates with 2-qubit gates and single qubit gates : details are in Wikipedia.
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
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
from qulacs import QuantumCircuit, QuantumState
from qulacs.gate import *
import random
from pwn import *
from Crypto.Cipher import AES
import base64
 
 
# CNOT gate
def ADD_CNOT(a, b):
    return "CNOT " + str(a) + "," + str(b) + "; "
 
# single qubit gates
def ADD_single(s, a):
    return s + " " + str(a) + "; "
 
# toffoli
def ADD_toffoli(a, b, c): 
    ret = ""
    ret += ADD_single("H", c)
    ret += ADD_CNOT(b, c)
 
    ret += ADD_single("TDAG", c)
    ret += ADD_CNOT(a, c)
 
    ret += ADD_single("T", c)
    ret += ADD_CNOT(b, c)
 
    ret += ADD_single("TDAG", c)
    ret += ADD_CNOT(a, c)
 
    ret += ADD_single("T", b)
    ret += ADD_single("T", c)
 
    ret += ADD_CNOT(a, b)
    ret += ADD_single("H", c)
 
    ret += ADD_single("T", a)
    ret += ADD_single("TDAG", b)
 
    ret += ADD_CNOT(a, b)
 
    return ret 
 
 
 
= 128
xi, xip = 0.980.98
= (xi * (1 + xi))**0.5 - xi
Np = int(N * (1 + 2*xi + 2*(xi*(1+xi))**0.5 + xip))
 
# shor error correction
decoder = ""
decoder += ADD_CNOT(01)
decoder += ADD_CNOT(34)
decoder += ADD_CNOT(67)
decoder += ADD_CNOT(02)
decoder += ADD_CNOT(35)
decoder += ADD_CNOT(68)
decoder += ADD_toffoli(120)
decoder += ADD_toffoli(453)
decoder += ADD_toffoli(786)
decoder += ADD_single("H"0)
decoder += ADD_single("H"3)
decoder += ADD_single("H"6)
decoder += ADD_CNOT(03)
decoder += ADD_CNOT(06)
decoder += ADD_toffoli(360)
 
= remote('others.ctf.zer0pts.com''11099')
 
def get_bin():
    s = r.recvline()
    s = s.split()[-1].decode()
    return int(s, 2)
 
r.recvline()
for i in range(0860):
    r.sendline(decoder)
 
measure = get_bin()
bb = (1 << 400- 1
r.sendline(bin(bb))
 
ba = get_bin()
xa = get_bin()
= get_bin()
 
= []
for i in range(Np):
    if (ba >> i) & 1 == (bb >> i) & 1 == 0:
        l.append(i)
 
chosen = []
for i in range(Np):
    if (m & (1 << i)) > 0:
        chosen.append(i)
 
= 0
for i in sorted(list(set(l) - set(chosen))):
    k = (k << 1| ((measure >> i) & 1)
 
key = int.to_bytes(k, N // 8'big')
= r.recvline()
vv = s.split()[-1]
 
fin = base64.b64decode(vv)
iv = fin[:16]
ct = fin[16:]
 
cipher = AES.new(key, AES.MODE_CBC, iv)
print(cipher.decrypt(ct))
 
cs

'CTF' 카테고리의 다른 글

zer0lfsr 2021 Quick Solutions/Thoughts  (0) 2021.07.05
LINE CTF Crypto Writeups  (0) 2021.03.22
AeroCTF 2021 Crypto Writeups  (0) 2021.02.28
UnionCTF Crypto Writeups  (2) 2021.02.21
CryptoHack All Solve, Again  (10) 2021.01.31
  • 연구 : 공부는 확실히 꽤 했고, 글도 블로그에 많이 썼다. 랩에서 받은 태스크 하나를 잘 수행했다. 연구 자체는 아직 어렵다...
  • 대수적 정수론 : 많이 못 봤다. 3장까지라도 방학 기간에 빠르게 읽을까..
  • Kaggle : 이건 딥러닝의 기초를 들으면서 병행하는 걸로 결정했다.
  • 위상 : 4단원까지 읽었고, 실해석학과 병행하면서 읽으면 될 것 같다.
  • CTF : 2월의 마지막 두 대회를 우승하면서 기분 좋게 마무리
  • 암호학 책 리딩 : 2장, 18장, 19장을 읽고 특히 19장은 정리했다. 공부하고 싶은 건 일단 다 한듯.
  • Rust 입문 : 당연히 아직 많이 부족하지만, 코드포스 Div 2를 순조롭게 칠 정도로 편해졌다.
  • 봉사활동 : 다 채웠고 문제 없이 장학금을 받을 수 있을 것으로 기대한다.
  • 건강 : 중간에 한 번 수치가 이상하게 나와서 정신적으로 충격을 받았는데, false alarm이었던 걸로.. 대신 몸에 더 신경쓰는 삶을 살기로 했다. 이 문제로 정신 없어서 일주일 정도 아무것도 못했다. 아무튼 폭식하지 않고 술도 가급적 안 마시는 걸로. 체중도 조금 줄일 생각...

2월에는 나름 잘 쓴 것 같은 좋은 블로그 글을 많이 쓴 것 같아서 기분이 좋다 :)

건강하고 재밌는 한 학기를 보낼 수 있었으면 좋겠다. 부담없이 과목을 수강하니, 재밌는 공부를 더 하고 싶다.

I participated in AeroCTF as a member of Super Guesser. We reached 1st place :)

We were able to grab first bloods on all three cryptography challenges.

Phoenix was solved by rbtree, and I solved horcrux and boggart.

 

Phoenix

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
#!/usr/bin/env sage
 
= 2
= GF(p)
P.<x> = PolynomialRing(F)
 
 
class Cipher:
    def __init__(self, size, params):
        self.size = size
        self.params = params
 
    def sequence(self, key):
        while True:
            key = key * self.params[0]
            yield key + self.params[1]
 
    def encrypt(self, key, data, strength):
        for value, pbit in zip(self.sequence(key), data):
            xbit = sum(value[i] for i in range(0, strength, 2))
            ybit = mul(value[i] for i in range(1, strength, 2))
            
            yield int(pbit) ^^ int(xbit) ^^ int(ybit)
 
 
def main():
    size = 256
    length = 1024
    strength = 10
 
    q = P.irreducible_element(size, 'minimal_weight')
    R.<x> = P.quo(q)
 
    key, a, b = [R.random_element() for _ in range(3)]
 
    with open('flag.txt''rb'as file:
        flag = file.read().strip()
 
    message = int.from_bytes(flag, 'big')
    assert message.bit_length() < size
    plaintext = list(map(int, bin(message)[2:]))
    padding = [0* (length - len(plaintext))
 
    cipher = Cipher(size, [a, b])
    ciphertext = list(cipher.encrypt(key, padding + plaintext, strength))
    result = int(''.join(map(str, ciphertext)), 2)
 
    print(a)
    print(b)
    print(result)
 
 
if __name__ == '__main__':
    main()
 
cs

 

Solution (due to rbtree)

Let's look at the encryption scheme. 

 

We have a key $k \in GF(2^{256})$ and known parameters $a, b \in GF(2^{256})$.

To encrypt a plaintext $m$, it is first converted to binary then padded at the front with zeros to have length 1024. 

Then, we take the first 10 bits of the polynomial representation of $a^i k +b$ for each $i$, and calculate $$xbit_i = \sum \text{(even index bits)}$$ $$ybit_i = \prod \text{(odd index bits)}$$ where calculations are done in $GF(2)$. 

Finally, we encrypt our text with a stream cipher, i.e. XOR our plaintext with $\{xbit_i \oplus ybit_i\}_{i \ge 1}$. 

 

We now make some observations.

  • It's very likely that $ybit_i$ is $0$, since it is a product of 5 bits. Heuristically, it has a $31/32$ chance of being $0$. 
  • The message bit length is guaranteed to be less than 256. This means the length of the padding is at least 768.
  • If we write each 256 coefficients of $k$ as variables, then $xbit_i$ can be written as a linear combination of them (and $0, 1$)

 

Now we have our solution idea.

  • If we assume that $ybit_i = 0$ and know the ptxt/ctxt bits, we can recover $xbit_i$. This gives us a linear equation over our 256 variables.
  • Of course, we know 768 bits of the plaintext and all of the ciphertext, so we can do this for 768 different values of $i$.
  • We take 256 equations out of those 768 and hope that the system of equations have a unique solution.
  • If there is a unique solution, we can recover the key and perform the decryption process to get a candidate for our flag.
  • If our recovered flag has "Aero{" in it, we can finish the challenge. If not, we try another set of 256 equations.

 

We need to have $ybit = 0$ for all 256 bits/equations we selected. This has a probability of $$(31/32)^{256} \approx 0.0003$$ of succeeding, and we also need our system to be uniquely solvable. Thankfully, our algorithm terminates in about 2000 trials, and each trial doesn't take too long. It seems like SageMath is very efficient when it comes to matrices over $GF(2)$ :) 

 

Here's the solution script, by rbtree. It runs in a few minutes.

 

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
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
import random
from Crypto.Util.number import *
 
= 2
= GF(p)
P.<x> = PolynomialRing(F)
 
size = 256
length = 1024
strength = 10
 
= P.irreducible_element(size, 'minimal_weight')
R.<x> = P.quo(q)
 
class Cipher:
    def __init__(self, size, params):
        self.size = size
        self.params = params
 
    def sequence(self, key):
        while True:
            key = key * self.params[0]
            yield key + self.params[1]
 
    def encrypt(self, key, data, strength):
        for value, pbit in zip(self.sequence(key), data):
            xbit = sum(value[i] for i in range(0, strength, 2))
            ybit = mul(value[i] for i in range(1, strength, 2))
            
            yield int(pbit) ^^ int(xbit) ^^ int(ybit)
 
enc = 69824286833704501471834043923417254326103912707315595840737453739249974863266259092449058810542265536810346421685955365128856715192808287450464619418781355923155781710833586631897182535937891456025282049302526058466298304955387306232279075295308862156912873485647349272079984781574084434511227361370780842056
enc = list(map(int, bin(enc)[2:]))
enc = [0* (length - len(enc)) + enc
 
= x^255 + x^252 + x^246 + x^245 + x^240 + x^239 + x^236 + x^234 + x^233 + x^232 + x^231 + x^229 + x^228 + x^227 + x^224 + x^223 + x^218 + x^217 + x^216 + x^215 + x^210 + x^209 + x^204 + x^202 + x^200 + x^198 + x^197 + x^196 + x^192 + x^188 + x^187 + x^182 + x^181 + x^180 + x^178 + x^177 + x^176 + x^174 + x^173 + x^172 + x^167 + x^166 + x^161 + x^160 + x^157 + x^155 + x^154 + x^151 + x^150 + x^149 + x^148 + x^147 + x^146 + x^144 + x^140 + x^137 + x^135 + x^133 + x^132 + x^130 + x^129 + x^126 + x^122 + x^119 + x^118 + x^115 + x^112 + x^111 + x^109 + x^107 + x^106 + x^105 + x^104 + x^101 + x^100 + x^99 + x^97 + x^96 + x^94 + x^92 + x^87 + x^86 + x^84 + x^83 + x^81 + x^79 + x^75 + x^71 + x^69 + x^68 + x^67 + x^66 + x^65 + x^63 + x^62 + x^61 + x^56 + x^55 + x^53 + x^52 + x^50 + x^46 + x^44 + x^43 + x^41 + x^39 + x^38 + x^37 + x^36 + x^35 + x^34 + x^33 + x^32 + x^30 + x^29 + x^27 + x^24 + x^21 + x^17 + x^16 + x^14 + x^13 + x^12 + x^11 + x^10 + x^9 + x^5 + x^4 + x^3 + x + 1
= x^255 + x^254 + x^250 + x^247 + x^243 + x^242 + x^241 + x^238 + x^235 + x^232 + x^229 + x^227 + x^222 + x^221 + x^219 + x^218 + x^217 + x^216 + x^215 + x^211 + x^207 + x^206 + x^204 + x^202 + x^201 + x^197 + x^195 + x^193 + x^192 + x^190 + x^189 + x^188 + x^186 + x^184 + x^181 + x^180 + x^179 + x^178 + x^176 + x^173 + x^172 + x^169 + x^167 + x^165 + x^161 + x^160 + x^158 + x^149 + x^147 + x^146 + x^145 + x^140 + x^138 + x^137 + x^134 + x^133 + x^132 + x^130 + x^129 + x^128 + x^126 + x^125 + x^124 + x^121 + x^120 + x^118 + x^117 + x^114 + x^112 + x^111 + x^110 + x^109 + x^108 + x^107 + x^106 + x^105 + x^101 + x^96 + x^95 + x^94 + x^93 + x^92 + x^90 + x^89 + x^88 + x^86 + x^85 + x^84 + x^83 + x^81 + x^80 + x^79 + x^78 + x^77 + x^76 + x^71 + x^70 + x^69 + x^68 + x^67 + x^64 + x^63 + x^59 + x^56 + x^55 + x^53 + x^50 + x^46 + x^43 + x^42 + x^40 + x^38 + x^37 + x^35 + x^34 + x^33 + x^25 + x^23 + x^22 + x^21 + x^18 + x^16 + x^14 + x^13 + x^12 + x^11 + x^10 + x^8 + x^3 + x^2 + x
= sum(b[i] for i in range(0102))
 
vs = [ x^i for i in range(size) ]
vecs = []
for idx in range(length - size):
    for i in range(size):
        vs[i] = vs[i] * a
    
    arr = [int( sum(vs[j][i] for i in range(0102)) ) for j in range(size)]
    val = int(t) ^^ enc[idx]
 
    vecs.append( (arr, val) )
 
 
cipher = Cipher(size, [a, b])
cnt = 0
while True:
    cnt += 1
    random.shuffle(vecs)
 
    mat = Matrix(F, [vecs[i][0for i in range(size)])
    vec = Matrix(F, [vecs[i][1for i in range(size)]).transpose()
    try:
        vec = mat.inverse() * vec
 
        key = 0
        for i in range(size):
            if int(vec[i][0]):
                key += x^i
 
        plaintext = list(cipher.encrypt(key, enc, strength))
        plaintext = int(''.join(map(str, plaintext)), 2)
        plaintext = long_to_bytes(plaintext)
 
        if b'Aero{' in plaintext:
            while True:
                print(plaintext)
        else:
            print("FAIL2", cnt)
 
    except:
        print("FAIL", cnt)
        pass
cs

 

 

Horcrux

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
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
#!/usr/bin/env python3.8
 
from os import urandom
from gmpy2 import next_prime
from random import randrange, getrandbits
from Crypto.Cipher import AES
from fastecdsa.curve import Curve
 
 
def bytes_to_long(data):
    return int.from_bytes(data, 'big')
 
 
def generate_random_point(p):
    while True:
        a, x, y = (randrange(0, p) for _ in range(3))
        b = (pow(y, 2, p) - pow(x, 3, p) - a * x) % p
 
        if (4 * pow(a, 3, p) + 27 * pow(b, 2, p)) % p != 0:
            break
 
    return Curve(None, p, a, b, None, x, y).G
 
 
class DarkWizard:
    def __init__(self, age):
        self.power = int(next_prime(getrandbits(age)))
        self.magic = generate_random_point(self.power)
        self.soul = randrange(0self.power)
 
    def create_horcrux(self, location, weakness):
        # committing a murder
        murder = self.cast_spell(b'AVADA KEDAVRA')
 
        # splitting the soul in half
        self.soul = self.soul * pow(2-1self.power) % self.power
 
        # making a horcrux
        horcrux = (self.soul + murder) * self.magic
 
        # nobody should know location and weakness of the horcrux
        horcrux.x ^= location
        horcrux.y ^= weakness
 
        return horcrux
 
    def cast_spell(self, spell_name):
        spell = bytes_to_long(spell_name)
 
        return spell %~ spell
 
 
def encrypt(key, plaintext):
    cipher = AES.new(key=key, mode=AES.MODE_ECB)
    padding = b'\x00' * (AES.block_size - len(plaintext) % AES.block_size)
 
    return cipher.encrypt(plaintext + padding)
 
 
def main():
    wizard_age = 3000
    horcruxes_count = 2
 
    wizard = DarkWizard(wizard_age)
    print(f'Wizard\'s power:\n{hex(wizard.power)}\n')
    print(f'Wizard\'s magic:\n{wizard.magic}\n')
 
    key = urandom(AES.key_size[0])
    horcrux_length = len(key) // horcruxes_count
 
    for i in range(horcruxes_count):
        key_part = key[i * horcrux_length:(i + 1* horcrux_length]
 
        horcrux_location = bytes_to_long(key_part[:horcrux_length // 2])
        horcrux_weakness = bytes_to_long(key_part[horcrux_length // 2:])
 
        horcrux = wizard.create_horcrux(horcrux_location, horcrux_weakness)
        print(f'Horcrux #{i + 1}:\n{horcrux}\n')
 
    with open('flag.txt''rb'as file:
        flag = file.read().strip()
 
    ciphertext = encrypt(key, flag)
    print(f'Ciphertext:\n{ciphertext.hex()}')
 
 
if __name__ == '__main__':
    main()
 
cs

 

Solution

We first take a look at what's going on.

  • We have a 3000 bit prime $p$, which we know.
  • We have a point $G$ on some elliptic curve $ y^2 = x^3 + ax + b$ over $\mathbb{F}_p$ - we know $G$ but we don't know $a, b$.
  • We have a 128 bit AES key we have to recover : this is split into two "horcruxes" of 64 bits.
  • Each horcrux is split into two parts of 32 bits - they are XOR'ed to the coordinates of some point on the curve then given.

 

Writing this in mathematical equations, we have

  • We are given a 3000 bit prime $p$ and three coordinates $(x_0, y_0), (x_1, y_1), (x_2, y_2)$
  • 128 bit AES key is equal to $\alpha_1 || \beta_1 || \alpha_2 || \beta_2$ with each part being 32 bits, i.e. less than $2^{32}$.
  • $(x_0, y_0)$ lies on a elliptic curve $y^2 = x^3 + ax + b$ over $\mathbb{F}_p$, so $y_0^2 = x_0^3 + ax_0 + b$ in $\mathbb{F}_p$.
  • $(x_1 \oplus \alpha_1, y_1 \oplus \beta_1)$ lies on the curve, so $(y_1 \oplus \beta_1)^2 = (x_1 \oplus \alpha_1)^3 + a(x_1 \oplus \alpha_1) + b$
  • $(x_2 \oplus \alpha_2, y_1 \oplus \beta_2)$ lies on the curve, so $(y_2 \oplus \beta_2)^2 = (x_2 \oplus \alpha_2)^3 + a(x_2 \oplus \alpha_2) + b$

 

Since XOR is hard to deal with, we change it to normal addition - of course, it's clear that

  • $U, V \ge 0$, $V < 2^b$ implies $|(U \oplus V) - U| < 2^b$. 

Therefore, we can change each XOR to addition, by changing our constraints from $0 \le \alpha_i, \beta_i < 2^{32}$ to $|\alpha_i|, |\beta_i| < 2^{32}$.

 

Now we have to solve, given $p, x_0, y_0, x_1, y_1, x_2, y_2$, the system of equations $$\begin{equation*} \begin{aligned}  y_0^2 &= x_0^3 + ax_0 + b \\ (y_1 + \beta_1)^2 &= (x_1 + \alpha_1)^3 + a(x_1 + \alpha_1) + b \\ (y_2 + \beta_2)^2 &= (x_2 + \alpha_2)^3 + a(x_2 + \alpha_2) + b \end{aligned} \end{equation*}$$ over $\mathbb{F}_p$ and $|\alpha_i|, |\beta_i| < 2^{32}$. Now what do we do? There were a few ideas at this point...

 

Idea 1 : Lattices

Of course, "small solutions" imply that we have to use lattices. Maybe we can directly use them here?

For example, using github.com/rkm0959/Inequality_Solving_with_CVP we can solve a system of linear inequalities.

Turns out, it's quite difficult to do so, since our system has terms like $a \alpha_1$ and $a \alpha_2$. 

Because we don't know the value of $a$ and its size can be large, we can't deal with this term using this method.

 

To use the tool in the link, we need all terms of our equations to be either

  • linear in terms of our variables (i.e. variable multiplied by a known constant)
  • small (so that we can move this to the lower bound / upper bound vector for the CVP)

In this case, $a \alpha_1$ and $a \alpha_2$ term does not satisfy this requirement.

 

Idea 2 : Determinant

The failure of our Idea 1 leads to the key idea : removing $a$ in our system completely.

 

This can be done by realizing that $[1, -a, -b]^{T}$ is in the kernel of the matrix $$ \left[ \begin{matrix} y_0^2 - x_0^3 & x_0 & 1 \\ (y_1 + \beta_1)^2 - (x_1 + \alpha_1)^3 & (x_1 + \alpha_1) & 1 \\ (y_2 + \beta_2)^2 - (x_2 + \alpha_2)^3 & (x_2 + \alpha_2) & 1 \end{matrix} \right]$$ which means its determinant must be $0$ in $\mathbb{F}_p$. This gives us a polynomial equation in $\alpha_1, \alpha_2, \beta_1, \beta_2$. Now what?

 

The first idea was multivariable coppersmith - I tried to use defund's implementation, but the matrix size got way too large.

The second idea is, of course, lattices from Idea 1. If we expand our determinant, we will get a equation of the form $$\sum \text{(constant coefficient)} \cdot \text{(some monomial)} = 0$$ we can simply bound each monomials with $|\alpha_i|, |\beta_i| < 2^{32}$, then solve using the CVP toolkit above.

 

To be more precise : our equation will be in a form of $$\sum_{i=1}^k c_i m_i \equiv -c_{k+1} \pmod{p}$$ where $c_i$ is a known constant and $m_i$ is some non-constant monomial over $\alpha_i$, $\beta_i$. (Consider $m_{k+1}$ as $1$)

We solve this system by utilizing the following inequalities and plugging them into the CVP toolkit : 

$$-2^{32 \cdot \deg(m_i)} < m_i < 2^{32 \cdot \deg(m_i)}$$ $$ \sum_{i=1}^k c_i m_i+ pX = -c_{k+1}$$ where $\deg(m_i)$ is the degree of the monomial. Note that abuse of notation was used - $m_i$ is used both as a monomial, and as the actual value of the monomial with the true values of $\alpha_i$, $\beta_i$ plugged in. For further details, please take a look in the github repo and the solution.

 

UPD : Mystiz gives a good insight on entropy to explain why this line of thinking works. Link is here.

This is a very useful way to think about whether or not the above CVP toolkit will work, when it comes to the uniqueness of the solution.

 

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
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
from sage.modules.free_module_integer import IntegerLattice
 
# https://github.com/rkm0959/Inequality_Solving_with_CVP
def Babai_CVP(mat, target):
    M = IntegerLattice(mat, lll_reduce=True).reduced_basis
    G = M.gram_schmidt()[0]
    diff = target
    for i in reversed(range(G.nrows())):
        diff -=  M[i] * ((diff * G[i]) / (G[i] * G[i])).round()
    return target - diff
 
 
def solve(mat, lb, ub, weight = None):
    num_var  = mat.nrows()
    num_ineq = mat.ncols()
 
    max_element = 0 
    for i in range(num_var):
        for j in range(num_ineq):
            max_element = max(max_element, abs(mat[i, j]))
 
    if weight == None:
        weight = num_ineq * max_element
 
    if len(lb) != num_ineq:
        print("Fail: len(lb) != num_ineq")
        return
 
    if len(ub) != num_ineq:
        print("Fail: len(ub) != num_ineq")
        return
 
    for i in range(num_ineq):
        if lb[i] > ub[i]:
            print("Fail: lb[i] > ub[i] at index", i)
            return
 
    # scaling process begins
    max_diff = max([ub[i] - lb[i] for i in range(num_ineq)])
    applied_weights = []
 
    for i in range(num_ineq):
        ineq_weight = weight if lb[i] == ub[i] else max_diff // (ub[i] - lb[i])
        applied_weights.append(ineq_weight)
        for j in range(num_var):
            mat[j, i] *= ineq_weight
        lb[i] *= ineq_weight
        ub[i] *= ineq_weight
 
    # Solve CVP
    target = vector([(lb[i] + ub[i]) // 2 for i in range(num_ineq)])
    result = Babai_CVP(mat, target)
 
    for i in range(num_ineq):
        if (lb[i] <= result[i] <= ub[i]) == False:
            print("Fail : inequality does not hold after solving")
            break
 
    ## recover your result
    return result, applied_weights
 
= 0x27266a1284b761f793a529b9664693a6b1db36864a8664898d98d8010ec51afffaef2d1c79ab2078c1e0b289a1719d34b0a081ca325ba2b367017ee8e0824aaa9488409e76e923d0f0fc917ddfc1b0534e93a74246405dadfd1683f0dc31682eed0fa6b95fc235c845e16d2ef40463b7e746668dad82981fc4e05933aca410c65b36f89738f7d97502f6626c38f595338f3864638d8613fb74c16b63f3969a49ebd103ef354ed756c3cd5e67f1d2dbe5acdbc088bd6c1d503acef4ec59e4a09efac4729ca796ad25217fe74e7a0c7ef5a3e1fcd9eb9288fb89e842ef0b16642f7e84e27df4bcb623726e2c44ef46be07f9b5a5f92fe2c77d0de79fa6d46193b064207125d8935c2ff04f63e2f858e98d2518077dc58e13307f01d65ae953efd70980f3aeed320b7a6b66eb0c578dc3f05d426f412c4e3c7a9bcc68f27fe236fde41400371a39f53828824f5de3d5902cd3e7dcaee58b89c1a234188e391d412e7bc598d4f10b2bcb26aab7cd09194e80be046022ee8471
 
x0 = 0x4f69c16e493693a9342b7b9bb0ae42d0e7425996151c631a5620ad4d7ed372d285f04df82975a379080af08322fd927cf8ea9702f4533b5981351dd12358ca61c6e34af3f902b13d2981cae7cd2416e910bfe2c90f3b64c02bf933b1d11ff8ee384f5b507d9a0b1be38b15e7db3ba700ed32122a0163ec66530911fe142e098187771f5a248627a1709800069d402c1a61a17bd053ff6541f981898df6670420b0f3ea6ebd6d01d24f6fe9e3091a7c36df0c1ad7596e8ae090c54528385422911ff8e1dc201ab56844bc92acabc495c442104de5c34a5fe6661b7213500f4ccef6e76b94c34d60772fcf0c1250e66812d85efcda4b323250f740d0d57d40697fedae0bdbaceb0582cd7c82a27610e7f9fba5ac8f84e006d034dcf481f2dc9338acabd51c28afc1b4fc1d0cd7c1cd5ad21109e8708e9458e5301868c68920ab1aefb1e9184383002e6c893f1793443f8388f3f1b6e1f4ecaecb4cf000f57f677156efdacb20d60a35b8337ee416957aadffdc889dce487
y0 = 0x2393c955e5f44ab3ac052dafd83554222728393b8fd20630f3c4f122d8c86d872a7b692b3782f12a6e57352c46328f85fad2d73eaccfcbf7beb47f97da2148c4f3fc7d6751bf66569c97a64f732e7cc71767ba11b1419732035c0285ff6973fba230a3d315fba7820855208e07e6fc5bc4b2cc3868aacde3b8e9b2075bfe861b2ebccd9b6c836d85dc319290263961344d348fe8faa0aa3ebca76fe514a7b7ba313a8200727b22a8714f8bba9e5ec2c549e5bfb5857c050d1ff0471d4b01426c2ee583a7d4b6c30df5ad2d9a6902f574416f8d55ed192b29d521e5d23a54e5062400b539468dca9aa3dc558feac63c88fc696d42434ad9a83551e6860aeb6a4d84ca80713387fa8c56a1e473a82af63ec03a71202aba0ae46fb9a97132f4c92d332327e2c11b79008586b22d92d60c3155e88b5e1f9c193b363ca28f0990400afb6e8148458708b89c6023c0d5ebb746bcd754fe37a84ee4dea6baa273b2ef31a864e9586f01bae855cc0d6f6055b2546b2a918664bdb6
 
 
x1 = 0x2700ff7abe81679c770d3171c993c55da4a47cda3360e33d2b763e65517f307a39d401a256ee0634f3f1c6c244aad34422dccfe9ebd1803339972095eb2b0a57c82ee0db9ba94365cfb5270c4924c5c1ec00db2aff529aa923b113d2ca8ad6a6774fb7c655cd101ea63bc5a6ea0261f8da82d455219c7584d7de0757b2fda627c4684d3bac8f899c24178c7e0ecef6c226892b86043d3853cdb777889ce2901d8496bf0232dd000208eb2ce77e953c478551b1112ebf4b02f0086726210a50dc20ac08eb15d846084f9324b4f1f5ec73e3ef7e4a5207e04ebb866f731201e5f626084d1b61c158cfd0fdbaae8b8dd23ba599689c74a790933a3e77daa1c95fbde63b74381aff2b98f41cfebb7f1b220a4f2e8c3361734d7bd6648c720efc0a5f978917d8b84b4764e416762884f00104981e62d876d460bd8c1095cd755d8f31c5377e5ef935da77ff82e823b49d817a1f91bd4d155306173eb07efa4567d362e1cddbf0483873d5efc9286c36a58e5b3d995f3d01ca60
y1 = 0x1bf9a696ef976644dd903fb148892044fe65dd16bf12966a6b6e43be4b3b52aaa348a76ed5a53bbb400a366c59c96cf88a9c6a713273a722c0c8cd6f42b7c6f5e1911d2d323780bce65be80eef4d375dd3425a9ff832e4166765c7687aaf3c6f3c1ae7c561c46c49f0075bb68d48b95be5fd2c94996a42ba255eb638ae5ce449064cc81831f2239dbb93f4944693ee42a507dc2878988928dc9ed5bb5cb3b3eba9ac414b4171e505d285edaf02d13c0748e47b17a498cb0c15883977a11cde98995a022999a555fb3f1a3b9226dc4122f3db6c86036b1b0e6ada10b23b8c89ddb590f2186191c0f1148a04a8e35c905d4053943554c58e8f0e2ccac42c2307532e2b8f55b74fccb8ab24a977c557d10bd7c8621bef3e7f326d00c3ff28a52ee85ec61bb8aef14f67ac80da5885a93f2840a5e44d9800b0352163667ba851bb15c4083955e1fba5cd9d6e7bf103a2bbb0fbc868136cf2871815762514c691e6352af18324d1ccdc21da3e1c4c6aa771aa9fccf343ad64b8
 
x2 = 0xa439486ae1dbd96ca0796947609757b9dc068d8ed8287f98261c0c77e0940d29e25d27e16b97bc6983be505b8295886a5e880664ecb2d9161036f5beb1dbdac08da0cc2797d574e0feac67bbcb43275ab9e9d936aa17fd9a0a13a2f95e111b5e0893c4d4c3afae3bfa4a60fc5d0673215dd876b8bc960fc765401135dc708562cbeb572fcefaef1fef51782f6e0b495c0799cf89a9f47ffec11b0e780fb41e80d75681f7c9dd5fbf5e93e1dcda8091ca6a84de5b82b3abc9194913119e429781c12c11ebc6b6e80c01683344319879dec8fe59dd2f9c5b7b6e5e211153ae3b5fd161edae59777f3e76ecdeaef518495e512119451c6a97a8f471c4349ba7df8c3e9c1ba67f7a4fef57551d58b8c87607dc830c8074eec50841a1e365be90b528e4f39b8e19e7617191ac5118bb1abc44739d65384d915cb72e226122965ebae1581f55ecb12e523b0e904b3c0a1b26cf506ca030a68fe40d34c0912272277cd0d605ab057dbb5423cca28629661ad163232e766e80949
y2 = 0x29cfc760529c48d68bc086a5b4403f1d4446db04abe243c99baf659b5e67cd6cdac9a658f273f4c682b9e13dda72aaa1ede42c69c98640f2fc58eabeef143a65334e4a236a6e72a157d92ab1a541c9bcf7d953b386a40d68880312ed1e900f6ba481bf134515f5a4c245a0d924db0e5105fc44fae71fc991df85219403ba5d48f68d5d91151f8411b4cbb6971bd63030b523989fa7ec94c14136ac712d4e91eca4c930d79caab7c328043c762917c4b3f868ce037cae9f19a579272c6b13ca8c19d349f5777bf6ac9c078d8472c92582daf96b30d4f7b8fdc004a36b792c133e2b6511956480892636ea91a3361afd8a3afbe3a5bf889feb5a5dc143f5a917347591634218066f2f71b36afd5257f6637152ac9a0965daa881ddc86ca8a8545b255cbb14e7738297a55c7428d9a0e79dee37f801fc1d49d205be809e0cd1b8a1b9d1d5b1b1a9ae98e7c564718cb17d10a3dd01c2914d7bb96c45a4fc942c2ed628354882464407c3208938282471aa2b8016e46c998e4
 
 
= Integers(p)
 
# determinant calculation
# a, b = alpha_1, beta_1
# c, d = alpha_2, beta_2
P.<a, b, c, d> = PolynomialRing(R)
= 0
+= (y0 ** 2 - x0 ** 3* (x1 + a) * 1
+= (x0) * 1 * ((y2 + d) ** 2 - (x2 + c) ** 3)
+= 1 * ((y1 + b) ** 2 - (x1 + a) ** 3* (x2 + c)
-= (y0 ** 2 - x0 ** 3* 1 * (x2 + c)
-= x0 * ((y1 + b) ** 2 - (x1 + a) ** 3* (1)
-= 1 * (x1 + a) * ((y2 + d) ** 2 - (x2 + c) ** 3)
 
# coefficient
= f.coefficients()
 
# degree
print(f.monomials())
# [a^3*c, a*c^3, a^3, a^2*c, b^2*c, a*c^2, c^3, a*d^2, a^2, b^2, a*c, b*c, c^2, a*d, d^2, a, b, c, d, 1]
= [poly.degree() for poly in f.monomials()]
print(D)
# [4, 4, 3, 3, 3, 3, 3, 3, 2, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 0]
 
 
# CVP setup
print(len(C)) # 20
 
= Matrix(ZZ, 2020)
lb = [0* 20
ub = [0* 20
 
# monomial bounds
for i in range(019):
    M[i, i] = 1
    lb[i] = -(2 ** (32 * D[i]))
    ub[i] = (2 ** (32 * D[i]))
 
# the polynomial
for i in range(019):
    M[i, 19= C[i]
M[1919= p
lb[19= int((-C[19]) % p)
ub[19= int((-C[19]) % p)
 
# solve CVP
result, applied_weights = solve(M, lb, ub)
 
# recover variables
alpha_1 = result[15// applied_weights[15]
beta_1 = result[16// applied_weights[16]
alpha_2 = result[17// applied_weights[17]
beta_2 = result[18// applied_weights[18]
 
print(alpha_1, beta_1, alpha_2, beta_2)
 
# recover actual key
alpha_1 = x1 ^^ (x1 + alpha_1)
beta_1 = y1 ^^ (y1 + beta_1)
alpha_2 = x2 ^^ (x2 + alpha_2)
beta_2 = y2 ^^ (y2 + beta_2)
 
print(alpha_1, beta_1, alpha_2, beta_2)
# the rest is trivial decryption
cs

 

 

Boggart

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
63
64
65
66
67
68
69
70
71
72
73
#!/usr/bin/env python3.8
 
from gmpy import next_prime
from random import getrandbits
 
 
def bytes_to_long(data):
    return int.from_bytes(data, 'big')
 
 
class Wardrobe:
    @staticmethod
    def create_boggarts(fear, danger):
        # for each level of danger we're increasing fear
        while danger > 0:
            fear = next_prime(fear)
 
            if getrandbits(1):
                yield fear
                danger -= 1
 
 
class Wizard:
    def __init__(self, magic, year, experience):
        self.magic = magic
        self.knowledge = year - 1 # the wizard is currently studying the current year
        self.experience = experience
 
    def cast_riddikulus(self, boggart):
        # increasing the wizard's experience by casting the riddikulus charm
        knowledge, experience = self.knowledge, self.experience
 
        while boggart > 1:
            knowledge, experience = experience, (experience * self.experience - knowledge) % self.magic
            boggart -= 1
 
        self.experience = experience
 
 
def main():
    year = 3
    bits = 512
    boggart_fear = 31337
    boggart_danger = 16
 
    neutral_magic, light_magic, dark_magic = [getrandbits(bits) for _ in range(3)]
    magic = next_prime(neutral_magic | light_magic) * next_prime(neutral_magic | dark_magic)
 
    print('Hello. I am Professor Remus Lupin. Today I\'m going to show you how to deal with the boggart.')
 
    print(neutral_magic)
    print(magic)
 
    with open('flag.txt''rb'as file:
        flag = file.read().strip()
 
    # some young wizards without knowledge of the riddikulus charm
    harry_potter = Wizard(magic, year, bytes_to_long(b'the boy who lived'))
    you = Wizard(magic, year, bytes_to_long(flag))
 
    for boggart in Wardrobe.create_boggarts(boggart_fear, boggart_danger):
        # wizards should train to learn the riddikulus charm
        harry_potter.cast_riddikulus(boggart)
        you.cast_riddikulus(boggart)
 
    # wizard's experience should be increased
    print(harry_potter.experience)
    print(you.experience)
 
 
if __name__ == '__main__':
    main()
 
cs

 

Solution

There is a lot going on in this code, so let's take a look at it.

  • We have random 512 bit integer $M_N, M_L, M_D$ corresponding to neutral/light/dark magic
  • The value magic is created as $M = \text{next_prime}(M_N | M_L) \times \text{next_prime}(M_N | M_D)$.
  • We are given the values of $M_N$ and $M$ - so the factorization of $M$ is unknown.

Now, the text "the boy who lived" and the flag goes through some encryption process, and we know their ciphertexts.

  • The encryption process itself is quite random - it uses the cast_ridikkulus algorithm randomly.
  • To be more precise, it starts the parameter "fear" at 31337, and increases it to the next prime each time.
  • Then, with a 50% chance, the text gets encrypted with the said value of "fear" and the cast_ridikkulus algorithm.
  • If we have used the cast_ridikkulus algorithm 16 times, the encryption process terminates.

 

cast_riddikulus algorithm

 

Naturally, we have to look at the cast_riddikulus algorithm and how to reverse/decrypt it.

We note that the value of self.knowledge is fixed at $2$ at all times. Denote the self.experience as $x$, and the boggart as $n$.

 

The algorithm calculates a linear recurrence - it starts with $a_0=2$, $a_1 = x$, and moves up by $$a_{i} = xa_{i-1} - a_{i-2}$$ Let's forget about the modulos and solve this recurrence. Using standard techniques, we get $$a_n = \left( \frac{x - \sqrt{x^2 - 4}}{2} \right)^n + \left( \frac{x + \sqrt{x^2- 4}}{2} \right)^n $$ This looks kinda scary, but then we can see that this is just $$x = t + t^{-1} \implies a_n = t^n + t^{-n}$$ Note that this implication holds even when we consider the fact that we are working modulo $M$.

 

This has a lot of meaning - for example, the encryption function is commutative in a sense that encrypting a text (experience) with (boggart) $m$ then $n$ is equivalent to encrypting with $n$ then $m$, which is equivalent to encryption once with $mn=nm$. 

 

It also implies that under some "mild" conditions, we can reverse the encryption function as follows - 

  • Given $ctxt \equiv t^n +  t^{-n} \pmod{M}$, we solve for the value $t^n \pmod{M}$.
  • We solve for $t \pmod{M}$ using the value of $t^n \pmod{M}$.
  • We get our original value $ptxt = t + t^{-1} \pmod{M}$.

In the above algorithm, we needed

  • Factorization of $M$, and we need $t$ to be coprime to $M$
  • $t^n + t^{-n} \equiv ctxt \pmod{M}$ to be solvable, i.e. $ctxt^2 - 4$ being a quadratic residue
  • $n$ being coprime to $\phi(M)$, for the discrete root to be easy to solve 
  • (of course, we can use adlemann root extraction, but that's not very fast)

Let's review our requirements and how difficult they are

  • Factorization of $M$ : this is hard but must be done - and this is probably why we are given $M_N$.
  • $t$ being coprime to $M$ : this will probably hold even if we don't care about it
  • $ctxt^2 - 4$ being a quadratic residue : if not, we can just get our random values of $M_N, M_L, M_D, M$ again.
  • $n$ being coprime to $\phi(M)$ : in our context, we will use 16 primes larger than $30000$ as $n$, so it'll most likely be true.

Therefore, it's time to go back and try to factorize $M$ with our information of $M_N$.

 

factorization of $M$

 

The key idea is that the value of $M_N$ gives us known bits of each prime $p, q$ of $M$.

  • $M_N|M_L$ and $M_N|M_D$ obviously has all bits that $M_N$ has. 
  • Then, we apply the next_prime function. Due to the small prime gap, this is similar to adding a few thousand to each number.
  • Therefore, some lower bits will definitely change, and even higher bits can change due to carry

Consider $M'_N = M_N - (M_N \pmod{2^{16}})$, i.e. $M_N$ but changed the $16$ LSB to 0.

  • If no carry happens at the 16th bit, our primes will still have all the bits of $M'_N$.
  • If a carry happens at the 16th bit, we hope our primes will still have all the bits of $M'_N + 2^{16}$.

Therefore, we write $M''_N = M'_N \& (M'_N + 2^{16})$ and claim that $p, q$ all have the bits of $M''_N$.

 

It is known that with around 57% of the bits of $p, q$ known, we can factorize $N$ using the following algorithm -

  • At stage $k$, calculate all possible $(p \pmod{2^k}, q \pmod{2^k})$ such that $pq \equiv N \pmod{2^k}$.
  • Of course, we must respect the known bits of $p, q$ while calculating this.
  • We can change $p \pmod{2^k}$ to $p \pmod{2^{k+1}}$ or $p + 2^k \pmod{2^{k+1}}$ while moving to the next stage.
  • If we find a set of $(p, q)$ such that $pq = N$, the algorithm terminates. 
  • This can be written in a depth-first search style or a breadth-first search style. The latter works better on my machine.

We require that $M''_N$ has a lot of ones in its binary representation. We simply restart the connection until we get good parameters.

For the details of the factorization algorithm, please refer to the solution script below.

 

how did we encrypt this?

 

Now we know the factorization of $M$. We can now check if $ctxt^2 - 4$ is a quadratic residue modulo $M$.

We test this for the ciphertext of the "the boy who lived" and the ciphertext of the flag. If this check fails, we run our algorithm again with a new connection for a good set of parameters. We now assume $ctxt^2 - 4$ is a quadratic residue.

 

We reasonably guess that our encryption algorithm will have used 16 primes out of the 35 consecutive primes larger than 31337.

We also verify that $p-1$ and $q-1$ are all coprime to these 35 candidate primes. If not, we look for a good set of parameters.

 

Now we can encrypt and decrypt as we like. Let's find out which 16 primes were used for encryption.

 

To do this, we use meet-in-the-middle algorithm. Basically what we do is

  • Split the set of 35 primes into two pieces, and assume that 8 primes were used in each piece
  • Calculate the encryption of "the boy who lived" with each combination of 8 primes of the first piece
  • Calculate the decryption of the ciphertext with each combination of 8 primes of the second piece
  • Find a match, and if there is one, we have found our set of 16 primes.
  • If not, we can try a new split of the primes and continue on.

Now we know which primes were used for encryption. All we need to do is decrypt the ciphertext of the flag.

 

additional notes

 

Note 1 : speeding up cast_riddikulus algorithm

  • This is a linear recurrence, so we can speedup using matrix multiplication
  • However, this is also a Lucas sequence used in Williams's $p+1$ algorithm
  • Therefore, we can just copy the implementation from Wikipedia to achieve a speedup.

Note 2 : requirements relaxation

  • If we assume that the flag is small, we only need to work with mod $p$ or mod $q$, instead of mod $M$.
  • Therefore, we can relax our quadratic residue conditions to "QR modulo at least one of the primes"
  • Thankfully, I was very lucky with the parameters and got QR modulo $M$.

Note 3 : meet-in-the-middle, but not quite middle

  • If we try to run the code, we can see that encryption is faster than decryption.
  • Therefore, it's reasonable to give the encryption part more work than just "half" of the bunch.

The solution code is displayed in parts for reader's understanding. Feel free to comment and ask questions...

 

reading parameters

1
2
3
4
5
6
7
8
9
10
11
12
13
def recvint():
    T = r.recvline()
    return int(T[:-1].decode(), 10)
 
 
def get_vals():
    r.recvline()
    C = recvint()
    N = recvint()
    H = bytes_to_long(b'the boy who lived')
    H_AFTER = recvint()
    F_AFTER = recvint()
    return C, N, H_AFTER, F_AFTER
cs

 

factorization of $M$

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
def bsum(x):
    ret = 0
    while x > 0:
        ret += x % 2
        x //= 2
    return ret 
 
def factor(N, C, cand, k):
    if len(cand) == 0:
        return
    print(k, len(cand))
    # fix pq mod 2^{k+1} by adding 2^k or not
    ret = []
    if (C >> k) % 2 == 1# we know p, q has 2^k
        for p, q in cand:
            if p * q == N:
                print(p)
                print(q)
                return
            pp = p + (1 << k)
            qq = q + (1 << k)
            if (pp * qq - N) % (1 << (k + 1)) == 0 and pp * qq <= N:
                ret.append((pp, qq))
        factor(N, C, ret, k + 1)
    else:
        for p, q in cand:
            if p * q == N:
                print(p)
                print(q)
                return
            for i in range(02):
                for j in range(02):
                    pp = p + i * (1 << k)
                    qq = q + j * (1 << k)
                    if (pp * qq - N) % (1 << (k + 1)) == 0 and pp * qq <= N:
                        ret.append((pp, qq))
        factor(N, C, ret, k + 1)
 
# M_N to M''_N
-=  C % (1 << 16)
&= (C + (1 << 16))
# check number of 1s, if bsum(C) is small (<256), stop and rerun
print(C.bit_length(), bsum(C)) 
# factorize
factor(N, C, [(11)], 1)
cs

 

encryption and decryption

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
def mlucas(v, a, n): # encryption fast, from Wikipedia
    v1, v2 = v, (v**2 - 2) % n
    for bit in bin(a)[3:]:
        if bit == "0":
            v1, v2 = ((v1**2 - 2) % n, (v1*v2 - v) % n)
        else:
            v1, v2 = ((v1*v2 - v) % n, (v2**2 - 2) % n)
    return v1
 
def REV(val, ex, pr): # decryption modulo a prime
    # t^n + 1/t^n = val
    if val == 0# if invalid
        return 0 # return invalid
    cc = tonelli((val * val - 4) % pr, pr) # modular sqrt
    if cc == -1# failure
        return 0 # return failure
    tn = ((val + cc) * inverse(2, pr)) % pr
    t = pow(tn, inverse(ex, pr-1), pr)
    return ((t + inverse(t, pr))) % pr
cs

 

meet-in-the-middle

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
primes = [31337313573137931387313913139331397314693147731481314893151131513315173153131541315433154731567315733158331601316073162731643316493165731663316673168731699317213172331727317293174131751317693177131793]
 
= primes[15:35]
= primes[:15]
 
LV = []
RV = []
 
# generally, encryption is faster than decryption
# therefore, we leave the "heavy" work for the encryption part...
 
for combo_L in tqdm(itertools.combinations(L, 8)):
    CC = list(combo_L)
    cur = H
    # you can also just multiply all primes and encrypt at once
    for pr in combo_L:
        cur = mlucas(cur, pr, p)
    CC.append(cur)
    LV.append(CC)
 
print("LEFT DONE")
 
for combo_R in tqdm(itertools.combinations(R, 8)):
    CC = list(combo_R)
    cur = HAFTER
    # you can also just multiply all primes and decrypt at once
    for pr in combo_R:
        cur = REV(cur, pr, p)
    CC.append(cur)
    RV.append(CC)
 
print("RIGHT DONE")
 
for i in tqdm(range(len(LV))):
    for j in range(len(RV)):
        if LV[i][8== RV[j][8]:
            print("FOUND!!")
            print(LV[i])
            print(RV[j])
cs

 

the finish, and the parameters I got

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
= bytes_to_long(b'the boy who lived')
sys.setrecursionlimit(10 ** 6)
 
= 12124756848659098434025945489515506912896022954145117746118560512007665385702760439414990812257455576297349156226093149988609289245714223348281989890389750
= 157597985389833012337654133040126048344064622845161536236706459270343905778002470548499258715513540516431526518690532156245280894778788692043941237295606686168037171464988128988463706375526180496632421973522548093894845498612792150825707672843107252573999144787226703076358545319417530365329436368718460943493
HAFTER = 66161881147822169408519540711422900962287264738494143175834051626001922954586728648835878096124744364195826536091510407493007528877139856387261499433277826944946254511824024047480941829026088269865298686453128715170657018128276813244425143986311708022950785583195028647859774987948632731985531259912781472862
FAFTER = 149186530719822614329126547638374064715014252925601014676661223009475822460330945440469384214084001910035138025738722725987466200681944900264994344927428683388976167111544750466576538413516454786176229441173029050647235653998791477157269246962955063391947778663841551982999293815571149539542758304215156142104
 
-=  C % (1 << 16)
&= (C + (1 << 16))
 
# factor(N, C, [(1, 1)], 1)
 
= 12558711464274427739528720572494472142909592647353129013838950445222814801805965383430302364628487022743397586481672449715551542652546057434522020868473011
= 12548897698474048380978452887676419841595083766206501465313606366388795637681128899285066184154566275021196792336453455837893284460576050862858626214885863
 
# recovered set of primes
pr = [31543315673157331607316493166731687316993135731379313973146931477315133151731531]
hp = FAFTER % p 
hq = FAFTER % q 
 
# reverse
for i in pr:
    hp = REV(hp, i, p)
    hq = REV(hq, i, q)
 
flag, tt = CRT(hp, p, hq, q) # CRT
print(long_to_bytes(flag))
cs

'CTF' 카테고리의 다른 글

LINE CTF Crypto Writeups  (0) 2021.03.22
zer0pts CTF 2021 Crypto Writeups  (1) 2021.03.07
UnionCTF Crypto Writeups  (2) 2021.02.21
CryptoHack All Solve, Again  (10) 2021.01.31
0x41 CTF SPN Writeup  (0) 2021.01.31

이어서 간다. 이번 글에서 다룰 주제는

  • OSPEP를 이용한 DRS의 tight contraction factor
  • OSPEP를 이용한 optimal parameter selection
  • CVXPY를 이용한 위 주제들과 관련된 구현 몇 가지

 

Part 3. Tight Analytic Contraction Factors for DRS

결과부터 소개한다. 우리가 앞에서 준비했던, 논문의 정리 4.1.만 소개한다. 그래도 논문 링크를 열고 정리 4.1.~4.4.를 한 번 확인하는 것을 추천.

결과도 매우 중요하지만, 아무래도 일단 관심이 가는 것은 대체 이렇게 복잡한 것을 어떻게 증명했냐는 것이다. 논문의 말을 옮기면,

  • 기본적으로 Mathematica 같은 Computer Algebra System은 필요하다
  • OSPEP에서 얻은 Primal/Dual Problem을 활용한다. 이들을 풀기 위해서, 다음과 같은 접근을 사용한다.
  • Primal : 핵심 아이디어는 worst-case operator가 $\mathbb{R}^2$ 위에서 존재할 것이라는 추측이다. 이는 Primal Problem의 해 $G^{*}$의 rank가 2 이하라는 것과 동치인데, 이 추측의 근거는 complementary slackness라고 한다. 이러한 2차원 위에서의 worst-case analysis는 SDP 대신 non-convex quadratic programming 문제로 분석할 수 있고, 이 문제에 대응되는 Karush–Kuhn–Tucker condition을 이용하여 문제를 해결했다고 한다.
  • Dual : 역시 여기서도 Dual Problem의 해로 얻어지는 행렬 $S^{*}$가 rank 2 이하라는 가정을 한다. 근거는 Primal과 동일하다. 이 가정으로 $\rho^2$과 $\lambda_\mu^A$, $\lambda_\beta^B$ 사이의 관계식을 추가적으로 얻고, 이를 constraint로 하여 $\rho^2$을 최소화했다고 한다. 말이 쉽지만 정말 어려워 보인다..
  • 이렇게 Primal/Dual 문제를 (analytical) 푼 다음, 그 결과가 같음을 확인하면 자동으로 증명이 끝난다.

Appendix를 구경하면 이 결과의 엄밀한 증명을 직접 볼 수 있다. ㄹㅇ 가슴이 웅장해진다... 

 

Part 4 : Optimal Parameter Selection with OSPEP

여기서는 operator class $\mathcal{Q}_1$, $\mathcal{Q}_2$를 알고 있을 때 parameter $\alpha, \theta$를 optimal 하게 설정하는 방법을 다룬다.

이는 결국 $\rho^2_{*} (\alpha, \theta)$를 해당 $\alpha, \theta$에 대한 optimal contraction factor라고 할 때, $$\rho^2_{*} = \inf_{\alpha, \theta > 0} \rho^2_{*}(\alpha, \theta)$$와 대응되는 $\alpha_{*}$, $\theta_{*}$를 찾는 문제와 같다. 원래 $\theta \in (0, 2)$라는 조건이 있으나, 이 조건은 relax 해도 무방하다. 

 

$\rho^2_{*}(\alpha, \theta)$의 값 자체는 SDP를 풀어서 계산할 수 있으니, 이차원 위에서 미분 불가능한 함수를 최소화하는 일반적인 알고리즘이 있다면 이를 사용하여 문제를 해결할 수도 있을 것이다. 그러나, 조금 더 작업을 거치면 더욱 편하게 문제를 풀 수 있다. $\theta$를 최적화 문제 안에 집어넣자.

 

Dual Problem $$\begin{equation*} \begin{aligned} \text{minimize} \quad & \rho^2 \\ \text{subject to} \quad & \lambda_\mu^A, \lambda_\beta^B \ge 0 \\ \quad & S = -M_O - \lambda_\mu^A M_\mu^A - \lambda_\beta^B M_{\beta}^B + \rho^2 M_I \succeq 0 \end{aligned} \end{equation*}$$을 다시 생각해보자. 우리가 $\alpha, \theta$를 고정된 값이 아니라 최적화 문제의 변수로 두지 못하는 이유는, 행렬 $S$의 entry가 $\alpha, \theta$의 nonlinear term을 포함하고 있기 때문이다. $\alpha$는 아예 분모에 있고, $\theta$는 $\theta^2$ 항이 $S$에 있다. 여기서 $\theta$에 대한 quadratic term을 제거하기 위해, Schur Complement를 사용한다.

 

$$\begin{equation*} \begin{aligned} S &= -M_O - \lambda_\mu^A M_\mu^A - \lambda_\beta^B M_{\beta}^B + \rho^2 M_I \\ S' &= - \lambda_\mu^A M_\mu^A - \lambda_\beta^B M_{\beta}^B + \rho^2 M_I \end{aligned} \end{equation*} $$라고 정의하면, $$M_O = (1, \theta, -\theta)^T (1, \theta, -\theta)$$이므로, Schur Complement의 성질에서 $S \succeq 0$는 $v = (1, \theta, -\theta)^T$라 했을 때 $$\overline{S} = \left[ \begin{matrix} S' & v \\ v^T & 1 \end{matrix} \right] \succeq 0$$과 동치이다. 또한, $S'$은 각 entry에 $\theta$에 대한 non-linear term이 없으므로, $\theta$까지 SDP에 넣을 수 있다.

 

즉, 새로운 문제를 $$\begin{equation*} \begin{aligned} \text{minimize} \quad & \rho^2 \\ \text{subject to} \quad & \lambda_\mu^A, \lambda_\beta^B, \theta \ge 0 \\ \quad & S' = - \lambda_\mu^A M_\mu^A - \lambda_\beta^B M_{\beta}^B + \rho^2 M_I \\ \quad & v = (1, \theta, -\theta)^T \\ \quad & \overline{S} = \left[ \begin{matrix} S' & v \\ v^T & 1 \end{matrix} \right] \succeq 0 \end{aligned} \end{equation*}$$로 바꾸면, $\theta$에 대한 최적화까지 자동으로 된다. 이제 남은 것은 $\alpha$에 대한 최적화다.

 

이렇게 해서 얻은 주어진 $\alpha$에 대한 optimal contraction factor를 $\rho^2_{*}(\alpha)$라고 하자. 이를 최소화하면 된다.

변수 하나에 대한 (미분 불가능한) 함수의 minimization은 다양한 방법이 있다. 그런데, 실험 결과에 따르면

  •  함수 $\rho^2_{*}(\alpha)$는 연속함수이다. (별로 놀랍지는 않은 사실)
  • 함수 $\rho^2_{*}(\alpha)$는 unimodal 함수이다. (놀라운 사실)

그러므로, 이를 최소화하기 위해서 (PS를 했던 사람의 입장에서 쉽게 생각하면) 삼분탐색을 사용해도 좋다.

사실 scipy나 matlab 등에서 이미 derivative-free optimization이 구현이 되어있기 때문에, 이걸 가져다 쓰면 된다.

 

Part 5 : Programming OSPEP with CVXPY

이건 논문에는 없고 (정확하게는 CVXPY로 하지 않고 MATLAB 등으로 한 것 같다) 그냥 내가 심심해서 하는 것이다.

논문의 저자들이 사용한 코드를 보고 싶다면, github.com/AdrienTaylor/OperatorSplittingPerformanceEstimation로 가자. 

내가 여기서 짠 코드들은 github.com/rkm0959/rkm0959_implements에 가면 있다. 

 

Task 1 : Verifying Theorem 4.1. Numerically

CVXPY를 이용하여 Theorem 4.1.을 검증해보자. 방법은 다음과 같다.

  • $\alpha, \theta, \mu, \beta$를 적당히 랜덤하게 선정한다.  
  • 대응되는 Primal/Dual Problem을 CVXPY로 해결한다. SDP를 풀어야하니, blas+lapack을 설치해야 한다.
  • 또한, Theorem 4.1.을 이용하여 optimal contraction factor를 계산한다. 
  • Theorem 4.1.에는 Case가 총 5가지가 있다. 이때, 어떤 Case에서 계산이 되었는지도 기록한다.
  • Primal Problem, Dual Problem, Theorem 4.1.로 얻은 답이 일치하는지 확인한다.
  • 위 과정을 여러 번 반복한 후, 5가지의 Case를 각각 몇 번 다루었는지를 확인한다.

이를 구현한 코드는 다음과 같다. 

 

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
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
import cvxpy as cp 
import matplotlib.pyplot as plt 
import numpy as np
import math
 
def THM_4_1_scaled(theta, mu, beta):
    assert 0 < theta < 2 and mu > 0 and beta > 0
 
    if mu * beta - mu + beta < 0 and theta <= 2 * (beta + 1* (mu - beta - mu * beta) / (mu + mu * beta - beta - beta * beta - 2 * mu * beta * beta):
        return 1, abs(1 - theta * beta / (beta + 1))
    if mu * beta - mu - beta > 0 and theta <= 2 * (mu * mu + beta * beta + mu * beta + mu + beta - mu * mu * beta * beta) / (mu * mu + beta * beta + mu * mu * beta + mu * beta * beta + mu + beta - 2 * mu * mu * beta * beta):
        return 2, abs(1 - theta * (1 + mu * beta) / (1 + mu) / (1 + beta))
    if theta >= 2 * (mu * beta + mu + beta) / (2 * mu * beta + mu + beta):
        return 3, abs(1 - theta)
    if mu * beta + mu - beta < 0 and theta <= 2 * (mu + 1* (beta - mu - mu * beta) / (beta + mu * beta - mu - mu * mu - 2 * mu * mu * beta):
        return 4, abs(1 - theta * mu / (mu + 1))
    
    ret = (2 - theta) / 4
    ret = ret * ((2 - theta) * mu * (beta + 1+ theta * beta * (1 - mu))
    ret = ret * ((2 - theta) * beta * (mu + 1+ theta * mu * (1 - beta))
    ret = ret / mu / beta
    ret = ret / (2 * mu * beta * (1 - theta) + (2 - theta) * (mu + beta + 1))
    return 5, math.sqrt(ret)
 
def THM_4_1_unscaled(alpha, theta, mu, beta):
    return THM_4_1_scaled(theta, mu * alpha, beta / alpha)
 
def gen_M_O(theta):
    ret = [[1, theta, -theta], 
           [theta, theta * theta, -theta * theta], 
           [-theta, -theta * theta, theta * theta]]
    return np.array(ret)
 
def gen_M_I():
    ret = np.zeros((33))
    ret[00= 1
    return ret
 
def gen_M_mu_A(alpha, mu):
    ret = [[0-1/20], 
           [-1/2-(1+alpha*mu), 1],
           [010]]
    return np.array(ret)
 
def gen_M_beta_B(alpha, beta):
    ret = [[-beta/alpha, 0, beta/alpha + 1/2], 
           [000],
           [beta/alpha + 1/20-beta/alpha-1]]
    return np.array(ret)
 
 
def Primal(alpha, theta, mu, beta):
    G = cp.Variable((33), symmetric = True)
    
    M_O = gen_M_O(theta)
    M_I = gen_M_I()
    M_mu_A = gen_M_mu_A(alpha, mu)
    M_beta_B = gen_M_beta_B(alpha, beta)
 
    constraints = []
    constraints.append(G >> 0)
    constraints.append(cp.trace(M_I @ G) == 1)
    constraints.append(cp.trace(M_mu_A @ G) >= 0)
    constraints.append(cp.trace(M_beta_B @ G) >= 0)
 
    objective = cp.Maximize(cp.trace(M_O @ G))
 
    problem = cp.Problem(objective, constraints)
    problem.solve()
 
    return math.sqrt(problem.value)
 
def Dual(alpha, theta, mu, beta):
    lambda_mu_A = cp.Variable()
    lambda_beta_B = cp.Variable()
    rho_sq = cp.Variable()
 
    M_O = gen_M_O(theta)
    M_I = gen_M_I()
    M_mu_A = gen_M_mu_A(alpha, mu)
    M_beta_B = gen_M_beta_B(alpha, beta)
 
    constraints = []
    constraints.append(lambda_mu_A >= 0)
    constraints.append(lambda_beta_B >= 0)
    constraints.append(-M_O - lambda_mu_A * M_mu_A - lambda_beta_B * M_beta_B + rho_sq * M_I >> 0)
 
    objective = cp.Minimize(rho_sq)
 
    problem = cp.Problem(objective, constraints)
    problem.solve()
 
    return math.sqrt(problem.value)
 
eps = 0.0005
checked = [00000]
 
for _ in range(300):
    alpha = np.random.uniform(low = 0.5, high = 2.0)
    theta = np.random.uniform(low = 0.05, high = 1.95)
    mu = np.random.uniform(low = 0.1, high = 3.9)
    beta = np.random.uniform(low = 0.1, high = 3.9)
 
    assert alpha > 0 and theta > 0
    assert mu > 0 and beta > 0 
 
    val_1 = Primal(alpha, theta, mu, beta)
    val_2 = Dual(alpha, theta, mu, beta)
    whi, val_3 = THM_4_1_unscaled(alpha, theta, mu, beta)
 
    checked[whi - 1+= 1
 
    if abs(val_2 - val_1) > eps or abs(val_3 - val_1) > eps or abs(val_3 - val_2) > eps:
        print(val_1, val_2, val_3)
        print("Failed at Case", whi)
        print("Parameters:", alpha, theta, mu, beta)
        break
 
print("Checked Cases")
for i in range(5):
    print("Case #{} Verified {} Times!".format(i + 1, checked[i]))
 
print("Finished!")
cs

 

이렇게 하여 Theorem 4.1.의 결과를 검수할 수 있다. SDP solver가 문제를 잘 해결할 수 있도록 극단적인 값들은 넣지 않았다.

각 코드의 의미는 거의 자명하니 (식을 그대로 옮긴 것에 불과하다) 자세한 설명은 생략하도록 하겠다. Task 1은 이렇게 끝.

 

Task 2 : Optimal Parameter Selection with OSPEP

이번에는 Part 4에서 다루었던 내용을 구현해보자. 우리가 할 일은

  • 주어진 $\mu, \beta$에 대하여 optimal parameter $\alpha, \theta$를 계산
  • $\rho^2_{*} (\alpha)$의 그래프를 그리고, 실제로 unimodal function인지 눈으로 확인
  • 물론, 이러한 내용 및 그래프는 논문에도 있으며, DRS가 아니라 DYS로 설명되어 있다.

Part 4의 내용을 그대로 CVXPY로 옮기면, 다음과 같은 코드를 얻을 수 있다.

 

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
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
import cvxpy as cp 
import matplotlib.pyplot as plt 
import numpy as np
import math
 
def THM_4_1_scaled(theta, mu, beta):
    assert 0 < theta < 2 and mu > 0 and beta > 0
 
    if mu * beta - mu + beta < 0 and theta <= 2 * (beta + 1* (mu - beta - mu * beta) / (mu + mu * beta - beta - beta * beta - 2 * mu * beta * beta):
        return 1, abs(1 - theta * beta / (beta + 1))
    if mu * beta - mu - beta > 0 and theta <= 2 * (mu * mu + beta * beta + mu * beta + mu + beta - mu * mu * beta * beta) / (mu * mu + beta * beta + mu * mu * beta + mu * beta * beta + mu + beta - 2 * mu * mu * beta * beta):
        return 2, abs(1 - theta * (1 + mu * beta) / (1 + mu) / (1 + beta))
    if theta >= 2 * (mu * beta + mu + beta) / (2 * mu * beta + mu + beta):
        return 3, abs(1 - theta)
    if mu * beta + mu - beta < 0 and theta <= 2 * (mu + 1* (beta - mu - mu * beta) / (beta + mu * beta - mu - mu * mu - 2 * mu * mu * beta):
        return 4, abs(1 - theta * mu / (mu + 1))
    
    ret = (2 - theta) / 4
    ret = ret * ((2 - theta) * mu * (beta + 1+ theta * beta * (1 - mu))
    ret = ret * ((2 - theta) * beta * (mu + 1+ theta * mu * (1 - beta))
    ret = ret / mu / beta
    ret = ret / (2 * mu * beta * (1 - theta) + (2 - theta) * (mu + beta + 1))
    return 5, math.sqrt(ret)
 
def THM_4_1_unscaled(alpha, theta, mu, beta):
    return THM_4_1_scaled(theta, mu * alpha, beta / alpha)
 
def gen_Base():
    ret = np.zeros((44))
    ret[03= 1
    ret[30= 1
    ret[33= 1
    return ret
 
def gen_M_theta():
    ret = np.zeros((44))
    ret[13= 1
    ret[23= -1
    ret[31= 1
    ret[32= -1
    return ret
 
def gen_M_I():
    ret = np.zeros((44))
    ret[00= 1
    return ret
 
def gen_M_mu_A(alpha, mu):
    ret = [[0-1/200], 
           [-1/2-(1+alpha*mu), 10],
           [0100],
           [0000]]
    return np.array(ret)
 
def gen_M_beta_B(alpha, beta):
    ret = [[-beta/alpha, 0, beta/alpha + 1/20], 
           [0000],
           [beta/alpha + 1/20-beta/alpha-10], 
           [0000]]
    return np.array(ret)
 
 
def opt_val(alpha, mu, beta, retrieve = False):
    lambda_mu_A = cp.Variable()
    lambda_beta_B = cp.Variable()
    theta = cp.Variable()
    rho_sq = cp.Variable()
 
    Base = gen_Base()
    M_theta = gen_M_theta()
    M_I = gen_M_I()
    M_mu_A = gen_M_mu_A(alpha, mu)
    M_beta_B = gen_M_beta_B(alpha, beta)
 
    constraints = []
    constraints.append(lambda_mu_A >= 0)
    constraints.append(lambda_beta_B >= 0)
    constraints.append(theta >= 0)
    constraints.append(Base + theta * M_theta - lambda_mu_A * M_mu_A - lambda_beta_B * M_beta_B + rho_sq * M_I >> 0)
 
    objective = cp.Minimize(rho_sq)
 
    problem = cp.Problem(objective, constraints)
    problem.solve()
 
    if retrieve == True:
        return math.sqrt(problem.value), theta.value
    else:
        return math.sqrt(problem.value)
 
mu = 0.53
beta = 1.35
 
alpha_L = 0.05
alpha_R = 3.95
 
print("Begin Ternary Search")
while alpha_R - alpha_L >= 0.0001:
    left = alpha_L + (alpha_R - alpha_L) / 3
    right = alpha_R - (alpha_R - alpha_L) / 3
    if opt_val(left, mu, beta) < opt_val(right, mu, beta):
        alpha_R = right
    else:
        alpha_L = left
 
opt_alpha = (alpha_L + alpha_R) / 2
opt_rho, opt_theta = opt_val(opt_alpha, mu, beta, True)
 
print("Optimal alpha, theta", opt_alpha, opt_theta)
print("Optimal contraction factor", opt_rho)
print("Check with Theorem 4.1.", THM_4_1_unscaled(opt_alpha, opt_theta, mu, beta))
 
print("Begin Graphing")
= np.linspace(0.253.7550)
= np.array([opt_val(alpha, mu, beta) for alpha in X])
 
plt.plot(X, Y)
plt.show()
 
plt.savefig("Graph.PNG")
 
print("Finished!")
cs

 

여기서 얻는 결론은 $\mu = 0.53$, $\beta = 1.35$일 때, optimal parameter는 $$\alpha_{*} = 1.5949751, \quad \theta_{*} = 1.4244099$$이고, 이때 얻어지는 optimal contraction factor는 $$\rho_{*} = 0.5010598$$이라는 것이다. 또한, 이 경우에 대한 $\rho^2_{*} (\alpha)$의 그래프를 그려서 unimodal function임을 확인할 수 있다. 그래프는 아래와 같다.

당연하지만, $x$축은 $\alpha$의 값, $y$축은 $\rho^2_{*}(\alpha)$의 값이다. 이렇게 Task 2와 논문 정리도 끝 :)

이제 본격적으로 논문을 읽어보자. 논문에서는 예시 문제를 Davis-Yin Splitting (DYS)로 두고 문제를 해결한다.

하지만 이 글에서는 예시 문제를 Douglas-Rachford Splitting (DRS)로 두고 해결하고자 한다. 이유는

  • 일단 나도 손계산 연습을 해보고 싶고, operator가 3개인 DYS보다는 2개인 DRS가 간단함
  • 앞선 글에서 DRS를 어느 정도 설명했으니 독자가 이해하기도 편할 것
  • 이 논문의 Theorem 4에서 DRS에 대한 정리를 하는데, 직접 코드를 짜서 verify 해보고 싶음
  • 그리고 솔직히 논문 내용 그대로 옮기면 재미없음 ㄹㅇㅋㅋ 대신 흐름은 그대로 유지할 것

이번 글에서는 다음 주제를 다룬다.

  • Operator Splitting Performance Estimation Problem을 SDP로 변환하는 방법
  • OSPEP로 얻어진 SDP의 Dual Problem과 그 의미 및 활용법

 

Part 1 : Operator Splitting Performance Estimation Problem to SDP

먼저 DRS를 다시 볼 필요가 있다. 우리는 $R_{\alpha A}R_{\alpha B}$가 nonexpansive operator이므로, 이를 $I$와 함께 섞어서 $$\frac{1}{2} I + \frac{1}{2} R_{\alpha A} R_{\alpha B}$$로 fixed point iteration을 돌렸다. 하지만 굳이 이렇게 섞을 필요는 당연히 없고, $\theta \in (0, 2)$에 대하여 $$\frac{2 - \theta}{2} I + \frac{\theta}{2} R_{\alpha A} R_{\alpha B}$$로 섞어도 된다. 여기에 $R_{\alpha A} = 2J_{\alpha A} - I$, $R_{\alpha B} = 2J_{\alpha B} - I$를 대입하면 $$I - \theta J_{\alpha B} + \theta J_{\alpha A} (2J_{\alpha B} - I)$$를 우리의 operator로 얻는다. 이제 목표는 이 operator가 "얼마나 contractive"한지 파악하는 것이다.

 

우리는 Performance Estimation Problem의 작은 문제에서 $\mathcal{F}_{\mu ,L}$에서 문제를 생각했다.

이번에 다루어 볼 operator class를 알아보자. 앞에서 정의한 네 가지 operator class들과 이들의 교집합을 생각한다.

  • maximal monotone operator의 집합 $\mathcal{M}$
  • $\mu$-strongly monotone이면서 maximal monotone인 operator의 집합 $\mathcal{M}_\mu$
  • $L$-Lipschitz operator의 집합 $\mathcal{L}_L$
  • $\beta$-cocoercive operator의 집합 $\mathcal{C}_{\beta}$

 

이제 본격적으로 문제를 서술할 준비가 되었다. $\mathcal{Q}_1, \mathcal{Q}_2$가 각각 operator class고, (단, 이들은 $\mathcal{M}$의 subset) $$T(z ; A, B, \alpha, \theta) = I - \theta J_{\alpha B} + \theta J_{\alpha A} (2J_{\alpha B} - I)$$이라 하자. 우리가 풀고자 하는 최적화 문제는, 고정된 $\alpha, \theta > 0$에 대하여, $$\begin{equation*} \begin{aligned} \text{maximize} \quad & \frac{ ||T(z;A,B,\alpha, \theta) - T(z';A,B,\alpha,\theta)||^2}{||z-z'||^2} \\ \text{subject to} \quad & A \in \mathcal{Q}_1, B \in \mathcal{Q}_2, z \neq z' \end{aligned} \end{equation*}$$이다. 즉, 최악의 경우에서도 $T$를 $L$-Lipschitz라고 부를 수 있는 가장 최적의 $L$을 구하는 것이다.

 

이제 문제를 조금씩 변형해가며 SDP로 만들어보자. Performance Estimation Problem과 크게 다르지 않다.

 

Step 1 : $T(z; A, B, \alpha, \theta)$를 표현하기 위해서, 다음과 같이 변수를 설정한다. $$z_B = J_{\alpha B} z, \quad z_A = J_{\alpha A}(2z_B - z), \quad T(z;A,B,\alpha, \theta) = z - \theta (z_B - z_A)$$

 

Step 2 : 분모에 있는 $||z-z'||^2$가 불편하다. 이를 일반성을 잃지 않고 $1$로 강제하고 싶은 마음이 자연스럽다.

다행히도 이는 가능하다. Operator class $\mathcal{Q}$가 임의의 $\gamma > 0$에 대하여 $$A \in \mathcal{Q} \iff (\gamma^{-1}I)A(\gamma I) \in \mathcal{Q}$$라면 $\mathcal{Q}$를 homogeneous 하다고 한다. 우리가 다루는 operator class는 전부 homogeneous.

이제 $\gamma = ||z-z'||$라고 하면, $z$들을 전부 $\gamma^{-1}z$로 교체하고 $A, B$를 $(\gamma^{-1} I)A(\gamma I)$, $(\gamma^{-1} I)B(\gamma I)$로 교체하자.

이래도 문제 없음은 homogeneous 성질에서 나온다. 또한, Step 1에서 얻은 식 정리를 그대로 사용해도 괜찮다는 것을 열심히 계산하면 알 수 있다.

결국 $||z-z'||=1$인 경우에서만 문제를 해결해도 무방하다는 사실을 얻는다. 제약 조건을 하나 추가하고, 분모를 날리자.

 

중간 점검을 하자면, 결국 우리가 풀 문제는 $$\begin{equation*} \begin{aligned} \text{maximize} \quad & ||z - \theta (z_B - z_A) - z' + \theta(z'_B - z'_A)||^2 \\ \text{subject to} \quad & A \in \mathcal{Q}_1, B \in \mathcal{Q}_2, ||z-z'||^2 = 1 \\ \quad & z_B = J_{\alpha B} z, z_A = J_{\alpha A}(2z_B - z) \\ \quad & z'_B = J_{\alpha B} z', z'_A = J_{\alpha A}(2z'_B- z') \end{aligned} \end{equation*}$$

 

Step 3 : 이제 슬슬 operator class $\mathcal{Q}$가 눈에 거슬린다. Interpolation의 아이디어를 적용하면 좋겠다.

뒤에서 Interpolation을 적용할 때 알게 되겠지만, 필요한 것은 두 점에 대한 Interpolation 조건이다.

즉, Graph에 $(x_1, q_1), (x_2, q_2)$가 속하는 operator $T \in \mathcal{Q}$가 존재할 조건만이 필요하다. 결론을 제시하면,

  • $\mathcal{M}$ : $\langle q_1 - q_2, x_1 - x_2 \rangle \ge 0$
  • $\mathcal{M}_\mu$ : $\langle q_1 - q_2 , x_1 - x_2 \rangle \ge \mu ||x_1 - x_2||^2$
  • $\mathcal{C}_{\beta}$ : $\langle q_1 - q_2, x_1 - x_2 \rangle \ge \beta || q_1 - q_2||^2$
  • $\mathcal{L}_L$ : $||q_1- q_2||^2 \le L^2 ||x_1 - x_2||^2$
  • 단, $L$-Lipschitz 조건을 넣고 싶다면 무조건 $\mathcal{M} \cap \mathcal{L}_L$ 형태로 operator class를 잡아야 한다
  • $\mathcal{Q}$가 위 operator class의 교집합이라면, 단순히 대응되는 조건을 모두 가져오면 된다.
  • 즉, 두 점이 operator class의 조건을 만족시키면 무조건 interpolation도 된다는 것이 결론. 점의 개수가 많아지면 성립하지 않는다.
  • TMI : monotone operator는 일반적인 경우에서도 maximal monotone operator로 interpolate 할 수 있다 (Zorn's Lemma)

논문의 Theorem 4.1을 직접 검증해보기 위해서, 대응되는 세팅인 $\mathcal{Q}_1 = \mathcal{M}_\mu$, $\mathcal{Q}_2 = \mathcal{M} \cap \mathcal{C}_\beta$를 생각해보자.

  • $B$가 지나야 하는 점은 $(z_B, (z-z_B)/\alpha)$와 $(z'_B, (z'-z'_B)/\alpha)$
  • 그러니 얻는 식은 $\langle ((z-z')-(z_B - z'_B))/\alpha, z_B - z'_B \rangle \ge \beta ||((z-z')-(z_B - z'_B))/\alpha||^2$
  • $A$가 지나야 하는 점은 $(z_A, (2z_B - z-z_A)/\alpha)$와 $(z'_A, (2z'_B-z'-z'_A)/\alpha)$
  • 그러니 얻는 식은 $\langle (2(z_B - z'_B) - (z-z') - (z_A - z'_A))/\alpha, z_A- z'_A \rangle \ge \mu ||z_A - z'_A||^2$

추가적으로, $||z-z'||^2 = 1$이라는 조건도 가지고 있다. 눈치챘겠지만, $z-z'$, $z_A-z_A'$, $z_B - z'_B$과 같이 두 변수의 차로 식이 표현된다.

 

이를 치환하고 정리하면, 우리가 원하는 것은 결국 $$\begin{equation*} \begin{aligned} \text{maximize} \quad & ||z - \theta (z_B - z_A)||^2 \\ \text{subject to} \quad & ||z||^2=1 \\ \quad & \langle z-z_B, z_B \rangle \ge (\beta/\alpha) ||z-z_B||^2 \\ \quad & \langle 2z_B - z , z_A \rangle \ge (1+\alpha \mu)||z_A||^2 \end{aligned} \end{equation*}$$를 얻는다. 이제 SDP를 꺼낼 준비가 완료되었다. 마찬가지로 Gram Matrix $$G = \left[ \begin{matrix} \langle z, z \rangle & \langle z_A, z \rangle & \langle z_B, z \rangle \\ \langle z, z_A \rangle & \langle z_A, z_A \rangle & \langle z_B, z_A \rangle \\ \langle z, z_B \rangle & \langle z_A, z_B \rangle & \langle z_B, z_B \rangle \end{matrix} \right]$$를 꺼내자. 작업하고 있는 공간의 차원이 3 이상이라는 가정하에서, 우리가 아는 것은 $G$가 Positive Semidefinite이라는 사실 뿐이다.

 

이제 다음과 같이 행렬을 정의하자. $$ M_O = \left[ \begin{matrix} 1 & \theta & -\theta \\ \theta & \theta^2 & -\theta^2 \\ -\theta & -\theta^2 & \theta^2 \end{matrix} \right], \quad M_I = \left[ \begin{matrix} 1 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{matrix} \right]$$ $$ M_\beta^B = \left[ \begin{matrix} -\beta/\alpha & 0 & \beta/\alpha + 1/2 \\ 0 & 0 & 0 \\ \beta/\alpha + 1/2 & 0 & -\beta/\alpha - 1 \end{matrix} \right], \quad M_\mu^A = \left[ \begin{matrix} 0 & -1/2 & 0 \\ -1/2 & -(1+\alpha \mu) & 1 \\ 0 & 1 & 0 \end{matrix} \right]$$ 그러면 최종적으로 우리가 원하는 SDP 형태를 얻는다. $$\begin{equation*} \begin{aligned} \text{maximize} \quad & Tr(M_OG) \\ \text{subject to} \quad & G \succeq 0, Tr(M_IG) = 1 \\ \quad & Tr(M_\mu^A G), Tr(M_\beta^BG) \ge 0 \end{aligned} \end{equation*}$$

 

특히, 이 SDP formulation의 해를 가지고 "worst case"에 대응되는 operator $A, B$를 복원할 수 있다.

$Tr(M_O G)$가 최대가 되는 경우의 $G$를 찾았다면, Gram Matrix의 성질을 이용하여 $z, z_A, z_B$를 복원한 다음 interpolate 하면 끝이다.

 

Part 2 : The Dual Problem and its Interpretations/Usage

이제 얻은 SDP의 dual problem을 생각해보자. 전형적인 계산을 하면 결국 $$\begin{equation*} \begin{aligned} \text{minimize} \quad & \rho^2 \\ \text{subject to} \quad & \lambda_\mu^A, \lambda_\beta^B \ge 0 \\ \quad & S = -M_O - \lambda_\mu^A M_\mu^A - \lambda_\beta^B M_{\beta}^B + \rho^2 M_I \succeq 0 \end{aligned} \end{equation*}$$를 얻는다. 또한, "웬만하면" (regularity condition) Strong Duality가 성립함은 Slater's Condition으로 보일 수 있다.

이는 논문의 정리 3.3.으로, 나름 전형적인 증명이고 테크니컬하므로 생략. 중요한 것은 이 Dual Problem을 어떻게 바라보느냐, 그리고 어떻게 써먹느냐다.

 

눈치를 챘겠지만, notation을 $M_{\mu}^A$, $M_{\beta}^B$로 둔 이유는 사실 이 행렬들, 정확하게는 이 행렬들과 관련된 trace에 대한 조건이 $A$의 $\mu$-strongly monotone, $B$의 $\beta$-cocoercivity를 encode 한 것이기 때문이다. 이와 같이 operator class에 $A, B$가 속하게 되면 이로부터 얻어지는 부등식들이 있고, 우리가 증명할 수 있는 것은 기본적으로 이러한 부등식들의 nonnegative combination이다. 이 생각을 이어가면 결국 얻을 수 있는 결론은

  • dual variable $\lambda_\mu^A$, $\lambda_\beta^B$은 $A, B$에 대응되는 부등식들에 대한 가중치이다. (그러니 $0$ 이상이다)
  • feasible point에서의 $\rho$ 값은 해당 가중치를 사용하였을 때 얻어지는 contraction factor이며, 이는 optimal 값에 대한 upper bound다.

즉, "가중치를 이렇게 넣었더니 contraction factor가 이렇게 나왔다"는 사실 "dual problem의 feasible point를 찾았다"와 같다.

 

강조하자면, Primal/Dual Problem의 직관적인 의미는 각각

  • Primal : 변수로 들어가는 것은 "최악의 경우에 해당하는 operator", 얻어지는 것은 "이런 contraction factor는 기대할 수 없다"는 결론
  • Dual : 변수로 들어가는 것은 "증명에 사용될 가중치", 얻어지는 것은 "이런 contraction factor는 얻어질 수 있다"는 결론
  • 특히, Duality Gap이 "웬만하면" (Slater's Condition) 0이므로, 정말 아무런 걱정도 할 필요가 없음

 

이제 Dual Problem의 의미를 알았으니, 활용하는 법도 파악할 수가 있다.

  • Contraction Factor에 대한 Upper Bound를 확보할 수 있다. 애초에 dual problem에서 항상 나오는 결론이다.
  • Dual Problem을 풀면 "증명에 사용될 가중치"를 얻을 수 있다! 그러니 여러 parameter에 대해서 문제를 수치적으로 풀어보고, 이를 통해서 어떤 가중치가 사용되는지 확인한 후, 이에 대한 analytical 공식을 때려맞춘뒤 엄밀한 증명을 시작할 수 (정확하게는, 시도할 수) 있다.
  • Convex Optimization에서 수렴성 증명 중 상당 부분은 정말 단순 계산인데, 그 계산각을 재는 게 어렵다. (후에 Acceleration을 다루면 Lyapunov Function에서 또 이 이야기가 나올 것이다) 부등식들을 섞어먹을 가중치를 선정하는 것도 난이도가 상당하다. Dual Problem은 여기서 도움을 준다.

글 3개에 걸쳐서 arxiv.org/pdf/1812.00146.pdf 논문을 간단하게 정리한다.

저번 글과 비슷한 제목과 비슷한 아이디어를 갖고 있으나, 그 context가 약간 다르다.

 

저번 글에서는 Gradient Descent라는 거의 모두가 아는 context가 있어서 설명이 쉬웠으나, 이번에는 약간 다르다. 

일단 Operator Splitting이 무엇인지를 알아야 context를 이해할 수 있다. 그러니 이번 글에서는 이를 빠르게 설명한다.

  • Disclaimer : 사실 아래에서 $\mathbb{R}^n$으로 말하는 것은 Real Hilbert Space로 보는 것이 맞는 것으로 알고 있다. 굳이 $\mathbb{R}^n$으로 쓰는 이유는 내가 실해석도 모르고 함수해석도 모르고 여러모로 모르는 게 많기 때문. 올해 공부를 하면 어떻게 해결할 수 있지 않을까 기대중이다.
  • 아래 내용의 대부분은 www.math.snu.ac.kr/~ernestryu/courses/optimization.html의 6, 7강에서 공부할 수 있다.

 

set-valued operator는 $\mathbb{R}^n$의 한 점을 $\mathbb{R}^n$의 subset으로 보내는 operator이다.

만약 set-valued operator의 값이 empty set이거나 크기가 1인 집합이라면, 이를 single-valued operator라고 한다.

 

Example (Subgradient) : $f : \mathbb{R}^n \rightarrow \mathbb{R} \cup \{\pm \infty\}$가 convex 함수이고 모든 $y \in \mathbb{R}^n$에 대하여 $$f(y) \ge f(x) + \langle g, y-x \rangle$$이 성립하면 $g$를 $f$의 subgradient at $x$라고 부르며, subgradient의 모음을 $\partial f(x)$라고 쓴다. 이는 set-valued operator.

  • $0 \in \partial f(x)$라면 $x$가 $f$의 값을 최소화하며, 그 역도 성립한다. 
  • $\partial (f+g)(x) \supseteq \partial f(x) + \partial g(x)$임은 자명. 단, 집합끼리의 덧셈은 Minkowski Sum.
  • $\text{dom} f \cap \text{int } \text{dom} g \neq \emptyset$이면 $\partial (f+g)(x) = \partial f(x) + \partial g(x)$. 증명은 여기. 단 domain은 $f(x) < \infty$인 점들의 집합.

 

함수에 대하여 Convex, $L$-smooth, $\mu$-strongly convex 등을 정의하는 것을 operator에도 비슷하게 할 수 있다.

operator $A$가 있다. 각 $x, y \in \mathbb{R}^n$와 $u \in Ax$, $v \in Ay$에 대하여 ($A(x), A(y)$를 $Ax, Ay$라고 간단하게 표기한다)

  • $\langle u-v, x-y \rangle \ge 0$이면, $A$를 monotone operator라 한다.
  • $\langle u-v, x-y \rangle \ge \mu ||x-y||^2$이면, $A$를 $\mu$-strongly monotone operator라 한다.
  • $\langle u-v, x-y \rangle \ge \beta||u-v||^2$이면, $A$를 $\beta$-cocoercive라 한다.
  • $||u-v|| \le L||x-y||$라면 $A$를 $L$-Lipschitz라 한다. 물론, $\mu, \beta, L > 0$.
  • 특히, $1$-Lipschitz인 operator를 nonexpansive operator라고 한다.

 

operator $A$에 대응되는 graph를 $\{(x, u) : u \in Ax\}$라고 정의하자.

또한, operator $A$가 다음 조건을 만족하면, $A$를 maximal monotone 하다고 정의한다.

  • $A$는 monotone operator
  • operator $B$가 다음 조건을 만족하면, $A = B$
  • $B$는 monotone operator이며, $B$의 graph는 $A$의 graph를 포함한다.

 

operator $A$의 inverse를 $\{(u, x) : u \in Ax\}$를 graph로 갖는 operator라고 정의하자.

$A$의 resolvent operator를 $J_A = (I+A)^{-1}$로 정의한다. 물론, $I$는 identity operator이다.

또한, reflected resolvent operator를 $R_A = 2J_A - I$로 정의한다. resolvent는 뒤에서 중요한 역할을 한다.

 

 

최적화 문제와의 관련성을 본격적으로 논의해보자. 기본적으로, 다음과 같은 세팅을 거친다.

  • 최적화 문제를 operator에 대한 문제로 바꾼다. 예를 들어, $f$의 최솟값을 구하는 문제는 $\partial f$의 zero를 찾는 문제가 된다.
  • 이를 operator $T$의 fixed point를 찾는 문제로, 즉 $z \in Tz$인 $z$를 찾는 문제로 바꾼다. 이를 fixed point encoding이라고 부른다.
  • operator $T$의 fixed point를 찾기 위해서, $z^{k+1} = Tz^k$를 반복한다. 이것이 수렴하면 성공이다.

operator $T$가 identity operator $I$와 nonexpansive operator $N$의 convex combination이라면 (단, $I$에 대응되는 비율이 nonzero) $T$를 averaged operator라고 부른다. 이 경우, $T$가 fixed point를 가진다면 $z^{k+1} = Tz^k$가 잘 수렴함을 증명할 수 있다. 증명 방식은 나중에 다루도록 하겠다.

특히, $T$가 $L$-Lipschitz고 $L<1$이라면, 수렴 속도가 기하급수적으로 빠르다는 것은 쉽게 증명할 수 있는 사실이다.

 

이제 fixed point encoding의 예시를 몇 개 보이는 것으로 이론적인 준비를 마치도록 하겠다.

 

Proximal Point Method : $A$가 maximal monotone이면 다음이 성립한다.

  • $J_A$가 averaged operator : $R_A = 2J_A - I$가 nonexpansive 함을 쉽게 보일 수 있다. 특히, $J_A$는 single-valued.
  • Minty Surjectivity Theorem : $J_A$의 domain은 $\mathbb{R}^n$이다. 즉, 모든 $x \in \mathbb{R}^n$에 대하여 $J_Ax$는 nonempty.

또한, $f$가 convex, closed, proper (CCP) function이라면, $\partial f$가 maximal monotone operator임을 보일 수 있다.

  • Convex : $f(\theta x + (1-\theta)y) \le \theta f(x) + (1-\theta)f(y)$ for $\theta \in [0, 1]$
  • Closed : $\{(x, \alpha) \in \mathbb{R}^n \times \mathbb{R} : f(x) \le \alpha\}$ is closed
  • Proper : $f(x) \neq -\infty$ always, and $f(x) < \infty$ somewhere

이제 준비를 할 수 있다. $f$가 CCP function이라고 하면, 

  • $f$의 최솟값을 구하는 문제를 $\partial f$의 zero를 찾는 문제로 바꿀 수 있다.
  • 이를 $J_{\alpha \partial f} = (I + \alpha \partial f)^{-1}$의 fixed point를 찾는 문제로 바꿀 수 있다. 단, $\alpha > 0$.
  • 그러니 $x^{k+1} = J_{\alpha \partial f}(x^k)$를 반복하여 문제를 해결할 수 있다. 이것이 Proximal Point Method.

 

Douglas-Rachford Splitting (DRS) : 이번에는 maximal monotone operator $A, B$가 있다고 가정하고, $$ 0 \in (A+B)x$$인 $x$를 찾는 문제를 생각하자. 열심히 계산하면, 이는 결국 $$x = J_{\alpha B}z, \quad R_{\alpha A}R_{\alpha B} z = z$$인 $(x, z)$를 찾는 문제와 같음을 알 수 있다. ($\alpha > 0$) 그러니 fixed point encoding을 얻는다.

그런데 $R_{\alpha A} R_{\alpha B}$는 nonexpansive 하지만 averaged operator는 아닐 수 있다. 그러므로, 이 문제를 해소하기 위하여 $$z^{k+1} = \left( \frac{1}{2} + \frac{1}{2} R_{\alpha A} R_{\alpha B} \right) z^k$$라고 하면 이는 averaged operator를 사용하였으니 수렴이 보장되는 fixed point iteration이 된다. 이를 통해서 문제를 해결할 수 있다.

이처럼 operator 여러 개가 섞인 문제를 쪼개서 푸는 것을 operator splitting이라고 한다.

 

예고 : 목표는 Operator Splitting Method에 대해서 Performance Estimation을 하는 것이다. 예를 들어, DRS에서는 iteration을 averaged operator $$ \frac{1}{2} I + \frac{1}{2} R_{\alpha A} R_{\alpha B}$$를 이용하여 했었다. $A, B$가 단순히 maximal monotone이 아니라 추가적인 성질을 가지고 있다면, 이 operator가 $L$-Lipschitz가 ($L<1$) 될 수도 있을 것이다. 이때, 우리가 보장할 수 있는 최선의 $L$ 값은 얼마일까? 그리고 이러한 문제를 어떻게 응용할 수 있을까? 다음 글에서 알아보도록 하자.

이 글에서는 논문 www.di.ens.fr/~ataylor/share/PESTO_CDC_2017.pdf를 정리한다.

(지금 보니까 저 논문은 toolkit 소개고, 아이디어 자체는 arxiv.org/pdf/1502.05666.pdf를 보는 게 더 좋을 것 같다)

 

보통 우리는 문제가 주어지면

  • 이를 해결하는 알고리즘을 설계하고
  • 최악의 경우에 어떻게 작동하는지 분석

한다. 만약 PS를 하는 사람이라면 답을 정확하게 구하는 알고리즘을 설계한 후, 최악의 경우의 시간복잡도를 분석할 것이다.

비슷하게, 최적화에서는 보통 답을 iterative 하게 근사하는 알고리즘을 설계한 후, 최악의 경우의 수렴성을 분석한다.

이 논문은 "최악의 경우의 수렴성"을 분석하는 자동화된 toolkit을 제시하는 논문이다. toolkit 링크는 여기다. Matlab으로 작성되어 있다.

 

toolkit의 기본 작동 원리를 알기 위해서, 작은 문제를 하나 풀어보자.

 

문제 (Informal) : $\mu$-strongly convex, $L$-smooth function $F$가 있다. $F$를 최소화하는 문제를 풀고자 한다. Initial point $x_0$에서 Gradient Descent를 적용하여 $x_1 = x_0 - \gamma \nabla F(x_0)$를 얻었을 때, $||\nabla F(x_1)||$의 크기를 얼마나 잘 줄일 수 있을까? 최악의 경우를 비교하여, 최선의 $\gamma$를 구할 수 있을까?

  • $\mu$-strongly convex, $L$-smooth의 의미는 이 글이 글에서 매우 잘 정리되어 있다. 
  • $\mu$-strongly convex, $L$-smooth function의 집합을 $\mathcal{F}_{\mu, L}$이라고 쓴다.

이 문제를 엄밀하게 서술하려면, 이를 최적화 문제 형태로 서술해야 한다. $$\begin{equation*} \begin{aligned} \text{maximize}  \quad  & ||\nabla F(x_1)||^2 \\  \text{subject to} \quad & F \in \mathcal{F}_{\mu, L} \\ \quad & x_1 = x_0 - \gamma \nabla F(x_0) \end{aligned} \end{equation*}$$ 그런데 이렇게만 쓰면 답이 $\infty$인 것이 뻔하므로, 추가적인 제약조건이 필요하다. 예를 들어, $|| \nabla F(x_0) ||^2$의 크기에 제약조건을 걸 수 있다. 그러면 $$\begin{equation*} \begin{aligned} \text{maximize}  \quad  & ||\nabla F(x_1)||^2 \\  \text{subject to} \quad & F \in \mathcal{F}_{\mu, L} \\ \quad & x_1 = x_0 - \gamma \nabla F(x_0) \\ \quad & ||\nabla F(x_0)||^2 \le R^2 \end{aligned} \end{equation*}$$ 라는 최적화 문제를 얻는다. 그런데 $F \in \mathcal{F}_{\mu, L}$가 변수라는 점이 문제 해결을 어렵게 한다. 

 

여기서 멋진 아이디어가 사용된다. 먼저, 우리가 $F$라는 함수에 접근하는 점은 $x_0, x_1$ 두 곳이다. 

그러니 사실 우리는 $f_i = F(x_i)$, $g_i = \nabla F(x_i)$의 값에만 관심을 가져도 충분하다. 이제 최적화 문제를 쓰면 $$\begin{equation*} \begin{aligned} \text{maximize}  \quad  & ||g_1||^2 \\ \text{subject to} \quad & \exists F \in \mathcal{F}_{\mu, L} \text{  s.t.  } F(x_i) = f_i, \nabla F(x_i) = g_i \\ \quad & x_1 = x_0 - \gamma g_0, \text{  } ||g_0||^2 \le R^2 \end{aligned} \end{equation*}$$의 형태가 된다. 그러니, 우리가 관심을 갖게 되는 문제는

  • 언제 $F(x_i) = f_i$, $\nabla F(x_i) = g_i$를 만족하는 $F \in \mathcal{F}_{\mu, L}$이 존재할까?

가 된다. 이런 문제는 Lagrange Interpolation의 향기가 나고, 실제로 Interpolation이라고 불린다. 

그리고 $F \in \mathcal{F}_{\mu, L}$로 interpolation이 되는 조건은 이 논문에서 이미 증명이 되었다.

 

결론은 각 $i \neq j$에 대하여 부등식 $$f_i \ge f_j + \langle g_j, x_i - x_j \rangle + \frac{1}{2(1-\mu / L)} \left( \frac{1}{L} ||g_i - g_j||^2 + \mu ||x_i - x_j||^2  - 2 \frac{\mu}{L} \langle g_j - g_i, x_j - x_i \rangle \right)$$이 성립하는 것이, 원하는 interpolation 조건을 만족하는 $F \in \mathcal{F}_{\mu, L}$이 존재할 필요충분조건이라는 것이다.

 

이제 최적화 문제의 constraint 첫 줄을 interpolation 조건으로 교체할 수 있다. 또한, $x_1 = x_0 - \gamma g_0$를 대입하여 $x_1$ 변수를 제거할 수 있다.

 

남은 문제를 바라보면 다루고 있는 변수가

  • $f_0$, $f_1$ : 이 값들은 constraint에서 그냥 혼자 있는 값들이다.
  • $g_0, g_1, x_0$ : 이 값들은 objective나 constraint나 무조건 다른 값과 내적이 되어있다. 

이 성질 때문에, 문제를 Semidefinite Programming (SDP) 형태로 바꿀 수 있다. Gram Matrix $$G = \left[ \begin{matrix} \langle x_0, x_0 \rangle & \langle g_0, x_0 \rangle & \langle g_1, x_0 \rangle \\ \langle x_0, g_0 \rangle & \langle g_0, g_0 \rangle & \langle g_1, g_0 \rangle \\ \langle x_0, g_1 \rangle & \langle g_0, g_1 \rangle & \langle g_1, g_1 \rangle \end{matrix} \right]$$를 생각하자. Gram Matrix의 성질에 의하여, 다음이 성립한다.

  • $G$는 Positive Semidefinite이며, 특히 풀고 있는 문제가 $\mathbb{R}^d$ 위의 문제면 rank가 최대 $d$.
  • Positive Semidefinite이며 rank가 최대 $d$인 행렬이 주어지면, 이를 Gram Matrix로 보고 vector $x_0, g_0, g_1$을 복원할 수 있음.

그러므로 우리의 모든 문제를 $f_0$, $f_1$과 $G$에 대한 문제로 볼 수 있으며, $G$에 대해 명시할 조건은 $G$가 Positive Semidefinite이라는 점과 rank가 $d$ 이하라는 점이다. 하지만 우리의 관심사는 일반적인 $d$에 대한 것이므로, rank 조건은 제거해도 무방하다. (당장 $d \ge 3$이기만 해도 rank 조건이 무의미하다)

 

이제 남은 문제를 보면 진짜 단순한 Semidefinite Programming 문제고, 이는 standard solver로 해결할 수 있다.

 

논문에 의하면, 이렇게 SDP로 문제를 변형하여 해결할 수 있는 경우는 다음을 포함한다.

  • Class of Functions : Smooth Convex, Strongly Convex, Bounded Subgradients, Bounded Domains, Non-Convex Smooth 등
  • Methods : First-Order Methods : Subgradient, Proximal Operator 등
  • Performance Measure (objective 설정) : $||x_n - x_{*}||^2$, $||\nabla F(x_n)||^2$, $F(x_n) - F(x_{*})$ 등
  • Initial Condition : Performance Measure와 동일함. 물론, $x_n$ 대신 initial point $x_0$에 대한 조건.

그리고 이 toolkit은 이러한 경우 대부분을 지원한다고 한다. 갖고 놀기 좋은 것 같은데, 아직 쓰는 법은 모른다 ㅎㅎ;

 

궁금한 점은 이 toolkit으로 답을 $\mu$, $L$에 대한 수식으로 얻으려면 수식이 변수로 들어가있는 Symbolic한 SDP를 해결해야 할텐데, 이건 대체 어떻게 한 건지 궁금하다. toolkit 보면 수식으로 답을 얻어주는 것 같은데... 관련 논문 읽어봐도 아이디어는 알겠으나 이 점이 계속 헷갈린다. 더 공부하는 걸로.

 

UPD : 직접 돌려본 결과, 수치적으로 계산해준다. 그냥 예시에 답이 이렇게 나와야 한다고 주석으로 써있는 것..

 

arxiv.org/pdf/1502.05666.pdf 논문의 추가적인 내용은

  • Interpolation 관련 엄밀한 증명
  • SDP로 문제를 변환하는 것과, 해당 문제의 성질에 대한 엄밀한 증명
  • SDP로 변환된 문제의 Lagrangian Dual을 생각하여, Upper Bound를 계산
  • 특히, Slater Condition을 이용하여 "웬만하면" Duality Gap이 없음을 증명
  • Performance Estimation Problem의 여러 응용 : 이걸로 문제를 많이 해결/제시하는 것 같다

I participated in UnionCTF as a member of Super Guesser. We reached 1st place :)

 

Mordell Primes

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
from Crypto.Util.number import bytes_to_long
from secrets import k, FLAG
 
assert k < 2^128
assert FLAG.startswith(b'union{')
 
= EllipticCurve(QQ,[0,1,0,78,-16])
= E(1,8)
= k*P
= (k+1)*P
 
= Q[0].numerator()
= R[0].numerator()
 
assert is_prime(p)
assert is_prime(q)
 
= 0x10001
= p*q
= bytes_to_long(FLAG)
= pow(m,e,N)
 
print(f'N = {N}')
print(f'c = {c}')
 
cs

Solution

The main intuition is that the numerator of $kP$ will grow very fast as $k$ increases. This can be verified by simply checking for small $k$.

Therefore, we can try to simply brute force the value of $k$ and see if the numerators multiply to $N$.

 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
= EllipticCurve(QQ,[0,1,0,78,-16])
= E(1,8)
 
for i in range(240):
    Q = i * P
    cc = Q[0].numerator()
    if N % cc == 0:
        print(cc) # p, q
 
= cc
= N / cc
= 65537
phi = (p-1* (q-1)
= inverse(e, phi)
print(long_to_bytes(pow(c, d, N)))
cs

 

Human Server

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
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
import os, random, hashlib, textwrap, json
from Crypto.Cipher import AES
from Crypto.Util.Padding import pad, unpad
from Crypto.Util.number import getPrime, long_to_bytes
 
 
from fastecdsa.curve import secp256k1
from fastecdsa.point import Point
 
FLAG = b'union{XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX}'
 
CURVE = secp256k1
ORDER = CURVE.q
= CURVE.G
 
class EllipticCurveKeyExchange():
    def __init__(self):
        self.private = random.randint(0,ORDER)
        self.public = self.get_public_key()
        self.recieved = None
        self.nonce = None
        self.key = None
 
    def get_public_key(self):
        A = G * self.private
        return A
 
    def send_public(self):
        return print(json.dumps({"Px" : self.public.x, "Py" : self.public.y}))
 
    def receive_public(self, data):
        """
        Remember to include the nonce for ultra-secure key exchange!
        """
        Px = int(data["Px"])
        Py = int(data["Py"])
        self.recieved = Point(Px, Py, curve=secp256k1)
        self.nonce = int(data['nonce'])
 
    def get_shared_secret(self):
        """
        Generates the ultra secure secret with added nonce randomness
        """
        assert self.nonce.bit_length() > 64
        self.key = (self.recieved * self.private).x ^ self.nonce
 
    def check_fingerprint(self, h2: str):
        """
        If this is failing, remember that you must send the SAME
        nonce to both Alice and Bob for the shared secret to match
        """
        h1 = hashlib.sha256(long_to_bytes(self.key)).hexdigest()
        return h1 == h2
 
    def send_fingerprint(self):
        return hashlib.sha256(long_to_bytes(self.key)).hexdigest()
 
def print_header(title: str):
    print('\n\n'+'*'*64+'\n'+'*'+title.center(62)+'*\n'+'*'*64+'\n\n')
 
def input_json(prompt: str):
    data = input(prompt)
    try:
        return json.loads(data)
    except:
        print({"error""Input must be sent as a JSON object"})
        exit()
 
def encrypt_flag(shared_secret: int):
    iv = os.urandom(16)
    key = hashlib.sha1(long_to_bytes(shared_secret)).digest()[:16]
    cipher = AES.new(key, AES.MODE_CBC, iv)
    ciphertext = cipher.encrypt(pad(FLAG, 16))
 
    data = {}
    data['iv'= iv.hex()
    data['encrypted_flag'= ciphertext.hex()
    return print(json.dumps(data))
 
 
Alice = EllipticCurveKeyExchange()
Bob = EllipticCurveKeyExchange()
 
print_header('Welcome!'
message = "Hello! Thanks so much for jumping in to help. Ever since everyone left WhatsApp, we've had a hard time keeping up with communications. We're hoping by outsourcing the message exchange to some CTF players we'll keep the load down on our servers... All messages are end-to-end encrypted so there's no privacy issues at all, we've even rolling out our new ultra-secure key exchange with enhanced randomness! Again, we really appreciate the help, feel free to add this experience to your CV!"
welcome = textwrap.fill(message, width=64)          
print(welcome)
 
print_header('Alice sends public key')
Alice.send_public()
 
print_header("Please forward Alice's key to Bob")
alice_to_bob = input_json('Send to Bob: ')
Bob.receive_public(alice_to_bob)
 
print_header('Bob sends public key')
Bob.send_public()
 
print_header("Please forward Bob's key to Alice")
bob_to_alice = input_json('Send to Alice: ')
Alice.receive_public(bob_to_alice)
            
Alice.get_shared_secret()
Bob.get_shared_secret()
 
print_header('Key verification in progress')
alice_happy = Alice.check_fingerprint(Bob.send_fingerprint())
bob_happy = Bob.check_fingerprint(Alice.send_fingerprint())
if not alice_happy or not bob_happy:
    print({"error""Alice and Bob panicked: Potential MITM attack in progress!!"})
    exit()
 
print_header('Alice sends encrypted flag to Bob')
encrypt_flag(Alice.key)
 
 
cs

Solution

We have to do a MITM attack, but there are some twists - Alice and Bob do a check later to see if they have the same shared secret.

Suppose Alice's secret key, public key is $d_1$, $d_1G$ and Bob's secret key, public key is $d_2$, $d_2G$.

 

Now let's forward Alice's key as $e_1G$ and give the nonce $r_1$ to Bob, where we know the value $e_1$.

This will make Bob's shared secret equal to $(e_1d_2G).x \oplus r_1$. Since $d_2G$ is a known value and $e_1, r_1$ are known, we can calculate this.

 

Similarly, let's forward Bob's key as $e_2G$ and give the nonce $r_2$ to Alice, where we know the value $e_2$.

This will make Alice's shared secret equal to $(e_2d_1G).x \oplus r_2$. We have to make this equal to Bob's shared secret. 

 

This can be done by selecting $r_2$ as $$r_2 = (e_1d_2G).x \oplus (e_2d_1G).x \oplus r_1$$ and we can pass the check between Alice and Bob. We now know the shared secret, so the challenge is solved.

 

Of course, we can just choose $e_1 = e_2 = 1$ and some random value for $r_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
63
64
65
66
67
68
69
70
71
# curve parameter
= 0x79BE667EF9DCBBAC55A06295CE870B07029BFCDB2DCE28D959F2815B16F81798
= 0x483ADA7726A3C4655DA4FBFC0E1108A8FD17B448A68554199C47D08FFB10D4B8
# random nonce
nonce = 0x7918768917649876918679818976981769817691968917698769
 
= remote('134.122.111.232'54321)
r.recvuntil("Alice sends public key")
r.recvline()
r.recvline()
r.recvline()
 
AliceKey = json.loads(r.recvline())
AX = AliceKey["Px"]
AY = AliceKey["Py"]
 
r.recvuntil("Please forward Alice's key to Bob")
r.recvline()
r.recvline()
r.recvline()
 
data = {
    "Px" : X,
    "Py" : Y,
    "nonce" : nonce
}
 
r.send(json.dumps(data))
 
r.recvuntil("Bob sends public key")
r.recvline()
r.recvline()
r.recvline()
 
BobKey = json.loads(r.recvline())
BX = BobKey["Px"]
BY = BobKey["Py"]
 
r.recvuntil("Please forward Bob's key to Alice")
r.recvline()
r.recvline()
r.recvline()
 
 
nonce2 = nonce ^ AX ^ BX
 
data = {
    "Px" : X,
    "Py" : Y,
    "nonce" : nonce2
}
 
r.send(json.dumps(data))
 
shared_secret = BX ^ nonce
 
r.recvuntil("Alice sends encrypted flag to Bob")
r.recvline()
r.recvline()
r.recvline()
 
fin = json.loads(r.recvline())
 
iv = bytes.fromhex(fin["iv"])
enc = bytes.fromhex(fin["encrypted_flag"])
 
key = hashlib.sha1(long_to_bytes(shared_secret)).digest()[:16]
cipher = AES.new(key, AES.MODE_CBC, iv)
flag = cipher.decrypt(enc)
 
print(flag)
cs

 

Cr0wn St3rling

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
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
#!/usr/bin/env python3
 
from secrets import flag, musical_key
from Crypto.Util.number import isPrime
import math
 
 
def sieve_for_primes_to(n):
    # Copyright Eratosthenes, 204 BC
    size = n//2
    sieve = [1]*size
    limit = int(n**0.5)
    for i in range(1, limit):
        if sieve[i]:
            val = 2*i+1
            tmp = ((size-1- i)//val
            sieve[i+val::val] = [0]*tmp
    return [2+ [i*2+1 for i, v in enumerate(sieve) if v and i > 0]
 
 
def is_quasi_prime(n, primes):
    # novel class of semi-prime numbers
    # https://arxiv.org/pdf/1903.08570.pdf
    p2 = 0
    for p1 in primes:
        if n % p1 == 0:
            p2 = n//p1
            break
    if isPrime(p2) and not p1 in [23and not p2 in [23]:
        return True
    return False
 
 
def bbp_pi(n):
    # Bailey-Borwein-Plouffe Formula
    # sounds almost as cool as Blum-Blum-Shub
    # nth hex digit of pi
    def S(j, n):
        s = 0.0
        k = 0
        while k <= n:
            r = 8*k+j
            s = (s + pow(16, n-k, r) / r) % 1.0
            k += 1
        t = 0.0
        k = n + 1
        while 1:
            newt = t + pow(16, n-k) / (8*k+j)
            if t == newt:
                break
            else:
                t = newt
            k += 1
        return s + t
 
    n -= 1
    x = (4*S(1, n) - 2*S(4, n) - S(5, n) - S(6, n)) % 1.0
    return "%02x" % int(x * 16**2)
 
 
def digital_root(n):
    # reveals Icositetragon modalities when applied to Fibonacci sequence
    return (n - 1) % 9 + 1 if n else 0
 
 
def fibonacci(n):
    # Nature's divine proportion gives high-speed oscillations of infinite
    # wave values of irrational numbers
    assert(n >= 0)
    if n < digital_root(2):
        return n
    else:
        return fibonacci(n - 1+ fibonacci(n - 2)
 
 
def is_valid_music(music):
    # Leverage music's infinite variability
    assert(all(c in "ABCDEFG" for c in music))
 
 
def is_valid_number(D):
    # Checks if input symbolizes the digital root of oxygen
    assert(8==D)
 
 
def get_key(motif):
    is_valid_music(motif)
    is_valid_number(len(motif))
    # transpose music onto transcendental frequencies
    indexes = [(ord(c)-0x40)**for i, c in enumerate(motif)]
    size = sum(indexes)
    assert(size < 75000# we will go larger when we have quantum
    return indexes, size
 
 
def get_q_grid(size):
    return [i for i in range(size) if is_quasi_prime(i, sieve_for_primes_to(math.floor(math.sqrt(size))+1))]
 
 
if __name__ == "__main__":
    print("[+] Oscillating the key")
    key_indexes, size = get_key(musical_key)
    print("[+] Generating quasi-prime grid")
    q_grid = get_q_grid(size)
    # print(f"indexes: {key_indexes}  size: {size}  len(q_grid): {len(q_grid)}")
 
    out = []
    for i, p in enumerate(flag):
        print(f"[+] Entangling key and plaintext at position {i}")
        index = key_indexes[i % len(key_indexes)] * fibonacci(i)
        q = q_grid[index % len(q_grid)]
        key_byte_hex = bbp_pi(q)
        # print(f"index: {index:10}  fib: {fibonacci(i):10}  q-prime: {q:10}  keybyte: {key_byte_hex:10}")
        out.append(ord(p) ^ int(key_byte_hex, 16))
 
    print(f"[+] Encrypted: {bytes(out).hex()}")
 
cs

Solution

We start by making some quick observations -

  • The length of key_indexes must be 8 due to the check is_valid_number.
  • There are only $7^8$ keys since there are 7 possible choices for each entry of the musical key.
  • Actually, the first entry of key_indexes is guaranteed to be 1. Make that $7^7$ valid key_indexes.
  • Fibonacci function is implemented pretty badly, and we can change that to simple dynamic programming
  • There's nothing we need to do about the $\pi$ function
  • The q_grid function is suboptimal, but it runs in decent time so no need to worry about that.

We will now decrease the number of key_index candidates by using the known plaintext "union".

  • Each entry of key_indexes works in a separable way, so we can work with each entry separately
  • There are 7 possible choices for a specific key_indexes entry
  • The problem is we do not know the length of q_grid (we do not know the variable $size$ in advance)
  • For $i < 5$, the maximum value of the index is $7^4 \cdot 3$, which is small compared to the length of the q_grid for $size = 75000$.
  • Therefore, we can reasonably assume that the length of q_grid is larger than the index.
  • Now we can try out each candidate for a specific key_indexes entry - this is sufficient to get the first 5 entries of key_indexes.

We now have only $7^3$ candidates to work with. We can try every single one of them.

  • The only slow part of the brute force is generating the q_grid and calculating its length.
  • This can be done very fast by calculating q_grid for $size = 75000$ and using binary search to calculate the "actual length of the q_grid" (i.e. number of entries that are less than the actual "size" variable)
  • However, since we have a lot of time, we can just calculate q_grid for $size = 75000$ then simply calculate the actual length of the q_grid with brute force. It doesn't take too long for us to get the flag anyway..
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
q_grid = get_q_grid(75005
enc = bytes.fromhex('76f64667220717784affa07cf6b8be52c7d8348d778a41615efa9e53f2566b27fd96eb984c08')
fib = [01]
for i in range(250):
    fib.append(fib[i-1+ fib[i-2])
 
# the first note of the musical key does not matter since the key_index will become 1 anyway
 
# known plaintext
ptxt = "union"
for i in range(05):
    if i == 0# first note doesn't matter
        continue
    for j in range(18): # try all 7 possible notes
        idx = (j ** i) * fib[i]
        q = q_grid[idx] # we assume idx < actual length of q_grid
        key_byte_hex = bbp_pi(q)
        out = enc[i] ^ int(key_byte_hex, 16)
        if out == ord(ptxt[i]): # match
            print(chr(64 + j)) # output the note
 
# results in C D A D -> set the first five notes as A C D A D
 
 
# brute force 7^3
for i in range(07):
    for j in range(07):
        for k in range(07):
            T = ['A''C''D''A''D', chr(65+i), chr(65+j), chr(65+k)]
            key_indexes , size = get_key(T)
            # remove the assertion for size in get_key
            if size >= 75000:
                continue
            # compute the "actual" q_grid
            q_qgrid = []
            for val in q_grid:
                if val < size:
                    q_qgrid.append(val)
            # compute the plaintext
            out = []
            for ii in range(0len(enc)):
                idx = key_indexes[ii % 8* fib[ii]
                q = q_qgrid[idx % len(q_qgrid)]
                key_byte_hex = bbp_pi(q)
                out.append(enc[ii] ^ int(key_byte_hex, 16))
            print(bytes(out)) # wait for it..
cs

 

Neo-Classical Key Exchange

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
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
import os
from hashlib import sha1
from random import randint
from Crypto.Cipher import AES
from Crypto.Util.Padding import pad, unpad
 
FLAG = b"union{XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX}"
 
def list_valid(l):
    x = l // 2
    checked = set([x])
    while x * x != l:
        x = (x + (l // x)) // 2
        if x in checked: return False
        checked.add(x)
    return True
 
def list_iter(n):
    x = n
    y = (x + 1// 2
    while y < x:
        x = y
        y = (x + n // x) // 2
    return x
 
def list_mul(l1,l2,p):
    X, Y = len(l1), len(l2)
    Z = list_iter(X)
    assert X == Y
    assert list_valid(X)
    l3 = []
    for i in range(Z):
        for j in range(Z):
            prod_list = [x*for x,y in zip(l1[Z*i:Z*(i+1)], l2[j::Z])]
            sum_list = sum(prod_list) % p
            l3.append(sum_list)
    return l3
 
def list_exp(l0, e, p):
    exp = bin(e)[3::]
    l = l0
    for i in exp:
        l = list_mul(l,l,p)
        if i == '1':
            l = list_mul(l,l0,p)
    return l
 
def gen_public_key(G,p):
    k = randint(2,p-1)
    B = list_exp(G,k,p)
    return B,k
 
def gen_shared_secret(M,k,p):
    S = list_exp(M,k,p)
    return S[0]
 
def encrypt_flag(shared_secret: int):
    # Derive AES key from shared secret
    key = sha1(str(shared_secret).encode('ascii')).digest()[:16]
    # Encrypt flag
    iv = os.urandom(16)
    cipher = AES.new(key, AES.MODE_CBC, iv)
    ciphertext = cipher.encrypt(pad(FLAG, 16))
    # Prepare data to send
    data = {}
    data['iv'= iv.hex()
    data['encrypted_flag'= ciphertext.hex()
    return data
 
= 64050696188665199345192377656931194086566536936726816377438460361325379667067
= [37474442957545178764106324981526765864975539603703225974060597893616967420393,59548952493843765553320545295586414418025029050337357927081996502641013504519311002066525512164709938000874013049550644788296268367056724529039089424037491386031482454287572407012381137953191558164465623529992046661815621863204773420708638990322428536520731257757756431087939910637434308755686013682215836263249525491465214495369731073552931306211582961157162030422899032923981311376221021836681925641294064263844659958138617889034069800460358141630174638641532727035735045369266322629011656427579578656066165031820538678953220132827396471587929444138998790449514672948945562632375967133202143205396916253265051473730587605323370564700860148975988622662724284710157566957213620913591119857266366110422436201592848917393008725709237548443793017124298122562856326649394385871891424162512325919431373810173511592710316040978823523358078859249102260718794394402264910240234942692440221176187631522440611503354694959423849000390378955527116778192120808910199353602982871650771627510964301491382871751987924260652304319543941114891709993329929124030812683307477971502992612959253926958823705349101783144068766123926443603026261539055007453105405205925131925190516128252182445043410788021004743874404334678085306701603781467743100927869431963764735143299058921864701886615585381428010877330553242342651823130483453772716228097418145796292233111277774478016673520810772503991055566747628629443375207256050745127045919163601367018956550763591458462169205918826786898398213162402878653481728846096771599941966230969939621034765177430371547059243127032356850437797415676110660436413346535063433156355547532408592015995190002391616368774565349584890853755466839699622482020499285870283811476739960099513665661150287720594400034444826365313288645670526357669076978338398633256587,23887025917289715287437926862183042001010845671403682948840305587666551788353]
A,a = gen_public_key(G,p)
B,b = gen_public_key(G,p)
assert gen_shared_secret(A,b,p) == gen_shared_secret(B,a,p)
 
shared_secret = gen_shared_secret(B,a,p)
encrypted_flag = encrypt_flag(shared_secret)
 
print(f"Alice's public key: {A}"
print(f"Bob's public key: {B}")
print(f"Encrypted flag: {encrypted_flag}")
 
 
 
 
cs

Solution

The first thing we notice is that list_valid checks if the length of the list is a square, and list_iter simply returns the square root.

Now we look at the main function, the list multiplication. If we stare at it for a while (and notice the "length is a square" stuff) we can find out that this list multiplication is just a matrix multiplication, where the list is the matrix entries in order.

 

Now we have to solve a Diffie-Hellman problem with matrices over $\mathbb{F}_p$. We will do it by solving discrete logarithm problem.

  • Idea 1 : The range of the secret keys
  • The range of the secret keys are from $2$ to $p$ - this implies that we do not need the entire discrete logarithm, just some part of it.
  • Idea 2 : The determinant
  • This idea can give us the required information of the (matrix) discrete logarithm if we calculate the respective discrete logarithm (over finite field $\mathbb{F}_p$) on the determinants. However, the primes used here are quite large, which makes this approach infeasible.

Instead of $5 \times 5$ matrices that we actually have to deal with, we first look at $2 \times 2$ matrices over the reals.

The idea is that we want to look at how matrix powers behave in a more relaxed setting, and find some properties we can use in $\mathbb{F}_p$.

 

We try to solve $A^k = B$. Take the Jordan Canonical Form $J$ with $A = PJP^{-1}$. Then we want to solve $$J^k = P^{-1}BP$$ so we may simply assume that $A$ is in a Jordan Canonical Form. If $A$ is in the form of $$A = \left[ \begin{matrix} \lambda_1 & 0 \\ 0 & \lambda_2 \end{matrix} \right]$$ then nothing special happens - we still need a logarithm to compute $k$, which is infeasible in $\mathbb{F}_p$. However, if $A$ is in the form of $$A = \left[ \begin{matrix} \lambda & 1 \\ 0 & \lambda \end{matrix} \right]$$ then we have some interesting phenomenon : $$A^k = \left[ \begin{matrix} \lambda^k & k \lambda^{k-1} \\ 0 & \lambda^k \end{matrix} \right]$$ and now we can get $k$ as $(k\lambda^{k-1} \times \lambda) / \lambda^k$, all of which we know and can be done in $\mathbb{F}_p$.

 

Let's come back to $\mathbb{F}_p$. If we want to use the Jordan Canoncial Form trick, we need our characteristic polynomial to completely split in $\mathbb{F}_p$, and we also need a block of size $2$. Our matrix satisfies all of these properties, so we can just use the exact same technique here as well.

 

Of course, we can see that we have calculated the discrete logarithm modulo $p$, not the full discrete logarithm.

However, because of Idea 1 (secret key range is small) this is enough to solve the problem. Here are some additional notes.

  • The Jordan Canonical Form of our matrix $G$ has 4 blocks of size 1, 1, 1, 2. 
  • This actually shows that the order of $G$ is $p(p-1)$ - this can be seen as the $p-1$ term dealing with the diagonal entries (Fermat's Theorem) and the $p$ term dealing with the off-diagonal term in our block of size 2. 
  • If we want to solve for the discrete logarithm modulo $p-1$, we need discrete logarithm over $\mathbb{F}_p$ which we decided infeasible.
  • I think there are some similarities with this challenge and TetCTF 2021 unevaluated.
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
= 64050696188665199345192377656931194086566536936726816377438460361325379667067
 
def MAT(X):
    M = Matrix(GF(p), 55)
    cnt = 0 
    for i in range(05):
        for j in range(05):
            M[i, j] = X[cnt]
            cnt += 1
    return M
 
= 
= 
= 
 
GMAT = MAT(G)
AMAT = MAT(A)
BMAT = MAT(B)
 
print(GMAT.charpoly().factor())
print(AMAT.charpoly().factor())
print(BMAT.charpoly().factor())
 
 
J, P = GMAT.jordan_form(transformation = True)
AMAT = P^-1 * AMAT * P
BMAT = P^-1 * BMAT * P
print(J)
print("")
print(AMAT)
print("")
print(BMAT)
 
print(AMAT[34/ AMAT[33* J[33]) # dlog
print(BMAT[34/ BMAT[33* J[33]) # dlog
# the rest is relatively trivial (get shared secret, decrypt flag)
cs

 

Why is a raven

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
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
import hashlib
from base64 import b64encode
from Crypto.Cipher import AES
 
= 2^216*3^137 - 1
F.<i> = GF(p^2, modulus=x^2+1)
= EllipticCurve(F, [06010])
 
xQ20=0x000C7461738340EFCF09CE388F666EB38F7F3AFD42DC0B664D9F461F31AA2EDC6B4AB71BD42F4D7C058E13F64B237EF7DDD2ABC0DEB0C6C
xQ21=0x00025DE37157F50D75D320DD0682AB4A67E471586FBC2D31AA32E6957FA2B2614C4CD40A1E27283EAAF4272AE517847197432E2D61C85F5
yQ20=0x001D407B70B01E4AEE172EDF491F4EF32144F03F5E054CEF9FDE5A35EFA3642A11817905ED0D4F193F31124264924A5F64EFE14B6EC97E5
yQ21=0x000E7DEC8C32F50A4E735A839DCDB89FE0763A184C525F7B7D0EBC0E84E9D83E9AC53A572A25D19E1464B509D97272AE761657B4765B3D6
xP20=0x0003CCFC5E1F050030363E6920A0F7A4C6C71E63DE63A0E6475AF621995705F7C84500CB2BB61E950E19EAB8661D25C4A50ED279646CB48
xP21=0x001AD1C1CAE7840EDDA6D8A924520F60E573D3B9DFAC6D189941CB22326D284A8816CC4249410FE80D68047D823C97D705246F869E3EA50
yP20=0x001AB066B84949582E3F66688452B9255E72A017C45B148D719D9A63CDB7BE6F48C812E33B68161D5AB3A0A36906F04A6A6957E6F4FB2E0
yP21=0x000FD87F67EA576CE97FF65BF9F4F7688C4C752DCE9F8BD2B36AD66E04249AAF8337C01E6E4E1A844267BA1A1887B433729E1DD90C7DD2F
xQ30=0x0012E84D7652558E694BF84C1FBDAAF99B83B4266C32EC65B10457BCAF94C63EB063681E8B1E7398C0B241C19B9665FDB9E1406DA3D3846
xQ31=0x0000000
yQ30=0x0000000
yQ31=0x000EBAAA6C731271673BEECE467FD5ED9CC29AB564BDED7BDEAA86DD1E0FDDF399EDCC9B49C829EF53C7D7A35C3A0745D73C424FB4A5FD2
xP30=0x0008664865EA7D816F03B31E223C26D406A2C6CD0C3D667466056AAE85895EC37368BFC009DFAFCB3D97E639F65E9E45F46573B0637B7A9
xP31=0x0000000
yP30=0x0006AE515593E73976091978DFBD70BDA0DD6BCAEEBFDD4FB1E748DDD9ED3FDCF679726C67A3B2CC12B39805B32B612E058A4280764443B
yP31=0x0000000
 
Q_A = E(xQ20 + xQ21*i, yQ20 + yQ21*i)
P_A = E(xP20 + xP21*i, yP20 + yP21*i)
Q_B = E(xQ30 + xQ31*i, yQ30 + yQ31*i)
P_B = E(xP30 + xP31*i, yP30 + yP31*i)
 
# Computes an l^e-isogeny out of E from a set Ss of kernel generators
# as a composition of e l-isogenies
def computeIsogeny(E, Ss, l, e):
    S_tmps = Ss
    E_tmp = E
    ϕ = None
    for k in range(e):
        R_tmps = S_tmps
        for _ in range(e-k-1):
            R_tmps = [ l*R_tmp for R_tmp in R_tmps ]
        ϕ_k = E_tmp.isogeny(kernel=R_tmps)
 
        S_tmps = [ ϕ_k(S_tmp) for S_tmp in S_tmps ]
        E_tmp = ϕ_k.codomain()
        if ϕ is None:
            ϕ = ϕ_k
        else:
            ϕ = ϕ_k * ϕ
    return ϕ
 
k_A = randint(02^216-1)
S_A = P_A + k_A*Q_A
ϕ_A = computeIsogeny(E, [S_A], 2216)
Alice = (ϕ_A.codomain().a_invariants(), ϕ_A(P_B).xy(), ϕ_A(Q_B).xy(), ϕ_A(P_A).xy(), ϕ_A(Q_A).xy())
print(f"Alice = {Alice}")
 
k_B = randint(03^137-1
S_B = P_B + k_B*Q_B
ϕ_B = computeIsogeny(E, [S_B], 3137)
Bob = (ϕ_B.codomain().a_invariants(), ϕ_B(P_A).xy(), ϕ_B(Q_A).xy())
print(f"Bob = {Bob}")
 
(E_B, B_P_A, B_Q_A) = Bob
E_B = EllipticCurve(F, E_B)
B_P_A = E_B(B_P_A)
B_Q_A = E_B(B_Q_A)
B_S_A = B_P_A + k_A*B_Q_A
jAlice = computeIsogeny(E_B, [B_S_A], 2216).codomain().j_invariant()
 
(E_A, A_P_B, A_Q_B, _, _) = Alice
E_A = EllipticCurve(F, E_A)
A_P_B = E_A(A_P_B)
A_Q_B = E_A(A_Q_B)
A_S_B = A_P_B + k_B*A_Q_B
jBob = computeIsogeny(E_A, [A_S_B], 3137).codomain().j_invariant()
 
assert jAlice == jBob
 
flag = open("flag.txt""rb").read().strip()
assert len(flag) == 32
 
sk = hashlib.sha256(str(jAlice).encode('ascii')).digest()[:16]
cipher = AES.new(sk, AES.MODE_CBC)
print(f"iv = {b64encode(cipher.iv)}")
print(f"ct = {b64encode(cipher.encrypt(flag))}")
 
cs

 

Solution

This is Supersingular Isogeny Diffie-Hellman protocol. Looking at the Wikipedia article and back, we see that we have extra information, $$\phi_A(P_A), \phi_A(Q_A)$$ Since $\phi_A$ has $S_A = P_A + k_AQ_A$ as it's kernel, we see that $$ \phi_A(P_A + k_AQ_A) = \phi_A(P_A) + k_A \phi_A(Q_A) = 0$$ Also, since the values here are in a subgroup of order $2^{216}$, we can calculate $k_A$ (i.e. discrete logarithm) efficiently with Pohlig-Hellman.

 

Now that we know $k_A$, one of the secret keys, we can just follow our algorithm to calculate the shared secret.

 

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
63
64
65
66
67
68
69
70
71
72
73
74
75
= 2^216*3^137 - 1
F.<i> = GF(p^2, modulus=x^2+1)
= EllipticCurve(F, [06010])
 
def computeIsogeny(E, Ss, l, e):
    S_tmps = Ss
    E_tmp = E
    ϕ = None
    for k in range(e):
        R_tmps = S_tmps
        for _ in range(e-k-1):
            R_tmps = [ l*R_tmp for R_tmp in R_tmps ]
        ϕ_k = E_tmp.isogeny(kernel=R_tmps)
 
        S_tmps = [ ϕ_k(S_tmp) for S_tmp in S_tmps ]
        E_tmp = ϕ_k.codomain()
        if ϕ is None:
            ϕ = ϕ_k
        else:
            ϕ = ϕ_k * ϕ
    return ϕ
 
xQ20=0x000C7461738340EFCF09CE388F666EB38F7F3AFD42DC0B664D9F461F31AA2EDC6B4AB71BD42F4D7C058E13F64B237EF7DDD2ABC0DEB0C6C
xQ21=0x00025DE37157F50D75D320DD0682AB4A67E471586FBC2D31AA32E6957FA2B2614C4CD40A1E27283EAAF4272AE517847197432E2D61C85F5
yQ20=0x001D407B70B01E4AEE172EDF491F4EF32144F03F5E054CEF9FDE5A35EFA3642A11817905ED0D4F193F31124264924A5F64EFE14B6EC97E5
yQ21=0x000E7DEC8C32F50A4E735A839DCDB89FE0763A184C525F7B7D0EBC0E84E9D83E9AC53A572A25D19E1464B509D97272AE761657B4765B3D6
xP20=0x0003CCFC5E1F050030363E6920A0F7A4C6C71E63DE63A0E6475AF621995705F7C84500CB2BB61E950E19EAB8661D25C4A50ED279646CB48
xP21=0x001AD1C1CAE7840EDDA6D8A924520F60E573D3B9DFAC6D189941CB22326D284A8816CC4249410FE80D68047D823C97D705246F869E3EA50
yP20=0x001AB066B84949582E3F66688452B9255E72A017C45B148D719D9A63CDB7BE6F48C812E33B68161D5AB3A0A36906F04A6A6957E6F4FB2E0
yP21=0x000FD87F67EA576CE97FF65BF9F4F7688C4C752DCE9F8BD2B36AD66E04249AAF8337C01E6E4E1A844267BA1A1887B433729E1DD90C7DD2F
xQ30=0x0012E84D7652558E694BF84C1FBDAAF99B83B4266C32EC65B10457BCAF94C63EB063681E8B1E7398C0B241C19B9665FDB9E1406DA3D3846
xQ31=0x0000000
yQ30=0x0000000
yQ31=0x000EBAAA6C731271673BEECE467FD5ED9CC29AB564BDED7BDEAA86DD1E0FDDF399EDCC9B49C829EF53C7D7A35C3A0745D73C424FB4A5FD2
xP30=0x0008664865EA7D816F03B31E223C26D406A2C6CD0C3D667466056AAE85895EC37368BFC009DFAFCB3D97E639F65E9E45F46573B0637B7A9
xP31=0x0000000
yP30=0x0006AE515593E73976091978DFBD70BDA0DD6BCAEEBFDD4FB1E748DDD9ED3FDCF679726C67A3B2CC12B39805B32B612E058A4280764443B
yP31=0x0000000
 
Q_A = E(xQ20 + xQ21*i, yQ20 + yQ21*i)
P_A = E(xP20 + xP21*i, yP20 + yP21*i)
Q_B = E(xQ30 + xQ31*i, yQ30 + yQ31*i)
P_B = E(xP30 + xP31*i, yP30 + yP31*i)
 
Alice = ((0606001602663961206370988750155271268843249113105575064100544244329723627508642651491029799456469448718421085928092642609765039188242*+ 1712256002728528379082564430864971406334041478918891587209542135490151021806041516369731041317073413033514408943949475062553807198022241966002746626067354539975418232642475646213518284443509300453664841759671344583904727702765870170961054029730219838438930886730*+ 6295409890219768050656401214128285134680958664621604435525993857295024885942522198730338761351461854361939143813589229958450834937), (20693617845250531673791079572257479372406496374051176010221583150895284635664420984163027961195027146723306399733543575074371571546*+ 1257944050482622779039989846116832908011645379538188503193888759983069361986464587598537959460710680527206312828714104047432447257917003339040898697587060167865681315428279965204095022676751761236717662173320135824191474549296911745414760875386583097246892970743*+ 2450227209495560745928392442008559789430024267104386893781343959329588604681368476319376824183170770268590193199446339985032818433), (24196999226902675779045571226331502647086872832197399777068255320010293863359017283213324144431537822661506383681353187559191999771*+ 1403187277399873350729873144059600531393910895713791231342921226797472498403919424333885862617451889202534903916737871843637458172210067956801857468023514754660550671095579147019588599811273848235704360993890497575424172688000039308192770149207724142524545451074*+ 9704956586624341710808951764565861341114293792082683257320510639095567055930904001585984811139137392398321200917960531515122425604), (21482597851178884503353076937735588452957177768806776910812379517549572253101759233236102347370343956671258496673645283509995572850*+ 1409607890280707835592859895681613004561924579015938475117693274564675323421118094150575882783331409169038834693561947366544280925913679392986650554551011681934303695650088628896811869083453967966676089303335417699532232752393700725181942165609165070072195990421*+ 22303973329492411565669001646989446651767650420180177485066511378012378529261175557912535448570067170663048114739927772127080694786), (5031508630808576087782892353323275460174142373365249589846782383660521445945988018058115279743582819518492860550820586178080959929*+ 203618640430880223098324249541882883341295202367378909200012993621765252931980356906286991355846680733796871308320906361027504960035326896702997677262072220524322872052674185107431056651996898750306495168544570294686542579294013185895403025686718275379409582021*+ 7018087072590614715963128743406699936749123349867893045580774389322712378049822434865393438816413618294542843639485193888664986503))
Bob = ((0602510377837984569005668272963938229152759212776314258952545654072169410901298850712372306096701112903052487282410712205434980682770*+ 153361659901846951954864105534225425250790425873535004318638204601924663903808975912970767591961237416790729864700484204884248322513335813103700469160284801960413086144549776993448017319107340684719947118153850729369660724596130930055245047262733708054423015655*+ 17338799868192494405204548102094725343306200019012693785648257092061088427620739800739315465276207804594142456129224673452195357714), (2771195673566835722887342815479686444463738278075014443830431110033775598812266459191044593910730473384513927831673567602258951977*+ 1469564256406938038163605778762348165866304542092994743998859506798641754551769151744125414548886984617946375481700338419227462646318564301293503451778592169644157610474379393936432543000343986957900909392616771402521075243703340903344745798060095728893961976940*+ 19628333731314344646926186187873836263784148023213378217245128148326949516664862760385029489345376891188072162779569669305066964933), (22650214990480384681860370457554162699319780968761610136003430194938807060051700030135993836240770510069117319449743388225281444184*+ 337715599698804126803907287393594153129537868872259108204032694805267651900616891555563288402498128590848003390275824058761569305417806681788782120983952163360625445316482360798557639190977860032884873427321883793075291472918577432082376117267831746467121303568*+ 21533180999838449404329422084189008931697041812999894837076670480949795196804840300494281304611360768987802541355649881398701758313))
 
(E_A, A_P_B, A_Q_B, A_P_A, A_Q_A) = Alice
E_A = EllipticCurve(F, E_A)
A_P_B = E_A(A_P_B)
A_Q_B = E_A(A_Q_B)
A_P_A = E_A(A_P_A)
A_Q_A = E_A(A_Q_A)
 
# print(discrete_log(A_P_A, A_Q_A, ord=2^216, operation='+'))
k_A = 2^216 - 77755655179930472801709066068873804442103726663917450704829188611
(E_B, B_P_A, B_Q_A) = Bob
E_B = EllipticCurve(F, E_B)
B_P_A = E_B(B_P_A)
B_Q_A = E_B(B_Q_A)
B_S_A = B_P_A + k_A*B_Q_A
jAlice = computeIsogeny(E_B, [B_S_A], 2216).codomain().j_invariant()
 
print(jAlice)
 
jAlice = "11056265280320404779614673772059927218754555738553445716393833134024824685344573644810486464601209854921719950024169587404070224902*i + 13483412480684215684949244771895365050495321597976589093887140592592966008044279045153836163302049958922591278477269370463266951423"
sk = hashlib.sha256(jAlice.encode('ascii')).digest()[:16]
 
iv = b'XSglnu+2ZwFuHGE8ddIuJQ=='
iv = base64.b64decode(iv)
ct = b'4VR9ty+lFW6oQoWTVHiDE7A9uKw0KbQzpnCWOGVQXGo='
ct = base64.b64decode(ct)
 
cipher = AES.new(sk, AES.MODE_CBC, iv)
print(cipher.decrypt(ct))
cs

 

'CTF' 카테고리의 다른 글

zer0pts CTF 2021 Crypto Writeups  (1) 2021.03.07
AeroCTF 2021 Crypto Writeups  (0) 2021.02.28
CryptoHack All Solve, Again  (10) 2021.01.31
0x41 CTF SPN Writeup  (0) 2021.01.31
TetCTF 2021 Crypto Writeups  (1) 2021.01.03