-
Notifications
You must be signed in to change notification settings - Fork 0
/
ENUM.cpp
322 lines (296 loc) · 15.8 KB
/
ENUM.cpp
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
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
#include<iostream>
#include <vector>
#include <cmath>
#include<random>
#include<time.h>
#include<tuple>
/*
~$ g++ -o ENUM ENUM.cpp
~$ ./ENUM
*/
/*
参考文献：安田雅哉、青野良範著『格子暗号解読のための数学的基礎』(2019)
*/
using namespace std;
//単位行列の作成
vector<vector<double>> identity_mat(int n){
vector<vector<double>> A(n, vector<double>(n));
for(int i = 0; i < n; i++){
A.at(i).at(i) = 1.0;
}
return A;
}
//行列の表示
void print_mat(vector<vector<double>> A){
cout << "[\n";
for(int i = 0; i < A.size(); i++){
cout << "[";
for(int j = 0; j < A.at(0).size(); j++) cout << A.at(i).at(j) << " ";
cout << "]" << endl;
}
cout << "]\n";
}
//vectorの表示
void print_vec(vector<double> v){
cout << "[";
for(int i = 0; i < v.size(); i++) cout << v.at(i) << " ";
cout << "]\n";
}
//zero-vectorか否かの判定
bool zero_test(vector<double> x){
int counter = 0;
for(int i = 0; i < x.size(); i++){
if(x.at(i) == 0) counter++;
}
if(counter == x.size()){return true;}else{return false;}
}
//内積
double dot(vector<double> x, vector<double> y)
{
double z = 0;
for(int i = 0; i < x.size(); i++){
z += x.at(i) * y.at(i);
}
return z;
}
//転置行列
vector<vector<double>> trans(vector<vector<double>> A)
{
vector<vector<double>> B(A.at(0).size(), vector<double>(A.size()));
for(int i = 0; i < A.at(0).size(); i++){
for(int j = 0; j < A.size(); j++) B.at(i).at(j) = A.at(j).at(i);
}
return B;
}
//行列式
double det(vector<vector<double>> A){
int n = A.size();
double det = 1.0, buf;
for(int i = 0; i < n; i++){
for(int j = 0; j < n; j++){
if(i < j){
buf = A.at(j).at(i) / A.at(i).at(i);
for(int k = 0; k < n; k++) A.at(j).at(k) -= A.at(i).at(k) * buf;
}
}
}
for(int i = 0; i < n; i++) det *= A.at(i).at(i);
return det;
}
//行列の積
vector<vector<double>> mul_mat(vector<vector<double>> A, vector<vector<double>> B)
{
int n = A.size(), m = B.at(0).size();
vector<vector<double>> C(n, vector<double>(m));
for(int i = 0; i < n; i++){
for(int j = 0; j < m; j++){
for(int k = 0; k < B.size(); k++) C.at(i).at(j) += A.at(i).at(k) * B.at(k).at(j);
}
}
return C;
}
//Gram_Schmidtの直交化
tuple<vector<vector<double>>, vector<vector<double>>> Gram_Schmidt(vector<vector<double>> b)
{
int n = b.size(), m = b.at(0).size();
vector<vector<double>> GSOb(n, vector<double>(m));
vector<vector<double>> mu = identity_mat(n);
for(int i = 0; i < n; i++){
for(int j = 0; j < m; j++) GSOb.at(i).at(j) = b.at(i).at(j);
for(int j = 0; j < i; j++){
mu.at(i).at(j) = dot(b.at(i), GSOb.at(j)) / dot(GSOb.at(j), GSOb.at(j));
for(int k = 0; k < m; k++) GSOb.at(i).at(k) -= mu.at(i).at(j) * GSOb.at(j).at(k);
}
}
return forward_as_tuple(GSOb, mu);
//if(flag){return GSOb;}else{return mu;}
}
//部分size基簡約
tuple<vector<vector<double>>, vector<vector<double>>> size_reduce(vector<vector<double>> b, int i, int j){
int q;
vector<vector<double>> GSOb, mu;
tie(GSOb, mu) = Gram_Schmidt(b);
if(fabs(mu.at(i).at(j)) > 0.5){
q = round(mu.at(i).at(j));
for(int k = 0; k < b.size(); k++) b.at(i).at(k) -= (double)q * b.at(j).at(k);
for(int k = 0; k < j + 1; j++){
mu.at(i).at(k) -= (double)q * mu.at(j).at(k);
mu.at(j).at(j) = 1.0;
}
}
return forward_as_tuple(b, mu);
}
//二乗normの作成
vector<double> make_square_norm(vector<vector<double>> GSOb){
vector<double> B(GSOb.size());
for(int i = 0; i < GSOb.size(); i++) B.at(i) = dot(GSOb.at(i), GSOb.at(i));
return B;
}
//体積
double vol(vector<vector<double>> b){
double v = 1.0;
vector<vector<double>> GSOb, mu;
tie(GSOb, mu) = Gram_Schmidt(b);
for(int i = 0; i < b.size(); i++) v *= sqrt(dot(GSOb.at(i), GSOb.at(i)));
return v;
}
//Randomな正方正則行列を作成する。
vector<vector<double>> random_mat(int n){
vector<vector<double>> A(n, vector<double>(n));
while(1){
for(int i = 0; i < n; i++){
for(int j = 0; j < n; j++) A.at(i).at(j) = (2 * (rand() % 2) - 1) * (rand() % 100);
}
if(det(A) != 0) return A;
}
}
//係数vectorを格子に変換
vector<double> coef2lat(vector<double> v, vector<vector<double>> b){
vector<double> x(b.at(0).size());
for(int i = 0; i < b.size(); i++){
for(int j = 0; j < b.at(0).size(); j++) x.at(j) += v.at(i) * b.at(i).at(j);
}
return x;
}
//射影格子の基底の出力
vector<vector<double>> project_basis(int k, int l, vector<vector<double>> b){
int n = b.size(), m = b.at(0).size();
vector<vector<double>> GSOb, mu, pi_b(l - k + 1, vector<double>(m));
tie(GSOb, mu) = Gram_Schmidt(b);
for(int i = k; i <= l; i++){
for(int j = k; j < n; j++){
for(int h = 0; h < m; h++)pi_b.at(i - k).at(h) += (dot(b.at(i), GSOb.at(j))) / dot(GSOb.at(j), GSOb.at(j)) * GSOb.at(j).at(h);
}
}
return pi_b;
}
//与えられた上階より小さいvectorの数え上げ
vector<double> ENUM(vector<vector<double>> mu, vector<double> B, double R){
int n = B.size();
vector<vector<double>> sigma(n + 2, vector<double>(n + 1));
vector<int> r(n + 1), w(n + 1);
vector<double> c(n + 1), rho(n + 2), v(n + 1), x(n);
int k = 1, last_nonzero = 1;
v[1] = 1;
for(int i = 0; i < n + 1; i++) r.at(i) = i;
while(1){
rho.at(k) = rho.at(k + 1) + (v.at(k) - c.at(k)) * (v.at(k) - c.at(k)) * B.at(k - 1);
if(rho.at(k) <= R * R){
if(k == 1){
for(int i = 0; i < n; i++) x.at(i) = v.at(i + 1);
return x;
}
k--;
r.at(k - 1) = max(r.at(k - 1), r.at(k));
for(int i = r.at(k) ; i >= k + 1; i--){
sigma.at(i).at(k) = sigma.at(i + 1).at(k) + mu.at(i - 1).at(k - 1) * v.at(i);
}
c.at(k) = -sigma.at(k + 1).at(k);
v.at(k) = round(c.at(k));
w.at(k) = 1;
}else{
k ++;
if(k == n + 1){
v.clear();
return v;
}
r.at(k - 1) = k;
if(k >= last_nonzero){
last_nonzero = k;
v.at(k)++;
}else{
if(v.at(k) > c.at(k)){v.at(k) -= w.at(k);}else{v.at(k) += w.at(k);}
w.at(k)++;
}
}
}
}
//最短vectorの数え上げ
vector<double> ENUM_all(vector<vector<double>> mu, vector<double> B, double R, double eps, vector<vector<double>> b){
int n = B.size();
vector<double> ENUM_v(n), pre_ENUM_v(n), x(n), v;
while(1){
for(int i = 0; i < n; i++) pre_ENUM_v.at(i) = ENUM_v.at(i);
ENUM_v = ENUM(mu, B, R);
if(ENUM_v.empty()) return pre_ENUM_v;
v = coef2lat(ENUM_v, b);
R = sqrt(dot(v, v)) * eps;
//cout << R <<endl;
}
}
int main(){
vector<vector<double>> b;
vector<vector<double>> GSOb, mu, a;
bool random;
int n;//格子次元
cout << "Random?" <<endl;
cin >> random;
if(random){
cout << "lattice dimension：" << endl;
cin >> n;
b = random_mat(n);
}else{
//SVP-challenge
//40次元：OK(2.64586秒)
b = {{159, 243, 117, 687, -47, -186, 255, -185, 255, -27, 12, 104, -45, -165, -322, 353, -225, -29, 177, -394, 599, -19, 183, 13, 544, 150, 84, -193, 192, -725, -378, 79, -192, 302, 610, -194, 377, 149, 224, 201},
{-136, -915, 96, -714, 572, -273, 31, -45, 136, 290, -142, -111, 243, 126, -21, 361, 27, 139, 3, -109, 524, -339, 87, 16, -364, -25, 96, -623, 97, 392, 123, -51, -23, -304, -108, 111, -278, -29, -323, -253},
{68, 610, -48, 1065, 405, -119, -52, -115, 181, 504, 263, -246, -248, -156, -349, -264, -689, -118, -370, -2, -74, -128, -22, 55, 403, 455, 112, 341, 115, 136, -217, -118, -419, 194, 272, -91, -97, -225, 152, 259},
{8, -713, 136, -558, -294, -217, 562, -54, 44, -163, -201, 321, -142, 649, 282, 197, 34, 38, -579, 76, 712, -231, 603, -551, -26, -117, 10, 258, -473, -301, 345, -447, 173, -161, -257, -148, -141, 565, -320, -259},
{248, -213, 249, -410, 319, 58, 585, -335, 332, 11, -337, 59, 31, -107, 154, 430, 364, 304, -243, 364, 358, -452, -18, -722, 379, -178, -131, 284, 23, 643, 73, -340, 345, -167, -665, -240, -425, 379, -316, -245},
{-331, -276, -153, -23, -87, -49, -31, -235, -780, -204, 144, -131, -85, -48, 869, -374, -363, -47, -297, 59, 265, 174, -449, 308, -289, -482, 276, -10, -234, 501, 22, 89, 570, 85, 400, 272, 98, -478, 8, -12},
{-62, -922, -114, -368, -305, -534, -324, -20, -710, -740, -15, -171, -477, -285, 283, -368, 83, 85, 51, 304, 627, 13, -463, 393, -278, -84, 253, -233, 32, -303, 173, -141, 127, 103, -182, 326, 524, 168, 47, 1},
{63, 164, 426, 21, -365, 2, -72, 113, 193, 596, 26, -286, 132, 384, 11, -177, -206, 299, -584, 29, -712, -409, 428, -275, -536, -239, 158, 229, -2, 497, 532, 464, 17, -2, 232, -75, -216, -223, 127, -292},
{52, -85, 2, 106, -455, -520, -117, -96, -422, -748, 398, 325, -184, 346, -447, -192, 0, 494, 85, -530, -231, -25, 411, 295, 471, 108, -280, 279, 30, -697, 387, -74, -786, -2, 14, -105, 135, 334, 47, 257},
{-147, -137, -240, -38, 97, 330, -100, -603, -495, -388, 600, -342, 104, 558, 426, 365, -224, 467, 7, -277, 658, -37, -285, 24, 304, 135, 33, -121, -495, -308, -40, -251, -74, 79, 486, -24, -300, -133, -262, 15},
{137, 193, 65, -380, -330, -78, 157, -313, -117, -27, -121, -16, 389, 619, 484, 443, -306, -11, -292, -124, 892, -57, 360, 444, -292, 11, 257, -15, -547, -495, -119, -199, 107, -144, 595, -601, 324, 327, 125, -354},
{123, -290, -123, 303, -284, 702, -403, -97, 206, -186, 38, -197, -146, -6, 194, -117, 292, -345, 400, -39, -221, -22, -393, 190, -252, -289, 317, 57, 149, -606, 482, 12, 267, 154, 182, 224, 324, -531, -145, 26},
{113, 168, 1, 261, 1088, -210, -228, 201, 114, -14, 327, -522, -136, -678, -539, 299, 420, 347, 490, 48, 432, -86, 22, 62, -168, 536, 135, 48, 351, -135, -575, -296, -16, 432, 83, 578, -278, -227, 254, 260},
{-359, -322, -258, -244, 239, 165, 94, 216, -90, -116, 261, 190, 397, -72, 337, 816, 34, 135, 510, -511, 79, 220, 593, -38, -322, 90, 93, -40, 63, -255, -335, -515, -605, 126, -215, 239, -296, -107, 275, 124},
{-92, -178, -354, -708, 456, -827, -84, 384, 9, -274, 4, 864, 313, 513, -451, -65, 205, -81, -418, 69, 427, 59, 585, 331, -91, 446, 284, -327, -216, -162, -126, -184, 364, 80, 47, 179, 25, 370, -15, -23},
{363, 709, 605, 198, -430, 438, 108, -326, 530, 319, -157, -191, -348, -215, -564, -335, 505, -366, 295, 87, -18, -409, 143, -76, 57, 104, -69, 538, -92, -502, -302, -256, -17, 43, -303, -600, 505, -201, 244, -98},
{-414, 134, -531, -677, 540, 241, -187, 132, -199, 603, 168, 22, 381, 796, 287, -13, -232, 65, -224, -95, -795, 461, -697, 146, -74, 598, 221, -579, -266, 222, 121, 320, -426, 16, 240, 29, -404, 68, -271, 73},
{-619, -618, 156, -373, 325, 204, 147, -161, 399, 914, 201, 38, -429, 654, -106, 59, -372, -784, -204, -42, 368, -201, -554, -521, 36, 244, 96, 89, -733, -269, 363, -65, 217, -256, -144, -127, -179, -381, -825, -109},
{-255, -110, -78, -105, -279, 295, 315, -195, 410, 207, 151, 415, -131, 263, -152, 195, -95, 234, -377, 308, -640, -583, -138, -863, -71, 348, 141, 24, 178, 70, 720, 133, -943, 66, -182, -258, -173, -112, -124, 61},
{-4, -512, -613, -318, -330, 358, -219, 39, -319, -420, 263, -135, 0, 18, 591, 61, -191, -677, -44, 143, 305, 433, -176, -403, 155, -30, -212, 228, -418, -446, 702, -124, -92, -204, -361, 381, -292, 30, 13, -378},
{-503, -225, 212, -607, -687, 123, -452, -55, 170, -189, -240, 213, 210, 277, -225, 740, 324, 924, 36, -310, -600, -436, 136, -377, 353, -239, 207, -434, -220, 91, -203, 432, -360, -1, -861, -103, -53, -5, -145, -125},
{-221, 24, -66, -65, 538, -94, -514, -712, -149, 163, 656, -530, 949, 5, 40, -90, -500, -150, -219, -431, 566, -252, 302, 678, 312, 498, -102, -775, 127, 145, -461, -164, 160, 22, 466, 56, -106, -258, 182, -51},
{150, 585, -446, 422, 297, -473, 319, -104, -221, 78, 275, 551, -316, 151, -150, 79, -40, -172, 139, 63, 193, 114, -223, -128, 655, 222, -233, -164, 19, 386, 181, -264, 169, 99, -380, 31, -66, 289, -139, -144},
{221, 370, -334, -102, 34, 164, -9, 134, -155, -40, 180, -211, -364, 390, -349, -990, -298, -291, 15, 274, -326, -184, -769, -265, 646, 114, -36, 399, -84, 621, -119, -243, 409, -118, -775, -227, -122, 192, -201, 424},
{-151, 96, -201, -92, 33, 137, -97, 115, 556, 18, 177, 85, -319, -155, 113, -310, -194, -559, -229, 211, -461, -226, -98, -967, 622, 414, -59, 882, 42, 772, 124, -271, 307, -97, -318, -87, -393, -65, 34, 164},
{363, -70, 307, 27, 185, -4, 307, 429, -267, -564, 56, 291, 125, -233, 133, 349, 676, 384, 129, 400, -794, -74, -155, 127, -602, 161, -67, 192, 346, -147, 265, -218, -156, 25, -513, -76, -59, 153, 101, 483},
{-124, -291, -148, -271, 66, 166, 373, -757, -103, -292, -449, 159, -51, -166, -182, -510, 181, 237, 525, 43, 322, 120, -492, -197, 542, -319, -341, -201, 109, 36, -598, -83, -54, -401, -579, -123, 32, 231, 272, 235},
{47, 829, -14, 935, -69, 206, 439, -461, 375, 85, -284, 45, 324, 609, -514, -179, -224, -207, -509, -133, 99, -69, 286, 139, -245, 227, 214, 349, 5, -446, -88, -346, 366, -307, 226, -453, 330, -24, -171, 199},
{-115, -462, 198, -285, -481, 119, -160, 124, -125, -28, -50, 84, -449, 35, 120, 655, -177, 321, 3, 539, 658, -223, 279, -55, -614, 150, 476, -571, -297, -807, -145, 43, -453, -150, 348, 188, 256, 571, 71, -536},
{101, 603, -300, -481, 298, -44, -154, 456, 189, -698, 254, 407, -157, -79, -26, -986, 690, -325, 422, -452, -603, 308, 381, 393, -134, 133, -418, -26, 77, -541, 287, 246, 86, -15, -402, 425, -115, -358, 69, 321},
{160, -316, 397, 430, -283, 293, 236, -655, 527, 185, -457, -241, 353, -592, -90, 215, -188, 213, 2, 274, 435, -165, -321, -295, -46, 157, -17, -86, -211, -611, 91, -4, 142, 353, 171, 43, 185, -139, -27, 33},
{-410, 749, -628, 254, -302, 149, -322, 373, 124, 108, 275, -72, 389, 189, -7, 262, 67, -488, 220, -291, -201, -4, 235, 578, -61, 398, 145, 137, 662, -346, 396, 130, -627, 149, 76, -56, -2, -329, 104, 309},
{31, -179, -724, 312, -52, 515, -285, 195, 274, 190, 93, -64, 476, -214, -365, 69, 96, 7, -156, 648, -600, 323, 202, -165, -687, 266, 0, 222, 950, -149, -204, 169, -96, 416, -134, 269, -97, -363, 350, 197},
{-124, -713, -134, -199, -274, 235, -101, 83, -409, -197, 68, 379, -68, -232, 74, 201, -282, 179, 70, 122, 346, 131, -657, -172, -361, 188, 223, -893, -80, -178, 82, 690, 494, 215, -37, 131, 350, 77, 0, -600},
{95, 259, -118, -46, -119, 567, 174, 52, 347, 113, -298, 508, 12, 705, -310, -263, 96, 57, -958, 573, -686, 466, 53, -629, 627, -57, -524, 408, -202, -121, -33, -120, -52, -173, -439, 150, 30, 19, -434, 274},
{-11, -763, -57, -475, 346, -53, -84, 614, 115, 586, 558, -315, -317, -424, -234, -98, -49, -426, 14, -7, 180, -236, -107, -146, -503, 158, 72, 116, 5, -153, -125, -180, -804, 185, -214, 299, -114, 693, 115, 108},
{-161, -390, -66, -476, 321, -55, -900, 120, -3, 42, 447, -20, 252, -619, -233, -338, 213, 478, 48, -36, -298, -306, -338, 71, 238, 149, -349, -342, 419, 728, -260, -546, -739, -147, -475, -163, -413, -193, -94, -220},
{571, -274, 497, -102, -270, 530, -212, 3, 290, -15, -286, -105, 570, -455, 16, -299, 562, 559, 238, 739, -51, -494, -117, -175, -644, 135, -358, -21, -5, -158, -271, 19, 284, -330, -104, -96, 124, 303, -51, -502},
{468, 269, -343, -528, 627, 111, 225, -553, -651, 243, 58, -417, 446, -498, 64, -15, 344, 287, -219, -20, -30, 308, -287, 170, 108, -205, -252, -220, 48, 648, 348, -139, 574, -316, -89, 13, -209, 166, -215, -58},
{-544, 115, -183, -328, -330, -290, -71, 354, 379, -157, 272, 297, 10, 540, -141, 252, 82, -134, 55, 206, 252, -180, -215, 527, -481, -76, -591, -407, -144, -131, -101, 55, 620, -386, 479, -141, 423, 121, 241, -97}};
n = b.size();
}
cout << "input basis matrices：" << endl;
print_mat(b);
tie(GSOb, mu) = Gram_Schmidt(b);
vector<double> B = make_square_norm(GSOb), shortest_v;
double R = sqrt((double)n) * pow(vol(b), 1.0 / (double)n) , eps = 0.999;
clock_t start = clock();
shortest_v = ENUM_all(mu, B, R, eps, b);
//a = MLLL(b, 0.75);
clock_t end = clock();
print_vec(coef2lat(shortest_v, b));
cout << "norm=" << sqrt(dot(coef2lat(shortest_v, b), coef2lat(shortest_v, b))) << endl;
//print_mat(a);
cout << "Run time = " << (double)(end - start) / CLOCKS_PER_SEC << "[secs]" << endl;
b.shrink_to_fit();
GSOb.shrink_to_fit();
mu.shrink_to_fit();
B.shrink_to_fit();
shortest_v.shrink_to_fit();
return 0;
}