1 /* Copyright (C) 2008-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 "evaluator.h"
17 
18 FPLLL_BEGIN_NAMESPACE
19 
init_delta_def(int prec,double rho,bool withRoundingToEnumf)20 void ErrorBoundedEvaluator::init_delta_def(int prec, double rho, bool withRoundingToEnumf)
21 {
22   /* Computes error bounds on GSO
23       For all 0 <= i < d and 0 <= j <= i we have:
24       |r~_i - r_i| / r_i <= maxRelDR = d * rho ^ (i + 1) * 2 ^ (2 - prec)
25       |mu~_(i,j) - mu_(i,j)| <= d * rho ^ (i + 1) * 2 ^ (4 - prec)
26       The following formula is used to bound absolute error on r_i:
27       |r~_i - r_i| <= r_i * maxRelDR <= (r~_i + |r~_i - r_i| * maxRelDR)
28       ==> |r~_i - r_i| <= r~_i * maxRelDR / (1 - maxRelDR)  */
29   FP_NR<mpfr_t> rhoPow, maxRelDR, ftmp1;
30   input_error_defined = true;
31   for (int j = 0; j < d; j++)
32   {
33     rhoPow = rho;
34     rhoPow.pow_si(rhoPow, j + 1, GMP_RNDU);  // >= rho ^ (j + 1)
35     ftmp1 = d;
36     ftmp1.mul_2si(ftmp1, 2 - prec);         // = d * 2 ^ (2 - prec)
37     maxRelDR.mul(ftmp1, rhoPow, GMP_RNDU);  // >= |r~_j - r_j| / r_j
38     ftmp1 = 1.0;
39     ftmp1.sub(ftmp1, maxRelDR, GMP_RNDD);  // <= 1 - maxRelDR
40     ftmp1.div(maxRelDR, ftmp1, GMP_RNDU);  // >= maxRelDR/(1-maxRelDR)
41     ftmp1.mul(ftmp1, r(j, j));             // >= |r~_j - r_j|
42     max_dr_diag[j] = ftmp1;
43 
44     ftmp1 = d;
45     ftmp1.mul_2si(ftmp1, 4 - prec);      // = d * 2 ^ (4 - prec)
46     ftmp1.mul(ftmp1, rhoPow, GMP_RNDU);  // >= |mu~_(?,j) - mu_(?,j)|
47     max_dm_u[j] = ftmp1;
48   }
49   if (withRoundingToEnumf)
50   {
51     FP_NR<mpfr_t> halfULP;
52     halfULP = numeric_limits<double>::epsilon() * 0.5;
53     for (int j = 0; j < d; j++)
54     {
55       max_dr_diag[j].addmul(r(j, j), halfULP, GMP_RNDU);  // >= |double(r~_j) - r_j|
56       max_dm_u[j].add(max_dm_u[j], halfULP, GMP_RNDU);
57     }
58   }
59 }
60 
61 /* Main function for error evaluation
62    Returns maxDE such that in each loop of the enumeration, if
63      exact_squared_norm(vector) <= max_dist (when boundOnExactVal = true)
64    or
65      computed_squared_norm(vector) <= max_dist (when boundOnExactVal = false)
66    then
67      |exact_squared_norm(vector) - computed_squared_norm(vector)| <= maxDE
68    (computed_squared_norm(vector) is the variable 'new_dist' is enumerateLoop)
69 
70    In the following comments:
71    - The symbols +~ and *~ represent rounded fp addition and multiplication
72    - a,b,c,... represent exact values
73    - a~,b~,c~,... are approx. values used or computed by the fp algorithm */
74 
get_max_error_aux(const FP_NR<mpfr_t> & max_dist,bool boundOnExactVal,FP_NR<mpfr_t> & maxDE)75 bool ErrorBoundedEvaluator::get_max_error_aux(const FP_NR<mpfr_t> &max_dist, bool boundOnExactVal,
76                                               FP_NR<mpfr_t> &maxDE)
77 {
78 
79   FPLLL_CHECK(input_error_defined,
80               "Evaluator: error evaluation failed because the input error is undefined");
81 
82   FP_NR<mpfr_t> ulp, halfULP, K, tmp1, tmp2;
83   FP_NR<mpfr_t> rdiagTilde, minRDiag, maxRDiag, muTilde, maxMu, maxMuTildeX;
84   FP_NR<mpfr_t> maxC, maxCTilde, maxY, maxYTilde, maxY2Tilde, maxRY2Tilde;
85   FP_NR<mpfr_t> maxDC, maxDY, maxDY2, maxDRY2;
86   vector<FP_NR<mpfr_t>> maxX(d);
87 
88   ulp = numeric_limits<double>::epsilon();
89   halfULP.mul_2si(ulp, -1);
90   K = 1.0;
91   K.add(K, halfULP, GMP_RNDU);
92   maxDE = 0.0;
93 
94   for (int i = d - 1; i >= 0; i--)
95   {
96     maxC      = 0.0;
97     maxCTilde = 0.0;
98     maxDC     = 0.0;
99 
100     long rdiagExp = r(i, i).exponent();
101     tmp1.mul_2si(r(i, i), -rdiagExp);
102     tmp1 = tmp1.get_d();
103     rdiagTilde.mul_2si(tmp1, rdiagExp);  // = r~_i
104 
105     /* Computes bounds on:
106        C = mu(d-1,i) * x_(d-1) + ... + mu(j,i) * x_j
107        C~ = mu~(d-1,i) *~ x_(d-1) +~ ... +~ mu~(j,i) *~ x_j
108        DC = |C - C~|  */
109     for (int j = d - 1; j > i; j--)
110     {
111       maxMu.abs(mu(j, i));
112       maxMu.add(maxMu, max_dm_u[i], GMP_RNDU);  // >= |mu(j,i)|
113       maxC.addmul(maxMu, maxX[j], GMP_RNDU);
114       // now maxC >= |mu(d-1,i)| * |x_(d-1)| + ... + |mu(j,i)| * |x_j|
115       muTilde = fabs(mu(j, i).get_d());              // = |mu~(j,i)|
116       maxMuTildeX.mul(muTilde, maxX[j], GMP_RNDU);   // >= mu~(j,i) * x_j
117       maxDC.addmul(max_dm_u[i], maxX[j]);            // err1: initial error on mu
118       maxDC.addmul(maxMuTildeX, halfULP, GMP_RNDU);  // err2: rounding after *
119       maxMuTildeX.mul(maxMuTildeX, K, GMP_RNDU);     // >= mu~(j,i) *~ x_j
120       maxCTilde.addmul(maxMuTildeX, K, GMP_RNDU);
121       // now maxCTilde >= |mu(d-1,i)| *~ |x_(d-1)| +~ ... + |mu(j,i)| *~ |x_j|
122       maxDC.addmul(maxCTilde, halfULP, GMP_RNDU);  // err3: rounding after +
123       maxCTilde.mul(maxCTilde, K, GMP_RNDU);
124       // now maxCTilde >= |mu(d-1,i)| *~ |x_(d-1)| +~ ... +~ |mu(j,i)| *~ |x_j|
125     }
126 
127     if (boundOnExactVal)
128     {
129       // We have dist <= max_dist
130       minRDiag.sub(r(i, i), max_dr_diag[i], GMP_RNDD);  // <= r_i
131       if (minRDiag <= 0.0)
132         return false;
133       tmp1.div(max_dist, minRDiag, GMP_RNDU);  // >= dist / r_i
134       maxY.sqrt(tmp1, GMP_RNDU);               // >= |y_i|
135                                                /* DY = |y_i - y~_i|
136                                                      <= DC + |x_i - C~| * halfULP
137                                                      <= DC + (|x_i - C| + DC) * halfULP
138                                                      = DC * K + y * halfULP  */
139       maxDY.mul(maxY, halfULP, GMP_RNDU);
140       maxDY.addmul(maxDC, K, GMP_RNDU);      // >= |y_i - y~_i|
141       maxYTilde.add(maxY, maxDY, GMP_RNDU);  // >= |y~_i|
142       tmp1.add(maxY, maxC, GMP_RNDD);        // >= |x_i|
143       maxX[i].floor(tmp1);                   // x_i is an integer
144       tmp1 = maxY;
145     }
146     else
147     {
148       // We have dist~ <= max_dist
149       /* dist~ >= (y~_i *~ y~_i) *~ r~_i
150                >= (y~_i *~ y~_i) * r~_i - dist~ * halfULP
151          ==> y~_i *~ y~_i <= dist~ * K / r~_i  */
152       tmp1.mul(max_dist, K, GMP_RNDU);
153       tmp1.div(tmp1, rdiagTilde, GMP_RNDU);  // >= y~_i *~ y~_i
154                                              /* tmp1 >= y~_i *~ y~_i >= y~_i * y~_i - tmp1 * halfULP
155                                                 ==> y~_i <= sqrt(tmp1 * K)   */
156       tmp1.mul(tmp1, K, GMP_RNDU);
157       maxYTilde.sqrt(tmp1, GMP_RNDU);           // >= y~_i
158       maxDY.mul(maxYTilde, halfULP, GMP_RNDU);  // err2: rounding after +
159       maxDY.add(maxDY, maxDC, GMP_RNDU);        // err1: initial error on C
160                                                 // now maxDY >= |y_i - y~_i|
161       /* maxYTilde >= |C~_i -~ x_i| >= |C~_i - x_i| - maxYTilde * halfULP
162          ==> |x_i| <= C~_i + maxYTilde * K  */
163       tmp1 = maxCTilde;
164       tmp1.addmul(maxYTilde, K, GMP_RNDD);  // >= |x_i|
165       maxX[i].floor(tmp1);                  // x_i is an integer
166       tmp1 = maxYTilde;
167     }
168 
169     maxDY2.mul(maxDY, tmp1);
170     maxDY2.mul_2si(maxDY2, 1);
171     maxDY2.addmul(maxDY, maxDY, GMP_RNDU);
172     /* now, we have maxDY2 <= |y~_i * y~_i - y_i * y_i|
173        Case 1: tmp1 = maxY
174          |y~_i * y~_i - y_i * y_i| <= (maxY + maxDY) ^ 2 - maxY ^ 2
175        Case 2: tmp1 = maxY~
176          |y~_i * y~_i - y_i * y_i| <= (maxY~ + maxDY) ^ 2 - maxY~ ^ 2  */
177     maxY2Tilde.mul(maxYTilde, maxYTilde, GMP_RNDU);  // >= y~_i * y~_i
178     maxDY2.addmul(maxY2Tilde, halfULP, GMP_RNDU);
179     /* now, we have maxDY2 <= |y~_i *~ y~_i - y_i * y_i| */
180     maxY2Tilde.mul(maxY2Tilde, K, GMP_RNDU);            // >= y~_i *~ y~_i
181     maxRDiag.add(r(i, i), max_dr_diag[i], GMP_RNDU);    // >= r_i
182     maxRY2Tilde.mul(rdiagTilde, maxY2Tilde, GMP_RNDU);  // >= y~_i *~ y~_i * r~_i
183 
184     maxDRY2.mul(maxRDiag, maxDY2, GMP_RNDU);
185     maxDRY2.addmul(maxY2Tilde, max_dr_diag[i], GMP_RNDU);
186     maxDRY2.addmul(maxRY2Tilde, halfULP, GMP_RNDU);
187     /* now maxDRY2 >= |y~_i *~ y~_i - y_i * y_i| * r_i +
188                       |y~_i *~ y~_i| * |r~_i - r_i| +
189                       y~_i *~ y~_i * r~_i * halfULP
190                    >= |y~_i *~ y~_i * r~_i - y_i * y_i * r_i| +
191                       y~_i *~ y~_i * r~_i * halfULP  */
192 
193     /* Let u = dist_(i+1), u~ = dist~_(i+1), v = y_i * y_i * r_i and
194        v~ = y~_i *~ y~_i *~ r~_i. We have:
195          |dist_i - dist~_i| = |(u~ +~ v~) - (u + v)|
196                             <= |u~ +~ v~) - (u~ + v~)| + |(u~ + v~) - (u + v)|
197        Case 1: dist <= max_dist
198          |dist_i - dist~_i| <= (u~ + v~) * halfULP + |(u~ + v~) - (u + v)|
199                             <= (u + v) * halfULP + K * |(u~ + v~) - (u + v)|
200                             <= max_dist * halfULP + K * |(u~ + v~) - (u + v)|
201        Case 2: dist~ <= max_dist
202          |dist_i - dist~_i| <= (u~ +~ v~) * halfULP + |(u~ + v~) - (u + v)|
203                             <= (u~ +~ v~) * halfULP + K * |(u~ + v~) - (u + v)|
204                             <= max_dist * halfULP + K * |(u~ + v~) - (u + v)|  */
205     maxDE.add(maxDE, maxDRY2, GMP_RNDU);
206     maxDE.mul(maxDE, K, GMP_RNDU);
207     maxDE.addmul(max_dist, halfULP, GMP_RNDU);
208   }
209 
210   return true;
211 }
212 
eval_sol(const vector<FP_NR<mpfr_t>> & new_sol_coord,const enumf & new_partial_dist,enumf & max_dist)213 void FastErrorBoundedEvaluator::eval_sol(const vector<FP_NR<mpfr_t>> &new_sol_coord,
214                                          const enumf &new_partial_dist, enumf &max_dist)
215 {
216   // Assumes that the solution is valid
217   if (eval_mode == EVALMODE_SV)
218   {
219     FP_NR<mpfr_t> dist(new_partial_dist);
220     dist.mul_2si(dist, normExp);
221     this->process_sol(dist, new_sol_coord, max_dist);
222   }
223   else if (eval_mode == EVALMODE_PRINT)
224   {
225     std::cout << new_sol_coord << "\n";
226   }
227 }
228 
eval_sub_sol(int offset,const vector<FP_NR<mpfr_t>> & new_sub_sol_coord,const enumf & sub_dist)229 void FastErrorBoundedEvaluator::eval_sub_sol(int offset,
230                                              const vector<FP_NR<mpfr_t>> &new_sub_sol_coord,
231                                              const enumf &sub_dist)
232 {
233   FP_NR<mpfr_t> dist = sub_dist;
234   dist.mul_2si(dist, normExp);
235 
236   sub_solutions.resize(std::max(sub_solutions.size(), std::size_t(offset + 1)));
237 
238   if (sub_solutions[offset].second.empty() || dist < sub_solutions[offset].first)
239   {
240     sub_solutions[offset].first  = dist;
241     sub_solutions[offset].second = new_sub_sol_coord;
242     for (int i = 0; i < offset; ++i)
243       sub_solutions[offset].second[i] = 0.0;
244   }
245 }
246 
get_max_error(FP_NR<mpfr_t> & max_error,const FP_NR<mpfr_t> & last_partial_dist)247 bool FastErrorBoundedEvaluator::get_max_error(FP_NR<mpfr_t> &max_error,
248                                               const FP_NR<mpfr_t> &last_partial_dist)
249 {
250   FP_NR<mpfr_t> maxE, maxDE, maxOptDE, minOptE, one;
251 
252   if (solutions.empty() || !input_error_defined)
253     return false;
254   if (!get_max_error_aux(last_partial_dist, false, maxDE))
255     return false;
256 
257   // Exact norm of an optimal solution <= Exact norm of the result <= maxE
258   maxE.add(last_partial_dist, maxDE, GMP_RNDU);
259   // Error on the norm of an optimal solution <= maxOptDE
260   if (!get_max_error_aux(maxE, true, maxOptDE))
261     return false;
262   // Approx. norm of an optimal solution >= minOptE
263   minOptE.sub(last_partial_dist, maxOptDE, GMP_RNDD);
264 
265   one = 1.0;
266   max_error.div(maxE, minOptE, GMP_RNDU);
267   max_error.sub(max_error, one, GMP_RNDU);
268   return true;
269 }
270 
get_max_error(FP_NR<mpfr_t> & max_error,const FP_NR<mpfr_t> &)271 bool ExactErrorBoundedEvaluator::get_max_error(FP_NR<mpfr_t> &max_error, const FP_NR<mpfr_t> &)
272 {
273   max_error = 0.0;
274   return true;
275 }
276 
eval_sol(const vector<FP_NR<mpfr_t>> & new_sol_coord,const enumf &,enumf & max_dist)277 void ExactErrorBoundedEvaluator::eval_sol(const vector<FP_NR<mpfr_t>> &new_sol_coord, const enumf &,
278                                           enumf &max_dist)
279 {
280   int n = gso.get_cols_of_b();
281 
282   Z_NR<mpz_t> new_sol_dist;
283   vector<Z_NR<mpz_t>> new_sol, coord;
284 
285   gen_zero_vect(new_sol, n);
286   gen_zero_vect(coord, n);
287   new_sol_dist = 0;
288 
289   // Computes the distance between x and zero
290   for (int i = 0; i < d; i++)
291   {
292     coord[i].set_f(new_sol_coord[i]);
293   }
294   gso.sqnorm_coordinates(new_sol_dist, coord);
295 
296   if (int_max_dist < 0 || new_sol_dist <= int_max_dist)
297   {
298     if (eval_mode == EVALMODE_SV)
299     {
300       int_max_dist = new_sol_dist;
301 
302       this->process_sol(int_dist2Float(int_max_dist), new_sol_coord, max_dist);
303     }
304     else if (eval_mode == EVALMODE_PRINT)
305     {
306       cout << new_sol_coord << "\n";
307     }
308   }
309 }
310 
eval_sub_sol(int offset,const vector<FP_NR<mpfr_t>> & new_sub_sol_coord,const enumf &)311 void ExactErrorBoundedEvaluator::eval_sub_sol(int offset,
312                                               const vector<FP_NR<mpfr_t>> &new_sub_sol_coord,
313                                               const enumf &)
314 {
315   Z_NR<mpz_t> minusone;
316   minusone = -1;
317 
318   int n = gso.get_cols_of_b();
319   Z_NR<mpz_t> new_sol_dist;
320   vector<Z_NR<mpz_t>> new_sol, coord;
321 
322   gen_zero_vect(new_sol, n);
323   gen_zero_vect(coord, n);
324   new_sol_dist = 0;
325 
326   // Computes the distance between x[[offset,d)] and zero
327   for (int i = offset; i < d; i++)
328   {
329     coord[i].set_f(new_sub_sol_coord[i]);
330   }
331   gso.sqnorm_coordinates(new_sol_dist, coord);
332 
333   FP_NR<mpfr_t> subdist = int_dist2Float(new_sol_dist);
334 
335   sub_solutions.resize(std::max(sub_solutions.size(), std::size_t(offset + 1)));
336   if (sub_solutions[offset].second.empty() || subdist <= sub_solutions[offset].first)
337   {
338     sub_solutions[offset].first  = subdist;
339     sub_solutions[offset].second = new_sub_sol_coord;
340     for (int i = 0; i < offset; ++i)
341       sub_solutions[offset].second[i] = 0.0;
342   }
343 }
344 
int_dist2Float(Z_NR<mpz_t> int_dist)345 FP_NR<mpfr_t> ExactErrorBoundedEvaluator::int_dist2Float(Z_NR<mpz_t> int_dist)
346 {
347   FP_NR<mpfr_t> fMaxDist, maxDE;
348   fMaxDist.set_z(int_dist, GMP_RNDU);
349   bool result = get_max_error_aux(fMaxDist, true, maxDE);
350   FPLLL_CHECK(result, "ExactEvaluator: error cannot be bounded");
351   FPLLL_CHECK(maxDE <= r(0, 0), "ExactEvaluator: max error is too large");
352   fMaxDist.add(fMaxDist, maxDE);
353   //  fMaxDist.mul_2si(fMaxDist, -normExp);
354   return fMaxDist;
355 }
356 
357 FPLLL_END_NAMESPACE
358