Lord Riot's Blog
0%

摸鱼Writeup-3

发表于 更新于 分类于
虎符CTF2022部分赛题复现，因为是赛题复现，就用中文写了orz，虎符的题质量很高，笔者零散时间复现，持续更新ing

HCDSA

求阶

​ writeup中提到题中超椭曲线的阶可以通过这篇paper中提到的算法来计算，读这篇paper后发现这是一个利用格基规约寻找Ideal上生成元的方法，思路是以往没见过的。

论文中一些公式的应用都是源自其他论文的结论，笔者尚未有时间全部梳理，故只讨论论文中涉及的部分

​ 对于\(\mathbb{F}_p\) 上的超椭 \(E:Y^2 + Y = X^n\)，即\(h(X) = 1, f(X) = X^n\)，论文中讨论的情况为p较大，但n很小。设\(\zeta = e^{\frac{2\pi i}{n}}\)，即为复数域\(C\)上的1的一个n次根，根据一个确定的\(\mathbb{F}_p\)上非零非n次根元素\(\alpha\)可以定义一个multiplicative map \(\mathcal{X}(\alpha) = \zeta\)，multiplicative map即为具有\(\mathcal{X}(ab) = \mathcal{X}(a)\mathcal{X}(b)\)性质的映射。\(\mathbb{F}_p\) 上的元素可以表示为\(\alpha^k\)，即映射为\(\zeta^k\)。可以将这个映射推广到\(\mathbb{F}_p\)上，从而可以定义基于此映射的Jacobi sum \(J(\mathcal{X}, \mathcal{X}) = \sum_{y\in \mathbb{F}_p}\mathcal{X}(y)\mathcal{X}(1-y)\)。

​ 定义数域\(K= Q(\zeta)\)，即\(K\)为\(Q\)基于\(\zeta\)的代数扩张，可知\(K\)的自同构群的一组基为\(\sigma_i(\zeta) = \zeta^i, i = 1, 2, \dots, n-1\)。设\(T, N\)分别为\(K/Q\)扩张上的trace和norm，通过这篇paper(sci-hub上并没有资源，很尴尬)，结合trace和norm的定义，得到曲线\(E\)上点的个数\(M = p+1+T_{K/Q}(J(\mathcal{X}, \mathcal{X}))\)，通过zeta函数的转换得到曲线Jacobian上点的个数\(N = N_{K/Q}(J(\mathcal{X}, \mathcal{X})+1)\)，论文中该给出了求原曲线经过\(i, j\)扭曲后的阶。

​ 可以看到现在的关键即是求解\(J(\mathcal{X}, \mathcal{X})\)，求解了\(J(\mathcal{X}, \mathcal{X})\)即可计算出\(N\) 即得到原曲线的阶。考虑由\(p, \alpha^{\frac{p-1}{n}}\) 生成的\(\mathbb{Z}[\zeta]\) 上的素理想\(P\)，\(n<23\)时(包含了论文中的情况)，\(P\)是principal的，即可以只由一个生成元生成，假设得到了\(P\)的一个生成元为\(\beta\)，则可以计算出\(J(\mathcal{X}, \mathcal{X})\)。

​ 由\(\hat J = \prod_{i=1}^g \sigma_i^{-1}(\beta)\) 生成的理想与由\(J(\mathcal{X}, \mathcal{X})\)生成的理想是相同的，所以可以通过 $ $ 求得 \(J(\mathcal{X}, \mathcal{X})\)（为同一个Ideal的生成元，故其之间仅差了一个单位元1的根）。存在 \(r\in \{-1, 1\}, s \in Z\) 使得 \(r\zeta^s\hat J = J(\mathcal{X}, \mathcal{X})\)，将\(\hat J\) 表示为\(\zeta\) 多项式的形式 \(\sum_{j=0}^{n-2}a_j\zeta^j\)。对于\(J(\mathcal{X}, \mathcal{X})\)，有\(J(\mathcal{X}, \mathcal{X}) \equiv -1 \mod (\zeta -1)^2\)，将\(r\zeta^s\hat J = J(\mathcal{X}, \mathcal{X})\)代入，则有

\(\quad \quad r\zeta^s\hat J \equiv r(1+s(\zeta-1))\sum_{j=0}^{n-2}a_j(1+j(\zeta-1))\)

\(\quad \quad \equiv r(\sum_{j=0}^{n-2}a_j + (s\sum_{j=0}^{n-2}a_j+\sum_{j=0}^{n-2}ja_j)(\zeta-1))\mod (\zeta-1)^2\)

​ 容易发现当\(r = -\sum_{j=0}^{n-2}a_j, s= r \sum_{j=0}^{n-2}ja_j\)时，等式成立（代入即可），由此只要知道了\(\beta\)，可以计算出\(\hat J\)和\(r, s\)，也就可以得到\(J(\mathcal{X}, \mathcal{X})\)。

​ 接下来讨论如何求出Ideal \(P\)的生成元\(\beta\)。\(Z(\zeta)\)中的元素可以看作生成元\(\zeta\)的多项式\(x_{2g}\zeta^{2g} + \dots + x_1\zeta\)，将\(x_{2j-1}, x_{2j}\)映射为复数\(x_{2j-1} + ix_{2j}\)，则可以通过此法将\(P\)嵌入到\(\mathbb{C}^g\)上，再按照顺序展开为\((x_1, x_2, \dots, x_{2g})\)，由此映射到了\(R^{2g}\)上，构造了一个行列式维\(2^{-g}p\sqrt{D}\)的格\(L_P\)，其为\(Z(\zeta)\)映射到\(R^{2g}\)上格的子格。通过对格\(L_P\)进行LLL，根据规约后的结果\(\vec x = (x_1, x_2, \dots, x_{2g})\)，可以计算其在\(K/Q\)上的norm，如果为\(p\)，则为\(P\)的一个生成元，即可计算出\(\beta\)，如果为\(c\cdot p\)，则在Ideal \(P\) 上爆破遍历寻找norm为\(c\)的\(y\)，寻找生成元\(x/y\)。也是因为这个原因，当\(n\)较大时，\(c\) 较大，需要爆破的范围过大，所以无法使用此方法来寻找生成元，不过论文中提到规约后的向量其实相当一部分的norm均为\(p\)，实际的表现比想象中要好得多。可以通过格基规约来求得生成元的关键在于，Ideal \(P\) 上的元素映射到 \(R^{2g}\) 上后，向量的范数也就与本身元素的norm等价，而Ideal \(P\) 上的范数均为 \(cp, c\in Z\)，故只需要寻找范数最短的向量也就找到了norm最小的元素即生成元。

