1 //==============================================================================================
2 //
3 // This file is part of LiDIA --- a library for computational number theory
4 //
5 // Copyright (c) 1994--2001 the LiDIA Group. All rights reserved.
6 //
7 // See http://www.informatik.tu-darmstadt.de/TI/LiDIA/
8 //
9 //----------------------------------------------------------------------------------------------
10 //
11 // $Id$
12 //
13 // Author :Werner Backes (WB), Thorsten Lauer (TL)
14 // Changes : See CVS log
15 //
16 //==============================================================================================
17
18
19 #ifdef HAVE_CONFIG_H
20 # include "config.h"
21 #endif
22 #include "LiDIA/lattices/bf_lattice_gensys.h"
23
24
25
26 #ifdef LIDIA_NAMESPACE
27 namespace LiDIA {
28 #endif
29
30
31
32 //
33 // using doubles for approximation
34 //
35 //
36 // Schnorr - Euchner
37 // * Tr_lll_dbl()
38 // * Tr_lll_var_dbl();
39 // * Tr_lll_trans_dbl(T)
40 // * Tr_lll_trans_var_dbl(T)
41 // * Tr_lll_deep_insert_dbl() (not yet implemented)
42 // * Tr_lll_dbl_gensys()
43 // * Tr_lll_var_dbl_gensys()
44 // * Tr_lll_trans_dbl_gensys(T)
45 // * Tr_lll_trans_var_dbl_gensys(T)
46 //
47 // Modified lll
48 // * Tr_mlll_dbl(v) (not yet implemented, possible ?)
49 //
50
51
52 //
53 // Schnorr - Euchner version of lll optimized for bigfloats using doubles
54 // result : lll - reduced lattice for parameter y
55 //
Tr_lll_dbl()56 void bigfloat_lattice_gensys::Tr_lll_dbl()
57 {
58 debug_handler("bigfloat_lattice_gensys", "Tr_lll_dbl() [1]");
59 bool Fc;
60 bool korr;
61 bool repeat_loop;
62 lidia_size_t k;
63 sdigit r_prec;
64 sdigit n_prec;
65 sdigit old_prec;
66 double bound = std::exp(std::log(2.0)*52/2);
67 double invbound = 1/bound;
68 double tempdbl0;
69 double *mydel;
70 double **my;
71 double **lattice;
72 double *vect;
73 double *vect_sq;
74 bigint x_factor;
75 bigfloat conv;
76 bigfloat *tempvect1;
77
78 vectsize = columns;
79
80 //
81 // Compute needed precision for exact arithmetik
82 //
83 old_prec = bigfloat::get_precision();
84 r_prec = compute_read_precision();
85 n_prec = static_cast<sdigit>(r_prec/std::log(10.0))+1;
86 n_prec *= columns; // columns >= rows, basis !!!
87 n_prec = static_cast<sdigit>(n_prec*std::log(10.0))+1;
88 bigfloat::set_precision(n_prec);
89
90 //
91 // Allocating Memory for matrix
92 //
93 my = new double*[2*rows];
94 memory_handler(my, "bigfloat_lattice_gensys", "Tr_lll_dbl() :: "
95 "not enough memory !");
96 lattice = &my[rows];
97 mydel = new double[rows*(rows+columns+2)];
98 memory_handler(mydel, "bigfloat_lattice_gensys", "Tr_lll_dbl() :: "
99 "not enough memory !");
100 for (lidia_size_t i = 0; i < rows; i++) {
101 my[i] = &mydel[i*rows]; // rows x rows
102 lattice[i] = &mydel[(rows*rows)+(i*columns)]; // rows x columns
103 }
104
105 //
106 // Allocating Memory for vector
107 //
108
109 tempvect1 = new bigfloat[columns];
110 memory_handler(tempvect1, "bigfloat_lattice_gensys", "Tr_lll_dbl() :: "
111 "not enough memory !");
112 vect = &mydel[(rows+columns)*rows];
113 vect_sq = &vect[rows];
114
115 //
116 // Start
117 //
118 // Copy into floating point
119 //
120 for (lidia_size_t i = 0; i < rows; i++)
121 for (lidia_size_t j = 0; j < columns; j++)
122 value[i][j].doublify(lattice[i][j]);
123
124 k = 1;
125 reduction_steps = 0;
126 swapping_steps = 0;
127 correction_steps = 0;
128 vectsize = columns;
129 bfl_assign_zero_dbl(vect);
130 bfl_scalprod_dbl(vect[0], lattice[0], lattice[0]);
131
132 while (k < rows) {
133 Fc = true;
134 repeat_loop = false;
135 while (Fc) {
136 //
137 // begin of orthogonalization
138 //
139 debug_handler("bigfloat_lattice_gensys", "Tr_lll_dbl() [2]");
140 bfl_scalprod_dbl(vect[k], lattice[k], lattice[k]);
141 vect_sq[k] = std::sqrt(vect[k]);
142 for (lidia_size_t j = 0; j < k; ++j) {
143 bfl_scalprod_dbl(tempdbl0, lattice[k], lattice[j]);
144 //
145 // Check Precision, correct if nessesary
146 //
147 if (std::fabs(tempdbl0) < invbound*vect_sq[k]*vect_sq[j]) {
148 bfl_scalprod_bfl(tempmz0, value[k], value[j]);
149 tempmz0.doublify(my[j][k]);
150 }
151 else
152 my[j][k] = tempdbl0;
153
154 for (lidia_size_t i = 0; i < j; ++i)
155 my[j][k] -= my[i][j]*my[i][k]*vect[i];
156
157 tempdbl0 = my[j][k];
158 my[j][k] /= vect[j];
159 vect[k] -= my[j][k]*tempdbl0;
160 }
161 //
162 // end of orthogonalization
163 //
164 Fc = false;
165 korr = false;
166 //
167 // begin of reduction
168 //
169 debug_handler("bigfloat_lattice_gensys", "Tr_lll_dbl() [3]");
170 if (!repeat_loop) {
171 for (lidia_size_t j = k-1; j >= 0; j--)
172 if (std::fabs(my[j][k]) > 0.5) {
173 ++reduction_steps;
174 korr = true;
175 tempdbl0 = rint(my[j][k]);
176 if (std::fabs(tempdbl0) > bound)
177 Fc = true;
178
179 tempmz0.assign(tempdbl0);
180 bfl_scalmul_bfl(tempvect1, tempmz0, value[j]);
181 bfl_subtract_bfl(value[k], value[k], tempvect1);
182
183 for (lidia_size_t i = 0; i < j; ++i) {
184 my[i][k] -= tempdbl0*my[i][j];
185 }
186
187 my[j][k] -= tempdbl0;
188
189 }
190 }
191 //
192 // correction (not in paper)
193 //
194
195 if (korr) {
196 ++correction_steps;
197 for (lidia_size_t j = 0; j < columns; j++)
198 value[k][j].doublify(lattice[k][j]);
199 }
200
201 //
202 // end of reduction
203 //
204 if (Fc) {
205 if (k > 1)
206 --k;
207 }
208 else {
209 Fc = (!repeat_loop) && korr;
210 repeat_loop = !repeat_loop;
211 }
212 }
213
214 debug_handler("bigfloat_lattice_gensys", "Tr_lll_dbl() [4]");
215 if (y_par*vect[k-1] > my[k-1][k]*my[k-1][k]*vect[k-1]+vect[k]) {
216 ++swapping_steps;
217 bfl_swap_bfl(value[k], value[k-1]);
218 bfl_swap_dbl(lattice[k], lattice[k-1]);
219
220 if (k > 1)
221 --k;
222 if (k == 1) {
223 bfl_scalprod_dbl(vect[0], lattice[0], lattice[0]);
224 vect_sq[0] = std::sqrt(vect[0]);
225 }
226 }
227 else
228 ++k;
229 }
230 //
231 // Free allocated storage
232 // and restore precision
233 //
234 bigfloat::set_precision(old_prec);
235 delete[] tempvect1;
236 delete[] mydel;
237 delete[] my;
238 }
239
240 //
241 // Schnorr - Euchner version of lll optimized for bigints using doubles
242 // result : transformation lattice and reduced lattice
243 //
Tr_lll_trans_dbl(math_matrix<bigint> & Tr)244 void bigfloat_lattice_gensys::Tr_lll_trans_dbl(math_matrix< bigint > & Tr)
245 {
246 debug_handler("bigfloat_lattice_gensys", "Tr_lll_trans_dbl(Tr) [1]");
247 bool Fc;
248 bool korr;
249 bool repeat_loop;
250 lidia_size_t k;
251 sdigit r_prec;
252 sdigit n_prec;
253 sdigit old_prec;
254 double bound = std::exp(std::log(2.0)*52/2);
255 double invbound = 1/bound;
256 double tempdbl0;
257 double *mydel;
258 double **my;
259 double **lattice;
260 double *vect;
261 double *vect_sq;
262 bigint conv;
263 bigint x_factor;
264 bigint *tempvect2;
265 bigint *Tdel;
266 bigint **T;
267 bigint **TrAddr;
268 bigfloat *tempvect1;
269
270 vectsize = columns;
271
272 //
273 // Compute needed precision for exact arithmetik
274 //
275 old_prec = bigfloat::get_precision();
276 r_prec = compute_read_precision();
277 n_prec = static_cast<sdigit>(r_prec/std::log(10.0))+1;
278 n_prec *= columns; // columns >= rows, basis !!!
279 n_prec = static_cast<sdigit>(n_prec*std::log(10.0))+1;
280 bigfloat::set_precision(n_prec);
281
282 //
283 // Allocating Memory for matrix
284 //
285 my = new double*[2*rows];
286 memory_handler(my, "bigfloat_lattice_gensys", "Tr_lll_trans_dbl(Tr) :: "
287 "not enough memory !");
288 lattice = &my[rows];
289 T = new bigint*[rows];
290 memory_handler(T, "bigfloat_lattice_gensys", "Tr_lll_trans_dbl(Tr) :: "
291 "Not enough memory !");
292
293 mydel = new double[rows*(rows+columns+2)];
294 memory_handler(mydel, "bigfloat_lattice_gensys", "Tr_lll_trans_dbl(Tr) :: "
295 "not enough memory !");
296 Tdel = new bigint[rows*rows+rows];
297 memory_handler(Tdel, "bigfloat_lattice_gensys", "Tr_lll_trans_dbl(Tr) :: "
298 "Not enough memory !");
299 for (lidia_size_t i = 0; i < rows; i++) {
300 my[i] = &mydel[i*rows]; // rows x rows
301 lattice[i] = &mydel[(rows*rows)+(i*columns)]; // rows x columns
302 T[i] = &Tdel[i*rows];
303 T[i][i].assign_one();
304 }
305
306 //
307 // Allocating Memory for vector
308 //
309
310 tempvect1 = new bigfloat[columns];
311 memory_handler(tempvect1, "bigfloat_lattice_gensys", "Tr_lll_trans_dbl(Tr) :: "
312 "not enough memory !");
313 tempvect2 = &Tdel[rows*rows];
314 vect = &mydel[(rows+columns)*rows];
315 vect_sq = &vect[rows];
316
317 //
318 // Start
319 //
320 // Copy into floating point
321 //
322 for (lidia_size_t i = 0; i < rows; i++)
323 for (lidia_size_t j = 0; j < columns; j++)
324 value[i][j].doublify(lattice[i][j]);
325
326 k = 1;
327 reduction_steps = 0;
328 swapping_steps = 0;
329 correction_steps = 0;
330 vectsize = columns;
331 bfl_assign_zero_dbl(vect);
332 bfl_scalprod_dbl(vect[0], lattice[0], lattice[0]);
333
334 while (k < rows) {
335 Fc = true;
336 repeat_loop = false;
337 while (Fc) {
338 //
339 // begin of orthogonalization
340 //
341 debug_handler("bigfloat_lattice_gensys", "Tr_lll_trans_dbl(Tr) [2]");
342 bfl_scalprod_dbl(vect[k], lattice[k], lattice[k]);
343 vect_sq[k] = std::sqrt(vect[k]);
344 for (lidia_size_t j = 0; j < k; ++j) {
345 bfl_scalprod_dbl(tempdbl0, lattice[k], lattice[j]);
346 //
347 // Check Precision, correct if nessesary
348 //
349 if (std::fabs(tempdbl0) < invbound*vect_sq[k]*vect_sq[j]) {
350 bfl_scalprod_bfl(tempmz0, value[k], value[j]);
351 tempmz0.doublify(my[j][k]);
352 }
353 else
354 my[j][k] = tempdbl0;
355
356 for (lidia_size_t i = 0; i < j; ++i)
357 my[j][k] -= my[i][j]*my[i][k]*vect[i];
358
359 tempdbl0 = my[j][k];
360 my[j][k] /= vect[j];
361 vect[k] -= my[j][k]*tempdbl0;
362 }
363 //
364 // end of orthogonalization
365 //
366 Fc = false;
367 korr = false;
368 //
369 // begin of reduction
370 //
371 debug_handler("bigfloat_lattice_gensys", "Tr_lll_trans_dbl(Tr) [3]");
372 if (!repeat_loop) {
373 for (lidia_size_t j = k-1; j >= 0; j--)
374 if (std::fabs(my[j][k]) > 0.5) {
375 ++reduction_steps;
376 korr = true;
377 tempdbl0 = rint(my[j][k]);
378 if (std::fabs(tempdbl0) > bound)
379 Fc = true;
380
381 tempmz0.assign(tempdbl0);
382 tempmz0.bigintify(conv);
383 bfl_scalmul_bfl(tempvect1, tempmz0, value[j]);
384 bfl_subtract_bfl(value[k], value[k], tempvect1);
385 vectsize = rows;
386 bfl_scalmul_bin(tempvect2, conv, T[j]);
387 bfl_subtract_bin(T[k], T[k], tempvect2);
388 vectsize = columns;
389 for (lidia_size_t i = 0; i < j; ++i)
390 my[i][k] -= tempdbl0*my[i][j];
391
392 my[j][k] -= tempdbl0;
393 }
394 }
395
396 //
397 // correction (not in paper)
398 //
399
400 if (korr) {
401 ++correction_steps;
402 for (lidia_size_t j = 0; j < columns; j++)
403 value[k][j].doublify(lattice[k][j]);
404 }
405
406 //
407 // end of reduction
408 //
409 if (Fc) {
410 if (k > 1)
411 --k;
412 }
413 else {
414 Fc = (!repeat_loop) && korr;
415 repeat_loop = !repeat_loop;
416 }
417 }
418
419 debug_handler("bigfloat_lattice_gensys", "Tr_lll_trans_dbl(Tr) [4]");
420 if (y_par*vect[k-1] > my[k-1][k]*my[k-1][k]*vect[k-1]+vect[k]) {
421 ++swapping_steps;
422 bfl_swap_bfl(value[k], value[k-1]);
423 bfl_swap_dbl(lattice[k], lattice[k-1]);
424 bfl_swap_bin(T[k], T[k-1]);
425
426 if (k > 1)
427 --k;
428 if (k == 1) {
429 bfl_scalprod_dbl(vect[0], lattice[0], lattice[0]);
430 vect_sq[0] = std::sqrt(vect[0]);
431 }
432 }
433 else
434 ++k;
435 }
436 //
437 // Copy into math_matrix< bigint > Tr
438 //
439 Tr.set_no_of_rows(rows);
440 Tr.set_no_of_columns(rows);
441 TrAddr = Tr.get_data_address();
442 for (lidia_size_t i = 0; i < rows; i++)
443 for (lidia_size_t j = 0; j < rows; j++)
444 if (trans_flag)
445 LiDIA::swap(TrAddr[j][i], T[i][j]);
446 else
447 LiDIA::swap(TrAddr[j][i], T[i][j]);
448 //
449 // restore precision and
450 // free allocated storage
451 //
452 bigfloat::set_precision(old_prec);
453 delete[] tempvect1;
454 delete[] Tdel;
455 delete[] T;
456 delete[] mydel;
457 delete[] my;
458 }
459
460 //
461 // Schnorr - Euchner version of lll optimized for bigfloats using doubles
462 // result : lll - reduced lattice for parameter y
463 //
Tr_lll_dbl_gensys(lidia_size_t & rank)464 void bigfloat_lattice_gensys::Tr_lll_dbl_gensys(lidia_size_t& rank)
465 {
466 debug_handler("bigfloat_lattice_gensys", "Tr_lll_dbl_gensys(rank) [1]");
467 bool Fc;
468 bool korr;
469 bool break_flag;
470 bool repeat_loop;
471 bool vector_is_zero;
472 lidia_size_t k;
473 sdigit r_prec;
474 sdigit n_prec;
475 sdigit old_prec;
476 lidia_size_t mom_rows = rows;
477 double bound = std::exp(std::log(2.0)*52/2);
478 double invbound = 1/bound;
479 double u_zero = 0.0;
480 double tempdbl0;
481 double *mydel;
482 double **my;
483 double **lattice;
484 double *vect;
485 double *vect_sq;
486 bigint x_factor;
487 bigfloat u_zero_bfl;
488 bigfloat conv;
489 bigfloat *tempvect1;
490
491 vectsize = columns;
492
493 //
494 // Compute needed precision for exact arithmetik
495 //
496 old_prec = bigfloat::get_precision();
497 r_prec = compute_read_precision();
498 if (r_prec < 200)
499 u_zero = 1/std::exp(std::log(10.0)*r_prec);
500 else {
501 log(u_zero_bfl, 10);
502 LiDIA::multiply(u_zero_bfl, u_zero_bfl, r_prec);
503 exp(u_zero_bfl, u_zero_bfl);
504 LiDIA::divide(u_zero_bfl, 1, u_zero_bfl);
505 }
506 n_prec = static_cast<sdigit>(r_prec/std::log(10.0))+1;
507 n_prec *= (rows > columns)?rows:columns; // columns > rows ? da event. keine basis !!!
508 n_prec = static_cast<sdigit>(n_prec*std::log(10.0))+1;
509 bigfloat::set_precision(n_prec);
510
511 //
512 // Allocating Memory for matrix
513 //
514 my = new double*[2*rows];
515 memory_handler(my, "bigfloat_lattice_gensys", "Tr_lll_dbl_gensys(rank) :: "
516 "not enough memory !");
517 lattice = &my[rows];
518 mydel = new double[rows*(rows+columns+2)];
519 memory_handler(mydel, "bigfloat_lattice_gensys", "Tr_lll_dbl_gensys(rank) :: "
520 "not enough memory !");
521 for (lidia_size_t i = 0; i < rows; i++) {
522 my[i] = &mydel[i*rows]; // rows x rows
523 lattice[i] = &mydel[(rows*rows)+(i*columns)]; // rows x columns
524 }
525
526 //
527 // Allocating Memory for vector
528 //
529
530 tempvect1 = new bigfloat[columns];
531 memory_handler(tempvect1, "bigfloat_lattice_gensys", "Tr_lll_dbl_gensys(rank) :: "
532 "not enough memory !");
533
534 vect = &mydel[(rows+columns)*rows];
535 vect_sq = &vect[rows];
536
537 //
538 // Start
539 //
540 // Copy into floating point
541 //
542 for (lidia_size_t i = 0; i < rows; i++)
543 for (lidia_size_t j = 0; j < columns; j++)
544 value[i][j].doublify(lattice[i][j]);
545
546 k = 1;
547 reduction_steps = 0;
548 swapping_steps = 0;
549 correction_steps = 0;
550 vectsize = columns;
551 break_flag = false;
552 bfl_assign_zero_dbl(vect);
553 bfl_scalprod_dbl(vect[0], lattice[0], lattice[0]);
554
555 while (k < mom_rows) {
556 Fc = true;
557 repeat_loop = false;
558 while (Fc) {
559 //
560 // begin of orthogonalization
561 //
562 debug_handler("bigfloat_lattice_gensys", "Tr_lll_dbl_gensys(rank) [2]");
563 bfl_scalprod_dbl(vect[k], lattice[k], lattice[k]);
564 vect_sq[k] = std::sqrt(vect[k]);
565 for (lidia_size_t j = 0; j < k; ++j) {
566 bfl_scalprod_dbl(tempdbl0, lattice[k], lattice[j]);
567 //
568 // Correction Step I, only if needed, |lattice[k]*lattice[j]| to big
569 //
570 if (std::fabs(tempdbl0) < invbound*vect_sq[k]*vect_sq[j]) {
571 //
572 // Compute scalar product correct
573 //
574 bfl_scalprod_bfl(tempmz0, value[k], value[j]);
575 tempmz0.doublify(my[j][k]);
576 }
577 else
578 my[j][k] = tempdbl0;
579
580 for (lidia_size_t i = 0; i < j; ++i)
581 my[j][k] -= my[i][j]*my[i][k]*vect[i];
582
583 tempdbl0 = my[j][k];
584 my[j][k] /= vect[j];
585 vect[k] -= my[j][k]*tempdbl0;
586 }
587 //
588 // end of orthogonalization
589 //
590 Fc = false;
591 korr = false;
592 //
593 // begin of reduction
594 //
595 debug_handler("bigfloat_lattice_gensys", "Tr_lll_dbl_gensys(rank) [3]");
596 if (!repeat_loop) {
597 for (lidia_size_t j = k-1; j >= 0; j--)
598 if (std::fabs(my[j][k]) > 0.5) {
599 ++reduction_steps;
600 korr = true;
601 tempdbl0 = rint(my[j][k]);
602 if (std::fabs(tempdbl0) > bound)
603 Fc = true;
604
605 tempmz0.assign(tempdbl0);
606 bfl_scalmul_bfl(tempvect1, tempmz0, value[j]);
607 bfl_subtract_bfl(value[k], value[k], tempvect1);
608
609 for (lidia_size_t i = 0; i < j; ++i)
610 my[i][k] -= tempdbl0*my[i][j];
611
612 my[j][k] -= tempdbl0;
613 }
614 }
615
616 //
617 // Correction Step II,
618 // lattice[k]=(double )value[k]
619 //
620
621 if (korr) {
622 ++correction_steps;
623 for (lidia_size_t j = 0; j < columns; j++)
624 value[k][j].doublify(lattice[k][j]);
625 }
626 //
627 // Check if it`s 0 - vector
628 // delete then
629 //
630
631 if (r_prec < 200) {
632 bfl_scalprod_dbl(tempdbl0, lattice[k], lattice[k]);
633 tempdbl0 = std::sqrt(tempdbl0);
634 if (std::fabs(tempdbl0) < u_zero)
635 vector_is_zero = true;
636 else
637 vector_is_zero = false;
638 }
639 else {
640 bigfloat::set_precision(r_prec);
641 bfl_scalprod_bfl(tempmz0, value[k], value[k]);
642 sqrt(tempmz0, tempmz0);
643 if (abs(tempmz0) < u_zero_bfl)
644 vector_is_zero = true;
645 else
646 vector_is_zero = false;
647 bigfloat::set_precision(n_prec);
648 }
649
650 if (vector_is_zero) {
651 for (lidia_size_t j = k; j < mom_rows-1; j++) {
652 bfl_swap_dbl(lattice[j], lattice[j+1]);
653 bfl_swap_bfl(value[j], value[j+1]);
654 }
655 bfl_assign_zero_bfl(value[--mom_rows]);
656 bfl_assign_zero_dbl(lattice[mom_rows]);
657 k = 1;
658 bfl_assign_zero_dbl(vect);
659 bfl_scalprod_dbl(vect[0], lattice[0], lattice[0]);
660 vect_sq[0] = std::sqrt(vect[0]);
661 break_flag = true;
662 korr = false;
663 Fc = false;
664 break;
665 }
666
667 //
668 // end of reduction
669 //
670 if (Fc) {
671 if (k > 1)
672 --k;
673 }
674 else {
675 Fc = (!repeat_loop) && korr;
676 repeat_loop = !repeat_loop;
677 }
678 }
679
680 if (!break_flag) {
681 //
682 // Check lll - condition for parameter y_par
683 //
684 debug_handler("bigfloat_lattice_gensys", "Tr_lll_dbl_gensys(rank) [4]");
685 if (y_par*vect[k-1] > my[k-1][k]*my[k-1][k]*vect[k-1]+vect[k]) {
686 //
687 // ::swap vectors k and k-1
688 //
689 ++swapping_steps;
690 bfl_swap_bfl(value[k], value[k-1]);
691 bfl_swap_dbl(lattice[k], lattice[k-1]);
692
693 if (k > 1)
694 --k;
695 if (k == 1) {
696 bfl_scalprod_dbl(vect[0], lattice[0], lattice[0]);
697 vect_sq[0] = std::sqrt(vect[0]);
698 }
699 }
700 else
701 ++k;
702 }
703 break_flag = false;
704 }
705 rank = mom_rows;
706 //
707 // Free allocated storage and
708 // restore precision
709 //
710 bigfloat::set_precision(old_prec);
711 delete[] tempvect1;
712 delete[] mydel;
713 delete[] my;
714 }
715
716 //
717 // Schnorr - Euchner version of lll optimized for bigints using doubles
718 // result : transformation lattice and reduced lattice
719 //
Tr_lll_trans_dbl_gensys(math_matrix<bigint> & Tr,lidia_size_t & rank)720 void bigfloat_lattice_gensys::Tr_lll_trans_dbl_gensys(math_matrix< bigint > & Tr, lidia_size_t& rank)
721 {
722 debug_handler("bigfloat_lattice_gensys", "Tr_lll_trans_dbl_gensys(Tr, rank) [1]");
723 bool Fc;
724 bool break_flag;
725 bool korr;
726 bool repeat_loop;
727 bool vector_is_zero;
728 sdigit r_prec;
729 sdigit n_prec;
730 sdigit old_prec;
731 lidia_size_t k;
732 lidia_size_t mom_rows = rows;
733 double bound = std::exp(std::log(2.0)*52/2);
734 double invbound = 1/bound;
735 double u_zero = 0.0;
736 double tempdbl0;
737 double *mydel;
738 double **my;
739 double **lattice;
740 double *vect;
741 double *vect_sq;
742 bigint conv;
743 bigint x_factor;
744 bigint *tempvect2;
745 bigint *Tdel;
746 bigint **T;
747 bigint **TrAddr;
748 bigfloat u_zero_bfl;
749 bigfloat *tempvect1;
750
751 vectsize = columns;
752
753 //
754 // Compute needed precision for exact arithmetik
755 //
756 old_prec = bigfloat::get_precision();
757 r_prec = compute_read_precision();
758 if (r_prec < 200)
759 u_zero = 1/std::exp(std::log(10.0)*r_prec);
760 else {
761 log(u_zero_bfl, 10);
762 LiDIA::multiply(u_zero_bfl, u_zero_bfl, r_prec);
763 exp(u_zero_bfl, u_zero_bfl);
764 LiDIA::divide(u_zero_bfl, 1, u_zero_bfl);
765 }
766 n_prec = static_cast<sdigit>(r_prec/std::log(10.0))+1;
767 n_prec *= (rows > columns)?rows:columns; // columns > rows ? da event. keine basis !!!
768 n_prec = static_cast<sdigit>(n_prec*std::log(10.0))+1;
769 bigfloat::set_precision(n_prec);
770
771 //
772 // Allocating Memory for matrix
773 //
774 my = new double*[2*rows];
775 memory_handler(my, "bigfloat_lattice_gensys", "Tr_lll_trans_dbl(Tr) :: "
776 "not enough memory !");
777 lattice = &my[rows];
778 T = new bigint*[rows];
779 memory_handler(T, "bigfloat_lattice_gensys", "Tr_lll_trans_dbl(Tr) :: "
780 "Not enough memory !");
781
782 mydel = new double[rows*(rows+columns+2)];
783 memory_handler(mydel, "bigfloat_lattice_gensys", "Tr_lll_trans_dbl(Tr) :: "
784 "not enough memory !");
785 Tdel = new bigint[rows*rows+rows];
786 memory_handler(Tdel, "bigfloat_lattice_gensys", "Tr_lll_trans_dbl(Tr) :: "
787 "Not enough memory !");
788 for (lidia_size_t i = 0; i < rows; i++) {
789 my[i] = &mydel[i*rows]; // rows x rows
790 lattice[i] = &mydel[(rows*rows)+(i*columns)]; // rows x columns
791 T[i] = &Tdel[i*rows];
792 T[i][i].assign_one();
793 }
794
795 //
796 // Allocating Memory for vector
797 //
798
799 tempvect1 = new bigfloat[columns];
800 memory_handler(tempvect1, "bigfloat_lattice_gensys", "Tr_lll_trans_dbl(Tr) :: "
801 "not enough memory !");
802 tempvect2 = &Tdel[rows*rows];
803 vect = &mydel[(rows+columns)*rows];
804 vect_sq = &vect[rows];
805
806 //
807 // Start
808 //
809 // Copy into floating point
810 //
811 for (lidia_size_t i = 0; i < rows; i++)
812 for (lidia_size_t j = 0; j < columns; j++)
813 value[i][j].doublify(lattice[i][j]);
814
815 k = 1;
816 reduction_steps = 0;
817 swapping_steps = 0;
818 correction_steps = 0;
819 vectsize = columns;
820 break_flag = false;
821 bfl_assign_zero_dbl(vect);
822 bfl_scalprod_dbl(vect[0], lattice[0], lattice[0]);
823
824 while (k < mom_rows) {
825 Fc = true;
826 repeat_loop = false;
827 while (Fc) {
828 //
829 // begin of orthogonalization
830 //
831 debug_handler("bigfloat_lattice_gensys", "Tr_lll_trans_dbl_gensys(Tr, rank) [2]");
832 bfl_scalprod_dbl(vect[k], lattice[k], lattice[k]);
833 vect_sq[k] = std::sqrt(vect[k]);
834 for (lidia_size_t j = 0; j < k; ++j) {
835 bfl_scalprod_dbl(tempdbl0, lattice[k], lattice[j]);
836 //
837 // Correction Step I, only if needed, |lattice[k]*lattice[j]| to big
838 //
839 if (std::fabs(tempdbl0) < invbound*vect_sq[k]*vect_sq[j]) {
840 //
841 // Compute the correct scalar product
842 //
843 bfl_scalprod_bfl(tempmz0, value[k], value[j]);
844 tempmz0.doublify(my[j][k]);
845 }
846 else
847 my[j][k] = tempdbl0;
848
849 for (lidia_size_t i = 0; i < j; ++i)
850 my[j][k] -= my[i][j]*my[i][k]*vect[i];
851
852 tempdbl0 = my[j][k];
853 my[j][k] /= vect[j];
854 vect[k] -= my[j][k]*tempdbl0;
855 }
856 //
857 // end of orthogonalization
858 //
859 Fc = false;
860 korr = false;
861 //
862 // begin of reduction
863 //
864 debug_handler("bigfloat_lattice_gensys", "Tr_lll_trans_dbl(Tr) [3]");
865 if (!repeat_loop) {
866 for (lidia_size_t j = k-1; j >= 0; j--)
867 if (std::fabs(my[j][k]) > 0.5) {
868 ++reduction_steps;
869 korr = true;
870 tempdbl0 = rint(my[j][k]);
871 if (std::fabs(tempdbl0) > bound)
872 Fc = true;
873
874 tempmz0.assign(tempdbl0);
875 tempmz0.bigintify(conv);
876 bfl_scalmul_bfl(tempvect1, tempmz0, value[j]);
877 bfl_subtract_bfl(value[k], value[k], tempvect1);
878 vectsize = rows;
879 bfl_scalmul_bin(tempvect2, conv, T[j]);
880 bfl_subtract_bin(T[k], T[k], tempvect2);
881 vectsize = columns;
882 for (lidia_size_t i = 0; i < j; ++i)
883 my[i][k] -= tempdbl0*my[i][j];
884
885 my[j][k] -= tempdbl0;
886 }
887 }
888
889 //
890 // Corrections Step II
891 // lattice[k]=(double )value[k]
892 //
893
894 if (korr) {
895 ++correction_steps;
896 for (lidia_size_t j = 0; j < columns; j++)
897 value[k][j].doublify(lattice[k][j]);
898 }
899 //
900 // Check if it`s 0 - vector
901 // then delete it
902 //
903 if (r_prec < 200) {
904 bfl_scalprod_dbl(tempdbl0, lattice[k], lattice[k]);
905 tempdbl0 = std::sqrt(tempdbl0);
906 if (std::fabs(tempdbl0) < u_zero)
907 vector_is_zero = true;
908 else
909 vector_is_zero = false;
910 }
911 else {
912 bigfloat::set_precision(r_prec);
913 bfl_scalprod_bfl(tempmz0, value[k], value[k]);
914 sqrt(tempmz0, tempmz0);
915 if (abs(tempmz0) < u_zero_bfl)
916 vector_is_zero = true;
917 else
918 vector_is_zero = false;
919 bigfloat::set_precision(n_prec);
920 }
921
922 if (vector_is_zero) {
923 for (lidia_size_t j = 0; j < mom_rows-1; j++) {
924 bfl_swap_bfl(value[j], value[j+1]);
925 bfl_swap_dbl(lattice[j], lattice[j+1]);
926 bfl_swap_bin(T[j], T[j+1]);
927 }
928 bfl_assign_zero_bfl(value[--mom_rows]);
929 vectsize = rows;
930 bfl_assign_zero_bin(T[mom_rows]);
931 vectsize = columns;
932 k = 1;
933 bfl_assign_zero_dbl(vect);
934 bfl_scalprod_dbl(vect[0], lattice[0], lattice[0]);
935 vect_sq[0] = std::sqrt(vect[0]);
936 break_flag = true;
937 korr = false;
938 Fc = false;
939 break;
940 }
941
942 //
943 // end of reduction
944 //
945 if (Fc) {
946 if (k > 1)
947 --k;
948 }
949 else {
950 Fc = (!repeat_loop) && korr;
951 repeat_loop = !repeat_loop;
952 }
953 }
954
955 if (!break_flag) {
956 //
957 // Check lll - condition for parameter y_par
958 //
959 debug_handler("bigfloat_lattice_gensys", "Tr_lll_trans_dbl(Tr) [4]");
960 if (y_par*vect[k-1] > my[k-1][k]*my[k-1][k]*vect[k-1]+vect[k]) {
961 //
962 // swap vectors k and k-1
963 //
964 ++swapping_steps;
965 bfl_swap_bfl(value[k], value[k-1]);
966 bfl_swap_dbl(lattice[k], lattice[k-1]);
967 bfl_swap_bin(T[k], T[k-1]);
968
969 if (k > 1)
970 --k;
971 if (k == 1) {
972 bfl_scalprod_dbl(vect[0], lattice[0], lattice[0]);
973 vect_sq[0] = std::sqrt(vect[0]);
974 }
975 }
976 else
977 ++k;
978 }
979 break_flag = false;
980 }
981 rank = mom_rows;
982 //
983 // Copy into math_matrix< bigint > Tr
984 //
985 Tr.set_no_of_rows(rows);
986 Tr.set_no_of_columns(rows);
987 TrAddr = Tr.get_data_address();
988 for (lidia_size_t i = 0; i < rows; i++)
989 for (lidia_size_t j = 0; j < rows; j++)
990 if (trans_flag)
991 LiDIA::swap(TrAddr[j][i], T[i][j]);
992 else
993 LiDIA::swap(TrAddr[j][i], T[j][i]);
994 //
995 // Free allocated storage and
996 // restore precision
997 //
998 bigfloat::set_precision(old_prec);
999 delete[] tempvect1;
1000 delete[] Tdel;
1001 delete[] T;
1002 delete[] mydel;
1003 delete[] my;
1004 }
1005
1006
1007 //
1008 // Schnorr - Euchner version of lll optimized for bigfloats using doubles
1009 // result : lll - reduced lattice for parameter y
1010 //
Tr_lll_var_dbl(math_matrix<bigint> & Bi,sdigit bit_len)1011 void bigfloat_lattice_gensys::Tr_lll_var_dbl(math_matrix< bigint > & Bi, sdigit bit_len)
1012 {
1013 debug_handler("bigfloat_lattice_gensys", "Tr_lll_var_dbl(Bi, bit_len) [1]");
1014 bool Fc;
1015 bool korr;
1016 bool repeat_loop;
1017 lidia_size_t k;
1018 sdigit old_prec;
1019 sdigit bit_diff = 0;
1020 sdigit bi_bit_len;
1021 sdigit old_bit_prec;
1022 double bound = std::exp(std::log(2.0)*52/2);
1023 double invbound = 1/bound;
1024 double tempdbl0;
1025 double *mydel;
1026 double **my;
1027 double **lattice;
1028 double *vect;
1029 double *vect_sq;
1030 bigint tempbiz;
1031 bigint **BiAddr;
1032 bigint *tempvect1;
1033
1034 vectsize = columns;
1035 old_prec = bigfloat::get_precision();
1036 old_bit_prec = static_cast<sdigit>(static_cast<double>(old_prec)/std::log(2.0)*std::log(10.0)+1);
1037 bigfloat::set_precision(20);
1038 BiAddr = Bi.get_data_address();
1039 //
1040 // Allocating Memory for matrix
1041 //
1042 my = new double*[2*rows];
1043 memory_handler(my, "bigfloat_lattice_gensys", "Tr_lll_var_dbl(Bi) :: "
1044 "not enough memory !");
1045 lattice = &my[rows];
1046 mydel = new double[rows*(rows+columns+2)];
1047 memory_handler(mydel, "bigfloat_lattice_gensys", "Tr_lll_var_dbl(Bi) :: "
1048 "not enough memory !");
1049 for (lidia_size_t i = 0; i < rows; i++) {
1050 my[i] = &mydel[i*rows]; // rows x rows
1051 lattice[i] = &mydel[(rows*rows)+(i*columns)]; // rows x columns
1052 }
1053
1054 //
1055 // Allocating Memory for vector
1056 //
1057
1058 tempvect1 = new bigint[columns];
1059 memory_handler(tempvect1, "bigfloat_lattice_gensys", "Tr_lll_var_dbl(Bi) :: "
1060 "not enough memory !");
1061 vect = &mydel[(rows+columns)*rows];
1062 vect_sq = &vect[rows];
1063
1064 //
1065 // Start
1066 //
1067 // Copy into floating point
1068 //
1069 for (lidia_size_t i = 0; i < rows; i++)
1070 for (lidia_size_t j = 0; j < columns; j++) {
1071 bi_bit_len = BiAddr[i][j].bit_length();
1072 if (bi_bit_len > 8*SIZEOF_DOUBLE) {
1073 bit_diff = bi_bit_len-8*SIZEOF_DOUBLE;
1074 shift_right(tempbiz, BiAddr[i][j], bit_diff);
1075 tempmz0.assign(tempbiz);
1076 }
1077 else {
1078 bit_diff = 0;
1079 tempmz0.assign(BiAddr[i][j]);
1080 }
1081 shift_right(tempmz0, tempmz0, bit_len-bit_diff);
1082 tempmz0.doublify(lattice[i][j]);
1083 }
1084
1085 k = 1;
1086 reduction_steps = 0;
1087 swapping_steps = 0;
1088 correction_steps = 0;
1089 vectsize = columns;
1090 bfl_assign_zero_dbl(vect);
1091 bfl_scalprod_dbl(vect[0], lattice[0], lattice[0]);
1092
1093 while (k < rows) {
1094 Fc = true;
1095 repeat_loop = false;
1096 while (Fc) {
1097 //
1098 // begin of orthogonalization
1099 //
1100 debug_handler("bigfloat_lattice_gensys", "Tr_lll_var_dbl(Bi, Factor) [2]");
1101 bfl_scalprod_dbl(vect[k], lattice[k], lattice[k]);
1102 vect_sq[k] = std::sqrt(vect[k]);
1103 for (lidia_size_t j = 0; j < k; ++j) {
1104 bfl_scalprod_dbl(tempdbl0, lattice[k], lattice[j]);
1105 //
1106 // Check Precision, correct if nessesary
1107 //
1108 if (std::fabs(tempdbl0) < invbound*vect_sq[k]*vect_sq[j]) {
1109 bfl_scalprod_bin(tempbiz, BiAddr[k], BiAddr[j]);
1110 bi_bit_len = tempbiz.bit_length();
1111 if (bi_bit_len > 8*SIZEOF_DOUBLE) {
1112 bit_diff = bi_bit_len-8*SIZEOF_DOUBLE;
1113 shift_right(tempbiz, tempbiz, bit_diff);
1114 tempmz0.assign(tempbiz);
1115 }
1116 else {
1117 bit_diff = 0;
1118 tempmz0.assign(tempbiz);
1119 }
1120 shift_right(tempmz0, tempmz0, 2*bit_len-bit_diff);
1121 tempmz0.doublify(my[j][k]);
1122 }
1123 else
1124 my[j][k] = tempdbl0;
1125
1126 for (lidia_size_t i = 0; i < j; ++i)
1127 my[j][k] -= my[i][j]*my[i][k]*vect[i];
1128
1129 tempdbl0 = my[j][k];
1130 my[j][k] /= vect[j];
1131 vect[k] -= my[j][k]*tempdbl0;
1132 }
1133 //
1134 // end of orthogonalization
1135 //
1136
1137 Fc = false;
1138 korr = false;
1139 //
1140 // begin of reduction
1141 //
1142 debug_handler("bigfloat_lattice_gensys", "Tr_lll_var_dbl(Bi, Factor) [3]");
1143 if (!repeat_loop) {
1144 for (lidia_size_t j = k-1; j >= 0; j--)
1145 if (std::fabs(my[j][k]) > 0.5) {
1146 ++reduction_steps;
1147 korr = true;
1148 tempdbl0 = rint(my[j][k]);
1149 if (std::fabs(tempdbl0) > bound)
1150 Fc = true;
1151
1152 tempmz0.assign(tempdbl0);
1153 tempmz0.bigintify(tempbiz);
1154 bfl_scalmul_bin(tempvect1, tempbiz, BiAddr[j]);
1155 bfl_subtract_bin(BiAddr[k], BiAddr[k], tempvect1);
1156
1157 for (lidia_size_t i = 0; i < j; ++i)
1158 my[i][k] -= tempdbl0*my[i][j];
1159
1160 my[j][k] -= tempdbl0;
1161 }
1162 }
1163
1164 //
1165 // correction (not in paper)
1166 //
1167
1168 if (korr) {
1169 ++correction_steps;
1170 for (lidia_size_t j = 0; j < columns; j++) {
1171 bi_bit_len = BiAddr[k][j].bit_length();
1172 if (bi_bit_len > 8*SIZEOF_DOUBLE) {
1173 bit_diff = bi_bit_len-8*SIZEOF_DOUBLE;
1174 shift_right(tempbiz, BiAddr[k][j], bit_diff);
1175 tempmz0.assign(tempbiz);
1176 }
1177 else {
1178 bit_diff = 0;
1179 tempmz0.assign(BiAddr[k][j]);
1180 }
1181 shift_right(tempmz0, tempmz0, bit_len-bit_diff);
1182 tempmz0.doublify(lattice[k][j]);
1183 }
1184 }
1185
1186 //
1187 // end of reduction
1188 //
1189 if (Fc) {
1190 if (k > 1)
1191 --k;
1192 }
1193 else {
1194 Fc = (!repeat_loop) && korr;
1195 repeat_loop = !repeat_loop;
1196 }
1197 }
1198
1199
1200 debug_handler("bigfloat_lattice_gensys", "Tr_lll_var_dbl(Bi, Factor) [4]");
1201 if (y_par*vect[k-1] > my[k-1][k]*my[k-1][k]*vect[k-1]+vect[k]) {
1202 ++swapping_steps;
1203 bfl_swap_bin(BiAddr[k], BiAddr[k-1]);
1204 bfl_swap_dbl(lattice[k], lattice[k-1]);
1205
1206 if (k > 1)
1207 --k;
1208 if (k == 1) {
1209 bfl_scalprod_dbl(vect[0], lattice[0], lattice[0]);
1210 vect_sq[0] = std::sqrt(vect[0]);
1211 }
1212 }
1213 else
1214 ++k;
1215 }
1216 bigfloat::set_precision(old_prec);
1217
1218 for (lidia_size_t i = 0; i < rows; i++)
1219 for (lidia_size_t j = 0; j < columns; j++) {
1220 bi_bit_len = BiAddr[i][j].bit_length();
1221 if (bi_bit_len > old_bit_prec) {
1222 bit_diff = bi_bit_len-old_bit_prec;
1223 shift_right(tempbiz, BiAddr[i][j], bit_diff);
1224 tempmz0.assign(tempbiz);
1225 }
1226 else {
1227 bit_diff = 0;
1228 tempmz0.assign(BiAddr[i][j]);
1229 }
1230 shift_right(value[i][j], tempmz0, bit_len-bit_diff);
1231 }
1232
1233 //
1234 // Free allocated storage
1235 // and restore precision
1236 //
1237 delete[] tempvect1;
1238 delete[] mydel;
1239 delete[] my;
1240 }
1241
1242
1243 //
1244 // Schnorr - Euchner version of lll optimized for bigfloats using doubles
1245 // result : lll - reduced lattice for parameter y
1246 //
Tr_lll_trans_var_dbl(math_matrix<bigint> & Bi,math_matrix<bigint> & Tr,sdigit bit_len)1247 void bigfloat_lattice_gensys::Tr_lll_trans_var_dbl(math_matrix< bigint > & Bi,
1248 math_matrix< bigint > & Tr,
1249 sdigit bit_len)
1250 {
1251 debug_handler("bigfloat_lattice_gensys", "Tr_lll_trans_var_dbl(Bi, Tr, Factor, bit_len) [1]");
1252 bool Fc;
1253 bool korr;
1254 bool repeat_loop;
1255 lidia_size_t k;
1256 sdigit old_prec;
1257 sdigit old_bit_prec;
1258 sdigit bit_diff;
1259 sdigit bi_bit_len;
1260 double bound = std::exp(std::log(2.0)*52/2);
1261 double invbound = 1/bound;
1262 double tempdbl0;
1263 double *mydel;
1264 double **my;
1265 double **lattice;
1266 double *vect;
1267 double *vect_sq;
1268 bigint tempbiz;
1269 bigint *Tdel;
1270 bigint **T;
1271 bigint **TrAddr;
1272 bigint **BiAddr;
1273 bigint *tempvect1;
1274
1275 vectsize = columns;
1276 old_prec = bigfloat::get_precision();
1277 old_bit_prec = static_cast<sdigit>(static_cast<double>(old_prec)/std::log(2.0)*std::log(10.0)+1);
1278 bigfloat::set_precision(20);
1279 BiAddr = Bi.get_data_address();
1280 //
1281 // Allocating Memory for matrix
1282 //
1283 my = new double*[2*rows];
1284 memory_handler(my, "bigfloat_lattice_gensys", "Tr_lll_trans_var_dbl(Bi, Tr, Factor) :: "
1285 "not enough memory !");
1286 lattice = &my[rows];
1287 T = new bigint*[rows];
1288 memory_handler(T, "bigfloat_lattice_gensys", "Tr_lll_trans_dbl(Tr) :: "
1289 "Not enough memory !");
1290
1291 mydel = new double[rows*(rows+columns+2)];
1292 memory_handler(mydel, "bigfloat_lattice_gensys", "Tr_lll_trans_var_dbl(Bi, Tr, Factor) :: "
1293 "not enough memory !");
1294 Tdel = new bigint[rows*rows+columns];
1295 memory_handler(Tdel, "bigfloat_lattice_gensys", "Tr_lll_trans_dbl(Tr) :: "
1296 "Not enough memory !");
1297
1298 for (lidia_size_t i = 0; i < rows; i++) {
1299 my[i] = &mydel[i*rows]; // rows x rows
1300 T[i] = &Tdel[i*rows];
1301 T[i][i].assign_one();
1302 lattice[i] = &mydel[(rows*rows)+(i*columns)]; // rows x columns
1303 }
1304
1305 //
1306 // Allocating Memory for vector
1307 //
1308
1309 tempvect1 = &Tdel[rows*rows];
1310 vect = &mydel[(rows+columns)*rows];
1311 vect_sq = &vect[rows];
1312
1313 //
1314 // Start
1315 //
1316 // Copy into floating point
1317 //
1318 for (lidia_size_t i = 0; i < rows; i++)
1319 for (lidia_size_t j = 0; j < columns; j++) {
1320 bi_bit_len = BiAddr[i][j].bit_length();
1321 if (bi_bit_len > 8*SIZEOF_DOUBLE) {
1322 bit_diff = bi_bit_len-8*SIZEOF_DOUBLE;
1323 shift_right(tempbiz, BiAddr[i][j], bit_diff);
1324 tempmz0.assign(tempbiz);
1325 }
1326 else {
1327 bit_diff = 0;
1328 tempmz0.assign(BiAddr[i][j]);
1329 }
1330 shift_right(tempmz0, tempmz0, bit_len-bit_diff);
1331 tempmz0.doublify(lattice[i][j]);
1332 }
1333
1334 k = 1;
1335 reduction_steps = 0;
1336 swapping_steps = 0;
1337 correction_steps = 0;
1338 vectsize = columns;
1339 Fc = false;
1340 bfl_assign_zero_dbl(vect);
1341 bfl_scalprod_dbl(vect[0], lattice[0], lattice[0]);
1342
1343 while (k < rows) {
1344 Fc = true;
1345 repeat_loop = false;
1346 while (Fc) {
1347 //
1348 // begin of orthogonalization
1349 //
1350 debug_handler("bigfloat_lattice_gensys", "Tr_lll_trans_var_dbl(Bi, Tr, Factor) [2]");
1351 bfl_scalprod_dbl(vect[k], lattice[k], lattice[k]);
1352 vect_sq[k] = std::sqrt(vect[k]);
1353 for (lidia_size_t j = 0; j < k; ++j) {
1354 bfl_scalprod_dbl(tempdbl0, lattice[k], lattice[j]);
1355 //
1356 // Check Precision, correct if nessesary
1357 //
1358 if (std::fabs(tempdbl0) < invbound*vect_sq[k]*vect_sq[j]) {
1359 bfl_scalprod_bin(tempbiz, BiAddr[k], BiAddr[j]);
1360 bi_bit_len = tempbiz.bit_length();
1361 if (bi_bit_len > 8*SIZEOF_DOUBLE) {
1362 bit_diff = bi_bit_len-8*SIZEOF_DOUBLE;
1363 shift_right(tempbiz, tempbiz, bit_diff);
1364 tempmz0.assign(tempbiz);
1365 }
1366 else {
1367 bit_diff = 0;
1368 tempmz0.assign(tempbiz);
1369 }
1370 shift_right(tempmz0, tempmz0, 2*bit_len-bit_diff);
1371 tempmz0.doublify(my[j][k]);
1372 }
1373 else
1374 my[j][k] = tempdbl0;
1375
1376 for (lidia_size_t i = 0; i < j; ++i)
1377 my[j][k] -= my[i][j]*my[i][k]*vect[i];
1378
1379 tempdbl0 = my[j][k];
1380 my[j][k] /= vect[j];
1381 vect[k] -= my[j][k]*tempdbl0;
1382 }
1383 //
1384 // end of orthogonalization
1385 //
1386 Fc = false;
1387 korr = false;
1388 //
1389 // begin of reduction
1390 //
1391 debug_handler("bigfloat_lattice_gensys", "Tr_lll_trans_var_dbl(Bi, Tr, Factor) [3]");
1392 if (!repeat_loop) {
1393 for (lidia_size_t j = k-1; j >= 0; j--)
1394 if (std::fabs(my[j][k]) > 0.5) {
1395 ++reduction_steps;
1396 korr = true;
1397 tempdbl0 = rint(my[j][k]);
1398 if (std::fabs(tempdbl0) > bound)
1399 Fc = true;
1400
1401 tempmz0.assign(tempdbl0);
1402 tempmz0.bigintify(tempbiz);
1403 bfl_scalmul_bin(tempvect1, tempbiz, BiAddr[j]);
1404 bfl_subtract_bin(BiAddr[k], BiAddr[k], tempvect1);
1405 vectsize = rows;
1406 bfl_scalmul_bin(tempvect1, tempbiz, T[j]);
1407 bfl_subtract_bin(T[k], T[k], tempvect1);
1408 vectsize = columns;
1409
1410 for (lidia_size_t i = 0; i < j; ++i)
1411 my[i][k] -= tempdbl0*my[i][j];
1412
1413 my[j][k] -= tempdbl0;
1414 }
1415 }
1416
1417 //
1418 // correction (not in paper)
1419 //
1420
1421 if (korr) {
1422 ++correction_steps;
1423 for (lidia_size_t j = 0; j < columns; j++) {
1424 bi_bit_len = BiAddr[k][j].bit_length();
1425 if (bi_bit_len > 8*SIZEOF_DOUBLE) {
1426 bit_diff = bi_bit_len-8*SIZEOF_DOUBLE;
1427 shift_right(tempbiz, BiAddr[k][j], bit_diff);
1428 tempmz0.assign(tempbiz);
1429 }
1430 else {
1431 bit_diff = 0;
1432 tempmz0.assign(BiAddr[k][j]);
1433 }
1434 shift_right(tempmz0, tempmz0, bit_len-bit_diff);
1435 tempmz0.doublify(lattice[k][j]);
1436 }
1437 }
1438
1439 //
1440 // end of reduction
1441 //
1442 if (Fc) {
1443 if (k > 1)
1444 --k;
1445 }
1446 else {
1447 Fc = (!repeat_loop) && korr;
1448 repeat_loop = !repeat_loop;
1449 }
1450 }
1451
1452 debug_handler("bigfloat_lattice_gensys", "Tr_lll_trans_var_dbl(Bi, Tr, Factor) [4]");
1453 if (y_par*vect[k-1] > my[k-1][k]*my[k-1][k]*vect[k-1]+vect[k]) {
1454 ++swapping_steps;
1455 bfl_swap_bin(BiAddr[k], BiAddr[k-1]);
1456 bfl_swap_dbl(lattice[k], lattice[k-1]);
1457 bfl_swap_bin(T[k], T[k-1]);
1458
1459 if (k > 1)
1460 --k;
1461 if (k == 1) {
1462 bfl_scalprod_dbl(vect[0], lattice[0], lattice[0]);
1463 vect_sq[0] = std::sqrt(vect[0]);
1464 }
1465 }
1466 else
1467 ++k;
1468 }
1469 //
1470 // restore precision and
1471 // generate bigfloats
1472 //
1473 bigfloat::set_precision(old_prec);
1474 for (lidia_size_t i = 0; i < rows; i++)
1475 for (lidia_size_t j = 0; j < columns; j++) {
1476 bi_bit_len = BiAddr[i][j].bit_length();
1477 if (bi_bit_len > old_bit_prec) {
1478 bit_diff = bi_bit_len-old_bit_prec;
1479 shift_right(tempbiz, BiAddr[i][j], bit_diff);
1480 tempmz0.assign(tempbiz);
1481 }
1482 else {
1483 bit_diff = 0;
1484 tempmz0.assign(BiAddr[i][j]);
1485 }
1486 shift_right(value[i][j], tempmz0, bit_len-bit_diff);
1487 }
1488
1489 //
1490 // Copy into math_matrix< bigint > Tr
1491 //
1492 Tr.set_no_of_rows(rows);
1493 Tr.set_no_of_columns(rows);
1494 TrAddr = Tr.get_data_address();
1495 for (lidia_size_t i = 0; i < rows; i++)
1496 for (lidia_size_t j = 0; j < rows; j++)
1497 if (trans_flag)
1498 LiDIA::swap(TrAddr[j][i], T[i][j]);
1499 else
1500 LiDIA::swap(TrAddr[j][i], T[j][i]);
1501
1502 //
1503 // Free allocated storage
1504 //
1505 delete[] Tdel;
1506 delete[] T;
1507 delete[] mydel;
1508 delete[] my;
1509 }
1510
1511 //
1512 // Schnorr - Euchner version of lll optimized for bigfloats using doubles
1513 // result : lll - reduced lattice for parameter y
1514 //
Tr_lll_var_dbl_gensys(math_matrix<bigint> & Bi,sdigit bit_len,lidia_size_t & rank)1515 void bigfloat_lattice_gensys::Tr_lll_var_dbl_gensys(math_matrix< bigint > & Bi,
1516 sdigit bit_len,
1517 lidia_size_t& rank)
1518 {
1519 debug_handler("bigfloat_lattice_gensys", "Tr_lll_var_dbl_gensys(Bi, Factor, bit_len, rank) [1]");
1520 bool Fc;
1521 bool korr;
1522 bool break_flag;
1523 bool vector_is_zero;
1524 bool repeat_loop;
1525 lidia_size_t k;
1526 lidia_size_t mom_rows = rows;
1527 sdigit old_prec;
1528 sdigit old_bit_prec;
1529 sdigit bit_diff;
1530 sdigit bi_bit_len;
1531 double bound = std::exp(std::log(2.0)*52/2);
1532 double invbound = 1/bound;
1533 double tempdbl0;
1534 double *mydel;
1535 double **my;
1536 double **lattice;
1537 double *vect;
1538 double *vect_sq;
1539 bigint tempbiz;
1540 bigint **BiAddr;
1541 bigint *tempvect1;
1542
1543 vectsize = columns;
1544 old_prec = bigfloat::get_precision();
1545 old_bit_prec = static_cast<sdigit>(static_cast<double>(old_prec)/std::log(2.0)*std::log(10.0)+1);
1546 bigfloat::set_precision(20);
1547 BiAddr = Bi.get_data_address();
1548 //
1549 // Allocating Memory for matrix
1550 //
1551 my = new double*[2*rows];
1552 memory_handler(my, "bigfloat_lattice_gensys", "Tr_lll_var_dbl_gensys(Bi, Factor) :: "
1553 "not enough memory !");
1554 lattice = &my[rows];
1555 mydel = new double[rows*(rows+columns+2)];
1556 memory_handler(mydel, "bigfloat_lattice_gensys", "Tr_lll_var_dbl_gensys(Bi, Factor) :: "
1557 "not enough memory !");
1558 for (lidia_size_t i = 0; i < rows; i++) {
1559 my[i] = &mydel[i*rows]; // rows x rows
1560 lattice[i] = &mydel[(rows*rows)+(i*columns)]; // rows x columns
1561 }
1562
1563 //
1564 // Allocating Memory for vector
1565 //
1566
1567 tempvect1 = new bigint[columns];
1568 memory_handler(tempvect1, "bigfloat_lattice_gensys", "Tr_lll_var_dbl_gensys(Bi, Factor) :: "
1569 "not enough memory !");
1570 vect = &mydel[(rows+columns)*rows];
1571 vect_sq = &vect[rows];
1572
1573 //
1574 // Start
1575 //
1576 // Copy into floating point
1577 //
1578 for (lidia_size_t i = 0; i < rows; i++)
1579 for (lidia_size_t j = 0; j < columns; j++) {
1580 bi_bit_len = BiAddr[i][j].bit_length();
1581 if (bi_bit_len > 8*SIZEOF_DOUBLE) {
1582 bit_diff = bi_bit_len-8*SIZEOF_DOUBLE;
1583 shift_right(tempbiz, BiAddr[i][j], bit_diff);
1584 tempmz0.assign(tempbiz);
1585 }
1586 else {
1587 bit_diff = 0;
1588 tempmz0.assign(BiAddr[i][j]);
1589 }
1590 shift_right(tempmz0, tempmz0, bit_len-bit_diff);
1591 tempmz0.doublify(lattice[i][j]);
1592 }
1593
1594 k = 1;
1595 reduction_steps = 0;
1596 swapping_steps = 0;
1597 correction_steps = 0;
1598 vectsize = columns;
1599 break_flag = false;
1600 vector_is_zero = false;
1601 bfl_assign_zero_dbl(vect);
1602 bfl_scalprod_dbl(vect[0], lattice[0], lattice[0]);
1603
1604 while (k < mom_rows) {
1605 Fc = true;
1606 repeat_loop = false;
1607 while (Fc) {
1608 //
1609 // begin of orthogonalization
1610 //
1611 debug_handler("bigfloat_lattice_gensys", "Tr_lll_var_dbl_gensys(Bi, Factor) [2]");
1612 bfl_scalprod_dbl(vect[k], lattice[k], lattice[k]);
1613 vect_sq[k] = std::sqrt(vect[k]);
1614 for (lidia_size_t j = 0; j < k; ++j) {
1615 bfl_scalprod_dbl(tempdbl0, lattice[k], lattice[j]);
1616 //
1617 // Check Precision, correct if nessesary
1618 //
1619 if (std::fabs(tempdbl0) < invbound*vect_sq[k]*vect_sq[j]) {
1620 bfl_scalprod_bin(tempbiz, BiAddr[k], BiAddr[j]);
1621 bi_bit_len = tempbiz.bit_length();
1622 if (bi_bit_len > 8*SIZEOF_DOUBLE) {
1623 bit_diff = bi_bit_len-8*SIZEOF_DOUBLE;
1624 shift_right(tempbiz, tempbiz, bit_diff);
1625 tempmz0.assign(tempbiz);
1626 }
1627 else {
1628 bit_diff = 0;
1629 tempmz0.assign(tempbiz);
1630 }
1631 shift_right(tempmz0, tempmz0, 2*bit_len-bit_diff);
1632 tempmz0.doublify(my[j][k]);
1633 }
1634 else
1635 my[j][k] = tempdbl0;
1636
1637 for (lidia_size_t i = 0; i < j; ++i)
1638 my[j][k] -= my[i][j]*my[i][k]*vect[i];
1639
1640 tempdbl0 = my[j][k];
1641 my[j][k] /= vect[j];
1642 vect[k] -= my[j][k]*tempdbl0;
1643 }
1644 //
1645 // end of orthogonalization
1646 //
1647 Fc = false;
1648 korr = false;
1649 //
1650 // begin of reduction
1651 //
1652 debug_handler("bigfloat_lattice_gensys", "Tr_lll_var_dbl_gensys(Bi, Factor) [3]");
1653 if (!repeat_loop) {
1654 for (lidia_size_t j = k-1; j >= 0; j--)
1655 if (std::fabs(my[j][k]) > 0.5) {
1656 ++reduction_steps;
1657 korr = true;
1658 tempdbl0 = rint(my[j][k]);
1659 if (std::fabs(tempdbl0) > bound)
1660 Fc = true;
1661
1662 tempmz0.assign(tempdbl0);
1663 tempmz0.bigintify(tempbiz);
1664 bfl_scalmul_bin(tempvect1, tempbiz, BiAddr[j]);
1665 bfl_subtract_bin(BiAddr[k], BiAddr[k], tempvect1);
1666
1667 for (lidia_size_t i = 0; i < j; ++i) {
1668 my[i][k] -= tempdbl0*my[i][j];
1669 }
1670
1671 my[j][k] -= tempdbl0;
1672
1673 }
1674 }
1675
1676 //
1677 // correction (not in paper)
1678 //
1679
1680 if (korr) {
1681 ++correction_steps;
1682 for (lidia_size_t j = 0; j < columns; j++) {
1683 if (!(BiAddr[k][j].is_zero())) {
1684 bi_bit_len = BiAddr[k][j].bit_length();
1685 if (bi_bit_len > 8*SIZEOF_DOUBLE) {
1686 bit_diff = bi_bit_len-8*SIZEOF_DOUBLE;
1687 shift_right(tempbiz, BiAddr[k][j], bit_diff);
1688 tempmz0.assign(tempbiz);
1689 }
1690 else {
1691 bit_diff = 0;
1692 tempmz0.assign(BiAddr[k][j]);
1693 }
1694 shift_right(tempmz0, tempmz0, bit_len-bit_diff);
1695 tempmz0.doublify(lattice[k][j]);
1696 }
1697 else
1698 lattice[k][j] = 0.0;
1699 }
1700 bfl_scalprod_dbl(vect[k], lattice[k], lattice[k]);
1701 vect_sq[k] = std::sqrt(vect[k]);
1702 }
1703 //
1704 // Check if it`s 0 - vector
1705 // delete then
1706 //
1707
1708 vector_is_zero = true;
1709 for (lidia_size_t j = 0; j < columns; j++)
1710 if (!(BiAddr[k][j].is_zero())) {
1711 vector_is_zero = false;
1712 break;
1713 }
1714
1715 if (vector_is_zero) {
1716 for (lidia_size_t j = k; j < mom_rows-1; j++) {
1717 bfl_swap_dbl(lattice[j], lattice[j+1]);
1718 bfl_swap_bin(BiAddr[j], BiAddr[j+1]);
1719 }
1720 bfl_assign_zero_bin(BiAddr[--mom_rows]);
1721 bfl_assign_zero_dbl(lattice[mom_rows]);
1722 k = 1;
1723 bfl_scalprod_dbl(vect[0], lattice[0], lattice[0]);
1724 vect_sq[0] = std::sqrt(vect[0]);
1725 vector_is_zero = false;
1726 break_flag = true;
1727 korr = false;
1728 Fc = false;
1729 break;
1730 }
1731
1732 //
1733 // end of reduction
1734 //
1735 if (Fc) {
1736 if (k > 1)
1737 --k;
1738 }
1739 else {
1740 Fc = (!repeat_loop) && korr;
1741 repeat_loop = !repeat_loop;
1742 }
1743 }
1744 if (!break_flag) {
1745 debug_handler("bigfloat_lattice_gensys", "Tr_lll_var_dbl_gensys(Bi, Factor) [4]");
1746 if (y_par*vect[k-1] > my[k-1][k]*my[k-1][k]*vect[k-1]+vect[k]) {
1747 ++swapping_steps;
1748 bfl_swap_bin(BiAddr[k], BiAddr[k-1]);
1749 bfl_swap_dbl(lattice[k], lattice[k-1]);
1750
1751 if (k > 1)
1752 --k;
1753 if (k == 1) {
1754 bfl_scalprod_dbl(vect[0], lattice[0], lattice[0]);
1755 vect_sq[0] = std::sqrt(vect[0]);
1756 }
1757 }
1758 else
1759 ++k;
1760 }
1761 break_flag = false;
1762 }
1763 rank = mom_rows;
1764 bigfloat::set_precision(old_prec);
1765 for (lidia_size_t i = 0; i < rows; i++)
1766 for (lidia_size_t j = 0; j < columns; j++) {
1767 bi_bit_len = BiAddr[i][j].bit_length();
1768 if (bi_bit_len > old_bit_prec) {
1769 bit_diff = bi_bit_len-old_bit_prec;
1770 shift_right(tempbiz, BiAddr[i][j], bit_diff);
1771 tempmz0.assign(tempbiz);
1772 }
1773 else {
1774 bit_diff = 0;
1775 tempmz0.assign(BiAddr[i][j]);
1776 }
1777 shift_right(value[i][j], tempmz0, bit_len-bit_diff);
1778 }
1779
1780 //
1781 // Free allocated storage
1782 // and restore precision
1783 //
1784 delete[] tempvect1;
1785 delete[] mydel;
1786 delete[] my;
1787 }
1788
1789
1790 //
1791 // Schnorr - Euchner version of lll optimized for bigfloats using doubles
1792 // result : lll - reduced lattice for parameter y
1793 //
Tr_lll_trans_var_dbl_gensys(math_matrix<bigint> & Bi,math_matrix<bigint> & Tr,sdigit bit_len,lidia_size_t & rank)1794 void bigfloat_lattice_gensys::Tr_lll_trans_var_dbl_gensys(math_matrix< bigint > & Bi,
1795 math_matrix< bigint > & Tr,
1796 sdigit bit_len,
1797 lidia_size_t& rank)
1798 {
1799 debug_handler("bigfloat_lattice_gensys", "Tr_lll_trans_var_dbl_gensys(Bi, Tr, Factor, bit_len, rank) [1]");
1800 bool Fc;
1801 bool korr;
1802 bool break_flag;
1803 bool vector_is_zero;
1804 bool repeat_loop;
1805 lidia_size_t k;
1806 lidia_size_t mom_rows = rows;
1807 sdigit old_prec;
1808 sdigit old_bit_prec;
1809 sdigit bit_diff;
1810 sdigit bi_bit_len;
1811 double bound = std::exp(std::log(2.0)*52/2);
1812 double invbound = 1/bound;
1813 double tempdbl0;
1814 double *mydel;
1815 double **my;
1816 double **lattice;
1817 double *vect;
1818 double *vect_sq;
1819 bigint tempbiz;
1820 bigint *Tdel;
1821 bigint **T;
1822 bigint **TrAddr;
1823 bigint **BiAddr;
1824 bigint *tempvect1;
1825
1826 vectsize = columns;
1827 old_prec = bigfloat::get_precision();
1828 old_bit_prec = static_cast<sdigit>(static_cast<double>(old_prec)/std::log(2.0)*std::log(10.0)+1);
1829 bigfloat::set_precision(20);
1830 BiAddr = Bi.get_data_address();
1831 //
1832 // Allocating Memory for matrix
1833 //
1834 my = new double*[2*rows];
1835 memory_handler(my, "bigfloat_lattice_gensys", "Tr_lll_trans_var_dbl_gensys(Bi, Tr, Factor, rank) :: "
1836 "not enough memory !");
1837 lattice = &my[rows];
1838 T = new bigint*[rows];
1839 memory_handler(T, "bigfloat_lattice_gensys", "Tr_lll_trans_var_dbl_gensys(Bi, Tr, Factor, rank) :: "
1840 "Not enough memory !");
1841
1842 mydel = new double[rows*(rows+columns+2)];
1843 memory_handler(mydel, "bigfloat_lattice_gensys", "Tr_lll_trans_var_dbl_gensys(Bi, Tr, Factor, rank) :: "
1844 "not enough memory !");
1845 Tdel = new bigint[rows*rows+((rows > columns)?rows:columns)];
1846 memory_handler(Tdel, "bigfloat_lattice_gensys", "Tr_lll_trans_var_dbl_gensys(Bi, Tr, Factor, rank) :: "
1847 "Not enough memory !");
1848
1849 for (lidia_size_t i = 0; i < rows; i++) {
1850 my[i] = &mydel[i*rows]; // rows x rows
1851 T[i] = &Tdel[i*rows];
1852 T[i][i].assign_one();
1853 lattice[i] = &mydel[(rows*rows)+(i*columns)]; // rows x columns
1854 }
1855
1856 //
1857 // Allocating Memory for vector
1858 //
1859
1860 tempvect1 = &Tdel[rows*rows];
1861 vect = &mydel[(rows+columns)*rows];
1862 vect_sq = &vect[rows];
1863
1864 //
1865 // Start
1866 //
1867 // Copy into floating point
1868 //
1869 for (lidia_size_t i = 0; i < rows; i++)
1870 for (lidia_size_t j = 0; j < columns; j++) {
1871 bi_bit_len = BiAddr[i][j].bit_length();
1872 if (bi_bit_len > 8*SIZEOF_DOUBLE) {
1873 bit_diff = bi_bit_len-8*SIZEOF_DOUBLE;
1874 shift_right(tempbiz, BiAddr[i][j], bit_diff);
1875 tempmz0.assign(tempbiz);
1876 }
1877 else {
1878 bit_diff = 0;
1879 tempmz0.assign(BiAddr[i][j]);
1880 }
1881 shift_right(tempmz0, tempmz0, bit_len-bit_diff);
1882 tempmz0.doublify(lattice[i][j]);
1883 }
1884
1885 k = 1;
1886 reduction_steps = 0;
1887 swapping_steps = 0;
1888 correction_steps = 0;
1889 vectsize = columns;
1890 break_flag = false;
1891 vector_is_zero = false;
1892 bfl_assign_zero_dbl(vect);
1893 bfl_scalprod_dbl(vect[0], lattice[0], lattice[0]);
1894
1895 while (k < mom_rows) {
1896 Fc = true;
1897 repeat_loop = false;
1898 while (Fc) {
1899 //
1900 // begin of orthogonalization
1901 //
1902 debug_handler("bigfloat_lattice_gensys", "Tr_lll_trans_var_dbl_gensys(Bi, Tr, Factor, rank) [2]");
1903 bfl_scalprod_dbl(vect[k], lattice[k], lattice[k]);
1904 vect_sq[k] = std::sqrt(vect[k]);
1905 for (lidia_size_t j = 0; j < k; ++j) {
1906 bfl_scalprod_dbl(tempdbl0, lattice[k], lattice[j]);
1907 //
1908 // Check Precision, correct if nessesary
1909 //
1910 if (std::fabs(tempdbl0) < invbound*vect_sq[k]*vect_sq[j]) {
1911 bfl_scalprod_bin(tempbiz, BiAddr[k], BiAddr[j]);
1912 bi_bit_len = tempbiz.bit_length();
1913 if (bi_bit_len > 8*SIZEOF_DOUBLE) {
1914 bit_diff = bi_bit_len-8*SIZEOF_DOUBLE;
1915 shift_right(tempbiz, tempbiz, bit_diff);
1916 tempmz0.assign(tempbiz);
1917 }
1918 else {
1919 bit_diff = 0;
1920 tempmz0.assign(tempbiz);
1921 }
1922 shift_right(tempmz0, tempmz0, 2*bit_len-bit_diff);
1923 tempmz0.doublify(my[j][k]);
1924 }
1925 else
1926 my[j][k] = tempdbl0;
1927
1928 for (lidia_size_t i = 0; i < j; ++i)
1929 my[j][k] -= my[i][j]*my[i][k]*vect[i];
1930
1931 tempdbl0 = my[j][k];
1932 my[j][k] /= vect[j];
1933 vect[k] -= my[j][k]*tempdbl0;
1934 }
1935 //
1936 // end of orthogonalization
1937 //
1938 Fc = false;
1939 korr = false;
1940 //
1941 // begin of reduction
1942 //
1943 debug_handler("bigfloat_lattice_gensys", "Tr_lll_trans_var_dbl_gensys(Bi, Tr, Factor, rank) [3]");
1944 if (!repeat_loop) {
1945 for (lidia_size_t j = k-1; j >= 0; j--)
1946 if (std::fabs(my[j][k]) > 0.5) {
1947 ++reduction_steps;
1948 korr = true;
1949 tempdbl0 = rint(my[j][k]);
1950 if (std::fabs(tempdbl0) > bound)
1951 Fc = true;
1952
1953 tempmz0.assign(tempdbl0);
1954 tempmz0.bigintify(tempbiz);
1955 bfl_scalmul_bin(tempvect1, tempbiz, BiAddr[j]);
1956 bfl_subtract_bin(BiAddr[k], BiAddr[k], tempvect1);
1957 vectsize = rows;
1958 bfl_scalmul_bin(tempvect1, tempbiz, T[j]);
1959 bfl_subtract_bin(T[k], T[k], tempvect1);
1960 vectsize = columns;
1961
1962 for (lidia_size_t i = 0; i < j; ++i)
1963 my[i][k] -= tempdbl0*my[i][j];
1964
1965 my[j][k] -= tempdbl0;
1966 }
1967 }
1968
1969 //
1970 // correction (not in paper)
1971 //
1972
1973 if (korr) {
1974 ++correction_steps;
1975 for (lidia_size_t j = 0; j < columns; j++) {
1976 if (!(BiAddr[k][j].is_zero())) {
1977 bi_bit_len = BiAddr[k][j].bit_length();
1978 if (bi_bit_len > 8*SIZEOF_DOUBLE) {
1979 bit_diff = bi_bit_len-8*SIZEOF_DOUBLE;
1980 shift_right(tempbiz, BiAddr[k][j], bit_diff);
1981 tempmz0.assign(tempbiz);
1982 }
1983 else {
1984 bit_diff = 0;
1985 tempmz0.assign(BiAddr[k][j]);
1986 }
1987 shift_right(tempmz0, tempmz0, bit_len-bit_diff);
1988 tempmz0.doublify(lattice[k][j]);
1989 }
1990 else
1991 lattice[k][j] = 0.0;
1992 }
1993 }
1994 //
1995 // Check if it`s 0 - vector
1996 // delete then
1997 //
1998 vector_is_zero = true;
1999 for (lidia_size_t j = 0; j < columns; j++)
2000 if (!(BiAddr[k][j].is_zero())) {
2001 vector_is_zero = false;
2002 break;
2003 }
2004
2005 if (vector_is_zero) {
2006 for (lidia_size_t j = k; j < mom_rows-1; j++) {
2007 bfl_swap_dbl(lattice[j], lattice[j+1]);
2008 bfl_swap_bin(BiAddr[j], BiAddr[j+1]);
2009 bfl_swap_bin(T[j], T[j+1]);
2010 }
2011 bfl_assign_zero_bin(BiAddr[--mom_rows]);
2012 bfl_assign_zero_dbl(lattice[mom_rows]);
2013 vectsize = rows;
2014 bfl_assign_zero_bin(T[mom_rows]);
2015 vectsize = columns;
2016 k = 1;
2017 bfl_scalprod_dbl(vect[0], lattice[0], lattice[0]);
2018 vect_sq[k] = std::sqrt(vect[0]);
2019 vector_is_zero = false;
2020 break_flag = true;
2021 Fc = false;
2022 korr = false;
2023 break;
2024 }
2025
2026 //
2027 // end of reduction
2028 //
2029 if (Fc) {
2030 if (k > 1)
2031 --k;
2032 }
2033 else {
2034 Fc = (!repeat_loop) && korr;
2035 repeat_loop = !repeat_loop;
2036 }
2037 }
2038 if (!break_flag) {
2039 debug_handler("bigfloat_lattice_gensys", "Tr_lll_trans_var_dbl_gensys(Bi, Tr, Factor, rank) [4]");
2040 if (y_par*vect[k-1] > my[k-1][k]*my[k-1][k]*vect[k-1]+vect[k]) {
2041 ++swapping_steps;
2042 bfl_swap_bin(BiAddr[k], BiAddr[k-1]);
2043 bfl_swap_dbl(lattice[k], lattice[k-1]);
2044 bfl_swap_bin(T[k], T[k-1]);
2045
2046 if (k > 1)
2047 --k;
2048 if (k == 1) {
2049 bfl_scalprod_dbl(vect[0], lattice[0], lattice[0]);
2050 vect_sq[0] = std::sqrt(vect[0]);
2051 }
2052 }
2053 else
2054 ++k;
2055 }
2056 break_flag = false;
2057 }
2058 rank = mom_rows;
2059 //
2060 // restore precision and
2061 // generate bigfloats
2062 //
2063 bigfloat::set_precision(old_prec);
2064 for (lidia_size_t i = 0; i < rows; i++)
2065 for (lidia_size_t j = 0; j < columns; j++) {
2066 bi_bit_len = BiAddr[i][j].bit_length();
2067 if (bi_bit_len > old_bit_prec) {
2068 bit_diff = bi_bit_len-old_bit_prec;
2069 shift_right(tempbiz, BiAddr[i][j], bit_diff);
2070 tempmz0.assign(tempbiz);
2071 }
2072 else {
2073 bit_diff = 0;
2074 tempmz0.assign(BiAddr[i][j]);
2075 }
2076 shift_right(value[i][j], tempmz0, bit_len-bit_diff);
2077 }
2078
2079 //
2080 // Copy into math_matrix< bigint > Tr
2081 //
2082 Tr.set_no_of_rows(rows);
2083 Tr.set_no_of_columns(rows);
2084 TrAddr = Tr.get_data_address();
2085 for (lidia_size_t i = 0; i < rows; i++)
2086 for (lidia_size_t j = 0; j < rows; j++)
2087 if (trans_flag)
2088 LiDIA::swap(TrAddr[j][i], T[i][j]);
2089 else
2090 LiDIA::swap(TrAddr[j][i], T[j][i]);
2091
2092 //
2093 // Free allocated storage
2094 //
2095 delete[] Tdel;
2096 delete[] T;
2097 delete[] mydel;
2098 delete[] my;
2099 }
2100
2101
2102
2103 #ifdef LIDIA_NAMESPACE
2104 } // end of namespace LiDIA
2105 #endif
2106