1 /* Copyright (C) 2011 Xavier Pujol.
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 "util.h"
17
18 #ifdef DEBUG
19 int debug_depth = 0;
20 #endif
21
22 FPLLL_BEGIN_NAMESPACE
23
24 enum MinPrecAlgo
25 {
26 MINPREC_GSO,
27 MINPREC_L2
28 };
29
30 /* State of LDConvHelper (declared in nr.h, must be defined in exactly one
31 source file) */
32 #ifdef FPLLL_WITH_LONG_DOUBLE
33 mpfr_t LDConvHelper::temp;
34 bool LDConvHelper::temp_initialized = false;
35 #endif
36
37 /* State of the random generator (declared in nr.h, must be defined in exactly
38 one source file) */
39 bool RandGen::initialized = false;
40 gmp_randstate_t RandGen::gmp_state;
41
compute_min_prec(double & rho,int d,double delta,double eta,double epsilon,MinPrecAlgo algo)42 static int compute_min_prec(double &rho, int d, double delta, double eta, double epsilon,
43 MinPrecAlgo algo)
44 {
45 int old_prec = FP_NR<mpfr_t>::set_prec(53);
46 FP_NR<mpfr_t> f_minprec, f_rho, f_d, f_eta, f_delta, f_epsilon, tmp1, tmp2;
47
48 // These four conversions are exact
49 f_d = static_cast<double>(d);
50 f_eta = eta;
51 f_delta = delta;
52 f_epsilon = epsilon;
53 if (algo == MINPREC_L2)
54 {
55 // eta - 0.5 is an exact fp operation
56 if (f_epsilon > eta - 0.5)
57 f_epsilon = eta - 0.5;
58 tmp1 = 1.0;
59 tmp1.sub(tmp1, f_delta, GMP_RNDD);
60 if (f_epsilon > tmp1)
61 f_epsilon = tmp1;
62 // now fEpsilon <= min(epsilon, eta - 0.5, 1 - delta);
63 }
64 // Computes tmp1 >= (1 + eta) ^ 2 + epsilon
65 tmp1 = 1.0; // exact
66 tmp1.add(f_eta, tmp1, GMP_RNDU); // >= 1 + eta
67 tmp1.mul(tmp1, tmp1, GMP_RNDU); // >= (1 + eta) ^ 2
68 tmp1.add(tmp1, f_epsilon, GMP_RNDU);
69 // Computes tmp2 <= delta - eta ^ 2
70 tmp2.mul(f_eta, f_eta, GMP_RNDU);
71 tmp2.sub(f_delta, tmp2, GMP_RNDD);
72 FPLLL_CHECK(tmp2 > 0, "invalid LLL parameters, eta must be < sqrt(delta)");
73 // Computes rho >= ((1 + eta) ^ 2 + epsilon) / (delta - eta ^ 2)
74 f_rho.div(tmp1, tmp2, GMP_RNDU);
75 rho = f_rho.get_d(GMP_RNDU);
76
77 /* Computes minprec >= constant + 2 * log2(d) - log2(epsilon) + d * log2(rho)
78 (constant = 5 for GSO, 10 for LLL) */
79 tmp1.log(f_d, GMP_RNDU); // >= log(d)
80 tmp1.mul_2si(tmp1, 1); // >= 2 * log(d)
81 tmp2.log(f_epsilon, GMP_RNDD); // <= log(epsilon) (<= 0)
82 tmp1.sub(tmp1, tmp2, GMP_RNDU); // >= 2 * log(d) - log(epsilon)
83 tmp2.log(f_rho, GMP_RNDU); // >= log(rho)
84 tmp2.mul(f_d, tmp2, GMP_RNDU); // >= d * log(rho)
85 tmp1.add(tmp1, tmp2, GMP_RNDU); // >= 2*log(d)-log(epsilon)+d*log(rho)
86 tmp2 = 2.0; // exact
87 tmp2.log(tmp2, GMP_RNDD); // <= log(2)
88 tmp1.div(tmp1, tmp2, GMP_RNDU); // >= 2*log2(d)-log2(epsilon)+d*log2(rho)
89 tmp2 = (algo == MINPREC_L2) ? 10.0 : 5.0;
90 f_minprec.add(tmp1, tmp2, GMP_RNDU);
91 int minprec = static_cast<int>(ceil(f_minprec.get_d(GMP_RNDU)));
92 mpfr_free_cache();
93 FP_NR<mpfr_t>::set_prec(old_prec);
94 return minprec;
95 }
96
gso_min_prec(double & rho,int d,double delta,double eta,double epsilon)97 int gso_min_prec(double &rho, int d, double delta, double eta, double epsilon)
98 {
99 return compute_min_prec(rho, d, delta, eta, epsilon, MINPREC_GSO);
100 }
101
l2_min_prec(int d,double delta,double eta,double epsilon)102 int l2_min_prec(int d, double delta, double eta, double epsilon)
103 {
104 double rho;
105 return compute_min_prec(rho, d, delta, eta, epsilon, MINPREC_L2);
106 }
107
hlll_min_prec(int d_i,int n_i,double delta_d,double eta_d,double theta_d,double c_d)108 int hlll_min_prec(int d_i, int n_i, double delta_d, double eta_d, double theta_d, double c_d)
109 {
110 FPLLL_CHECK(delta_d < 1.0 && delta_d >= 0.25, "delta must be in [1/4, 1).");
111 FPLLL_CHECK(theta_d >= 0.0, "theta must be positive.");
112 FPLLL_CHECK(eta_d >= 0.5, "theta must be larger than or equal to 0.5.");
113 FPLLL_CHECK(eta_d - theta_d > 0.5, "eta - theta must be larger than 0.5.");
114
115 int old_prec = FP_NR<mpfr_t>::set_prec(53);
116 FP_NR<mpfr_t> d, n, delta, eta, theta, c, alpha, c0, c1, rho, phi, p0, p;
117 FP_NR<mpfr_t> ftmp0, ftmp1, ftmp2, ftmp3, ftmp4;
118
119 d = d_i;
120 n = n_i;
121 delta = delta_d;
122 eta = eta_d;
123 theta = theta_d;
124 c = c_d;
125
126 // ftmp0 = (1 + theta^2) * delta - eta^2
127 ftmp0 = (1.0 + theta * theta) * delta - eta * eta;
128 // ftmp0 = sqrt((1 + theta^2) * delta - eta^2)
129 ftmp0.sqrt(ftmp0);
130
131 // alpha = theta * eta + sqrt((1 + theta^2) * delta - eta^2) / (delta - eta^2)
132 alpha = (theta * eta + ftmp0) / (delta - eta * eta);
133
134 // ftmp0 = 3 / 2
135 ftmp0 = 3.0 / 2.0;
136 // ftmp0 = sqrt(3 / 2)
137 ftmp0.sqrt(ftmp0);
138 // ftmp1 = 1 - eta - theta
139 ftmp1 = 1.0 - eta - theta;
140 // ftmp1 = |1 - eta - theta|
141 ftmp1.abs(ftmp1);
142 ftmp2 = 6.0;
143 // ftmp2 = sqrt(6)
144 ftmp2.sqrt(ftmp2);
145 // ftmp3 = 1 + d * eta^2
146 ftmp3 = 1.0 + d * eta * eta;
147 // ftmp3 = sqrt(1 + d * eta^2)
148 ftmp3.sqrt(ftmp3);
149 // ftmp4 = sqrt(d)
150 ftmp4.sqrt(d);
151 // ftmp0 = 1 + |1 - eta - theta| * alpha / ((eta + theta) * (-1 + sqrt(3/2)))
152 ftmp0 = (1.0 + ftmp1 * alpha) / ((eta + theta) * (-1.0 + ftmp0));
153 // ftmp1 = 4 * sqrt(6) / (1 + eta) * sqrt(1 + d * eta^2)
154 ftmp1 = 4.0 * ftmp2 / (1.0 + eta) * ftmp3;
155 // ftmp0 = max(1 + |1 - eta - theta| * alpha / ((eta + theta) * (-1 + sqrt(3/2))),
156 // 4 * sqrt(6) / (1 + eta) * sqrt(1 + d * eta^2))
157 ftmp0.max_f(ftmp1);
158 // c0 = max(...) * n * sqrt(d)
159 c0 = ftmp0 * n * ftmp4;
160
161 // c1 = 8 * d * (n + 9) * c0
162 c1 = 8.0 * d * (n + 9.0) * c0;
163
164 // rho = (1 + eta + theta) * alpha
165 rho = (1.0 + eta + theta) * alpha;
166
167 // ftmp0 = rho^(d + 1) (since we want to compute phi(d))
168 ftmp0.pow_si(rho, d_i + 1);
169 // phi(d) = c1 * (1 + 1 / theta) * ftmp0
170 phi = c1 * (1.0 + 1.0 / theta) * ftmp0;
171
172 // ftmp0 = alpha^d
173 ftmp0.pow_si(alpha, d_i);
174 // ftmp0 = log(d^3 * phi(d) * alpha^d / theta)
175 ftmp0.log(d * d * d * phi * ftmp0 / theta);
176 // ftmp1 = log(2)
177 ftmp1.log(2);
178 // ftmp0 = log(d^3 * phi(d) * alpha^d / theta) / log(2)
179 ftmp0 = ftmp0 / ftmp1;
180 // p0 = log2(d^3 * phi(d) * alpha^d / theta) + 16 + c * d / 2
181 p0 = ftmp0 + 16.0 + c * d / 2.0;
182
183 // ftmp0 = log(1 - delta)
184 ftmp0.log(1.0 - delta);
185 // ftmp0 = log(1 - delta) / log(2)
186 ftmp0 = ftmp0 / ftmp1;
187 // ftmp2 = log(eta - theta - 1/2)
188 ftmp2.log(eta - theta - 0.5);
189 // ftmp2 = log(eta - theta - 1/2) / log(2)
190 ftmp2 = ftmp2 / ftmp1;
191 // p = p0 + 1 - log2(1 - delta) - log2(eta - theta - 1 / 2)
192 p = p0 + 1.0 - ftmp0 - ftmp2;
193
194 // Convert p in int
195 int p_i = static_cast<int>(ceil(p.get_d(GMP_RNDU)));
196
197 FP_NR<mpfr_t>::set_prec(old_prec);
198
199 return p_i;
200 }
201
202 /**
203 * Computes the volume of a d-dimensional hypersphere of radius 1.
204 */
sphere_volume(FP_NR<mpfr_t> & volume,int d)205 void sphere_volume(FP_NR<mpfr_t> &volume, int d)
206 {
207 FP_NR<mpfr_t> rtmp1;
208 volume = pow(M_PI, (double)(d / 2));
209
210 if (d % 2 == 0)
211 for (int i = 1; i <= d / 2; i++)
212 {
213 rtmp1 = (double)i;
214 volume.div(volume, rtmp1);
215 }
216 else
217 for (int i = 0; i <= d / 2; i++)
218 {
219 rtmp1 = 2.0 / (double)(2 * i + 1);
220 volume.mul(volume, rtmp1);
221 }
222 }
223
224 /**
225 * Estimates the cost of the enumeration for SVP.
226 */
cost_estimate(FP_NR<mpfr_t> & cost,const FP_NR<mpfr_t> & bound,const Matrix<FP_NR<mpfr_t>> & r,int dimMax)227 void cost_estimate(FP_NR<mpfr_t> &cost, const FP_NR<mpfr_t> &bound, const Matrix<FP_NR<mpfr_t>> &r,
228 int dimMax)
229 {
230 FP_NR<mpfr_t> det, level_cost, tmp1;
231 det = 1.0;
232 cost = 0.0;
233
234 for (int i = dimMax - 1; i >= 0; i--)
235 {
236 tmp1.div(bound, r(i, i));
237 det.mul(det, tmp1);
238
239 level_cost.sqrt(det);
240 sphere_volume(tmp1, dimMax - i);
241 level_cost.mul(level_cost, tmp1);
242
243 cost.add(cost, level_cost);
244 }
245 }
246
get_red_status_str(int status)247 const char *get_red_status_str(int status)
248 {
249 if (status >= 0 && status < RED_STATUS_MAX)
250 return RED_STATUS_STR[status];
251 else
252 return "unknown error";
253 }
254
zeros_first(ZZ_mat<ZT> & b,ZZ_mat<ZT> & u,ZZ_mat<ZT> & u_inv_t)255 template <class ZT> void zeros_first(ZZ_mat<ZT> &b, ZZ_mat<ZT> &u, ZZ_mat<ZT> &u_inv_t)
256 {
257 int i, d = b.get_rows();
258 for (i = d; i > 0 && b[i - 1].is_zero(); i--)
259 {
260 }
261 if (i > 0 && i < d)
262 {
263 b.rotate(0, i, d - 1);
264 if (!u.empty())
265 u.rotate(0, i, d - 1);
266 if (!u_inv_t.empty())
267 u_inv_t.rotate(0, i, d - 1);
268 }
269 }
270
zeros_last(ZZ_mat<ZT> & b,ZZ_mat<ZT> & u,ZZ_mat<ZT> & u_inv_t)271 template <class ZT> void zeros_last(ZZ_mat<ZT> &b, ZZ_mat<ZT> &u, ZZ_mat<ZT> &u_inv_t)
272 {
273 int i, d = b.get_rows();
274 for (i = 0; i < d && b[i].is_zero(); i++)
275 {
276 }
277 if (i > 0 && i < d)
278 {
279 b.rotate(0, i, d - 1);
280 if (!u.empty())
281 u.rotate(0, i, d - 1);
282 if (!u_inv_t.empty())
283 u_inv_t.rotate(0, i, d - 1);
284 }
285 }
286
287 template void zeros_first<mpz_t>(ZZ_mat<mpz_t> &, ZZ_mat<mpz_t> &, ZZ_mat<mpz_t> &);
288 template void zeros_last<mpz_t>(ZZ_mat<mpz_t> &, ZZ_mat<mpz_t> &, ZZ_mat<mpz_t> &);
289
290 #ifdef FPLLL_WITH_ZLONG
291 template void zeros_first<long>(ZZ_mat<long> &, ZZ_mat<long> &, ZZ_mat<long> &);
292 template void zeros_last<long>(ZZ_mat<long> &, ZZ_mat<long> &, ZZ_mat<long> &);
293 #endif
294
295 #ifdef FPLLL_WITH_ZDOUBLE
296 template void zeros_first<double>(ZZ_mat<double> &, ZZ_mat<double> &, ZZ_mat<double> &);
297 template void zeros_last<double>(ZZ_mat<double> &, ZZ_mat<double> &, ZZ_mat<double> &);
298 #endif
299
300 FPLLL_END_NAMESPACE
301