​ 梳理一下算法流程：

  1. 使用Ideal \(P\)的基\(\{p, \zeta-a, \zeta^2-a_2, \dots, \zeta^{n-2}-a_{n-2}\}\) 构造格\(L_P\)
  2. 对格\(L_P\)进行规约，检测短向量的norm是否为\(p\)，为\(p\)则其对应\(Z(\zeta)\)中元素即为Ideal \(P\)的一个生成元\(\beta\)，若短向量norm为\(c\cdot p\)，则爆破所有norm为\(c\)的\(y\)，通过\(x/y\)寻找生成元。
  3. 得到生成元后，可以计算出\(\hat J\)，将\(\hat J\)展开为\(\zeta\)的多项式，得到\(\hat J = \sum_{j=0}^{n-2}a_j\zeta^j\)，求出\(r \equiv -\sum_{j=0}^{n-2}a_j \mod n, s \equiv r \sum_{j=0}^{n-2}ja_j \mod n\)，即得到$ J(, ) = r^sJ $
  4. 将\(J(\mathcal{X}, \mathcal{X})\)代入原曲线的Jacobian上点数量的公式\(N = N_{K/Q}(J(\mathcal{X}, \mathcal{X})+1)\)，即得到了原曲线的阶。

​ 构造格\(L_P\)规约时注意，我一开始使用sage的CC作为复数域，但是其精度过低，之后使用欧拉公式转换到RR上进行运算，再使用一个放大因子即可。然而在最后使用\(\hat J\) 求 \(J(\mathcal{X}, \mathcal{X})\) 的过程中，我发现 \(s \equiv r \sum_{j=0}^{n-2}ja_j \mod n\) 并不能求出正确的\(J(\mathcal{X}, \mathcal{X})\)，反而在使用\(s \equiv r \sum_{j=0}^{n-2}(n-2-j)a_j \mod n\) 时有大概率正确。考虑与论文中不一致的地方在于，题中曲线为\(y^2 + y = 114514x^7\)，\(x^n\)的系数并不为1，所以需要在最后需要利用扭曲后曲线的求阶公式\(N_{i, j} = N_{K/Q}(J(\mathcal{X}, \mathcal{X}) + (-1)^i\zeta^{j})\)，扭曲公式为\(v^2+v+(1-\eta^i)/4 = \eta^i \alpha^j u^n\)，此处即为\(i = 0, \alpha^j \equiv 114514 \mod p\)，这里只需要求\(j \mod n\)，

求私钥

​ 求了阶之后，HCDSA的部分就比较显然了，构造HNP求解即可。我们知道\(S_i = r + H\cdot (d\otimes f_i)\)，故可以使用 \(\frac{S_m-S_n}{S_i-Sj} \equiv \frac{d\otimes f_m -d\otimes f_n}{d\otimes f_i - d \otimes f_j} \mod q\) 消掉未知的\(r, H\)。由于\(d\otimes f_i = d + f_i - 2 \cdot (d \wedge f_i)\)，由此可以写为\(\frac{f_m-f_n+2(d\wedge f_n - d\wedge f_m)}{d\otimes f_i - d \otimes f_j}\)。考虑二进制展开\(d = \sum_k 2^kd_k, f_m = \sum_k2^kf_{m,k}, f_n = \sum_k2^kf_{n,k}\)，则\(d \wedge f_m - d\wedge f_n = \sum_k 2^k d_k\cdot (f_{n,k}-f_{m,k})\)，故对于不同的\(i, j, m, n\)，均有关系式：

\((S_m-S_n)(f_i-f_j+\sum_k 2^{k+1} d_k\cdot (f_{j,k}-f_{i,k})) \equiv (S_i-S_j)(f_m-f_n+\sum_k 2^{k+1} d_k\cdot (f_{n,k}-f_{m,k})) \mod q\)

\(\Rightarrow \sum_k2^k (S_m-S_n) (f_{i,k}-f_{j, k}+2 d_k(f_{j, k}-f_{i, k})) \equiv \sum_k2^k (S_i-S_j) (f_{m,k}-f_{n, k}+2 d_k(f_{n, k}-f_{m, k})) \mod q\)

\(\Rightarrow (S_m-S_n) (f_{i}-f_{j}) - (S_i-S_j) (f_{m}-f_{n}) + \sum_k 2^{k+1}d_k ((S_m-S_n)(f_{j, k}-f_{i, k})- (S_i-S_j)(f_{n, k}-f_{m, k}))\equiv 0 \mod q\)

\(\Rightarrow \sum_k 2^{k}(2d_k-1) ((S_m-S_n)(f_{j, k}-f_{i, k})- (S_i-S_j)(f_{n, k}-f_{m, k}))\equiv 0 \mod q\)

设\(\beta_{i, j, m, n, k} = 2^k((S_m-S_n)(f_{j, k}-f_{i, k})- (S_i-S_j)(f_{n, k}-f_{m, k}))\)， 则有\(\sum_k (2d_k-1)\beta_{i, j, m, n, k} \equiv 0 \mod q\)

