Sage_10_3_Setup

Posted on In Word count in article: 12k Reading time ≈ 11 mins.

Sagemath 10.3 版本包含了很多令人惊喜的新功能，我在学习 flatter 的过程中顺便把以前的 Sagemath 替换成了新版。

0x01 安装 Sagemath 10.3

首先需要下载：

  • Miniconda3
  • Mambaforce: 孩子们我没意见 给 conda 提速
  • SageMath 10.3 的源码
1
2
3
wget https://repo.anaconda.com/miniconda/Miniconda3-latest-Linux-x86_64.sh
curl -L -O https://github.com/conda-forge/miniforge/releases/latest/download/Mambaforge-$(uname)-$(uname -m).sh
wget https://mirrors.aliyun.com/sagemath/src/sage-10.3.tar.gz

前两个都是安装脚本，运行后默认装在 ~ 路径下 (最好别改成自己的)，安装完成后激活 bashrc ， miniconda3 ：

1
2
3
sh Miniconda3-latest-Linux-x86_64.sh
source ~/miniconda3/bin/activate
conda --version

mambaforge ：

1
2
3
sh Mambaforge-Linux-x86_64.sh
source ~/mambaforge/bin/activate
mamba --version

如果想要去掉前面的 base ，可以：

1
2
echo "conda deactivate" >> ~/.bashrc
source ~/.bashrc

这样就可以 mamba activate/deactivate 直接控制打开/关闭了。

解压 Sagemath 源码压缩包：

1
2
tar xf ./sage-10.3.tar.gz
cd ./sage-10.3

在 build 之前先给 conda 换源：

1
2
3
4
5
6
channels:
  - https://mirrors.tuna.tsinghua.edu.cn/anaconda/pkgs/main/
  - https://mirrors.tuna.tsinghua.edu.cn/anaconda/pkgs/free/
  - https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/conda-forge/
  - https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/pytorch/
ssl_verify: true

最后按官网的教程做就行了，build 的时候需要在 sage 目录里执行：

1
2
3
4
5
6
7
8
9
10
11
12
export SAGE_NUM_THREADS=16
./bootstrap-conda
mamba env create --file src/environment-dev-3.11-linux.yml --name sage-dev
mamba activate sage-dev
./bootstrap
./configure --with-python=$CONDA_PREFIX/bin/python \
              --prefix=$CONDA_PREFIX
make
./bootstrap
pip install --no-build-isolation -v -v --editable ./pkgs/sage-conf_conda ./pkgs/sage-setup
pip install --no-build-isolation --config-settings editable_mode=compat -v -v --editable ./src
sage

一个惊喜的点是：environment-dev-3.11-linux.yml 包含了至关重要的 fplll 工具。有多重要呢？就处在 flatter 工具链的上游。换言之，Sage 会不会在将来的版本中直接加入 flatter 作为新的特性呢？

不管怎么说，Sage 到这里就可以开始用了。同时注意：sage -n jupyter --allow-roots 命令也可以联动 vscode 编程。

0x02 flatter & msolve

flatter: 一款用来大大改善 LLL 花费时间的工具。老老实实按 Github 上的教程走就行了：

1
2
3
4
5
6
7
$ sudo apt install libgmp-dev libmpfr-dev fplll-tools \
	libfplll-dev libeigen3-dev
$ mkdir build && cd ./build
$ cmake ..
$ make
$ sudo make install
$ sudo ldconfig

msolve : Multivariate polynomial system solver，直译为多元多项式系统求解器。为什么要用他呢？因为在 sage 的新特性中：

  • 使用 msolve 计算 Gröbner 基
1
2
3
4
5
6
7
8
9
10
11
12
13
14
from sage.rings.polynomial.msolve import groebner_basis_degrevlex

R.<a,b,c> = PolynomialRing(GF(101), 3, order='lex')
I = sage.rings.ideal.Katsura(R,3)
gb = groebner_basis_degrevlex(I); gb 
# [a + 2*b + 2*c - 1, b*c - 19*c^2 + 10*b + 40*c,
# b^2 - 41*c^2 + 20*b - 20*c, c^3 + 28*c^2 - 37*b + 13*c]
gb.universe() is R 
# False
gb.universe().term_order() 
# Degree reverse lexicographic term order
ideal(gb).transformed_basis(other_ring=R) 
# [c^4 + 38*c^3 - 6*c^2 - 6*c, 30*c^3 + 32*c^2 + b - 14*c,
# a + 2*b + 2*c - 1]
  • Gröbner 基，proof = False
