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