则可以构建格： \[ M = \begin{bmatrix} K\cdot \beta_{0, 1, 2, 3, 0} & \cdots & K\cdot \beta_{7, 8, 9, 10, 0}& 1 & \cdots & 0\\ \vdots & \ddots & \vdots & \vdots & \ddots & \vdots \\ K\cdot \beta_{0, 1, 2, 3, 767}& \cdots & K\cdot \beta_{7, 8, 9, 10, 767}& 0 & \cdots & 1 \\ K\cdot q & \cdots & 0 & 0 & \cdots & 0 \\ \vdots & \ddots & \vdots & \vdots & \ddots & \vdots\\ 0 & \cdots & K\cdot q & 0 & \cdots & 0 \ \end{bmatrix} \] 其中\(K\)为放大系数，这里取1即可，放缩优化习惯了orz。格上存在短向量\((0, \dots, 0, 2d_0-1, 2d_1-1, \dots, 2d_{767}-1)\)，规约矩阵\(M\)，这里遍历\(i, j, m, n\) 可以拿到几千组数据。可以选其中几十组构造格，但由于\(\sum_i f_i\)有2个零位，下标为250，767，所选数据至少也会产生两个0位，所以d的这两位是无法获得的，并且格上会存在形如\((0, \dots, 1, \dots, 0)\) 的短向量，作为优化可以把这两行剔除，剔除后整理起来麻烦些，exp里为不剔除的情况。还有一种情况为会存在固定的\(k_0, k_1\)使得对于所有\(i, j\) 有 \(f_{j, k_0} - f_{i, k_0} + f_{j, k_1} - f_{i, k_1} = 0\)，当\(k_1-k_0=1\) 时，格上会存在形如\((0, \dots, 2,1, \dots, 0)\) 的短向量，题中\(k_0 = 497, k_1=498\)时有此情况，此时如果不剔除有可能产生错误比特对 \((-1, 0)-(1, 1)\)。

具体选取 \(\beta\) 时其实可以少选些，只要产生的零位最少即可，大概需要40-50组数据即可恢复，有时候可能出现预期向量不在规约后格的情况，多跑几次(shuffle)即可，完整exp如下：

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
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
from sage.all import *
from Crypto.Util.number import *
from hashlib import sha256, sha512
import time


def H(_m):
    if not isinstance(_m, bytes):
        _m = str(_m).encode()
    return bytes_to_long(sha256(_m).digest() + sha512(_m).digest())


