1 /* Copyright (C) 2005-2008 Damien Stehle.
2    Copyright (C) 2007 David Cade.
3    Copyright (C) 2011 Xavier Pujol.
4    Copyright (C) 2019 Koen de Boer
5 
6    This file is part of fplll. fplll is free software: you
7    can redistribute it and/or modify it under the terms of the GNU Lesser
8    General Public License as published by the Free Software Foundation,
9    either version 2.1 of the License, or (at your option) any later version.
10 
11    fplll is distributed in the hope that it will be useful,
12    but WITHOUT ANY WARRANTY; without even the implied warranty of
13    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14    GNU Lesser General Public License for more details.
15 
16    You should have received a copy of the GNU Lesser General Public License
17    along with fplll. If not, see <http://www.gnu.org/licenses/>. */
18 
19 /* Template source file */
20 
21 #include "gso_interface.h"
22 
23 FPLLL_BEGIN_NAMESPACE
24 
25 template <class ZT, class FT>
invalidate_gso_row(int i,int new_valid_cols)26 inline void MatGSOInterface<ZT, FT>::invalidate_gso_row(int i, int new_valid_cols)
27 {
28   FPLLL_DEBUG_CHECK(i >= 0 && i < n_known_rows && new_valid_cols >= 0 && new_valid_cols <= i + 1);
29   gso_valid_cols[i] = min(gso_valid_cols[i], new_valid_cols);
30 }
31 
row_op_end(int first,int last)32 template <class ZT, class FT> void MatGSOInterface<ZT, FT>::row_op_end(int first, int last)
33 {
34 #ifdef DEBUG
35   FPLLL_DEBUG_CHECK(row_op_first == first && row_op_last == last);
36   row_op_first = row_op_last = -1;
37 #endif
38   for (int i = first; i < last; i++)
39   {
40     if (!enable_int_gram)
41     {
42       update_bf(i);
43       invalidate_gram_row(i);
44       for (int j = i + 1; j < n_known_rows; j++)
45         gf(j, i).set_nan();
46     }
47     invalidate_gso_row(i, 0);
48   }
49   for (int i = last; i < n_known_rows; i++)
50   {
51     invalidate_gso_row(i, first);
52   }
53 }
54 
get_max_gram()55 template <class ZT, class FT> inline ZT MatGSOInterface<ZT, FT>::get_max_gram()
56 {
57   ZT tmp;
58   if (enable_int_gram)
59   {
60     if (gptr == nullptr)
61     {
62       throw std::runtime_error("Error: gptr is equal to the nullpointer.");
63     }
64     Matrix<ZT> gr = *gptr;
65     tmp           = gr(0, 0);
66     for (int i = 0; i < n_known_rows; i++)
67       tmp = tmp.max_z(gr(i, i));
68   }
69   else
70   {
71     FT tmp1 = gf(0, 0);
72     for (int i = 0; i < n_known_rows; i++)
73       tmp1 = tmp1.max_f(gf(i, i));
74     tmp.set_f(tmp1);
75   }
76   return tmp;
77 }
78 
get_max_bstar()79 template <class ZT, class FT> inline FT MatGSOInterface<ZT, FT>::get_max_bstar()
80 {
81   FT tmp;
82   tmp = r(0, 0);
83   for (int i = 0; i < n_known_rows; i++)
84     tmp = tmp.max_f(r(i, i));
85   return tmp;
86 }
87 
get_max_mu_exp(int i,int n_columns)88 template <class ZT, class FT> long MatGSOInterface<ZT, FT>::get_max_mu_exp(int i, int n_columns)
89 {
90   FPLLL_DEBUG_CHECK(i >= 0 && i < n_known_rows && gso_valid_cols[i] >= n_columns);
91   long max_expo = LONG_MIN, expo;
92   for (int j = 0; j < n_columns; j++)
93   {
94     long expo2 = get_mu_exp(i, j, expo).exponent();
95     max_expo   = max(max_expo, expo + expo2);
96   }
97   return max_expo;
98 }
99 // This function *can* be included in the interface; it has no base-class specific functions.
100 // If this function is going to be included, it must be removed from gso.* and gsogram.*
101 // and also added to gso_interface.h as non-virtual.
102 /*
103 template <class ZT, class FT>
104 void MatGSOInterface<ZT, FT>::row_addmul_we(int i, int j, const FT &x, long expo_add)
105 {
106   FPLLL_DEBUG_CHECK(j >= 0 && i < n_known_rows && j < n_source_rows);
107   long expo;
108   long lx = x.get_si_exp_we(expo, expo_add);
109 
110   if (expo == 0)
111   {
112     if (lx == 1)
113       row_add(i, j);
114     else if (lx == -1)
115       row_sub(i, j);
116     else if (lx != 0)
117       row_addmul_si(i, j, lx);
118   }
119   else if (row_op_force_long)
120   {
121     row_addmul_si_2exp(i, j, lx, expo);
122   }
123   else
124   {
125     x.get_z_exp_we(ztmp2, expo, expo_add);
126     row_addmul_2exp(i, j, ztmp2, expo);
127   }
128 }
129 */
130 
update_gso_row(int i,int last_j)131 template <class ZT, class FT> bool MatGSOInterface<ZT, FT>::update_gso_row(int i, int last_j)
132 {
133   // FPLLL_TRACE_IN("Updating GSO up to (" << i << ", " << last_j << ")");
134   // FPLLL_TRACE("n_known_rows=" << n_known_rows << " n_source_rows=" << n_source_rows);
135   if (i >= n_known_rows)
136   {
137     discover_row();
138   }
139   FPLLL_DEBUG_CHECK(i >= 0 && i < n_known_rows && last_j >= 0 && last_j < n_source_rows);
140 
141   int j = max(0, gso_valid_cols[i]);
142 
143   for (; j <= last_j; j++)
144   {
145     get_gram(ftmp1, i, j);
146     FPLLL_DEBUG_CHECK(j == i || gso_valid_cols[j] >= j);
147     for (int k = 0; k < j; k++)
148     {
149       ftmp2.mul(mu(j, k), r(i, k));
150       ftmp1.sub(ftmp1, ftmp2);
151     }
152     r(i, j) = ftmp1;
153     if (i > j)
154     {
155       mu(i, j).div(ftmp1, r(j, j));
156       if (!mu(i, j).is_finite())
157         return false;
158     }
159   }
160 
161   gso_valid_cols[i] = j;  // = max(0, gso_valid_cols[i], last_j + 1)
162   // FPLLL_TRACE_OUT("End of GSO update");
163   return true;
164 }
165 
lock_cols()166 template <class ZT, class FT> void MatGSOInterface<ZT, FT>::lock_cols() { cols_locked = true; }
167 
unlock_cols()168 template <class ZT, class FT> void MatGSOInterface<ZT, FT>::unlock_cols()
169 {
170   n_known_rows = n_source_rows;
171   cols_locked  = false;
172 }
173 
174 template <class ZT, class FT>
apply_transform(const Matrix<FT> & transform,int src_base,int target_base)175 void MatGSOInterface<ZT, FT>::apply_transform(const Matrix<FT> &transform, int src_base,
176                                               int target_base)
177 {
178   int target_size = transform.get_rows(), src_size = transform.get_cols();
179   int old_d = d;
180   create_rows(target_size);
181   for (int i = 0; i < target_size; i++)
182   {
183     for (int j = 0; j < src_size; j++)
184     {
185       row_addmul(old_d + i, src_base + j, transform(i, j));
186     }
187   }
188   row_op_begin(target_base, target_base + target_size);
189   for (int i = 0; i < target_size; i++)
190   {
191     row_swap(target_base + i, old_d + i);
192   }
193   row_op_end(target_base, target_base + target_size);
194   remove_last_rows(target_size);
195 }
196 
197 template <class ZT, class FT>
get_current_slope(int start_row,int stop_row)198 double MatGSOInterface<ZT, FT>::get_current_slope(int start_row, int stop_row)
199 {
200   FT f, log_f;
201   long expo;
202   vector<double> x;
203   x.resize(stop_row);
204   for (int i = start_row; i < stop_row; i++)
205   {
206     update_gso_row(i);
207     f = get_r_exp(i, i, expo);
208     log_f.log(f, GMP_RNDU);
209     x[i] = log_f.get_d() + expo * std::log(2.0);
210   }
211   int n         = stop_row - start_row;
212   double i_mean = (n - 1) * 0.5 + start_row, x_mean = 0, v1 = 0, v2 = 0;
213   for (int i = start_row; i < stop_row; i++)
214   {
215     x_mean += x[i];
216   }
217   x_mean /= n;
218   for (int i = start_row; i < stop_row; i++)
219   {
220     v1 += (i - i_mean) * (x[i] - x_mean);
221     v2 += (i - i_mean) * (i - i_mean);
222   }
223   return v1 / v2;
224 }
225 
get_root_det(int start_row,int end_row)226 template <class ZT, class FT> FT MatGSOInterface<ZT, FT>::get_root_det(int start_row, int end_row)
227 {
228   start_row   = max(0, start_row);
229   end_row     = min(d, end_row);
230   FT h        = (double)(end_row - start_row);
231   FT root_det = get_log_det(start_row, end_row) / h;
232   root_det.exponential(root_det);
233   return root_det;
234 }
235 
get_log_det(int start_row,int end_row)236 template <class ZT, class FT> FT MatGSOInterface<ZT, FT>::get_log_det(int start_row, int end_row)
237 {
238   FT log_det = 0.0;
239   start_row  = max(0, start_row);
240   end_row    = min(d, end_row);
241   FT h;
242   for (int i = start_row; i < end_row; ++i)
243   {
244     get_r(h, i, i);
245     log_det += log(h);
246   }
247   return log_det;
248 }
249 
250 template <class ZT, class FT>
get_slide_potential(int start_row,int end_row,int block_size)251 FT MatGSOInterface<ZT, FT>::get_slide_potential(int start_row, int end_row, int block_size)
252 {
253   FT potential = 0.0;
254   int p        = (end_row - start_row) / block_size;
255   if ((end_row - start_row) % block_size == 0)
256   {
257     --p;
258   }
259   for (int i = 0; i < p; ++i)
260   {
261     potential += (p - i) * get_log_det(i * block_size, (i + 1) * block_size);
262   }
263   return potential;
264 }
265 
266 template <class FT>
adjust_radius_to_gh_bound(FT & max_dist,long max_dist_expo,int block_size,const FT & root_det,double gh_factor)267 void adjust_radius_to_gh_bound(FT &max_dist, long max_dist_expo, int block_size, const FT &root_det,
268                                double gh_factor)
269 {
270   double t = (double)block_size / 2.0 + 1;
271   t        = lgamma(t);
272   t        = pow(M_E, t * 2.0 / (double)block_size);
273   t        = t / M_PI;
274   FT f     = t;
275   f        = f * root_det;
276   f.mul_2si(f, -max_dist_expo);
277   f = f * gh_factor;
278   if (f < max_dist)
279   {
280     max_dist = f;
281   }
282 }
283 
284 template class MatGSOInterface<Z_NR<long>, FP_NR<double>>;
285 template class MatGSOInterface<Z_NR<double>, FP_NR<double>>;
286 template class MatGSOInterface<Z_NR<mpz_t>, FP_NR<double>>;
287 template void adjust_radius_to_gh_bound<FP_NR<double>>(FP_NR<double> &max_dist, long max_dist_expo,
288                                                        int block_size,
289                                                        const FP_NR<double> &root_det,
290                                                        double gh_factor);
291 
292 #ifdef FPLLL_WITH_LONG_DOUBLE
293 template class MatGSOInterface<Z_NR<long>, FP_NR<long double>>;
294 template class MatGSOInterface<Z_NR<double>, FP_NR<long double>>;
295 template class MatGSOInterface<Z_NR<mpz_t>, FP_NR<long double>>;
296 template void adjust_radius_to_gh_bound<FP_NR<long double>>(FP_NR<long double> &max_dist,
297                                                             long max_dist_expo, int block_size,
298                                                             const FP_NR<long double> &root_det,
299                                                             double gh_factor);
300 
301 #endif
302 
303 #ifdef FPLLL_WITH_QD
304 template class MatGSOInterface<Z_NR<long>, FP_NR<dd_real>>;
305 template class MatGSOInterface<Z_NR<double>, FP_NR<dd_real>>;
306 template class MatGSOInterface<Z_NR<mpz_t>, FP_NR<dd_real>>;
307 template void adjust_radius_to_gh_bound<FP_NR<dd_real>>(FP_NR<dd_real> &max_dist,
308                                                         long max_dist_expo, int block_size,
309                                                         const FP_NR<dd_real> &root_det,
310                                                         double gh_factor);
311 template class MatGSOInterface<Z_NR<long>, FP_NR<qd_real>>;
312 template class MatGSOInterface<Z_NR<double>, FP_NR<qd_real>>;
313 template class MatGSOInterface<Z_NR<mpz_t>, FP_NR<qd_real>>;
314 template void adjust_radius_to_gh_bound<FP_NR<qd_real>>(FP_NR<qd_real> &max_dist,
315                                                         long max_dist_expo, int block_size,
316                                                         const FP_NR<qd_real> &root_det,
317                                                         double gh_factor);
318 #endif
319 
320 #ifdef FPLLL_WITH_DPE
321 template class MatGSOInterface<Z_NR<long>, FP_NR<dpe_t>>;
322 template class MatGSOInterface<Z_NR<double>, FP_NR<dpe_t>>;
323 template class MatGSOInterface<Z_NR<mpz_t>, FP_NR<dpe_t>>;
324 template void adjust_radius_to_gh_bound<FP_NR<dpe_t>>(FP_NR<dpe_t> &max_dist, long max_dist_expo,
325                                                       int block_size, const FP_NR<dpe_t> &root_det,
326                                                       double gh_factor);
327 #endif
328 
329 template class MatGSOInterface<Z_NR<long>, FP_NR<mpfr_t>>;
330 template class MatGSOInterface<Z_NR<double>, FP_NR<mpfr_t>>;
331 template class MatGSOInterface<Z_NR<mpz_t>, FP_NR<mpfr_t>>;
332 template void adjust_radius_to_gh_bound<FP_NR<mpfr_t>>(FP_NR<mpfr_t> &max_dist, long max_dist_expo,
333                                                        int block_size,
334                                                        const FP_NR<mpfr_t> &root_det,
335                                                        double gh_factor);
336 
337 FPLLL_END_NAMESPACE
338