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