def lattice_reduction_get_order(_p, _n, constant, x_n_coefficient):
    assert isPrime(_p) and (_p - 1) % _n == 0
    # generate K
    PRQ = PolynomialRing(QQ, 'y')
    y = PRQ.gens()[0]
    K = QQ.extension((y ** _n - 1).factor()[-1][0], 'z')
    z = K.gens()[0]
    g = (_n - 1) // 2

    # generate alpha, a
    while True:
        alpha = randrange(1, _p - 1)
        if pow(alpha, _n, _p) != 1 and pow(alpha, (_p - 1) // _n, _p) != 1:
            break
    a = int(pow(alpha, (_p - 1) // _n, _p))

    # build P lattice
    MUL = 10 ** 500
    LP = Matrix(ZZ, _n - 1, _n - 1)
    LP[0] = vector(ZZ, [0 if i % 2 else _p * MUL for i in range(_n - 1)])
    for i in range(1, _n - 1):
        tmp_vec = []
        for j in range(1, g + 1):
            tmp_vec.append(RR((2 / _n * i * j * pi).cos() * MUL).round() - int(pow(a, i, _p)) * MUL)
            tmp_vec.append(RR((2 / _n * i * j * pi).sin() * MUL).round())
        LP[i] = vector(ZZ, tmp_vec)

    # reduce the lattice
    LPL = LP.LLL()

    # get beta coefficients of Z-basis, then get beta
    coefficients = LP.transpose().solve_right(LPL[0])
    tmp = _p * coefficients[0]
    for i in range(1, _n - 1):
        tmp -= int(pow(a, i, _p)) * coefficients[i]
    beta = K([tmp] + coefficients[1:].list())
    assert beta.norm() % _p == 0
    if beta.norm() != _p:
        c = beta.norm() // _p
        # brute force for a y has norm c
        for i in range(3 ** g):
            y = [0] * g
            for j in range(g):
                y.append(i % 3)
                i //= 3
            if K(y).norm() == c:
                beta /= K(y)
                break
    assert beta.norm() == _p

    # calculate J(X, X)
    sigma = K.automorphisms()
    J_hat = 1
    for i in range((_n - 1) // 2):
        J_hat *= sigma[i].preimage(beta)
    coefficients = J_hat.polynomial().coefficients()
    r = -int(sum(coefficients)) % _n
    if r != 1:
        assert r == _n - 1
        r = -1
    s = int(r * sum([coefficients[j] * j for j in range(_n - 1)])) % _n
    J_X_X = r * z ** s * J_hat

    # calculate twist i, j and get order
    if constant != 0:
        i = 1
    else:
        i = 0
    for j in range(_n):
        if pow(a, j, _p) == pow((-1) ** i * x_n_coefficient, (_p - 1) // _n, _p):
            return int((J_X_X + (-1) ** i * z ** j).norm())


def lattice_reduction_HCDSA_HNP_attack(_Sf, _q):
    bits = _q.bit_length()
    num = len(_Sf)
    S = [sfi[1] for sfi in _Sf]
    f = [sfi[0] for sfi in _Sf]
    fb = [list(map(int, bin(f[i])[2:].rjust(bits, '0')))[::-1] for i in range(num)]
    check_zero = sum(vector(ZZ, fbi) for fbi in fb).list()
    _zero_index = []
    for i in range(bits):
        if check_zero[i] == 0:
            _zero_index.append(i)
    for k0 in range(bits-3):
        for k1 in range(k0+1, k0+3):
            is_zero_index = True
            for i in range(num):
                if not is_zero_index:
                    break
                for j in range(i+1, num):
                    if fb[j][k0] - fb[i][k0] + fb[j][k1] - fb[i][k1] != 0:
                        is_zero_index = False
                        break
            if is_zero_index:
                _zero_index.extend([k0, k1])
    _zero_index = list(set(_zero_index))

    # calculate beta
    beta = []
    for i in range(num - 3):
        for j in range(i + 1, num - 2):
            for m in range(j + 1, num - 1):
                for n in range(m + 1, num):
                    beta.append([])
                    for k in range(bits):
                        beta[-1].append(2 ** k *
                                        ((S[m] - S[n]) * (fb[j][k] - fb[i][k]) -
                                         (S[i] - S[j]) * (fb[n][k] - fb[m][k])))

    # calculate the least number of beta
    use_num = 45
    got = False
    while not got:
        use_num += 1
        for i in range(20):
            shuffle(beta)
            if sum(vector(ZZ, beta_i) for beta_i in beta[:use_num]).list().count(0) == len(_zero_index):
                got = True
                break

    # build lattice
    # K = 2**bits
    K = 1
    M = Matrix(ZZ, bits + use_num, bits + use_num)
    for i in range(use_num):
        M[:, i] = K * vector(ZZ, beta[i] + [0] * i + [q] + [0] * (use_num - i - 1))
    M[:bits, use_num:use_num + bits] = identity_matrix(ZZ, bits)
    # M.delete_rows(_zero_index)

    # reduce lattice and get d
    start = time.time()
    ML = M.LLL()
    end = time.time()
    print('use %d seconds' % int(end - start))
    for l in ML:
        if l[:use_num].is_zero():
            tmp_db = l[use_num:use_num + bits][::-1]
            if set(tmp_db) != {-1, 0, 1}:
                continue
            _db0 = ['1' if dbi == 1 else '0' for dbi in tmp_db]
            _db1 = ['0' if dbi == 1 else '1' for dbi in tmp_db]
            for i in range(bits):
                if tmp_db[i] == 0 and bits-1-i not in _zero_index:
                    _zero_index.append(bits-1-i)
            for i in _zero_index:
                _db0[bits-1-i] = '0'
                _db1[bits-1-i] = '0'
            return [int(''.join(_db0), 2), int(''.join(_db1), 2)], _zero_index


# get order
p = 109912987917332562152579558699996713895127464171381146718513726836854765335393
PRp = PolynomialRing(GF(p), 'x')
x = PRp.gens()[0]
C2 = HyperellipticCurve(114514 * x ** 7, 1)
J2 = C2.jacobian()
G = [x - 1195618977299087187547891775318892467576479852329251168026008209153499647043, 61722058527536369025641486069921773295057337558818637400938099789108164705978]
G = J2(*G)
q = lattice_reduction_get_order(p, 7, 0, 114514)
assert q * G == J2(0)

# get most bits of d
Sf = 
all_d, zero_index = lattice_reduction_HCDSA_HNP_attack(Sf, q)

# # check correctness and get flag
msg = 'Welcome to HUFU CTF 2022 Quals!'
for possible_error in range(2**len(zero_index)):
    error_bit_array = list(map(int, bin(possible_error)[2:].rjust(len(zero_index), '0')))
    error = 0
    for index in range(len(zero_index)):
        error += 2**zero_index[index] * error_bit_array[index]
    for possible_d in all_d:
        d = possible_d + error
        R = H(str(d) + msg) * G
        if (H(str(d) + msg) + H(str(R) + str(d*G) + msg) * (d ^^ Sf[0][0])) % q == Sf[0][1]:
            print('d =', d)
            print('Flag of "HCDSA": HFCTF{{{}}}'.format(sha256(str(d).encode()).hexdigest()))
            break

最终得到d及flag：

HCDH

歧路

没仔细审题导致读了一篇别的paper，最后发现并不符合，但是觉得还挺有意思，整体流程很适合出题，之后应该会出到别的比赛里，心路历程就不在这里写了。

归途

重新思考题目，这道题官方wp没给具体思路，观察后发现是个很标准的基于超椭的DH，所以攻击点肯定在于这个曲线本身，算了下，发现这个曲线的embedding inidex很小只有3，考虑可以转换到有限域上的DLP问题。类似于MOV attack，只需要将曲线映射到\(2^{113*3}\) 上，然后通过pairing映射到有限域\(GF(2^{113*3})\)上，直接在小因子上求DLP，通过CRT构造一个模比\(2^{113}\)大的即可解出私钥，注意这里是求\(GF(2^{339})\)上的DLP，需要特定的index calculus算法，使用pari的fflog计算。

还需要注意的是，\(2^{113}\)上曲线的lift_x函数貌似有bug，算y坐标需要在\(2^{113}\)上解一元二次方程，Jacobian上同样没提供lift_x函数，需要在\(\mathbb{F}[2^{113}]\)上解\(v^2 + v \equiv x^5 + x^3 + f_1x + f_0 \mod u\)，同样可以化简为\(2^{113}\)上的二次方程。由于Sagemath未实现超椭上的pairing，一开始想自己实现个Miller算法，但发现挺麻烦，调了调bug懒得调了，直接使用MAGMA求解即可，注意MAGMA的\(GF(2^{339})\) 默认的模与sage的不同，sage求DLP时需要调整一下，完整exp如下：

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
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
from sage.all import *
from Crypto.Util.number import *
from hashlib import sha256, md5
from tqdm import tqdm

GF2_ext = GF(2 ** 113, 'z')
z = GF2_ext.gens()[0]
PR2_ext = PolynomialRing(GF2_ext, 'x')
x = PR2_ext.gens()[0]
C1 = HyperellipticCurve(x ^ 5 + x ^ 3 + x, 1)
J1 = C1.jacobian()


def calculate_order(p, n, f0, f1, genus):
    # h(x) = 1
    # f(x) = y^5 + y^3 + f1 * y + f0

    # Calc M_i
    M = []
    for i in range(genus):
        GFpi = GF(p ** (i + 1))
        Ppn = PolynomialRing(GFpi, 'w')
        w = Ppn.gens()[0]
        if i == 0:
            _f0 = f0.integer_representation() % p
            _f1 = f1.integer_representation() % p
        else:
            _f0 = GFpi.fetch_int(f0.integer_representation() % p ** (i + 1))
            _f1 = GFpi.fetch_int(f1.integer_representation() % p ** (i + 1))
        C = HyperellipticCurve(w ^ 5 + w ^ 3 + _f1 * w + _f0, 1)
        M.append(C.cardinality())
    print("Calc M done")

    # Calc a_i
    a = [0]
    for i in range(genus):
        a.append((M[i] - 1 - p ^ (i + 1) + a[i] ^ (i + 1)) / (i + 1))
    a = [1] + a[1:]
    print("Calc a done")

    # Calc gamma_i
    var('X')
    function = -genus * p
    for i in range(genus + 1):
        function += a[i] * X ^ (genus - i)
    gamma = list(map(lambda y: y.rhs(), solve([function == 0], X)))
    print("Calc gamma done")

    # Calc alpha_i
    alpha = []
    for i in range(genus):
        function = X ^ 2 - gamma[i] * X + p
        alpha.append(list(map(lambda y: y.rhs(), solve([function == 0], X))))
    print("Calc alpha done")

    # Calc #J
    order = 1
    for i in range(genus):
        order *= abs(1 - alpha[i][0] ^ n) ^ 2
    return int(order)


def solve_quadratic_equation_GF2n(c, m):
    # solve v**2 + v = c on 2**m
    v = 0
    for i in range((m + 1) // 2):
        v += c ** (2 ** (2 * i))
    return GF2_ext(v)


def get_v_from_u(u, f0, f1):
    v2_v = (x ** 5 + x ** 3 + f1 * x + f0) % u
    _e = u[1]
    _f = u[0]
    _c = v2_v[1]
    _d = v2_v[0]
    # v = ax+b, a-a^2e = c, b^2+b-a^2f = d
    _a = solve_quadratic_equation_GF2n(-_c * _e, 113) / -_e
    _b = solve_quadratic_equation_GF2n(_d + _a ** 2 * _f, 113)
    return PR2_ext([_b, _a])


def parse(point, shift, bits):
    pad = bits - 2 * shift - 1
    assert pad > 0
    _x, _y = point
    n = (_x[0].integer_representation() << shift) + _x[1].integer_representation()
    n = (n << pad) + randint(1, 2 ^ pad)
    n = (n << 1) + int(_y[0].integer_representation() % 2)
    return n


def re_parse(n, shift, bits):
    pad = bits - 2 * shift - 1
    symbol = n % 2
    n >>= 1
    n >>= pad
    x1 = n % 2 ** shift
    x0 = n >> shift
    u = PR2_ext([GF2_ext.fetch_int(x0), GF2_ext.fetch_int(x1), 1])
    assert u.is_irreducible()
    v = get_v_from_u(u, 0, 1)
    if v[0].integer_representation() % 2 != symbol:
        v = PR2_ext([v[0] + 1, v[1]])
    return J1([u, PR2_ext(v)])


def GF2n_Ph(_m, _g, _h, _order):
    _a = []
    for mi in _m:
        tmp_order = _order // mi
        _a.append(pari.fflog(_h ** tmp_order, _g ** tmp_order, mi))
        print(mi, _a)
    return crt(_a, _m)


Alice_pk = [x ^ 2 + (z ^ 111 + z ^ 107 + z ^ 106 + z ^ 104 + z ^ 102 + z ^ 100 + z ^ 97 + z ^ 95 + z ^ 94 + z ^ 93 + z ^ 92 + z ^ 89 + z ^ 88 + z ^ 84 + z ^ 83 + z ^ 77 + z ^ 76 + z ^ 72 + z ^ 69 + z ^ 68 + z ^ 66 + z ^ 62 + z ^ 61 + z ^ 60 + z ^ 59 + z ^ 58 + z ^ 56 + z ^ 55 + z ^ 54 + z ^ 53 + z ^ 50 + z ^ 47 + z ^ 46 + z ^ 45 + z ^ 43 + z ^ 42 + z ^ 39 + z ^ 38 + z ^ 35 + z ^ 34 + z ^ 32 + z ^ 31 + z ^ 25 + z ^ 24 + z ^ 23 + z ^ 21 + z ^ 19 + z ^ 18 + z ^ 17 + z ^ 16 + z ^ 15 + z ^ 13 + z ^ 12 + z ^ 11 + z ^ 7 + z ^ 5 + z ^ 3) * x + z ^ 112 + z ^ 111 + z ^ 110 + z ^ 107 + z ^ 106 + z ^ 102 + z ^ 100 + z ^ 97 + z ^ 96 + z ^ 95 + z ^ 94 + z ^ 90 + z ^ 86 + z ^ 85 + z ^ 84 + z ^ 82 + z ^ 79 + z ^ 78 + z ^ 77 + z ^ 76 + z ^ 75 + z ^ 74 + z ^ 71 + z ^ 69 + z ^ 64 + z ^ 62 + z ^ 61 + z ^ 58 + z ^ 56 + z ^ 55 + z ^ 52 + z ^ 51 + z ^ 49 + z ^ 47 + z ^ 44 + z ^ 43 + z ^ 41 + z ^ 33 + z ^ 29 + z ^ 28 + z ^ 27 + z ^ 26 + z ^ 25 + z ^ 24 + z ^ 23 + z ^ 21 + z ^ 18 + z ^ 17 + z ^ 15 + z ^ 12 + z ^ 11 + z ^ 10 + z ^ 9 + z ^ 7 + z ^ 4 + z ^ 3 + z + 1,
            (z ^ 112 + z ^ 110 + z ^ 109 + z ^ 103 + z ^ 102 + z ^ 101 + z ^ 100 + z ^ 99 + z ^ 97 + z ^ 96 + z ^ 95 + z ^ 94 + z ^ 93 + z ^ 92 + z ^ 89 + z ^ 88 + z ^ 86 + z ^ 84 + z ^ 82 + z ^ 80 + z ^ 79 + z ^ 75 + z ^ 74 + z ^ 73 + z ^ 68 + z ^ 64 + z ^ 60 + z ^ 59 + z ^ 55 + z ^ 50 + z ^ 48 + z ^ 46 + z ^ 43 + z ^ 42 + z ^ 39 + z ^ 38 + z ^ 37 + z ^ 36 + z ^ 35 + z ^ 34 + z ^ 28 + z ^ 25 + z ^ 23 + z ^ 22 + z ^ 20 + z ^ 17 + z ^ 15 + z ^ 13 + z ^ 10 + z ^ 3 + z) * x + z ^ 111 + z ^ 110 + z ^ 109 + z ^ 106 + z ^ 105 + z ^ 102 + z ^ 93 + z ^ 90 + z ^ 89 + z ^ 88 + z ^ 86 + z ^ 84 + z ^ 83 + z ^ 81 + z ^ 80 + z ^ 77 + z ^ 76 + z ^ 74 + z ^ 73 + z ^ 72 + z ^ 70 + z ^ 68 + z ^ 66 + z ^ 64 + z ^ 61 + z ^ 60 + z ^ 59 + z ^ 58 + z ^ 57 + z ^ 50 + z ^ 49 + z ^ 47 + z ^ 46 + z ^ 45 + z ^ 44 + z ^ 43 + z ^ 42 + z ^ 39 + z ^ 38 + z ^ 36 + z ^ 33 + z ^ 32 + z ^ 31 + z ^ 29 + z ^ 28 + z ^ 26 + z ^ 25 + z ^ 24 + z ^ 17 + z ^ 15 + z ^ 13 + z ^ 10 + z ^ 9 + z ^ 8 + z ^ 5 + z ^ 4 + z ^ 3 + z ^ 2 + z]
Alice_pk = J1(*Alice_pk)
q = ZZ(calculate_order(2, 113, z - z, z - z + 1, 2))
assert q * Alice_pk == J1(0)
# q = 7 * 10113049 * 320021624768405574452943847 * 4760137992283599860814226997712217
for k in range(2, 12):
    if (2 ** (113 * k) - 1) % q == 0:
        break

d = 402045457200978546368835668950711933106933538070552823738108748139869490234501184225175369065216758931906216168485716846341035514282333269655203225744353550872749071863611510874494552527466229826353578228691696842853401202429961474
Shared_key = re_parse(d, 113, 767)
check_d = parse(Shared_key, 113, 767)
assert check_d >> 541 == d >> 541 and check_d % 2 == d % 2

# MAGMA
'''
K<z>:=GF(2^113);
R<x>:=PolynomialRing(K);
f:=x^5+x^3+x;
C:=HyperellipticCurve(f, 1);
J:=Jacobian(C);
q:=107839786668602559178668060348078533079142294759817546985407099437057;
k:=3;

ug:=x^2 + (z^111 + z^107 + z^106 + z^104 + z^102 + z^100 + z^97 + z^95 + z^94 + z^93 + z^92 + z^89 + z^88 + z^84 + z^83 + z^77 + z^76 + z^72 + z^69 + z^68 + z^66 + z^62 + z^61 + z^60 + z^59 + z^58 + z^56 + z^55 + z^54 + z^53 + z^50 + z^47 + z^46 + z^45 + z^43 + z^42 + z^39 + z^38 + z^35 + z^34 + z^32 + z^31 + z^25 + z^24 + z^23 + z^21 + z^19 + z^18 + z^17 + z^16 + z^15 + z^13 + z^12 + z^11 + z^7 + z^5 + z^3)*x + z^112 + z^111 + z^110 + z^107 + z^106 + z^102 + z^100 + z^97 + z^96 + z^95 + z^94 + z^90 + z^86 + z^85 + z^84 + z^82 + z^79 + z^78 + z^77 + z^76 + z^75 + z^74 + z^71 + z^69 + z^64 + z^62 + z^61 + z^58 + z^56 + z^55 + z^52 + z^51 + z^49 + z^47 + z^44 + z^43 + z^41 + z^33 + z^29 + z^28 + z^27 + z^26 + z^25 + z^24 + z^23 + z^21 + z^18 + z^17 + z^15 + z^12 + z^11 + z^10 + z^9 + z^7 + z^4 + z^3 + z + 1;
vg:=(z^112 + z^110 + z^109 + z^103 + z^102 + z^101 + z^100 + z^99 + z^97 + z^96 + z^95 + z^94 + z^93 + z^92 + z^89 + z^88 + z^86 + z^84 + z^82 + z^80 + z^79 + z^75 + z^74 + z^73 + z^68 + z^64 + z^60 + z^59 + z^55 + z^50 + z^48 + z^46 + z^43 + z^42 + z^39 + z^38 + z^37 + z^36 + z^35 + z^34 + z^28 + z^25 + z^23 + z^22 + z^20 + z^17 + z^15 + z^13 + z^10 + z^3 + z)*x + z^111 + z^110 + z^109 + z^106 + z^105 + z^102 + z^93 + z^90 + z^89 + z^88 + z^86 + z^84 + z^83 + z^81 + z^80 + z^77 + z^76 + z^74 + z^73 + z^72 + z^70 + z^68 + z^66 + z^64 + z^61 + z^60 + z^59 + z^58 + z^57 + z^50 + z^49 + z^47 + z^46 + z^45 + z^44 + z^43 + z^42 + z^39 + z^38 + z^36 + z^33 + z^32 + z^31 + z^29 + z^28 + z^26 + z^25 + z^24 + z^17 + z^15 + z^13 + z^10 + z^9 + z^8 + z^5 + z^4 + z^3 + z^2 + z;
uh:=x^2 + (z^112 + z^110 + z^109 + z^106 + z^105 + z^104 + z^102 + z^97 + z^93 + z^92 + z^90 + z^88 + z^87 + z^85 + z^84 + z^81 + z^80 + z^79 + z^78 + z^76 + z^75 + z^69 + z^67 + z^63 + z^62 + z^61 + z^60 + z^58 + z^55 + z^54 + z^53 + z^52 + z^50 + z^49 + z^48 + z^47 + z^46 + z^44 + z^43 + z^42 + z^38 + z^37 + z^36 + z^34 + z^30 + z^28 + z^27 + z^26 + z^25 + z^23 + z^22 + z^21 + z^20 + z^17 + z^15 + z^7 + z^6 + z^5 + z^2 + z + 1)*x + z^112 + z^107 + z^104 + z^101 + z^99 + z^98 + z^96 + z^95 + z^93 + z^92 + z^90 + z^89 + z^88 + z^87 + z^85 + z^80 + z^79 + z^76 + z^75 + z^74 + z^73 + z^72 + z^71 + z^70 + z^67 + z^66 + z^64 + z^63 + z^59 + z^58 + z^56 + z^55 + z^52 + z^49 + z^48 + z^42 + z^41 + z^40 + z^39 + z^37 + z^36 + z^35 + z^34 + z^33 + z^32 + z^31 + z^30 + z^28 + z^26 + z^25 + z^23 + z^22 + z^21 + z^19 + z^17 + z^15 + z^14 + z^11 + z^4 + z^3 + z;
vh:=(z^112 + z^111 + z^110 + z^108 + z^105 + z^103 + z^100 + z^99 + z^98 + z^97 + z^96 + z^95 + z^94 + z^92 + z^88 + z^87 + z^86 + z^85 + z^83 + z^81 + z^80 + z^79 + z^76 + z^74 + z^72 + z^71 + z^70 + z^69 + z^66 + z^65 + z^63 + z^57 + z^56 + z^53 + z^52 + z^51 + z^41 + z^35 + z^34 + z^32 + z^30 + z^29 + z^27 + z^25 + z^22 + z^14 + z^12 + z^11 + z^10 + z^8 + z^5 + z^3 + z)*x + z^112 + z^110 + z^108 + z^107 + z^103 + z^98 + z^95 + z^92 + z^90 + z^88 + z^83 + z^82 + z^81 + z^76 + z^75 + z^72 + z^71 + z^70 + z^68 + z^67 + z^65 + z^64 + z^58 + z^56 + z^52 + z^51 + z^49 + z^45 + z^44 + z^43 + z^42 + z^41 + z^35 + z^28 + z^27 + z^25 + z^24 + z^23 + z^21 + z^17 + z^16 + z^14 + z^13 + z^10 + z^9 + z^7 + z^2 + z;

g:=J![ug,vg];
h:=J![uh,vh];

Jext:=BaseExtend(J,3);
Kk<zk>:=GF(2^339);
g3:=Jext!g;
h3:=Jext!h;
q3_q:= 107839786668602559178668060348078501925361143550851775802429124116481;
R:= q3_q * Random(Jext);

ge_pairing:=WeilPairing(g3, R, q);
he_pairing:=WeilPairing(h3, R, q);
print ge_pairing;
print he_pairing;
'''

PR = PolynomialRing(GF(2), 'xw')
xw = PR.gens()[0]
GF2k_ext = GF(2 ** (113 * k), 'zk', modulus=xw ^ 339 + xw ^ 7 + xw ^ 5 + xw ^ 3 + xw ^ 2 + xw + 1)
zk = GF2k_ext.gens()[0]
PR2k_ext = PolynomialRing(GF2k_ext, 'xk')
xk = PR2k_ext.gens()[0]
C3 = HyperellipticCurve(xk ^ 5 + xk ^ 3 + xk, 1)
J3 = C3.jacobian()


g3_pairing = zk ^ 337 + zk ^ 336 + zk ^ 334 + zk ^ 333 + zk ^ 332 + zk ^ 331 + zk ^ 330 + zk ^ 323 + zk ^ 319 + zk ^ 318 + zk ^ 313 + zk ^ 312 + zk ^ 310 + zk ^ 306 + zk ^ 304 + zk ^ 303 + zk ^ 299 + zk ^ 297 + zk ^ 296 + zk ^ 295 + zk ^ 293 + zk ^ 289 + zk ^ 287 + zk ^ 286 + zk ^ 285 + zk ^ 284 + zk ^ 282 + zk ^ 279 + zk ^ 278 + zk ^ 277 + zk ^ 275 + zk ^ 273 + zk ^ 269 + zk ^ 267 + zk ^ 265 + zk ^ 264 + zk ^ 261 + zk ^ 260 + zk ^ 258 + zk ^ 257 + zk ^ 255 + zk ^ 254 + zk ^ 251 + zk ^ 249 + zk ^ 243 + zk ^ 242 + zk ^ 241 + zk ^ 239 + zk ^ 237 + zk ^ 235 + zk ^ 232 + zk ^ 231 + zk ^ 230 + zk ^ 229 + zk ^ 228 + zk ^ 227 + zk ^ 223 + zk ^ 221 + zk ^ 220 + zk ^ 219 + zk ^ 214 + zk ^ 212 + zk ^ 209 + zk ^ 208 + zk ^ 206 + zk ^ 204 + zk ^ 200 + zk ^ 196 + zk ^ 193 + zk ^ 192 + zk ^ 190 + zk ^ 189 + zk ^ 187 + zk ^ 185 + zk ^ 184 + zk ^ 183 + zk ^ 180 + zk ^ 178 + zk ^ 176 + zk ^ 175 + zk ^ 173 + zk ^ 171 + zk ^ 169 + zk ^ 168 + zk ^ 166 + zk ^ 165 + zk ^ 162 + zk ^ 158 + zk ^ 157 + zk ^ 156 + zk ^ 155 + zk ^ 154 + zk ^ 151 + zk ^ 149 + zk ^ 148 + zk ^ 142 + zk ^ 140 + zk ^ 139 + zk ^ 137 + zk ^ 135 + zk ^ 134 + zk ^ 133 + zk ^ 131 + zk ^ 126 + zk ^ 125 + zk ^ 123 + zk ^ 119 + zk ^ 117 + zk ^ 115 + zk ^ 113 + zk ^ 112 + zk ^ 111 + zk ^ 108 + zk ^ 107 + zk ^ 106 + zk ^ 105 + zk ^ 101 + zk ^ 100 + zk ^ 99 + zk ^ 96 + zk ^ 93 + zk ^ 92 + zk ^ 91 + zk ^ 89 + zk ^ 87 + zk ^ 84 + zk ^ 82 + zk ^ 81 + zk ^ 80 + zk ^ 79 + zk ^ 78 + zk ^ 77 + zk ^ 75 + zk ^ 72 + zk ^ 70 + zk ^ 69 + zk ^ 68 + zk ^ 66 + zk ^ 65 + zk ^ 63 + zk ^ 60 + zk ^ 58 + zk ^ 57 + zk ^ 55 + zk ^ 53 + zk ^ 52 + zk ^ 51 + zk ^ 49 + zk ^ 48 + zk ^ 47 + zk ^ 45 + zk ^ 43 + zk ^ 40 + zk ^ 39 + zk ^ 38 + zk ^ 37 + zk ^ 35 + zk ^ 32 + zk ^ 30 + zk ^ 27 + zk ^ 26 + zk ^ 25 + zk ^ 20 + zk ^ 19 + zk ^ 15 + zk ^ 14 + zk ^ 13 + zk ^ 12 + zk ^ 11 + zk ^ 10 + zk ^ 8 + zk ^ 7 + zk ^ 5 + zk ^ 4 + zk ^ 3
h3_pairing = zk ^ 338 + zk ^ 336 + zk ^ 334 + zk ^ 333 + zk ^ 331 + zk ^ 329 + zk ^ 328 + zk ^ 327 + zk ^ 326 + zk ^ 324 + zk ^ 323 + zk ^ 322 + zk ^ 321 + zk ^ 320 + zk ^ 317 + zk ^ 314 + zk ^ 313 + zk ^ 307 + zk ^ 306 + zk ^ 301 + zk ^ 300 + zk ^ 297 + zk ^ 296 + zk ^ 295 + zk ^ 293 + zk ^ 291 + zk ^ 285 + zk ^ 282 + zk ^ 279 + zk ^ 278 + zk ^ 274 + zk ^ 273 + zk ^ 272 + zk ^ 270 + zk ^ 267 + zk ^ 260 + zk ^ 258 + zk ^ 257 + zk ^ 255 + zk ^ 254 + zk ^ 253 + zk ^ 252 + zk ^ 251 + zk ^ 249 + zk ^ 248 + zk ^ 247 + zk ^ 245 + zk ^ 242 + zk ^ 240 + zk ^ 239 + zk ^ 236 + zk ^ 235 + zk ^ 234 + zk ^ 233 + zk ^ 232 + zk ^ 229 + zk ^ 228 + zk ^ 226 + zk ^ 222 + zk ^ 221 + zk ^ 220 + zk ^ 216 + zk ^ 215 + zk ^ 214 + zk ^ 213 + zk ^ 212 + zk ^ 211 + zk ^ 210 + zk ^ 208 + zk ^ 207 + zk ^ 206 + zk ^ 205 + zk ^ 202 + zk ^ 198 + zk ^ 195 + zk ^ 194 + zk ^ 192 + zk ^ 187 + zk ^ 181 + zk ^ 179 + zk ^ 178 + zk ^ 177 + zk ^ 176 + zk ^ 175 + zk ^ 174 + zk ^ 173 + zk ^ 170 + zk ^ 167 + zk ^ 166 + zk ^ 162 + zk ^ 161 + zk ^ 160 + zk ^ 157 + zk ^ 155 + zk ^ 153 + zk ^ 151 + zk ^ 150 + zk ^ 147 + zk ^ 145 + zk ^ 141 + zk ^ 139 + zk ^ 137 + zk ^ 136 + zk ^ 135 + zk ^ 133 + zk ^ 132 + zk ^ 130 + zk ^ 129 + zk ^ 128 + zk ^ 127 + zk ^ 126 + zk ^ 125 + zk ^ 123 + zk ^ 113 + zk ^ 111 + zk ^ 107 + zk ^ 104 + zk ^ 99 + zk ^ 98 + zk ^ 92 + zk ^ 90 + zk ^ 89 + zk ^ 88 + zk ^ 82 + zk ^ 81 + zk ^ 78 + zk ^ 77 + zk ^ 74 + zk ^ 72 + zk ^ 71 + zk ^ 70 + zk ^ 69 + zk ^ 68 + zk ^ 67 + zk ^ 65 + zk ^ 60 + zk ^ 58 + zk ^ 56 + zk ^ 55 + zk ^ 54 + zk ^ 53 + zk ^ 49 + zk ^ 46 + zk ^ 44 + zk ^ 42 + zk ^ 41 + zk ^ 39 + zk ^ 31 + zk ^ 30 + zk ^ 23 + zk ^ 20 + zk ^ 19 + zk ^ 17 + zk ^ 16 + zk ^ 15 + zk ^ 13 + zk ^ 12 + zk ^ 11 + zk ^ 10 + zk ^ 8 + zk ^ 6 + zk ^ 4 + zk ^ 2

m = [7, 10113049, 320021624768405574452943847]
Bob_sk = GF2n_Ph(m, g3_pairing, h3_pairing, q)

assert Bob_sk * Alice_pk == Shared_key
print('Flag of "HCDH": HFCTF{{{}}}'.format(sha256(str(Bob_sk).encode()).hexdigest()))