1 /* Copyright (C) 2015 Martin Albrecht
2
3 This file is part of fplll. fplll is free software: you
4 can redistribute it and/or modify it under the terms of the GNU Lesser
5 General Public License as published by the Free Software Foundation,
6 either version 2.1 of the License, or (at your option) any later version.
7
8 fplll is distributed in the hope that it will be useful,
9 but WITHOUT ANY WARRANTY; without even the implied warranty of
10 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11 GNU Lesser General Public License for more details.
12
13 You should have received a copy of the GNU Lesser General Public License
14 along with fplll. If not, see <http://www.gnu.org/licenses/>. */
15
16 #include <cstring>
17 #include <gso.h>
18 #include <gso_gram.h>
19 #include <gso_interface.h>
20 #include <householder.h>
21 #include <nr/matrix.h>
22 //#include <random>
23 #include <test_utils.h>
24
25 using namespace std;
26 using namespace fplll;
27
28 #ifndef TESTDATADIR
29 #define TESTDATADIR ".."
30 #endif
31
matrix_relative_difference(Matrix<FT> r1,Matrix<FT> r2)32 template <class ZT, class FT> Matrix<FT> matrix_relative_difference(Matrix<FT> r1, Matrix<FT> r2)
33 {
34 Matrix<FT> diff_matrix = Matrix<FT>(r1.get_rows(), r1.get_cols());
35 diff_matrix.fill(0.0);
36 FT relativation_factor = 0.0;
37 for (int i = 0; i < r1.get_rows(); i++)
38 {
39 for (int j = 0; j < i; j++)
40 { // j < i, because r is lower-triangular, and has only 1 on the diagonal.
41 relativation_factor = abs(r1[i][j]) + abs(r2[i][j]);
42 if (relativation_factor.is_zero())
43 {
44 diff_matrix[i][j] = abs(r1[i][j] - r2[i][j]);
45 }
46 else
47 {
48 diff_matrix[i][j] = abs(r1[i][j] - r2[i][j]) / relativation_factor;
49 }
50 }
51 }
52 return diff_matrix;
53 }
54
55 // Returns true when the r-matrices of M1 and M2 are entry-wise equal, up to an error 'error'.
rs_are_equal(MatGSO<ZT,FT> M1,MatGSOGram<ZT,FT> M2,FT error)56 template <class ZT, class FT> bool rs_are_equal(MatGSO<ZT, FT> M1, MatGSOGram<ZT, FT> M2, FT error)
57 {
58 Matrix<FT> r1 = M1.get_r_matrix();
59 Matrix<FT> r2 = M2.get_r_matrix();
60 Matrix<FT> diff = matrix_relative_difference<ZT, FT>(r1, r2);
61
62 FT max_entry = 0.0;
63 max_entry = diff.get_max();
64 if (max_entry > error)
65 {
66 diff.print(cerr);
67 cerr << endl << endl;
68 return false;
69 }
70 return true;
71 }
72
assert_diag_R(FT rhd,int i,int & status)73 template <class FT> void assert_diag_R(FT rhd, int i, int &status)
74 {
75 if (rhd.cmp(0.0) <= 0)
76 {
77 cerr << "R(" << i << ", " << i << ") must be positive." << endl;
78 status = 1;
79 }
80 }
81
assert_mu_r_householder(FT rh,FT rhd,FT mu,FT r,int & status)82 template <class FT> void assert_mu_r_householder(FT rh, FT rhd, FT mu, FT r, int &status)
83 {
84 if (abs(mu - rh / rhd) > 0.001)
85 {
86 cerr << "Error (mu): " << mu << " != " << rh / rhd << endl;
87 status = 1;
88 }
89 if (abs(r - rh * rhd) > 0.001)
90 {
91 cerr << "Error (r): " << r << " != " << rh * rhd << endl;
92 status = 1;
93 }
94 }
95
96 /**
97 * Test if:
98 * mu(i, j) = R(i, j) / R(j, j)
99 * r(i, j) = R(i, j) * R(j, j)
100 */
test_householder(ZZ_mat<ZT> & A)101 template <class ZT, class FT> int test_householder(ZZ_mat<ZT> &A)
102 {
103 int status = 0;
104 ZZ_mat<ZT> U;
105 ZZ_mat<ZT> UT;
106 MatGSO<Z_NR<ZT>, FP_NR<FT>> M(A, U, UT, GSO_INT_GRAM);
107 // Since HOUSEHOLDER_DEFAULT, the exponent returned by get_R(_naively) must be equal to zero.
108 MatHouseholder<Z_NR<ZT>, FP_NR<FT>> Mhouseholder(A, U, UT, HOUSEHOLDER_DEFAULT);
109 M.update_gso();
110 // Here, we just need to refresh R. However, refresh_R() does not modify n_known_rows, but
111 // refresh_R_bf() does.
112 Mhouseholder.refresh_R_bf();
113 Mhouseholder.update_R();
114 Mhouseholder.update_R_naively();
115
116 FP_NR<FT> r;
117 FP_NR<FT> rh;
118 FP_NR<FT> rhd;
119 FP_NR<FT> mu;
120 long expo;
121
122 for (int i = 0; i < A.get_rows(); i++)
123 {
124 Mhouseholder.get_R(rhd, i, i, expo);
125 FPLLL_CHECK(expo == 0, "expo must be equal to 0");
126 assert_diag_R(rhd, i, status);
127 Mhouseholder.get_R_naively(rhd, i, i, expo);
128 FPLLL_CHECK(expo == 0, "expo must be equal to 0");
129 assert_diag_R(rhd, i, status);
130 }
131
132 for (int i = 0; i < A.get_rows(); i++)
133 {
134 for (int j = 0; j < i; j++)
135 {
136 M.get_r(r, i, j);
137 M.get_mu(mu, i, j);
138 Mhouseholder.get_R(rh, i, j, expo);
139 FPLLL_CHECK(expo == 0, "expo must be equal to 0");
140 Mhouseholder.get_R(rhd, j, j, expo);
141 FPLLL_CHECK(expo == 0, "expo must be equal to 0");
142 assert_mu_r_householder(rh, rhd, mu, r, status);
143 Mhouseholder.get_R_naively(rh, i, j, expo);
144 FPLLL_CHECK(expo == 0, "expo must be equal to 0");
145 Mhouseholder.get_R_naively(rhd, j, j, expo);
146 FPLLL_CHECK(expo == 0, "expo must be equal to 0");
147 assert_mu_r_householder(rh, rhd, mu, r, status);
148 }
149 }
150
151 return status;
152 }
153
test_ggso(ZZ_mat<ZT> & A)154 template <class ZT, class FT> int test_ggso(ZZ_mat<ZT> &A)
155 {
156 // TEST A
157 // Method:
158 // - Apply 'normal' MatGSO to A.
159 // - Extract r-matrix of A
160 // - Compute G = A^T A.
161 // - Apply gram MatGSO to G.
162 // - Extract r-matrix of G
163 // -> The r-matrices should be equal.
164
165 // TEST B
166 // Apply some 'random' elementary operation on A and on G
167 // (of course, the operations on A and G are only 'abstractly' the same)
168 // check if their r-matrices are still equal.
169
170 // TEST A
171 // ----------------------------
172 int r = A.r;
173
174 ZZ_mat<ZT> U;
175 ZZ_mat<ZT> UT;
176
177 MatGSO<Z_NR<ZT>, FP_NR<FT>> Mbuf(A, U, UT, GSO_INT_GRAM);
178 Mbuf.update_gso();
179 Matrix<Z_NR<ZT>> G = Mbuf.get_g_matrix();
180
181 MatGSO<Z_NR<ZT>, FP_NR<FT>> M(A, U, UT, GSO_DEFAULT);
182 M.update_gso();
183 MatGSOGram<Z_NR<ZT>, FP_NR<FT>> M2(G, U, UT, 1);
184 M2.update_gso();
185
186 ZZ_mat<ZT> A1(A);
187 MatGSO<Z_NR<ZT>, FP_NR<FT>> M3(A1, U, UT, GSO_INT_GRAM);
188 M3.update_gso();
189
190 FP_NR<FT> err = .001;
191 bool retvalue1 = rs_are_equal(M, M2, err);
192
193 // TEST B
194 // ------------------------
195
196 for (int i = 0; i < rand() % 10 + 1; i++)
197 {
198 int k = rand() % r;
199 int j = rand() % r;
200 M.move_row(k, j);
201 M2.move_row(k, j);
202 M3.move_row(k, j);
203 }
204 M.update_gso();
205 M2.update_gso();
206 M3.update_gso();
207 bool retvalue2 = rs_are_equal(M, M2, err);
208
209 M.row_op_begin(0, r);
210 M2.row_op_begin(0, r);
211 M3.row_op_begin(0, r);
212 for (int i = 0; i < rand() % 10 + 1; i++)
213 {
214 int k = rand() % r;
215 int j = rand() % r;
216 M.row_add(k, j);
217 M2.row_add(k, j);
218 M3.row_add(k, j);
219 }
220 M.row_op_end(0, r);
221 M2.row_op_end(0, r);
222 M3.row_op_end(0, r);
223
224 M.update_gso();
225 M2.update_gso();
226 M3.update_gso();
227 bool retvalue3 = rs_are_equal(M, M2, err);
228 bool retvalue4 = rs_are_equal(M3, M2, err);
229
230 return (!retvalue1) * 1 + (!retvalue2) * 2 + (!retvalue3) * 4 + (!retvalue4) * 8;
231 }
232
test_filename(const char * input_filename)233 template <class ZT, class FT> int test_filename(const char *input_filename)
234 {
235 ZZ_mat<ZT> A;
236 int status = read_file(A, input_filename);
237 // if status == 1, read_file fails.
238 if (status == 1)
239 {
240 return 1;
241 }
242 int retvalue = test_ggso<ZT, FT>(A);
243 if (retvalue & 1)
244 {
245 cerr
246 << input_filename
247 << " shows different GSO-outputs for grammatrix representation and basis representation.\n";
248 }
249 if (retvalue & 2)
250 {
251 cerr << input_filename
252 << " shows different GSO-outputs for grammatrix representation and "
253 "basis representation after moving rows.\n";
254 }
255 if (retvalue & 4)
256 {
257 cerr << input_filename
258 << " shows different GSO-outputs for grammatrix representation and "
259 "basis representation after adding rows.\n";
260 }
261 retvalue |= test_householder<ZT, FT>(A);
262 if (retvalue > 0)
263 {
264 return 1;
265 }
266 else
267 {
268 return 0;
269 }
270 }
271
272 /**
273 @brief Construct d × (d+1) integer relations matrix with bit size b and test LLL.
274
275 @param d dimension
276 @param b bit size
277 @param method LLL method to test
278 @param float_type floating point type to test
279 @param flags flags to use
280 @param prec precision used for is_lll_reduced
281
282 @return zero on success
283 */
284
test_int_rel(int d,int b)285 template <class ZT, class FT> int test_int_rel(int d, int b)
286 {
287 ZZ_mat<ZT> A;
288 A.resize(d, d + 1);
289 A.gen_intrel(b);
290 int retvalue = test_ggso<ZT, FT>(A);
291 if (retvalue >= 1)
292 {
293 cerr
294 << "Integer relation matrix with parameters " << d << " and " << b
295 << " shows different GSO-outputs for grammatrix representation and basis representation.\n";
296 return 1;
297 }
298 retvalue |= test_householder<ZT, FT>(A);
299 if (retvalue > 0)
300 return 1;
301
302 return 0;
303 }
304
main(int,char **)305 int main(int /*argc*/, char ** /*argv*/)
306 {
307
308 int status = 0;
309
310 status |= test_filename<mpz_t, double>(TESTDATADIR "/tests/lattices/example2_in");
311 status |= test_filename<mpz_t, double>(TESTDATADIR "/tests/lattices/example3_in");
312 status |= test_filename<mpz_t, double>(TESTDATADIR "/tests/lattices/example_cvp_in_lattice");
313 status |= test_filename<mpz_t, double>(TESTDATADIR "/tests/lattices/example_cvp_in_lattice2");
314 status |= test_filename<mpz_t, double>(TESTDATADIR "/tests/lattices/example_cvp_in_lattice3");
315 status |= test_filename<mpz_t, double>(TESTDATADIR "/tests/lattices/example_cvp_in_lattice4");
316 status |= test_filename<mpz_t, double>(TESTDATADIR "/tests/lattices/example_cvp_in_lattice5");
317 status |= test_int_rel<mpz_t, double>(50, 20);
318 status |= test_int_rel<mpz_t, double>(40, 10);
319
320 status |= test_filename<mpz_t, mpfr_t>(TESTDATADIR "/tests/lattices/example2_in");
321 status |= test_filename<mpz_t, mpfr_t>(TESTDATADIR "/tests/lattices/example3_in");
322 status |= test_filename<mpz_t, mpfr_t>(TESTDATADIR "/tests/lattices/example_cvp_in_lattice");
323 status |= test_filename<mpz_t, mpfr_t>(TESTDATADIR "/tests/lattices/example_cvp_in_lattice2");
324 status |= test_filename<mpz_t, mpfr_t>(TESTDATADIR "/tests/lattices/example_cvp_in_lattice3");
325 status |= test_filename<mpz_t, mpfr_t>(TESTDATADIR "/tests/lattices/example_cvp_in_lattice4");
326 status |= test_filename<mpz_t, mpfr_t>(TESTDATADIR "/tests/lattices/example_cvp_in_lattice5");
327 status |= test_int_rel<mpz_t, mpfr_t>(50, 20);
328 status |= test_int_rel<mpz_t, mpfr_t>(40, 10);
329
330 #ifdef FPLLL_WITH_LONG_DOUBLE
331 status |= test_filename<mpz_t, long double>(TESTDATADIR "/tests/lattices/example2_in");
332 status |= test_filename<mpz_t, long double>(TESTDATADIR "/tests/lattices/example3_in");
333 status |= test_filename<mpz_t, long double>(TESTDATADIR "/tests/lattices/example_cvp_in_lattice");
334 status |=
335 test_filename<mpz_t, long double>(TESTDATADIR "/tests/lattices/example_cvp_in_lattice2");
336 status |=
337 test_filename<mpz_t, long double>(TESTDATADIR "/tests/lattices/example_cvp_in_lattice3");
338 status |=
339 test_filename<mpz_t, long double>(TESTDATADIR "/tests/lattices/example_cvp_in_lattice4");
340 status |=
341 test_filename<mpz_t, long double>(TESTDATADIR "/tests/lattices/example_cvp_in_lattice5");
342 status |= test_int_rel<mpz_t, long double>(50, 20);
343 status |= test_int_rel<mpz_t, long double>(40, 10);
344 #endif
345 #ifdef FPLLL_WITH_QD
346 status |= test_filename<mpz_t, dd_real>(TESTDATADIR "/tests/lattices/example2_in");
347 status |= test_filename<mpz_t, dd_real>(TESTDATADIR "/tests/lattices/example3_in");
348 status |= test_filename<mpz_t, dd_real>(TESTDATADIR "/tests/lattices/example_cvp_in_lattice");
349 status |= test_filename<mpz_t, dd_real>(TESTDATADIR "/tests/lattices/example_cvp_in_lattice2");
350 status |= test_filename<mpz_t, dd_real>(TESTDATADIR "/tests/lattices/example_cvp_in_lattice3");
351 status |= test_filename<mpz_t, dd_real>(TESTDATADIR "/tests/lattices/example_cvp_in_lattice4");
352 status |= test_filename<mpz_t, dd_real>(TESTDATADIR "/tests/lattices/example_cvp_in_lattice5");
353 status |= test_int_rel<mpz_t, dd_real>(50, 20);
354 status |= test_int_rel<mpz_t, dd_real>(40, 10);
355
356 status |= test_filename<mpz_t, qd_real>(TESTDATADIR "/tests/lattices/example2_in");
357 status |= test_filename<mpz_t, qd_real>(TESTDATADIR "/tests/lattices/example3_in");
358 status |= test_filename<mpz_t, qd_real>(TESTDATADIR "/tests/lattices/example_cvp_in_lattice");
359 status |= test_filename<mpz_t, qd_real>(TESTDATADIR "/tests/lattices/example_cvp_in_lattice2");
360 status |= test_filename<mpz_t, qd_real>(TESTDATADIR "/tests/lattices/example_cvp_in_lattice3");
361 status |= test_filename<mpz_t, qd_real>(TESTDATADIR "/tests/lattices/example_cvp_in_lattice4");
362 status |= test_filename<mpz_t, qd_real>(TESTDATADIR "/tests/lattices/example_cvp_in_lattice5");
363 status |= test_int_rel<mpz_t, qd_real>(50, 20);
364 status |= test_int_rel<mpz_t, qd_real>(40, 10);
365 #endif
366 #ifdef FPLLL_WITH_DPE
367 status |= test_filename<mpz_t, dpe_t>(TESTDATADIR "/tests/lattices/example2_in");
368 status |= test_filename<mpz_t, dpe_t>(TESTDATADIR "/tests/lattices/example3_in");
369 status |= test_filename<mpz_t, dpe_t>(TESTDATADIR "/tests/lattices/example_cvp_in_lattice");
370 status |= test_filename<mpz_t, dpe_t>(TESTDATADIR "/tests/lattices/example_cvp_in_lattice2");
371 status |= test_filename<mpz_t, dpe_t>(TESTDATADIR "/tests/lattices/example_cvp_in_lattice3");
372 status |= test_filename<mpz_t, dpe_t>(TESTDATADIR "/tests/lattices/example_cvp_in_lattice4");
373 status |= test_filename<mpz_t, dpe_t>(TESTDATADIR "/tests/lattices/example_cvp_in_lattice5");
374 status |= test_int_rel<mpz_t, dpe_t>(50, 20);
375 status |= test_int_rel<mpz_t, dpe_t>(40, 10);
376 #endif
377
378 if (status == 0)
379 {
380 cerr << "All tests passed." << endl;
381 return 0;
382 }
383 else
384 {
385 return -1;
386 }
387
388 return 0;
389 }
390