1
2
3
4
5
6
7
8
R.<x, y> = PolynomialRing(QQ, 2)
I = Ideal([ x*y - 1, (x-2)^2 + (y-1)^2 - 1])
I.groebner_basis(algorithm='msolve') 
# Traceback (most recent call last):
# ...
# ValueError: msolve relies on heuristics; please use proof=False
I.groebner_basis(algorithm='msolve', proof=False) 
# [x*y - 1, x^2 + y^2 - 4*x - 2*y + 4, y^3 - 2*y^2 + x + 4*y - 4]
  • variety
1
2
3
4
5
6
7
8
9
from sage.rings.polynomial.msolve import variety
p = 536870909
R.<x, y> = PolynomialRing(GF(p), 2, order='lex')
I = Ideal([ x*y - 1, (x-2)^2 + (y-1)^2 - 1])
sorted(variety(I, GF(p^2), proof=False), key=str) 
# [{x: 1, y: 1},
# {x: 254228855*z2 + 114981228, y: 232449571*z2 + 402714189},
# {x: 267525699, y: 473946006},
# {x: 282642054*z2 + 154363985, y: 304421338*z2 + 197081624}]

0x03 Power!

工具都装好后，就可以写一个程序直观的感受 flatter 的威力了：

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
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
# Copied from https://gist.github.com/maple3142/18cde45198304f844eb619caff67e179 and altered a little bit XD.
from subprocess import check_output, DEVNULL, CalledProcessError
import itertools
import IPython
import time


def to_fplll_format(M):
    m, n = M.dimensions()
    ret = ""
    s = "["
    for i in range(m):
        s += "["
        for j in range(n):
            s += str(M[i, j])
            if j < n - 1:
                s += " "
        s += "]"
        ret += s + "\n"
        s = ""
    ret += "]"
    return ret


def from_fplll_format(s):
    rows = []
    for line in s.splitlines():
        line = line.lstrip("[").rstrip("\n").rstrip("]")
        if len(line) == 0:
            break

        row = [int(x) for x in line.split(" ") if len(x) > 0 and x != "]"]
        rows += [row]
    m = len(rows)
    n = len(rows[0])
    for row in rows:
        assert len(row) == n

    L = Matrix(ZZ, m, n)
    for i in range(m):
        for j in range(n):
            L[i, j] = rows[i][j]
    return L


def sort_by_norm(M):
    # idk why .norm() is really slow
    M2 = [(sum([y ^ 2 for y in x]), x) for x in M]
    M2.sort()
    return matrix([y for x, y in M2])


def LLL(M):
    return M.LLL()


def flatter(M):
    # compile https://github.com/keeganryan/flatter and put it in $PATH
    s = to_fplll_format(M)
    try:
        ret = check_output(["flatter"], input=s.encode(), stderr=DEVNULL)
        return from_fplll_format(ret.decode())
    except CalledProcessError:
        print("Flatter error, matrix written to /tmp/flatter_error")
        with open("/tmp/flatter_error", "w") as f:
            f.write(s)
        return M.change_ring(ZZ)


def subset_sum():
    print("=" * 40)
    print("Subset sum")
    print("=" * 40)
    # (very) low density subset sum
    n = 128
    B = vector([ZZ(randrange(0, 1 << 2048)) for _ in range(n)])
    s = random_vector(Zmod(2), n).change_ring(ZZ)
    C = B * s
    density = n / log(max(B), 2)
    print("density =", density.n())

    M = Matrix(ZZ, n + 1, n + 1)
    for i in range(n):
        M[i, i] = 2
        M[i, n] = B[i]
        M[n, i] = -1
    M[n, n] = -C
    M[:, n] *= M.det()

    print("flatter")
    print(IPython.get_ipython().run_line_magic("time", "flatter(M)[0]"))
    print("LLL")
    print(IPython.get_ipython().run_line_magic("time", "M.LLL()[0]"))
    print(vector([2 * x - 1 for x in s]))


def small_roots(f, bounds, m=1, d=None, reduction=LLL):
    # modified from https://github.com/defund/coppersmith/blob/master/coppersmith.sage
    if not d:
        d = f.degree()

    R = f.base_ring()
    N = R.cardinality()

    f /= f.coefficients().pop(0)
    f = f.change_ring(ZZ)

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

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

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

    B = reduction(B.dense_matrix())

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

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

    return []


