/
bernstein.py
189 lines (180 loc) · 6.49 KB
/
bernstein.py
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
'''
This file was *autogenerated* from the file bernstein.sage
sage code taken from Daniel J. Bernstein. Understanding binary-Goppa decoding. Cryptology ePrint Archive, Paper 2022/473. https://eprint.iacr.org/2022/473. 2022.
'''
from sage.all_cmdline import * # import sage library
_sage_const_0 = Integer(0); _sage_const_1 = Integer(1); _sage_const_2 = Integer(2); _sage_const_100 = Integer(100); _sage_const_3 = Integer(3); _sage_const_10 = Integer(10)
def interpolator(n,k,a,r):
a,r = list(a),list(r)
assert k.is_field()
assert len(a) == n and len(set(a)) == n and len(r) == n
kpoly = k['x']; (x,) = kpoly._first_ngens(1)
A = kpoly(prod(x-a[j] for j in range(n)))
Aprime = A.derivative()
return kpoly(sum(((r[i]/Aprime(a[i]))*(A//(x-a[i]))) for i in range(n)))
def approximant(t,k,A,B):
assert t >= _sage_const_0 and A.base_ring() == k and B.base_ring() == k
kpoly,n = A.parent(),A.degree()
assert n > B.degree()
M = [ [ B[t+n-_sage_const_1 -i-j] for i in range(t+_sage_const_1 )] + [-A[t+n-_sage_const_1 -i-j] for i in range(t) ] for j in range(_sage_const_2 *t)]
M = matrix(k,_sage_const_2 *t,_sage_const_2 *t+_sage_const_1 ,M)
ab = list(M.right_kernel().gens()[_sage_const_0 ])
a,b = kpoly(ab[:t+_sage_const_1 ]),kpoly(ab[t+_sage_const_1 :])
d = gcd(a,b)
return a//d,b//d
def interpolator_with_errors(n,t,k,alpha,r):
alpha,r = list(alpha),list(r)
assert k.is_field()
assert len(alpha) == n and len(set(alpha)) == n and len(r) == n
B = interpolator(n,k,alpha,r)
kpoly = B.parent()
A = kpoly(prod(kpoly([-alpha[j],_sage_const_1 ]) for j in range(n)))
a,b = approximant(t,k,A,B)
if a.divides(A):
if a*B-b*A == _sage_const_0 or (a*B-b*A).degree() < n-_sage_const_2 *t+a.degree():
return B-b*A//a
def goppa_errors(n,t,k,alpha,g,r):
alpha,r = list(alpha),list(r)
assert k.is_field() and k.characteristic() == _sage_const_2
assert g.base_ring() == k and g.degree() == t and g.is_squarefree()
assert len(alpha) == n and len(set(alpha)) == n and len(r) == n
kpoly = g.parent()
A = kpoly(prod(kpoly([-alpha[j],_sage_const_1 ]) for j in range(n)))
Aprime = A.derivative()
rtwist = [r[i]*Aprime(alpha[i])/g(alpha[i])**_sage_const_2 for i in range(n)]
B = interpolator(n,k,alpha,rtwist)
a,b = approximant(t,k,A,B)
aprime = a.derivative()
if a.divides(A):
if a.divides(g**_sage_const_2 *b-aprime):
if a*B-b*A == _sage_const_0 or (a*B-b*A).degree() < n-_sage_const_2 *t+a.degree():
return [k(a(alpha[j]) == _sage_const_0 ) for j in range(n)]
def test_interpolator():
for q in range(_sage_const_100 ):
q = ZZ(q)
if not q.is_prime_power():
continue
print('interp %d' % q)
sys.stdout.flush()
k = GF(q)
for loop in range(_sage_const_100 ):
n = randrange(q+_sage_const_1 )
a = list(k)
shuffle(a)
a = a[:n]
r = [k.random_element() for j in range(n)]
phi = interpolator(n,k,a,r)
assert phi.degree() < n
assert all(phi(aj) == rj for aj,rj in zip(a,r))
kpoly = phi.parent()
assert phi == kpoly.lagrange_polynomial(zip(a,r))
def test_approximant():
for q in range(_sage_const_100 ):
q = ZZ(q)
if not q.is_prime_power():
continue
print('approximant %d' % q)
sys.stdout.flush()
k = GF(q)
kpoly = k['x']; (x,) = kpoly._first_ngens(1)
for loop in range(_sage_const_100 ):
Adeg = randrange(_sage_const_100 )
A = kpoly([k.random_element() for j in range(Adeg)]+[_sage_const_1 ])
if Adeg == _sage_const_0 :
B = kpoly(_sage_const_0 )
else:
Bdeg = randrange(Adeg)
B = kpoly([k.random_element() for j in range(Bdeg+_sage_const_1 )])
# note that B could actually have lower degree
t = randrange(Adeg+_sage_const_3 )
a,b = approximant(t,k,A,B)
assert gcd(a,b) == _sage_const_1
assert a.degree() <= t
assert b.degree() < t
assert a != _sage_const_0
assert a*B-b*A == _sage_const_0 or (a*B-b*A).degree() < A.degree()-t
def test_interpolator_with_errors():
for q in range(_sage_const_100 ):
q = ZZ(q)
if not q.is_prime_power():
continue
print('interpolator_with_errors %d' % q)
sys.stdout.flush()
k = GF(q)
kpoly = k['x']; (x,) = kpoly._first_ngens(1)
for loop in range(_sage_const_100 ):
n = randrange(q+_sage_const_1 )
t = randrange(_sage_const_3 +n//_sage_const_2 )
a = list(k)
shuffle(a)
a = a[:n]
for known in True,False:
if known:
f = kpoly([k.random_element() for j in range(n-_sage_const_2 *t)])
r = list(map(f,a))
e = [k.random_element() for j in range(t)]+[_sage_const_0 ]*(n-t)
shuffle(e)
assert len([ej for ej in e if ej != _sage_const_0 ]) <= t
for j in range(n):
r[j] += e[j]
else:
f = 'unknown' # cut off data flow from previous iteration
r = [k.random_element() for j in range(n)]
f2 = interpolator_with_errors(n,t,k,a,r)
if f2 == None:
assert not known
else:
assert f2 == _sage_const_0 or f2.degree() < n-_sage_const_2 *t
if known:
assert f2 == f
assert len([j for j in range(n) if f2(a[j]) != r[j]]) <= t
def test_goppa_errors():
for m in range(_sage_const_1 ,_sage_const_10 ):
q = _sage_const_2 **m
print('goppa_errors %d' % q)
sys.stdout.flush()
k = GF(q)
kpoly = k['x']; (x,) = kpoly._first_ngens(1)
for loop in range(_sage_const_100 ):
while True:
n = randrange(q+_sage_const_1 )
t = randrange(_sage_const_3 +n//m)
if t >= n:
t = n
a = list(k)
shuffle(a)
a = a[:n]
g = kpoly([k.random_element() for j in range(t)]+[_sage_const_1 ])
if g.is_squarefree():
if all(g(aj) != _sage_const_0 for aj in a):
break
assert g.degree() == t
A = kpoly(prod(x-aj for aj in a))
Aprime = A.derivative()
for aj in a:
assert Aprime(aj) != _sage_const_0
for known in True,False:
if known:
f = kpoly([k.random_element() for j in range(n-_sage_const_2 *t)])
r = [(f*g**_sage_const_2 )(aj)/Aprime(aj) for aj in a]
if randrange(_sage_const_2 ):
e = [_sage_const_1 ]*t+[_sage_const_0 ]*(n-t)
else:
actualweight = randrange(t+_sage_const_1 )
e = [_sage_const_1 ]*actualweight+[_sage_const_0 ]*(n-actualweight)
shuffle(e)
assert len([ej for ej in e if ej != _sage_const_0 ]) <= t
for j in range(n):
r[j] += e[j]
else:
e = 'unknown' # cut off data flow from previous iteration
r = [k.random_element() for j in range(n)]
e2 = goppa_errors(n,t,k,a,g,r)
if e2 == None:
assert not known
else:
assert len(e2) == n
if known:
assert e2 == e
assert len([ej for ej in e2 if ej != _sage_const_0 ]) <= t
assert g.divides(sum((r[i]-e2[i])*A//(x-a[i]) for i in range(n)))