def univariate():
    print("=" * 40)
    print("Univariate Coppersmith")
    print("=" * 40)
    p = random_prime(2 ^ 1024, proof=False)
    q = random_prime(2 ^ 1024, proof=False)
    N = p * q
    bounds = (floor(N ^ 0.315),)
    roots = tuple(randrange(bound) for bound in bounds)
    R = Integers(N)
    P = PolynomialRing(R, 1, "x")
    x = P.gen()
    monomials = [x, x ^ 2, x ^ 3]
    f = sum(randrange(N) * monomial for monomial in monomials)
    f -= f(*roots)
    print("flatter")
    print(
        IPython.get_ipython().run_line_magic(
            "time", "small_roots(f, bounds, m=12, reduction=flatter)"
        )
    )
    print("LLL")
    print(IPython.get_ipython().run_line_magic("time", "small_roots(f, bounds, m=12)"))


def bivariate():
    def flatter_echelon(M):
        st = time.time()
        M = M.echelon_form(algorithm="pari", include_zero_rows=False)
        print("Echelon form done", time.time() - st)
        return flatter(M)

    def flatter_echelon_sort_by_norm(M):
        st = time.time()
        M = M.echelon_form(algorithm="pari", include_zero_rows=False)
        M = sort_by_norm(M)
        print("Echelon form + sort done", time.time() - st)
        return flatter(M)

    print("=" * 40)
    print("Bivariate Coppersmith")
    print("=" * 40)
    p = random_prime(2 ^ 1024, proof=False)
    q = random_prime(2 ^ 1024, proof=False)
    N = p * q
    bounds = (floor(N ^ 0.12), floor(N ^ 0.12))
    roots = tuple(randrange(bound) for bound in bounds)
    R = Integers(N)
    P = PolynomialRing(R, "x, y")
    x, y = P.gens()
    monomials = [x, y, x * y, x ^ 2, y ^ 2]
    f = sum(randrange(N) * monomial for monomial in monomials)
    f -= f(*roots)
    print("flatter")
    print(
        IPython.get_ipython().run_line_magic(
            "time", "small_roots(f, bounds, m=4, d=3, reduction=flatter)"
        )
    )
    print("flatter_echelon")
    print(
        IPython.get_ipython().run_line_magic(
            "time", "small_roots(f, bounds, m=4, d=3, reduction=flatter_echelon)"
        )
    )
    print("flatter_echelon_sort_by_norm")
    print(
        IPython.get_ipython().run_line_magic(
            "time",
            "small_roots(f, bounds, m=4, d=3, reduction=flatter_echelon_sort_by_norm)",
        )
    )
    print("LLL")
    print(
        IPython.get_ipython().run_line_magic("time", "small_roots(f, bounds, m=4, d=3)")
    )


if __name__ == "__main__":
    # copy this into sage repl or IPython wouldn't be available
    subset_sum()
    univariate()
    bivariate()
'''
========================================
Subset sum
========================================
density = 0.0625000360302351
flatter
CPU times: user 13.8 ms, sys: 114 µs, total: 13.9 ms
Wall time: 9.1 s
(1, -1, 1, -1, -1, 1, -1, 1, -1, 1, -1, 1, 1, 1, -1, 1, -1, 1, -1, -1, -1, 1, 1, 1, 1, -1, -1, -1, -1, 1, 1, -1, 1, 1, -1, -1, -1, -1, 1, -1, -1, -1, -1, 1, -1, 1, -1, 1, 1, 1, -1, -1, 1, 1, 1, -1, -1, 1, 1, -1, -1, 1, -1, 1, 1, -1, -1, 1, -1, -1, 1, -1, 1, 1, 1, -1, 1, 1, -1, 1, -1, -1, -1, 1, 1, 1, -1, -1, -1, 1, -1, -1, 1, 1, -1, 1, 1, 1, 1, 1, -1, -1, -1, -1, -1, -1, 1, -1, 1, -1, 1, 1, -1, 1, -1, -1, -1, -1, -1, -1, 1, 1, 1, 1, -1, 1, 1, 1, 0)
LLL
CPU times: user 13.5 s, sys: 9.3 ms, total: 13.5 s
Wall time: 13.5 s
(-1, 1, -1, 1, 1, -1, 1, -1, 1, -1, 1, -1, -1, -1, 1, -1, 1, -1, 1, 1, 1, -1, -1, -1, -1, 1, 1, 1, 1, -1, -1, 1, -1, -1, 1, 1, 1, 1, -1, 1, 1, 1, 1, -1, 1, -1, 1, -1, -1, -1, 1, 1, -1, -1, -1, 1, 1, -1, -1, 1, 1, -1, 1, -1, -1, 1, 1, -1, 1, 1, -1, 1, -1, -1, -1, 1, -1, -1, 1, -1, 1, 1, 1, -1, -1, -1, 1, 1, 1, -1, 1, 1, -1, -1, 1, -1, -1, -1, -1, -1, 1, 1, 1, 1, 1, 1, -1, 1, -1, 1, -1, -1, 1, -1, 1, 1, 1, 1, 1, 1, -1, -1, -1, -1, 1, -1, -1, -1, 0)
(-1, 1, -1, 1, 1, -1, 1, -1, 1, -1, 1, -1, -1, -1, 1, -1, 1, -1, 1, 1, 1, -1, -1, -1, -1, 1, 1, 1, 1, -1, -1, 1, -1, -1, 1, 1, 1, 1, -1, 1, 1, 1, 1, -1, 1, -1, 1, -1, -1, -1, 1, 1, -1, -1, -1, 1, 1, -1, -1, 1, 1, -1, 1, -1, -1, 1, 1, -1, 1, 1, -1, 1, -1, -1, -1, 1, -1, -1, 1, -1, 1, 1, 1, -1, -1, -1, 1, 1, 1, -1, 1, 1, -1, -1, 1, -1, -1, -1, -1, -1, 1, 1, 1, 1, 1, 1, -1, 1, -1, 1, -1, -1, 1, -1, 1, 1, 1, 1, 1, 1, -1, -1, -1, -1, 1, -1, -1, -1)
========================================
Univariate Coppersmith
========================================
flatter
CPU times: user 577 ms, sys: 10.3 ms, total: 587 ms
Wall time: 3.02 s
[(52898868887146752972212802565290540830347003906693038177010989263599521640793523062140183403701304030725704274812929807988885129712177601140616553514257943382667480845507998269206545006670478278,)]
LLL
CPU times: user 21.6 s, sys: 0 ns, total: 21.6 s
Wall time: 22 s
[(52898868887146752972212802565290540830347003906693038177010989263599521640793523062140183403701304030725704274812929807988885129712177601140616553514257943382667480845507998269206545006670478278,)]
========================================
Bivariate Coppersmith
========================================
flatter
Flatter error, matrix written to /tmp/flatter_error
CPU times: user 52.5 ms, sys: 0 ns, total: 52.5 ms
Wall time: 795 ms
[(0, 0)]
flatter_echelon
Echelon form done 0.05237317085266113
CPU times: user 158 ms, sys: 0 ns, total: 158 ms
Wall time: 20.4 s
[(20370105518095839079344261190460931964311641178643539749741763646736449383, 54661297645350919577042112920463577161355888674536247152335056763963988510)]
flatter_echelon_sort_by_norm
Echelon form + sort done 0.06062889099121094
CPU times: user 160 ms, sys: 9.74 ms, total: 169 ms
Wall time: 6.31 s
[(20370105518095839079344261190460931964311641178643539749741763646736449383, 54661297645350919577042112920463577161355888674536247152335056763963988510)]
LLL
CPU times: user 1.14 s, sys: 7 µs, total: 1.14 s
Wall time: 1.25 s
[(20370105518095839079344261190460931964311641178643539749741763646736449383, 54661297645350919577042112920463577161355888674536247152335056763963988510)]

'''

0x04 Some annoying question

首先，直接把原生的 flatter 拼接到解多元 copper会报错，原因是 flatter 在对行数大于列数的格规约时会直接报错，虽然 maple 引进了 flatter_echelon 和 flatter_echelon_sort_by_norm 这两个函数缓解了这个问题，但没有彻底根除问题，这个 issue 到现在还是 open 的。

其次，flatter 的维数提高会导致栈溢出崩溃。可以通过设置 sage: pari.allocatemem(10^9) 的方式缓解，但这个貌似是受限于硬件，没法有效解决这个问题，这个 issue 也是 open 的。

总之，优化 LLL 任重道远，我在一边查找更好的板子，同时也在分析源码，看看能在哪里做进一步优化。