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 small bigfloats for approximation
34 //
35 //
36 // Schnorr - Euchner
37 // * Tr_lll_bfl()
38 // * Tr_lll_var_bfl();
39 // * Tr_lll_trans_bfl(T)
40 // * Tr_lll_trans_var_bfl(T);
41 // * Tr_lll_deep_insert_bfl() (not yet implemented)
42 // * Tr_lll_bfl_gensys()
43 // * Tr_lll_var_bfl_gensys()
44 // * Tr_lll_trans_bfl_gensys(T)
45 // * Tr_lll_trans_var_bfl_gensys(T)
46 //
47 // Modified lll
48 // * Tr_mlll_bfl(v)
49 //
50
51
52 //
53 // Schnorr - Euchner version of lll optimized for bigfloats using small bigfloats
54 // result : lll - reduced lattice for parameter y
55 //
56
57 //
58 // only high precision if needed
59 //
Tr_lll_bfl(sdigit bit_prec)60 void bigfloat_lattice_gensys::Tr_lll_bfl(sdigit bit_prec)
61 {
62 debug_handler("bigfloat_lattice_gensys", "Tr_lll_bfl(bit_prec) [1]");
63 bool Fc;
64 bool korr;
65 bool repeat_loop;
66 lidia_size_t k;
67 sdigit old_prec;
68 sdigit short_prec;
69 sdigit r_prec;
70 sdigit n_prec;
71 sdigit bi_bit_len;
72 sdigit bit_diff;
73 sdigit x_exponent;
74 bigint x_factor;
75 bigfloat *tempvect1;
76 bigfloat bound = std::exp(std::log(2.0)*bit_prec/2);
77 bigfloat invbound;
78 bigfloat *mydel;
79 bigfloat **my;
80 bigfloat **lattice;
81 bigfloat *vect;
82 bigfloat *vect_sq;
83 bigfloat conv;
84 bigfloat halb(0.5);
85
86 LiDIA::divide(invbound, 1, bound);
87
88 //
89 // Compute needed precision for exact arithmetik
90 //
91 old_prec = bigfloat::get_precision();
92 short_prec = static_cast<sdigit>(static_cast<double>(bit_prec)*std::log(2.0)/std::log(10.0));
93 r_prec = compute_read_precision();
94 n_prec = static_cast<sdigit>(r_prec/std::log(10.0))+1;
95 n_prec *= columns; // columns >= rows, basis !!!
96 n_prec = static_cast<sdigit>(n_prec*std::log(10.0))+1;
97
98 vectsize = columns;
99
100 //
101 // Allocating Memory for matrix
102 //
103 my = new bigfloat*[rows*2];
104 memory_handler(my, "bigfloat_lattice_gensys", "Tr_lll_bfl(bit_prec) :: "
105 "not enough memory !");
106 lattice = &my[rows];
107 //
108 // think again
109 //
110 mydel = new bigfloat[rows*(rows+columns+2)+columns];
111 memory_handler(mydel, "bigfloat_lattice_gensys", "Tr_lll_bfl(bit_prec) :: "
112 "not enough memory !");
113 for (lidia_size_t i = 0; i < rows; i++) {
114 my[i] = &mydel[i*rows]; // rows x rows
115 lattice[i] = &mydel[(rows*rows)+(i*columns)]; // rows x columns
116 }
117
118 //
119 // Allocating Memory for vector
120 //
121
122 //
123 // think again Part II
124 //
125 tempvect1 = new bigfloat[columns];
126 memory_handler(tempvect1, "bigfloat_lattice_gensys", "Tr_lll_bfl(bit_prec) :: "
127 "not enough memory !");
128
129 vect = &mydel[(rows+columns)*rows];
130 vect_sq = &vect[rows];
131
132
133 //
134 // Start
135 //
136 // Copy into floating point
137 //
138 bigfloat::set_precision(short_prec);
139 for (lidia_size_t i = 0; i < rows; i++)
140 for (lidia_size_t j = 0; j < columns; j++) {
141 x_factor = value[i][j].mantissa();
142 x_exponent = value[i][j].exponent();
143 bi_bit_len = x_factor.bit_length();
144 if (bi_bit_len > bit_prec) {
145 bit_diff = bi_bit_len-bit_prec;
146 shift_right(x_factor, x_factor, bit_diff);
147 lattice[i][j].assign(x_factor);
148 shift_right(lattice[i][j], lattice[i][j], bit_diff);
149 }
150 else
151 lattice[i][j].assign(x_factor);
152 shift_left(lattice[i][j], lattice[i][j], x_exponent);
153 }
154
155 k = 1;
156 reduction_steps = 0;
157 swapping_steps = 0;
158 correction_steps = 0;
159 vectsize = columns;
160 bfl_assign_zero_bfl(vect);
161 bfl_scalprod_bfl(vect[0], lattice[0], lattice[0]);
162
163 while (k < rows) {
164 Fc = true;
165 repeat_loop = false;
166 while (Fc) {
167 //
168 // begin of orthogonalization
169 //
170 debug_handler("bigfloat_lattice_gensys", "Tr_lll_bfl(bit_prec) [2]");
171 bfl_scalprod_bfl(vect[k], lattice[k], lattice[k]);
172 sqrt(vect_sq[k], vect[k]);
173 for (lidia_size_t j = 0; j < k; ++j) {
174 bfl_scalprod_bfl(tempmz0, lattice[k], lattice[j]);
175 LiDIA::multiply(tempmz1, invbound, vect_sq[k]);
176 LiDIA::multiply(tempmz1, tempmz1, vect_sq[j]);
177 //
178 // Correction Step I, only if needed, |lattice[k]*lattice[j]| to big
179 //
180 if (abs(tempmz0).compare(tempmz1) < 0) {
181 //
182 // Compute the correct scalar product
183 //
184 bigfloat::set_precision(n_prec);
185 bfl_scalprod_bfl(tempmz0, value[k], value[j]);
186 x_factor = tempmz0.mantissa();
187 x_exponent = tempmz0.exponent();
188 bi_bit_len = x_factor.bit_length();
189 if (bi_bit_len > bit_prec) {
190 bit_diff = bi_bit_len-bit_prec;
191 shift_right(x_factor, x_factor, bit_diff);
192 my[j][k].assign(x_factor);
193 shift_left(my[j][k], my[j][k], bit_diff);
194 }
195 else
196 my[j][k].assign(x_factor);
197 shift_left(my[j][k], my[j][k], x_exponent);
198 bigfloat::set_precision(short_prec);
199 }
200 else
201 my[j][k].assign(tempmz0);
202
203 for (lidia_size_t i = 0; i < j; ++i) {
204 LiDIA::multiply(tempmz1, my[i][j], my[i][k]);
205 LiDIA::multiply(tempmz1, tempmz1, vect[i]);
206 LiDIA::subtract(my[j][k], my[j][k], tempmz1);
207 }
208
209 tempmz0.assign(my[j][k]);
210 LiDIA::divide(my[j][k], my[j][k], vect[j]);
211 LiDIA::multiply(tempmz1, my[j][k], tempmz0);
212 LiDIA::subtract(vect[k], vect[k], tempmz1);
213 }
214 //
215 // end of orthogonalization
216 //
217 Fc = false;
218 korr = false;
219 //
220 // begin of reduction
221 //
222 debug_handler("bigfloat_lattice_gensys", "Tr_lll_bfl(bit_prec) [3]");
223 if (!repeat_loop) {
224 for (lidia_size_t j = k-1; j >= 0; j--)
225 if (abs(my[j][k]).compare(halb) > 0) {
226 ++reduction_steps;
227 korr = true;
228 round (tempmz0, my[j][k]);
229 if (abs(tempmz0).compare(bound) > 0)
230 Fc = true;
231
232 bigfloat::set_precision(n_prec);
233 bfl_scalmul_bfl(tempvect1, tempmz0, value[j]);
234 bfl_subtract_bfl(value[k], value[k], tempvect1);
235 bigfloat::set_precision(short_prec);
236
237 for (lidia_size_t i = 0; i < j; ++i) {
238 LiDIA::multiply(tempmz1, tempmz0, my[i][j]);
239 LiDIA::subtract(my[i][k], my[i][k], tempmz1);
240 }
241 LiDIA::subtract(my[j][k], my[j][k], tempmz0);
242
243 }
244 }
245
246 //
247 // Correction Step II
248 // lattice=(bigfloat )value
249 //
250 if (korr) {
251 ++correction_steps;
252 for (lidia_size_t j = 0; j < columns; j++) {
253 x_factor = value[k][j].mantissa();
254 x_exponent = value[k][j].exponent();
255 bi_bit_len = x_factor.bit_length();
256 if (bi_bit_len > bit_prec) {
257 bit_diff = bi_bit_len-bit_prec;
258 shift_right(x_factor, x_factor, bit_diff);
259 lattice[k][j].assign(x_factor);
260 shift_left(lattice[k][j], lattice[k][j], bit_diff);
261 }
262 else
263 lattice[k][j].assign(x_factor);
264 shift_left(lattice[k][j], lattice[k][j], x_exponent);
265 }
266 }
267
268 //
269 // end of reduction
270 //
271 if (Fc) {
272 if (k > 1)
273 --k;
274 }
275 else {
276 Fc = (!repeat_loop) && korr;
277 repeat_loop = !repeat_loop;
278 }
279 }
280
281 debug_handler("bigfloat_lattice_gensys", "Tr_lll_bfl(bit_prec) [4]");
282 //
283 // Check lll - condition for parameter y_par
284 //
285 tempmz0.assign(y_par);
286 LiDIA::multiply(tempmz0, tempmz0, vect[k-1]);
287 LiDIA::square(tempmz1, my[k-1][k]);
288 LiDIA::multiply(tempmz1, tempmz1, vect[k-1]);
289 LiDIA::add(tempmz1, tempmz1, vect[k]);
290
291 if (tempmz0.compare(tempmz1) > 0) {
292 //
293 // ::swap vectors k and k-1
294 //
295 ++swapping_steps;
296 bfl_swap_bfl(value[k], value[k-1]);
297 bfl_swap_bfl(lattice[k], lattice[k-1]);
298
299 if (k > 1)
300 --k;
301 if (k == 1) {
302 bfl_scalprod_bfl(vect[0], lattice[0], lattice[0]);
303 sqrt(vect_sq[0], vect[0]);
304 }
305 }
306 else
307 ++k;
308 }
309
310 //
311 // Free allocated storage and
312 // restore precision
313 //
314 bigfloat::set_precision(old_prec);
315 delete[] tempvect1;
316 delete[] mydel;
317 delete[] my;
318 }
319
320
321
322 //
323 // Schnorr - Euchner version of lll optimized for bigints usign bigfloats
324 // result : lll - reduced lattice for parameter y
325 //
326
327 //
328 // only high precision if needed
329 //
Tr_lll_trans_bfl(math_matrix<bigint> & Tr,sdigit bit_prec)330 void bigfloat_lattice_gensys::Tr_lll_trans_bfl(math_matrix< bigint > & Tr, sdigit bit_prec)
331 {
332 debug_handler("bigfloat_lattice_gensys", "Tr_lll_trans_bfl(Tr, bit_prec) [1]");
333 bool Fc;
334 bool korr;
335 bool repeat_loop;
336 lidia_size_t k;
337 sdigit old_prec;
338 sdigit short_prec = static_cast<sdigit>(static_cast<double>(bit_prec)*std::log(2.0)/std::log(10.0));
339 sdigit r_prec;
340 sdigit n_prec;
341 sdigit bi_bit_len;
342 sdigit bit_diff;
343 sdigit x_exponent;
344 bigint x_factor;
345 bigfloat *tempvect1;
346 bigint *tempvect2;
347 bigfloat bound = std::exp(std::log(2.0)*bit_prec/2);
348 bigfloat invbound;
349 bigfloat *mydel;
350 bigfloat **my;
351 bigfloat **lattice;
352 bigint *Tdel;
353 bigint **T;
354 bigint **TrAddr;
355 bigfloat *vect;
356 bigfloat *vect1;
357 bigfloat *vect_sq;
358 bigint conv;
359 bigfloat halb(0.5);
360
361 LiDIA::divide(invbound, 1, bound);
362 //
363 // Compute needed precision for exact arithmetik
364 //
365 old_prec = bigfloat::get_precision();
366 r_prec = compute_read_precision();
367 n_prec = static_cast<sdigit>(r_prec/std::log(10.0))+1;
368 n_prec *= columns; // columns >= rows, basis !!!
369 n_prec = static_cast<sdigit>(n_prec*std::log(10.0))+1;
370
371 vectsize = columns;
372
373 //
374 // Allocating Memory for matrix
375 //
376 my = new bigfloat*[rows*2];
377 memory_handler(my, "bigfloat_lattice_gensys", "Tr_lll_trans_bfl(Tr, bit_prec) :: "
378 "not enough memory !");
379 lattice = &my[rows];
380 //
381 // think again
382 //
383 mydel = new bigfloat[(rows+2)*(rows+columns)-columns];
384 memory_handler(mydel, "bigfloat_lattice_gensys", "Tr_lll_trans_bfl(Tr, bit_prec) :: "
385 "not enough memory !");
386 T = new bigint*[rows];
387 memory_handler(T, "bigfloat_lattice_gensys", "Tr_lll_trans_bfl(Tr, bit_prec) :: "
388 "Not enough memory !");
389
390 Tdel = new bigint[rows*rows+rows];
391 memory_handler(Tdel, "bigfloat_lattice_gensys", "Tr_lll_trans_bfl(Tr, bit_prec) :: "
392 "Not enough memory !");
393
394 for (lidia_size_t i = 0; i < rows; i++) {
395 my[i] = &mydel[i*rows]; // rows x rows
396 lattice[i] = &mydel[(rows*rows)+(i*columns)]; // rows x columns
397 T[i] = &Tdel[i*rows];
398 T[i][i].assign_one();
399 }
400
401 //
402 // Allocating Memory for vector
403 //
404
405 //
406 // think again Part II
407 //
408 tempvect1 = new bigfloat[columns];
409 memory_handler(tempvect1, "bigfloat_lattice_gensys", "Tr_lll_trans_bfl(Tr, bit_prec) :: "
410 "not enough memory !");
411
412 tempvect2 = &Tdel[rows*rows];
413 vect = &mydel[(rows+columns)*rows];
414 vect_sq = &vect[rows];
415 vect1 = &vect[2*rows];
416
417
418 //
419 // Start
420 //
421 // Copy into floating point
422 //
423 bigfloat::set_precision(short_prec);
424 for (lidia_size_t i = 0; i < rows; i++)
425 for (lidia_size_t j = 0; j < columns; j++) {
426 x_factor = value[i][j].mantissa();
427 x_exponent = value[i][j].exponent();
428 bi_bit_len = x_factor.bit_length();
429 if (bi_bit_len > bit_prec) {
430 bit_diff = bi_bit_len-bit_prec;
431 shift_right(x_factor, x_factor, bit_diff);
432 lattice[i][j].assign(x_factor);
433 shift_right(lattice[i][j], lattice[i][j], bit_diff);
434 }
435 else
436 lattice[i][j].assign(x_factor);
437 shift_left(lattice[i][j], lattice[i][j], x_exponent);
438 }
439
440 k = 1;
441 reduction_steps = 0;
442 swapping_steps = 0;
443 correction_steps = 0;
444 vectsize = columns;
445 bfl_assign_zero_bfl(vect);
446 bfl_scalprod_bfl(vect[0], lattice[0], lattice[0]);
447
448 while (k < rows) {
449 Fc = true;
450 repeat_loop = false;
451 while (Fc) {
452 //
453 // begin of orthogonalization
454 //
455 debug_handler("bigfloat_lattice_gensys", "Tr_lll_trans_bfl(Tr, bit_prec) [2]");
456 bfl_scalprod_bfl(vect[k], lattice[k], lattice[k]);
457 sqrt(vect_sq[k], vect[k]);
458 for (lidia_size_t j = 0; j < k; ++j) {
459 bfl_scalprod_bfl(tempmz0, lattice[k], lattice[j]);
460 LiDIA::multiply(tempmz1, invbound, vect_sq[k]);
461 LiDIA::multiply(tempmz1, tempmz1, vect_sq[j]);
462 //
463 // Correction Step I, only if needed, |lattice[k]*lattice[j]| to big
464 //
465 if (abs(tempmz0).compare(tempmz1) < 0) {
466 //
467 // Compute the correct scalar product
468 //
469 bigfloat::set_precision(n_prec);
470 bfl_scalprod_bfl(tempmz0, value[k], value[j]);
471 x_factor = tempmz0.mantissa();
472 x_exponent = tempmz0.exponent();
473 bi_bit_len = x_factor.bit_length();
474 if (bi_bit_len > bit_prec) {
475 bit_diff = bi_bit_len-bit_prec;
476 shift_right(x_factor, x_factor, bit_diff);
477 my[j][k].assign(x_factor);
478 shift_left(my[j][k], my[j][k], bit_diff);
479 }
480 else
481 my[j][k].assign(x_factor);
482 shift_left(my[j][k], my[j][k], x_exponent);
483 bigfloat::set_precision(short_prec);
484 }
485 else
486 my[j][k].assign(tempmz0);
487
488 for (lidia_size_t i = 0; i < j; ++i) {
489 LiDIA::multiply(tempmz1, my[i][j], my[i][k]);
490 LiDIA::multiply(tempmz1, tempmz1, vect[i]);
491 LiDIA::subtract(my[j][k], my[j][k], tempmz1);
492 }
493
494 tempmz0.assign(my[j][k]);
495 LiDIA::divide(my[j][k], my[j][k], vect[j]);
496 LiDIA::multiply(tempmz1, my[j][k], tempmz0);
497 LiDIA::subtract(vect[k], vect[k], tempmz1);
498 }
499 //
500 // end of orthogonalization
501 //
502 Fc = false;
503 korr = false;
504 //
505 // begin of reduction
506 //
507 debug_handler("bigfloat_lattice_gensys", "Tr_lll_trans_bfl(Tr, bit_prec) [3]");
508 if (!repeat_loop) {
509 for (lidia_size_t j = k-1; j >= 0; j--)
510 if (abs(my[j][k]).compare(halb) > 0) {
511 ++reduction_steps;
512 korr = true;
513 round (tempmz0, my[j][k]);
514 if (abs(tempmz0).compare(bound) > 0)
515 Fc = true;
516
517 bigfloat::set_precision(n_prec);
518 bfl_scalmul_bfl(tempvect1, tempmz0, value[j]);
519 bfl_subtract_bfl(value[k], value[k], tempvect1);
520 vectsize = rows;
521 tempmz0.bigintify(conv);
522 bfl_scalmul_bin(tempvect2, conv, T[j]);
523 bfl_subtract_bin(T[k], T[k], tempvect2);
524 vectsize = columns;
525 bigfloat::set_precision(short_prec);
526
527 for (lidia_size_t i = 0; i < j; ++i) {
528 LiDIA::multiply(tempmz1, tempmz0, my[i][j]);
529 LiDIA::subtract(my[i][k], my[i][k], tempmz1);
530 }
531 LiDIA::subtract(my[j][k], my[j][k], tempmz0);
532 }
533 }
534
535 //
536 // Correction Step II
537 // lattice[k]=(bigfloat )value[k]
538 //
539 if (korr) {
540 ++correction_steps;
541 for (lidia_size_t j = 0; j < columns; j++) {
542 x_factor = value[k][j].mantissa();
543 x_exponent = value[k][j].exponent();
544 bi_bit_len = x_factor.bit_length();
545 if (bi_bit_len > bit_prec) {
546 bit_diff = bi_bit_len-bit_prec;
547 shift_right(x_factor, x_factor, bit_diff);
548 lattice[k][j].assign(x_factor);
549 shift_left(lattice[k][j], lattice[k][j], bit_diff);
550 }
551 else
552 lattice[k][j].assign(x_factor);
553 shift_left(lattice[k][j], lattice[k][j], x_exponent);
554 }
555 }
556 //
557 // end of reduction
558 //
559 if (Fc) {
560 if (k > 1)
561 --k;
562 }
563 else {
564 Fc = (!repeat_loop) && korr;
565 repeat_loop = !repeat_loop;
566 }
567 }
568
569 debug_handler("bigfloat_lattice_gensys", "Tr_lll_trans_bfl(Tr, bit_prec) [4]");
570 //
571 // Check lll - condition for parameter y_par
572 //
573 tempmz0.assign(y_par);
574 LiDIA::multiply(tempmz0, tempmz0, vect[k-1]);
575 LiDIA::square(tempmz1, my[k-1][k]);
576 LiDIA::multiply(tempmz1, tempmz1, vect[k-1]);
577 LiDIA::add(tempmz1, tempmz1, vect[k]);
578
579 if (tempmz0.compare(tempmz1) > 0) {
580 //
581 // ::swap vectors k and k-1
582 //
583 ++swapping_steps;
584 bfl_swap_bfl(value[k], value[k-1]);
585 bfl_swap_bfl(lattice[k], lattice[k-1]);
586 bfl_swap_bin(T[k], T[k-1]);
587
588 if (k > 1)
589 --k;
590 if (k == 1) {
591 bfl_scalprod_bfl(vect[0], lattice[0], lattice[0]);
592 sqrt(vect_sq[0], vect[0]);
593 }
594 }
595 else
596 ++k;
597 }
598 //
599 // Copy into math_matrix< bigint > Tr
600 //
601 Tr.set_no_of_rows(rows);
602 Tr.set_no_of_columns(rows);
603 TrAddr = Tr.get_data_address();
604 for (lidia_size_t i = 0; i < rows; i++)
605 for (lidia_size_t j = 0; j < rows; j++)
606 if (trans_flag)
607 LiDIA::swap(TrAddr[j][i], T[i][j]);
608 else
609 LiDIA::swap(TrAddr[j][i], T[j][i]);
610 //
611 // Free allocated storage
612 // and restore old precision
613 //
614 bigfloat::set_precision(old_prec);
615 delete[] tempvect1;
616 delete[] Tdel;
617 delete[] T;
618 delete[] mydel;
619 delete[] my;
620 }
621
622
623
624 //
625 // Schnorr - Euchner version of lll optimized for bigfloats using small bigfloats
626 // result : lll - reduced lattice for parameter y
627 //
628
629 //
630 // only high precision if needed
631 //
Tr_lll_bfl_gensys(lidia_size_t & rank,sdigit bit_prec)632 void bigfloat_lattice_gensys::Tr_lll_bfl_gensys(lidia_size_t& rank, sdigit bit_prec)
633 {
634 debug_handler("bigfloat_lattice_gensys", "Tr_lll_bfl_gensys(rank, bit_prec) [1]");
635 bool Fc;
636 bool korr;
637 bool break_flag;
638 bool repeat_loop;
639 lidia_size_t k;
640 lidia_size_t mom_rows = rows;
641 sdigit old_prec;
642 sdigit short_prec = static_cast<sdigit>(static_cast<double>(bit_prec)*std::log(2.0)/std::log(10.0));
643 sdigit r_prec;
644 sdigit n_prec;
645 sdigit bi_bit_len;
646 sdigit bit_diff;
647 sdigit x_exponent;
648 bigint x_factor;
649 bigfloat *tempvect1;
650 bigfloat bound = std::exp(std::log(2.0)*bit_prec/2);
651 bigfloat invbound;
652 bigfloat *mydel;
653 bigfloat **my;
654 bigfloat **lattice;
655 bigfloat *vect;
656 bigfloat *vect1;
657 bigfloat *vect_sq;
658 bigfloat conv;
659 bigfloat u_zero;
660 bigfloat halb(0.5);
661
662 LiDIA::divide(invbound, 1, bound);
663 //
664 // Compute needed precision for exact arithmetik
665 //
666 old_prec = bigfloat::get_precision();
667 r_prec = compute_read_precision();
668 log(u_zero, 10);
669 LiDIA::multiply(u_zero, u_zero, r_prec);
670 exp(u_zero, u_zero);
671 LiDIA::divide(u_zero, 1, u_zero);
672 n_prec = static_cast<sdigit>(r_prec/std::log(10.0))+1;
673 n_prec *= columns; // columns >= rows, basis !!!
674 n_prec = static_cast<sdigit>(n_prec*std::log(10.0))+1;
675
676 vectsize = columns;
677
678 //
679 // Allocating Memory for matrix
680 //
681 my = new bigfloat*[rows*2];
682 memory_handler(my, "bigfloat_lattice_gensys", "Tr_lll_bfl_gensys(rank, bit_prec) :: "
683 "not enough memory !");
684 lattice = &my[rows];
685 //
686 // think again
687 //
688 mydel = new bigfloat[(rows+2)*(rows+columns)-columns];
689 memory_handler(mydel, "bigfloat_lattice_gensys", "Tr_lll_bfl_gensys(rank, bit_prec) :: "
690 "not enough memory !");
691 for (lidia_size_t i = 0; i < rows; i++) {
692 my[i] = &mydel[i*rows]; // rows x rows
693 lattice[i] = &mydel[(rows*rows)+(i*columns)]; // rows x columns
694 }
695
696 //
697 // Allocating Memory for vector
698 //
699
700 //
701 // think again Part II
702 //
703 tempvect1 = new bigfloat[columns];
704 memory_handler(tempvect1, "bigfloat_lattice_gensys", "Tr_lll_bfl_gensys(rank, bit_prec) :: "
705 "not enough memory !");
706
707 vect = &mydel[(rows+columns)*rows];
708 vect_sq = &vect[rows];
709 vect1 = &vect[2*rows];
710
711
712 //
713 // Start
714 //
715 // Copy into floating point
716 //
717 bigfloat::set_precision(short_prec);
718 for (lidia_size_t i = 0; i < rows; i++)
719 for (lidia_size_t j = 0; j < columns; j++) {
720 x_factor = value[i][j].mantissa();
721 x_exponent = value[i][j].exponent();
722 bi_bit_len = x_factor.bit_length();
723 if (bi_bit_len > bit_prec) {
724 bit_diff = bi_bit_len-bit_prec;
725 shift_right(x_factor, x_factor, bit_diff);
726 lattice[i][j].assign(x_factor);
727 shift_right(lattice[i][j], lattice[i][j], bit_diff);
728 }
729 else
730 lattice[i][j].assign(x_factor);
731 shift_left(lattice[i][j], lattice[i][j], x_exponent);
732 }
733
734 k = 1;
735 reduction_steps = 0;
736 swapping_steps = 0;
737 correction_steps = 0;
738 vectsize = columns;
739 break_flag = false;
740 bfl_assign_zero_bfl(vect);
741 bfl_scalprod_bfl(vect[0], lattice[0], lattice[0]);
742
743 while (k < mom_rows) {
744 Fc = true;
745 repeat_loop = false;
746 while (Fc) {
747 //
748 // begin of orthogonalization
749 //
750 debug_handler("bigfloat_lattice_gensys", "Tr_lll_bfl_gensys(rank, bit_prec) [2]");
751 bfl_scalprod_bfl(vect[k], lattice[k], lattice[k]);
752 sqrt(vect_sq[k], vect[k]);
753 for (lidia_size_t j = 0; j < k; ++j) {
754 bfl_scalprod_bfl(tempmz0, lattice[k], lattice[j]);
755 LiDIA::multiply(tempmz1, invbound, vect_sq[k]);
756 LiDIA::multiply(tempmz1, tempmz1, vect_sq[j]);
757 //
758 // Correction Step I, only if needed, |lattice[k]*lattice[j]| to big
759 //
760 if (abs(tempmz0).compare(tempmz1) < 0) {
761 //
762 // Compute the correct scalar product
763 //
764 bigfloat::set_precision(n_prec);
765 bfl_scalprod_bfl(tempmz0, value[k], value[j]);
766 x_factor = tempmz0.mantissa();
767 x_exponent = tempmz0.exponent();
768 bi_bit_len = x_factor.bit_length();
769 if (bi_bit_len > bit_prec) {
770 bit_diff = bi_bit_len-bit_prec;
771 shift_right(x_factor, x_factor, bit_diff);
772 my[j][k].assign(x_factor);
773 shift_left(my[j][k], my[j][k], bit_diff);
774 }
775 else
776 my[j][k].assign(x_factor);
777 shift_left(my[j][k], my[j][k], x_exponent);
778 bigfloat::set_precision(short_prec);
779 }
780 else
781 my[j][k].assign(tempmz0);
782
783 for (lidia_size_t i = 0; i < j; ++i) {
784 LiDIA::multiply(tempmz1, my[i][j], my[i][k]);
785 LiDIA::multiply(tempmz1, tempmz1, vect[i]);
786 LiDIA::subtract(my[j][k], my[j][k], tempmz1);
787 }
788
789 tempmz0.assign(my[j][k]);
790 LiDIA::divide(my[j][k], my[j][k], vect[j]);
791 LiDIA::multiply(tempmz1, my[j][k], tempmz0);
792 LiDIA::subtract(vect[k], vect[k], tempmz1);
793 }
794 //
795 // end of orthogonalization
796 //
797 Fc = false;
798 korr = false;
799 //
800 // begin of reduction
801 //
802 debug_handler("bigfloat_lattice_gensys", "Tr_lll_bfl_gensys(rank, bit_prec) [3]");
803 if (!repeat_loop) {
804 for (lidia_size_t j = k-1; j >= 0; j--)
805 if (abs(my[j][k]).compare(halb) > 0) {
806 ++reduction_steps;
807 korr = true;
808 round (tempmz0, my[j][k]);
809 if (abs(tempmz0).compare(bound) > 0)
810 Fc = true;
811
812 bigfloat::set_precision(n_prec);
813 bfl_scalmul_bfl(tempvect1, tempmz0, value[j]);
814 bfl_subtract_bfl(value[k], value[k], tempvect1);
815 bigfloat::set_precision(short_prec);
816
817 for (lidia_size_t i = 0; i < j; ++i) {
818 LiDIA::multiply(tempmz1, tempmz0, my[i][j]);
819 LiDIA::subtract(my[i][k], my[i][k], tempmz1);
820 }
821 LiDIA::subtract(my[j][k], my[j][k], tempmz0);
822
823 }
824 }
825
826 //
827 // Correction Step II
828 // lattice=(bigfloat )value
829 //
830 if (korr) {
831 ++correction_steps;
832 for (lidia_size_t j = 0; j < columns; j++) {
833 x_factor = value[k][j].mantissa();
834 x_exponent = value[k][j].exponent();
835 bi_bit_len = x_factor.bit_length();
836 if (bi_bit_len > bit_prec) {
837 bit_diff = bi_bit_len-bit_prec;
838 shift_right(x_factor, x_factor, bit_diff);
839 lattice[k][j].assign(x_factor);
840 shift_left(lattice[k][j], lattice[k][j], bit_diff);
841 }
842 else
843 lattice[k][j].assign(x_factor);
844 shift_left(lattice[k][j], lattice[k][j], x_exponent);
845 }
846 bfl_scalquad_bfl(vect[k], lattice[k]);
847 sqrt(vect_sq[k], vect_sq[k]);
848 }
849 //
850 // Check if it`s 0 - vector
851 // delete it
852 //
853 if (abs(vect_sq[k]) < u_zero) {
854 for (lidia_size_t j = 0; j < mom_rows-1; j++) {
855 bfl_swap_bfl(lattice[j], lattice[j+1]);
856 bfl_swap_bfl(value[j], value[j+1]);
857 }
858 bfl_assign_zero_bfl(value[--mom_rows]);
859 bfl_assign_zero_bfl(lattice[mom_rows]);
860 k = 1;
861 bfl_assign_zero_bfl(vect);
862 bfl_scalprod_bfl(vect[0], lattice[0], lattice[0]);
863 sqrt(vect_sq[0], vect[0]);
864 break_flag = true;
865 korr = false;
866 Fc = false;
867 break;
868 }
869
870 //
871 // end of reduction
872 //
873 if (Fc) {
874 if (k > 1)
875 --k;
876 }
877 else {
878 Fc = (!repeat_loop) && korr;
879 repeat_loop = !repeat_loop;
880 }
881 }
882
883 if (!break_flag) {
884 debug_handler("bigfloat_lattice_gensys", "Tr_lll_bfl_gensys(rank, bit_prec) [4]");
885 //
886 // Check lll - condition for parameter y_par
887 //
888 tempmz0.assign(y_par);
889 LiDIA::multiply(tempmz0, tempmz0, vect[k-1]);
890 LiDIA::square(tempmz1, my[k-1][k]);
891 LiDIA::multiply(tempmz1, tempmz1, vect[k-1]);
892 LiDIA::add(tempmz1, tempmz1, vect[k]);
893
894 if (tempmz0.compare(tempmz1) > 0) {
895 //
896 // ::swap vectors k and k-1
897 //
898 ++swapping_steps;
899 bfl_swap_bfl(value[k], value[k-1]);
900 bfl_swap_bfl(lattice[k], lattice[k-1]);
901
902 if (k > 1)
903 --k;
904 if (k == 1) {
905 bfl_scalprod_bfl(vect[0], lattice[0], lattice[0]);
906 sqrt(vect_sq[0], vect[0]);
907 }
908 }
909 else
910 ++k;
911 }
912 break_flag = false;
913 }
914 rank = mom_rows;
915 //
916 // Free allocated storage
917 // and restore precision
918 //
919 bigfloat::set_precision(old_prec);
920 delete[] tempvect1;
921 delete[] mydel;
922 delete[] my;
923 }
924
925
926
927 //
928 // Schnorr - Euchner version of lll optimized for bigints usign bigfloats
929 // result : lll - reduced lattice for parameter y
930 //
931
932 //
933 // only high precision if needed
934 //
Tr_lll_trans_bfl_gensys(math_matrix<bigint> & Tr,lidia_size_t & rank,sdigit bit_prec)935 void bigfloat_lattice_gensys::Tr_lll_trans_bfl_gensys(math_matrix< bigint > & Tr, lidia_size_t& rank,
936 sdigit bit_prec)
937 {
938 debug_handler("bigfloat_lattice_gensys", "Tr_lll_trans_bfl_gensys(Tr, rank, bit_prec) [1]");
939 bool Fc;
940 bool korr;
941 bool break_flag;
942 bool repeat_loop;
943 lidia_size_t k;
944 lidia_size_t mom_rows = rows;
945 sdigit old_prec;
946 sdigit short_prec = static_cast<sdigit>(static_cast<double>(bit_prec)*std::log(2.0)/std::log(10.0));
947 sdigit r_prec;
948 sdigit n_prec;
949 sdigit bi_bit_len;
950 sdigit bit_diff;
951 sdigit x_exponent;
952 bigint x_factor;
953 bigfloat *tempvect1;
954 bigint *tempvect2;
955 bigfloat bound = std::exp(std::log(2.0)*bit_prec/2);
956 bigfloat invbound;
957 bigfloat *mydel;
958 bigfloat **my;
959 bigfloat **lattice;
960 bigfloat u_zero;
961 bigint *Tdel;
962 bigint **T;
963 bigint **TrAddr;
964 bigfloat *vect;
965 bigfloat *vect_sq;
966 bigint conv;
967 bigfloat halb(0.5);
968
969 LiDIA::divide(invbound, 1, bound);
970 //
971 // Compute needed precision for exact arithmetik
972 //
973 old_prec = bigfloat::get_precision();
974 r_prec = compute_read_precision();
975 log(u_zero, 10);
976 LiDIA::multiply(u_zero, u_zero, r_prec);
977 exp(u_zero, u_zero);
978 LiDIA::divide(u_zero, 1, u_zero);
979 n_prec = static_cast<sdigit>(r_prec/std::log(10.0))+1;
980 n_prec *= columns; // columns >= rows, basis !!!
981 n_prec = static_cast<sdigit>(n_prec*std::log(10.0))+1;
982
983 vectsize = columns;
984
985 //
986 // Allocating Memory for matrix
987 //
988 my = new bigfloat*[rows*2];
989 memory_handler(my, "bigfloat_lattice_gensys", "Tr_lll_trans_bfl_gensys(Tr, rank, bit_prec) :: "
990 "not enough memory !");
991 lattice = &my[rows];
992 //
993 // think again
994 //
995 mydel = new bigfloat[(rows+2)*(rows+columns)-columns];
996 memory_handler(mydel, "bigfloat_lattice_gensys", "Tr_lll_trans_bfl_gensys(Tr, rank, bit_prec) :: "
997 "not enough memory !");
998 T = new bigint*[rows];
999 memory_handler(T, "bigfloat_lattice_gensys", "Tr_lll_trans_bfl_gensys(Tr, rank, bit_prec) :: "
1000 "Not enough memory !");
1001
1002 Tdel = new bigint[rows*rows+rows];
1003 memory_handler(Tdel, "bigfloat_lattice_gensys", "Tr_lll_trans_bfl_gensys(Tr, rank, bit_prec) :: "
1004 "Not enough memory !");
1005
1006 for (lidia_size_t i = 0; i < rows; i++) {
1007 my[i] = &mydel[i*rows]; // rows x rows
1008 lattice[i] = &mydel[(rows*rows)+(i*columns)]; // rows x columns
1009 T[i] = &Tdel[i*rows];
1010 T[i][i].assign_one();
1011 }
1012
1013 //
1014 // Allocating Memory for vector
1015 //
1016
1017 //
1018 // think again Part II
1019 //
1020 tempvect1 = new bigfloat[columns];
1021 memory_handler(tempvect1, "bigfloat_lattice_gensys", "Tr_lll_trans_bfl_gensys(Tr, rank, bit_prec) :: "
1022 "not enough memory !");
1023
1024 tempvect2 = &T[0][rows*rows];
1025 vect = &my[0][(rows+columns)*rows];
1026 vect_sq = &vect[rows];
1027
1028
1029 //
1030 // Start
1031 //
1032 // Copy into floating point
1033 //
1034 bigfloat::set_precision(short_prec);
1035 for (lidia_size_t i = 0; i < rows; i++)
1036 for (lidia_size_t j = 0; j < columns; j++) {
1037 x_factor = value[i][j].mantissa();
1038 x_exponent = value[i][j].exponent();
1039 bi_bit_len = x_factor.bit_length();
1040 if (bi_bit_len > bit_prec) {
1041 bit_diff = bi_bit_len-bit_prec;
1042 shift_right(x_factor, x_factor, bit_diff);
1043 lattice[i][j].assign(x_factor);
1044 shift_right(lattice[i][j], lattice[i][j], bit_diff);
1045 }
1046 else
1047 lattice[i][j].assign(x_factor);
1048 shift_left(lattice[i][j], lattice[i][j], x_exponent);
1049 }
1050
1051 k = 1;
1052 reduction_steps = 0;
1053 swapping_steps = 0;
1054 correction_steps = 0;
1055 vectsize = columns;
1056 break_flag = false;
1057 bfl_assign_zero_bfl(vect);
1058 bfl_scalprod_bfl(vect[0], lattice[0], lattice[0]);
1059
1060 while (k < mom_rows) {
1061 Fc = true;
1062 repeat_loop = false;
1063 while (Fc) {
1064 //
1065 // begin of orthogonalization
1066 //
1067 debug_handler("bigfloat_lattice_gensys", "Tr_lll_trans_bfl_gensys(Tr, rank, bit_prec) [2]");
1068 bfl_scalprod_bfl(vect[k], lattice[k], lattice[k]);
1069 sqrt(vect_sq[k], vect[k]);
1070 for (lidia_size_t j = 0; j < k; ++j) {
1071 bfl_scalprod_bfl(tempmz0, lattice[k], lattice[j]);
1072 LiDIA::multiply(tempmz1, invbound, vect_sq[k]);
1073 LiDIA::multiply(tempmz1, tempmz1, vect_sq[j]);
1074 //
1075 // Correction Step I, only if needed, |lattice[k]*lattice[j]| to big
1076 //
1077 if (abs(tempmz0).compare(tempmz1) < 0) {
1078 //
1079 // Compute the correct scalar product
1080 //
1081 bigfloat::set_precision(n_prec);
1082 bfl_scalprod_bfl(tempmz0, value[k], value[j]);
1083 x_factor = tempmz0.mantissa();
1084 x_exponent = tempmz0.exponent();
1085 bi_bit_len = x_factor.bit_length();
1086 if (bi_bit_len > bit_prec) {
1087 bit_diff = bi_bit_len-bit_prec;
1088 shift_right(x_factor, x_factor, bit_diff);
1089 my[j][k].assign(x_factor);
1090 shift_left(my[j][k], my[j][k], bit_diff);
1091 }
1092 else
1093 my[j][k].assign(x_factor);
1094 shift_left(my[j][k], my[j][k], x_exponent);
1095 bigfloat::set_precision(short_prec);
1096
1097 }
1098 else
1099 my[j][k].assign(tempmz0);
1100
1101 for (lidia_size_t i = 0; i < j; ++i) {
1102 LiDIA::multiply(tempmz1, my[i][j], my[i][k]);
1103 LiDIA::multiply(tempmz1, tempmz1, vect[i]);
1104 LiDIA::subtract(my[j][k], my[j][k], tempmz1);
1105 }
1106
1107 tempmz0.assign(my[j][k]);
1108 LiDIA::divide(my[j][k], my[j][k], vect[j]);
1109 LiDIA::multiply(tempmz1, my[j][k], tempmz0);
1110 LiDIA::subtract(vect[k], vect[k], tempmz1);
1111 }
1112 //
1113 // end of orthogonalization
1114 //
1115 Fc = false;
1116 korr = false;
1117 //
1118 // begin of reduction
1119 //
1120 debug_handler("bigfloat_lattice_gensys", "Tr_lll_trans_bfl_gensys(Tr, rank, bit_prec) [3]");
1121 if (!repeat_loop) {
1122 for (lidia_size_t j = k-1; j >= 0; j--)
1123 if (abs(my[j][k]).compare(halb) > 0) {
1124 ++reduction_steps;
1125 korr = true;
1126 round (tempmz0, my[j][k]);
1127 if (abs(tempmz0).compare(bound) > 0)
1128 Fc = true;
1129
1130 bigfloat::set_precision(n_prec);
1131 bfl_scalmul_bfl(tempvect1, tempmz0, value[j]);
1132 bfl_subtract_bfl(value[k], value[k], tempvect1);
1133 vectsize = rows;
1134 bfl_scalmul_bin(tempvect2, conv, T[j]);
1135 bfl_subtract_bin(T[k], T[k], tempvect2);
1136 vectsize = columns;
1137 bigfloat::set_precision(short_prec);
1138
1139 for (lidia_size_t i = 0; i < j; ++i) {
1140 LiDIA::multiply(tempmz1, tempmz0, my[i][j]);
1141 LiDIA::subtract(my[i][k], my[i][k], tempmz1);
1142 }
1143 LiDIA::subtract(my[j][k], my[j][k], tempmz0);
1144 }
1145 }
1146
1147 //
1148 // Correction Step II
1149 // lattice[k]=(bigfloat )value[k]
1150 //
1151 if (korr) {
1152 ++correction_steps;
1153 for (lidia_size_t j = 0; j < columns; j++) {
1154 x_factor = value[k][j].mantissa();
1155 x_exponent = value[k][j].exponent();
1156 bi_bit_len = x_factor.bit_length();
1157 if (bi_bit_len > bit_prec) {
1158 bit_diff = bi_bit_len-bit_prec;
1159 shift_right(x_factor, x_factor, bit_diff);
1160 lattice[k][j].assign(x_factor);
1161 shift_left(lattice[k][j], lattice[k][j], bit_diff);
1162 }
1163 else
1164 lattice[k][j].assign(x_factor);
1165 shift_left(lattice[k][j], lattice[k][j], x_exponent);
1166 }
1167
1168 bfl_scalquad_bfl(vect[k], lattice[k]);
1169 sqrt(vect_sq[k], vect_sq[k]);
1170 }
1171 //
1172 // Check if it`s 0 - vector
1173 // delete it
1174 //
1175 if (abs(vect_sq[k]) < u_zero) {
1176 for (lidia_size_t j = 0; j < mom_rows-1; j++) {
1177 bfl_swap_bfl(lattice[j], lattice[j+1]);
1178 bfl_swap_bfl(value[j], value[j+1]);
1179 bfl_swap_bin(T[j], T[j+1]);
1180 }
1181 bfl_assign_zero_bfl(value[--mom_rows]);
1182 bfl_assign_zero_bin(T[mom_rows]);
1183 k = 1;
1184 bfl_assign_zero_bfl(vect);
1185 bfl_scalprod_bfl(vect[0], lattice[0], lattice[0]);
1186 sqrt(vect_sq[0], vect[0]);
1187 break_flag = true;
1188 korr = false;
1189 Fc = false;
1190 break;
1191 }
1192 //
1193 // end of reduction
1194 //
1195 if (Fc) {
1196 if (k > 1)
1197 --k;
1198 }
1199 else {
1200 Fc = (!repeat_loop) && korr;
1201 repeat_loop = !repeat_loop;
1202 }
1203 }
1204
1205 if (!break_flag) {
1206 debug_handler("bigfloat_lattice_gensys", "Tr_lll_trans_bfl_gensys(Tr, rank, bit_prec) [4]");
1207 //
1208 // Check lll - condition for parameter y_par
1209 //
1210 tempmz0.assign(y_par);
1211 LiDIA::multiply(tempmz0, tempmz0, vect[k-1]);
1212 LiDIA::square(tempmz1, my[k-1][k]);
1213 LiDIA::multiply(tempmz1, tempmz1, vect[k-1]);
1214 LiDIA::add(tempmz1, tempmz1, vect[k]);
1215
1216 if (tempmz0.compare(tempmz1) > 0) {
1217 //
1218 // ::swap vectors k and k-1
1219 //
1220 ++swapping_steps;
1221 bfl_swap_bfl(value[k], value[k-1]);
1222 bfl_swap_bfl(lattice[k], lattice[k-1]);
1223 bfl_swap_bin(T[k], T[k-1]);
1224
1225 if (k > 1)
1226 --k;
1227 if (k == 1) {
1228 bfl_scalprod_bfl(vect[0], lattice[0], lattice[0]);
1229 sqrt(vect_sq[0], vect[0]);
1230 }
1231 }
1232 else
1233 ++k;
1234 }
1235 break_flag = false;
1236 }
1237 //
1238 // Copy into math_matrix< bigint > Tr
1239 //
1240 Tr.set_no_of_rows(rows);
1241 Tr.set_no_of_columns(rows);
1242 TrAddr = Tr.get_data_address();
1243 for (lidia_size_t i = 0; i < rows; i++)
1244 for (lidia_size_t j = 0; j < rows; j++)
1245 if (trans_flag)
1246 LiDIA::swap(TrAddr[j][i], T[i][j]);
1247 else
1248 LiDIA::swap(TrAddr[j][i], T[j][i]);
1249 rank = mom_rows;
1250 //
1251 // Free allocated storage
1252 // and restore old precsion
1253 //
1254 bigfloat::set_precision(old_prec);
1255 delete[] tempvect1;
1256 delete[] Tdel;
1257 delete[] T;
1258 delete[] mydel;
1259 delete[] my;
1260 }
1261
1262
1263
1264 //
1265 // Buchmann Kessler version of lll for generating systems
1266 // result : transformation lattice
1267 //
Tr_lin_gen_system(math_matrix<bigint> & T,lidia_size_t & rk)1268 void bigfloat_lattice_gensys::Tr_lin_gen_system(math_matrix< bigint > & T, lidia_size_t& rk)
1269 {
1270 debug_handler("bigfloat_lattice_gensys", "Tr_lin_gen_system(T, rk) [1]");
1271
1272 bigint_lattice_gensys Atilde(rows, rows+columns), // The approximate lattice
1273 Ahead(rows, columns);
1274 bigint *help = new bigint[columns];
1275 bigint *rel = new bigint[rows];
1276 bigint *temp = new bigint[columns];
1277
1278 bigfloat vor, rechts;
1279 bigfloat zweipotq, alpha;
1280 bigfloat norm1, norm2;
1281
1282 bigint bi_norm1, bi_norm2;
1283 bigint zero;
1284 bigint **TAddr;
1285
1286 sdigit n2 = rows;
1287 sdigit old_prec;
1288 sdigit prec;
1289 sdigit bit_prec;
1290
1291
1292 zero.assign_zero();
1293 old_prec = bigfloat::get_precision();
1294 prec = compute_precision();
1295 bigfloat::set_precision(prec);
1296 alpha_compute(alpha);
1297 zwei_pot_q_compute(zweipotq, n2, alpha);
1298
1299
1300 for (lidia_size_t i = 0; i < rows; ++i)
1301 for (lidia_size_t j = 0; j < columns; ++j) {
1302 LiDIA::multiply(tempmz0, zweipotq, value[i][j]);
1303 tempmz0.bigintify(Ahead.value[i][j]);
1304 }
1305 for (lidia_size_t i = 0; i < rows; ++i) {
1306 Atilde.value[i][i].assign_one(); // 1, wenn i = j, 0 sonst
1307 for (lidia_size_t j = rows; j < (rows+columns); ++j)
1308 Atilde.value[i][j].assign(Ahead.value[i][j-rows]);
1309 }
1310
1311 prec = Atilde.compute_read_precision();
1312 bit_prec = static_cast<sdigit>(static_cast<double>(prec)*(std::log(10.0)/std::log(2.0)));
1313 bit_prec = ((bit_prec/300)+1)*52;
1314
1315 debug_handler("bigfloat_lattice_gensys", "Tr_lin_gen_system(T, rk) [2]");
1316 Atilde.y_par = y_par;
1317 Atilde.Tr_lll_bfl(bit_prec);
1318 debug_handler("bigfloat_lattice_gensys", "Tr_lin_gen_system(T, rk) [3]");
1319 reduction_steps = Atilde.reduction_steps;
1320 swapping_steps = Atilde.swapping_steps;
1321 correction_steps = Atilde.correction_steps;
1322
1323
1324 lidia_size_t l = 0;
1325 do {
1326 vectsize = columns;
1327 bfl_assign_zero_bin(help); // Initializes help with the zero - vector
1328 for (lidia_size_t j = 0; j < rows; ++j) {
1329 rel[j].assign(Atilde.value[l][j]);
1330 bfl_scalmul_bin(temp, rel[j], Ahead.value[j]);
1331 bfl_add_bin(help, help, temp);
1332 }
1333 bfl_l2_norm_bin(bi_norm2, help);
1334 norm2.assign(bi_norm2);
1335 vectsize = rows;
1336 bfl_l1_norm_bin(bi_norm1, rel);
1337 norm1.assign(bi_norm1);
1338 sqrt(norm1, norm1);
1339 ++l;
1340 LiDIA::divide(vor, n2, rows);
1341 vor.multiply_by_2();
1342 sqrt(vor, vor);
1343 LiDIA::multiply(vor, vor, norm1);
1344 LiDIA::divide(rechts, bigfloat(std::sqrt(static_cast<double>(n2))), bigfloat(2.0));
1345 LiDIA::multiply(rechts, rechts, norm1);
1346 }
1347 while ((zweipotq.compare(vor) > 0) && (norm2.compare(rechts) <= 0) && (l < rows));
1348 debug_handler("bigfloat_lattice_gensys", "Tr_lin_gen_system(T, rk) [4]");
1349 if (l >= rows)
1350 warning_handler("bigfloat_lattice_gensys", "Tr_lin_gen_system(T, rk) :: "
1351 "lattice of dimension 1");
1352
1353 rk = rows-l+1; // rows is the dimension of the lattice
1354 T.set_no_of_rows(rows);
1355 T.set_no_of_columns(rows);
1356 TAddr = T.get_data_address();
1357 for (lidia_size_t i = 0; i < rows; i++)
1358 for (lidia_size_t j = 0; j < rows-rk; j++)
1359 if (trans_flag)
1360 TAddr[i][j].assign_zero();
1361 else
1362 TAddr[j][i].assign_zero();
1363 for (lidia_size_t i = 0; i < rows; ++i)
1364 for (lidia_size_t j = rows-rk; j < rows; ++j)
1365 if (trans_flag)
1366 LiDIA::swap(TAddr[i][j], Atilde.value[j][i]);
1367 else
1368 LiDIA::swap(TAddr[j][i], Atilde.value[j][i]);
1369
1370 //
1371 // Free allocated storage
1372 // and restore old precsion
1373 //
1374 bigfloat::set_precision(old_prec);
1375 delete[] help;
1376 delete[] rel;
1377 delete[] temp;
1378 }
1379
1380
1381
1382 //
1383 // the modified lll-algorithm
1384 //
Tr_mlll_bfl(bigint * v)1385 void bigfloat_lattice_gensys::Tr_mlll_bfl(bigint* v)
1386 {
1387 debug_handler("bigfloat_lattice_gensys", "Tr_mlll_bfl(v)");
1388 lidia_size_t j, k = columns, l, m;
1389 sdigit old_prec;
1390 sdigit prec;
1391 bigint *tempvect0;
1392 bigfloat *tempvect1;
1393 bigfloat *tempvect2;
1394 bigfloat *tempvect3;
1395 bigfloat *B;
1396 bigfloat halb(0.5);
1397 bigfloat dreiviertel(0.75);
1398 bigfloat Mu, Bz;
1399 bigint tempbin0;
1400 bigfloat *mudel;
1401 bigfloat **mu, **bs;
1402 bigint *Trdel;
1403 bigint **Tr;
1404
1405 //
1406 // Allocating needed storage
1407 //
1408 // Lattices
1409 //
1410
1411 Tr = new bigint*[rows];
1412 memory_handler(T, "bigfloat_lattice_gensys", "Tr_mll_bfl(v) :: "
1413 "not enough memory !");
1414 mu = new bigfloat*[2*rows];
1415 memory_handler(mu, "bigfloat_lattice_gensys", "Tr_mlll_bfl(v) :: "
1416 "not enough memory !");
1417
1418 bs = &mu[rows];
1419 mudel = new bigfloat[rows*(rows+columns)+3*rows+columns];
1420 memory_handler(mudel, "bigfloat_lattice_gensys", "Tr_mlll_bfl(v) :: "
1421 "not enough memory !");
1422 Trdel = new bigint[rows*rows+rows];
1423 memory_handler(Trdel, "bigfloat_lattice_gensys", "Tr_mlll_bfl(v) :: "
1424 "not enough memory !");
1425
1426
1427 old_prec = bigfloat::get_precision();
1428 prec = compute_read_precision();
1429 bigfloat::set_precision(rows*prec);
1430 reduction_steps = 0;
1431 swapping_steps = 0;
1432 correction_steps = 0;
1433 for (lidia_size_t i = 0; i < rows; i++) {
1434 Tr[i] = &Trdel[i*rows];
1435 Tr[i][i].assign_one();
1436 mu[i] = &mudel[i*rows];
1437 bs[i] = &mudel[rows*rows+i*columns];
1438 }
1439
1440
1441 //
1442 // Vectors
1443 //
1444 tempvect2 = &mudel[rows*(rows+columns)];
1445 tempvect3 = &tempvect2[rows];
1446 B = &tempvect2[2*rows];
1447
1448 tempvect0 = &Trdel[rows*rows];
1449 tempvect1 = &tempvect2[3*rows];
1450
1451 //
1452 // Start Timer
1453 //
1454 vectsize = columns;
1455
1456
1457 for (lidia_size_t i = 0; i < k+1; i++) {
1458 for (lidia_size_t h = 0; h < i; h++) {
1459 mu[h][i].assign_zero();
1460 bfl_scalprod_bfl(mu[h][i], bs[h], value[i]);
1461
1462 LiDIA::divide(mu[h][i], mu[h][i], B[h]);
1463 }
1464 bfl_assign_zero_bfl(tempvect3);
1465 for (lidia_size_t h = 0; h < i; h++) {
1466 bfl_scalmul_bfl(tempvect2, mu[h][i], bs[h]);
1467 bfl_add_bfl(tempvect3, tempvect3, tempvect2);
1468 }
1469 bfl_subtract_bfl(bs[i], value[i], tempvect3);
1470 bfl_scalprod_bfl(B[i], bs[i], bs[i]);
1471 }
1472
1473 m = 1;
1474 l = 0;
1475
1476 while (true) {
1477 //
1478 // Reduction Step
1479 //
1480 if (abs(mu[l][m]).compare(halb) > 0) {
1481 reduction_steps++;
1482 round(tempmz0, mu[l][m]);
1483 tempmz0.bigintify(tempbin0);
1484 bfl_scalmul_bfl(tempvect1, tempmz0, value[l]);
1485 bfl_subtract_bfl(value[m], value[m], tempvect1);
1486 vectsize = rows;
1487 bfl_scalmul_bin(tempvect0, tempbin0, Tr[l]);
1488 bfl_subtract_bin(Tr[m], Tr[m], tempvect0);
1489 vectsize = columns;
1490
1491 LiDIA::subtract(mu[l][m], mu[l][m], tempmz0);
1492 for (lidia_size_t h = 0; h < l; h++) {
1493 LiDIA::multiply(tempmz1, tempmz0, mu[h][l]);
1494 LiDIA::subtract(mu[h][m], mu[h][m], tempmz1);
1495 }
1496 }
1497
1498
1499 j = 0;
1500 while ((j < k) && (value[m][j].is_approx_zero())) {
1501 j++;
1502 }
1503
1504 //
1505 // Exiting - condition
1506 //
1507 if (j == k) {
1508 for (lidia_size_t i = m; i < k; i++)
1509 bfl_assign_bfl(value[i], value[i+1]);
1510 bfl_assign_zero_bfl(value[rows-1]);
1511 vectsize = rows;
1512 bfl_assign_bin(v, Tr[m]);
1513 vectsize = columns;
1514 bigfloat::set_precision(old_prec);
1515 delete[] mudel;
1516 delete[] mu;
1517 delete[] Trdel;
1518 delete[] Tr;
1519 return;
1520 }
1521
1522
1523 if (l >= m-1) {
1524 LiDIA::square(tempmz1, mu[m-1][m]);
1525 LiDIA::subtract(tempmz0, dreiviertel, tempmz1);
1526 LiDIA::multiply(tempmz0, tempmz0, B[m-1]);
1527
1528 if (B[m].compare(tempmz0) < 0) {
1529 Mu.assign(mu[m-1][m]);
1530 LiDIA::square(tempmz0, Mu);
1531 LiDIA::multiply(tempmz0, tempmz0, B[m-1]);
1532 LiDIA::add(Bz, B[m], tempmz0);
1533 if (!Bz.is_approx_zero()) {
1534 LiDIA::divide(tempmz0, B[m-1], Bz);
1535 LiDIA::multiply(mu[m-1][m], Mu, tempmz0);
1536 LiDIA::multiply(B[m], B[m], tempmz0);
1537 for (lidia_size_t i = m+1; i < k+1; i++) {
1538 LiDIA::multiply(tempmz0, Mu, mu[m][i]);
1539 LiDIA::subtract(tempmz1, mu[m-1][i], tempmz0);
1540 LiDIA::multiply(tempmz0, tempmz1, mu[m-1][m]);
1541 LiDIA::add(mu[m-1][i], mu[m][i], tempmz0);
1542 mu[m][i].assign(tempmz1);
1543 }
1544 }
1545 B[m-1].assign(Bz);
1546 swapping_steps++;
1547 bfl_swap_bfl(value[m-1], value[m]);
1548 bfl_swap_bin(Tr[m-1], Tr[m]);
1549 for (j = 0; j <= m-2; j++) {
1550 tempmz1.assign(mu[j][m-1]);
1551 mu[j][m-1].assign(mu[j][m]);
1552 mu[j][m].assign(tempmz1);
1553 }
1554 if (m > 1)
1555 m--;
1556 l = m-1;
1557 continue;
1558 }
1559 }
1560 l--;
1561 if (l < 0) {
1562 m++;
1563 l = m-1;
1564 }
1565 }
1566 }
1567
1568
1569
Tr_lll_var_bfl(math_matrix<bigint> & Bi,sdigit bit_len,sdigit bit_prec)1570 void bigfloat_lattice_gensys::Tr_lll_var_bfl(math_matrix< bigint > & Bi,
1571 sdigit bit_len, sdigit bit_prec)
1572 {
1573 debug_handler("bigfloat_lattice_gensys", "Tr_lll_var_bfl(Bi, bit_len, bit_prec) [1]");
1574 bool Fc;
1575 bool korr;
1576 bool repeat_loop;
1577 sdigit old_prec;
1578 sdigit old_bit_prec;
1579 sdigit bit_diff;
1580 sdigit bi_bit_len;
1581 sdigit short_prec;
1582 lidia_size_t k;
1583 bigint tempbiz;
1584 bigint *tempvect1;
1585 bigint **BiAddr;
1586 bigfloat bound = std::exp(std::log(2.0)*bit_prec/2);
1587 bigfloat invbound;
1588 bigfloat *mydel;
1589 bigfloat **my;
1590 bigfloat **lattice;
1591 bigfloat *vect;
1592 bigfloat *vect_sq;
1593 bigfloat halb(0.5);
1594
1595 LiDIA::divide(invbound, 1, bound);
1596
1597 //
1598 // Compute needed precision for exact arithmetik
1599 //
1600 old_prec = bigfloat::get_precision();
1601 old_bit_prec = static_cast<sdigit>(static_cast<double>(old_prec)/std::log(2.0)*std::log(10.0)+1);
1602 short_prec = static_cast<sdigit>(static_cast<double>(bit_prec)*std::log(2.0)/std::log(10.0));
1603 BiAddr = Bi.get_data_address();
1604
1605 vectsize = columns;
1606
1607 //
1608 // Allocating Memory for matrix
1609 //
1610 my = new bigfloat*[rows*2];
1611 memory_handler(my, "bigfloat_lattice_gensys", "Tr_lll_var_bfl(Bi, bit_prec) :: "
1612 "not enough memory !");
1613 lattice = &my[rows];
1614 //
1615 // think again
1616 //
1617 mydel = new bigfloat[rows*(rows+columns+2)+columns];
1618 memory_handler(mydel, "bigfloat_lattice_gensys", "Tr_lll_var_bfl(Bi, Factor, bit_prec) :: "
1619 "not enough memory !");
1620 for (lidia_size_t i = 0; i < rows; i++) {
1621 my[i] = &mydel[i*rows]; // rows x rows
1622 lattice[i] = &mydel[(rows*rows)+(i*columns)]; // rows x columns
1623 }
1624
1625 //
1626 // Allocating Memory for vector
1627 //
1628
1629 //
1630 // think again Part II
1631 //
1632 tempvect1 = new bigint[columns];
1633 memory_handler(tempvect1, "bigfloat_lattice_gensys", "Tr_lll_var_bfl(Bi, Factor, bit_prec) :: "
1634 "not enough memory !");
1635
1636 vect = &mydel[(rows+columns)*rows];
1637 vect_sq = &vect[rows];
1638
1639
1640 //
1641 // Start
1642 //
1643 // Copy into floating point
1644 //
1645 old_prec = bigfloat::get_precision();
1646 bigfloat::set_precision(short_prec);
1647 for (lidia_size_t i = 0; i < rows; i++)
1648 for (lidia_size_t j = 0; j < columns; j++) {
1649 bi_bit_len = BiAddr[i][j].bit_length();
1650 if (bi_bit_len > bit_prec) {
1651 bit_diff = bi_bit_len-bit_prec;
1652 shift_right(tempbiz, BiAddr[i][j], bit_diff);
1653 tempmz0.assign(tempbiz);
1654 }
1655 else {
1656 bit_diff = 0;
1657 tempmz0.assign(BiAddr[i][j]);
1658 }
1659 shift_right(lattice[i][j], tempmz0, bit_len-bit_diff);
1660 }
1661
1662 k = 1;
1663 reduction_steps = 0;
1664 swapping_steps = 0;
1665 correction_steps = 0;
1666 vectsize = columns;
1667 bfl_assign_zero_bfl(vect);
1668 bfl_scalprod_bfl(vect[0], lattice[0], lattice[0]);
1669
1670 while (k < rows) {
1671 Fc = true;
1672 repeat_loop = false;
1673 while (Fc) {
1674 //
1675 // begin of orthogonalization
1676 //
1677 debug_handler("bigfloat_lattice_gensys", "Tr_lll_var_bfl(Bi, Factor, bit_prec) [2]");
1678 bfl_scalprod_bfl(vect[k], lattice[k], lattice[k]);
1679 sqrt(vect_sq[k], vect[k]);
1680 for (lidia_size_t j = 0; j < k; ++j) {
1681 bfl_scalprod_bfl(tempmz0, lattice[k], lattice[j]);
1682 LiDIA::multiply(tempmz1, invbound, vect_sq[k]);
1683 LiDIA::multiply(tempmz1, tempmz1, vect_sq[j]);
1684 //
1685 // Correction Step I, only if needed, |lattice[k]*lattice[j]| to big
1686 //
1687 if (abs(tempmz0).compare(tempmz1) < 0) {
1688 //
1689 // Compute the correct scalar product
1690 //
1691 bfl_scalprod_bin(tempbiz, BiAddr[k], BiAddr[j]);
1692 bi_bit_len = tempbiz.bit_length();
1693 if (bi_bit_len > bit_prec) {
1694 bit_diff = bi_bit_len-bit_prec;
1695 shift_right(tempbiz, tempbiz, bit_diff);
1696 tempmz0.assign(tempbiz);
1697 }
1698 else {
1699 bit_diff = 0;
1700 tempmz0.assign(tempbiz);
1701 }
1702 shift_right(my[j][k], tempmz0, 2*bit_len-bit_diff);
1703 // ::divide(my[j][k], tempbiz, SqrFactor);
1704 }
1705 else
1706 my[j][k].assign(tempmz0);
1707
1708 for (lidia_size_t i = 0; i < j; ++i) {
1709 LiDIA::multiply(tempmz1, my[i][j], my[i][k]);
1710 LiDIA::multiply(tempmz1, tempmz1, vect[i]);
1711 LiDIA::subtract(my[j][k], my[j][k], tempmz1);
1712 }
1713
1714 tempmz0.assign(my[j][k]);
1715 LiDIA::divide(my[j][k], my[j][k], vect[j]);
1716 LiDIA::multiply(tempmz1, my[j][k], tempmz0);
1717 LiDIA::subtract(vect[k], vect[k], tempmz1);
1718 }
1719 //
1720 // end of orthogonalization
1721 //
1722 Fc = false;
1723 korr = false;
1724
1725 //
1726 // begin of reduction
1727 //
1728 debug_handler("bigfloat_lattice_gensys", "Tr_lll_var_bfl(Bi, Factor, bit_prec) [3]");
1729 if (!repeat_loop) {
1730 for (lidia_size_t j = k-1; j >= 0; j--)
1731 if (abs(my[j][k]).compare(halb) > 0) {
1732 ++reduction_steps;
1733 korr = true;
1734 round (tempmz0, my[j][k]);
1735 if (abs(tempmz0).compare(bound) > 0)
1736 Fc = true;
1737
1738 tempmz0.bigintify(tempbiz);
1739 bfl_scalmul_bin(tempvect1, tempbiz, BiAddr[j]);
1740 bfl_subtract_bin(BiAddr[k], BiAddr[k], tempvect1);
1741
1742 for (lidia_size_t i = 0; i < j; ++i) {
1743 LiDIA::multiply(tempmz1, tempmz0, my[i][j]);
1744 LiDIA::subtract(my[i][k], my[i][k], tempmz1);
1745 }
1746 LiDIA::subtract(my[j][k], my[j][k], tempmz0);
1747
1748 }
1749 }
1750
1751 //
1752 // Correction Step II
1753 // lattice=(bigfloat )value
1754 //
1755 if (korr) {
1756 ++correction_steps;
1757 for (lidia_size_t j = 0; j < columns; j++) {
1758 bi_bit_len = BiAddr[k][j].bit_length();
1759 if (bi_bit_len > bit_prec) {
1760 bit_diff = bi_bit_len-bit_prec;
1761 shift_right(tempbiz, BiAddr[k][j], bit_diff);
1762 tempmz0.assign(tempbiz);
1763 }
1764 else {
1765 bit_diff = 0;
1766 tempmz0.assign(BiAddr[k][j]);
1767 }
1768 shift_right(lattice[k][j], tempmz0, bit_len-bit_diff);
1769 }
1770 }
1771
1772 //
1773 // end of reduction
1774 //
1775 if (Fc) {
1776 if (k > 1)
1777 --k;
1778 }
1779 else {
1780 Fc = (!repeat_loop) && korr;
1781 repeat_loop = !repeat_loop;
1782 }
1783 }
1784
1785 debug_handler("bigfloat_lattice_gensys", "Tr_lll_var_bfl(Bi, Factor, bit_prec) [4]");
1786 //
1787 // Check lll - condition for parameter y_par
1788 //
1789 tempmz0.assign(y_par);
1790 LiDIA::multiply(tempmz0, tempmz0, vect[k-1]);
1791 LiDIA::square(tempmz1, my[k-1][k]);
1792 LiDIA::multiply(tempmz1, tempmz1, vect[k-1]);
1793 LiDIA::add(tempmz1, tempmz1, vect[k]);
1794
1795 if (tempmz0.compare(tempmz1) > 0) {
1796 //
1797 // ::swap vectors k and k-1
1798 //
1799 ++swapping_steps;
1800 bfl_swap_bin(BiAddr[k], BiAddr[k-1]);
1801 bfl_swap_bfl(lattice[k], lattice[k-1]);
1802
1803 if (k > 1)
1804 --k;
1805 if (k == 1) {
1806 bfl_scalprod_bfl(vect[0], lattice[0], lattice[0]);
1807 sqrt(vect_sq[0], vect[0]);
1808 }
1809 }
1810 else
1811 ++k;
1812 }
1813 bigfloat::set_precision(old_prec);
1814 for (lidia_size_t i = 0; i < rows; i++)
1815 for (lidia_size_t j = 0; j < columns; j++) {
1816 bi_bit_len = BiAddr[i][j].bit_length();
1817 if (bi_bit_len > old_bit_prec) {
1818 bit_diff = bi_bit_len-old_bit_prec;
1819 shift_right(tempbiz, BiAddr[i][j], bit_diff);
1820 tempmz0.assign(tempbiz);
1821 }
1822 else {
1823 bit_diff = 0;
1824 tempmz0.assign(BiAddr[i][j]);
1825 }
1826 shift_right(value[i][j], tempmz0, bit_len-bit_diff);
1827 }
1828
1829 //
1830 // Free allocated storage
1831 // and restore old precision
1832 //
1833 bigfloat::set_precision(old_prec);
1834 delete[] tempvect1;
1835 delete[] mydel;
1836 delete[] my;
1837 }
1838
1839
1840
Tr_lll_trans_var_bfl(math_matrix<bigint> & Bi,math_matrix<bigint> & Tr,sdigit bit_len,sdigit bit_prec)1841 void bigfloat_lattice_gensys::Tr_lll_trans_var_bfl(math_matrix< bigint > & Bi, math_matrix< bigint > & Tr,
1842 sdigit bit_len, sdigit bit_prec)
1843 {
1844 debug_handler("bigfloat_lattice_gensys", "Tr_lll_trans_var_bfl(Bi, Tr, Factor, bit_len, bit_prec) [1]");
1845 bool Fc;
1846 bool korr;
1847 bool repeat_loop;
1848 sdigit old_prec;
1849 sdigit old_bit_prec;
1850 sdigit bit_diff;
1851 sdigit bi_bit_len;
1852 sdigit short_prec;
1853 lidia_size_t k;
1854 bigint tempbiz;
1855 bigint *tempvect1;
1856 bigint **BiAddr;
1857 bigint *Tdel;
1858 bigint **T;
1859 bigint **TrAddr;
1860 bigfloat bound = std::exp(std::log(2.0)*bit_prec/2);
1861 bigfloat invbound;
1862 bigfloat *mydel;
1863 bigfloat **my;
1864 bigfloat **lattice;
1865 bigfloat *vect;
1866 bigfloat *vect_sq;
1867 bigfloat halb(0.5);
1868
1869 LiDIA::divide(invbound, 1, bound);
1870
1871 //
1872 // Compute needed precision for exact arithmetik
1873 //
1874 old_prec = bigfloat::get_precision();
1875 old_bit_prec = static_cast<sdigit>(static_cast<double>(old_prec)/std::log(2.0)*std::log(10.0)+1);
1876 short_prec = static_cast<sdigit>(static_cast<double>(bit_prec)*std::log(2.0)/std::log(10.0));
1877 BiAddr = Bi.get_data_address();
1878
1879 vectsize = columns;
1880
1881 //
1882 // Allocating Memory for matrix
1883 //
1884 my = new bigfloat*[rows*2];
1885 memory_handler(my, "bigfloat_lattice_gensys", "Tr_lll_trans_var_bfl(Bi, Tr, Factor, bit_prec) :: "
1886 "not enough memory !");
1887 lattice = &my[rows];
1888 T = new bigint*[rows];
1889 memory_handler(T, "bigfloat_lattice_gensys", "Tr_lll_trans_var_bfl(Bi, Tr, Factor, bit_prec) :: "
1890 "Not enough memory !");
1891
1892
1893 //
1894 // think again
1895 //
1896 mydel = new bigfloat[rows*(rows+columns+2)+columns];
1897 memory_handler(mydel, "bigfloat_lattice_gensys", "Tr_lll_trans_var_bfl(Bi, Tr, Factor, bit_prec) :: "
1898 "not enough memory !");
1899 Tdel = new bigint[rows*rows+columns];
1900 memory_handler(Tdel, "bigfloat_lattice_gensys", "Tr_lll_trans_var_bfl(Bi, Tr, Factor, bit_prec) :: "
1901 "Not enough memory !");
1902
1903 for (lidia_size_t i = 0; i < rows; i++) {
1904 my[i] = &mydel[i*rows]; // rows x rows
1905 T[i] = &Tdel[i*rows];
1906 T[i][i].assign_one();
1907 lattice[i] = &mydel[(rows*rows)+(i*columns)]; // rows x columns
1908 }
1909
1910 //
1911 // Allocating Memory for vector
1912 //
1913
1914 //
1915 // think again Part II
1916 //
1917
1918 tempvect1 = &Tdel[rows*rows];
1919 vect = &mydel[(rows+columns)*rows];
1920 vect_sq = &vect[rows];
1921
1922 //
1923 // Start
1924 //
1925 // Copy into floating point
1926 //
1927 old_prec = bigfloat::get_precision();
1928 bigfloat::set_precision(short_prec);
1929 for (lidia_size_t i = 0; i < rows; i++)
1930 for (lidia_size_t j = 0; j < columns; j++) {
1931 bi_bit_len = BiAddr[i][j].bit_length();
1932 if (bi_bit_len > bit_prec) {
1933 bit_diff = bi_bit_len-bit_prec;
1934 shift_right(tempbiz, BiAddr[i][j], bit_diff);
1935 tempmz0.assign(tempbiz);
1936 }
1937 else {
1938 bit_diff = 0;
1939 tempmz0.assign(BiAddr[i][j]);
1940 }
1941 shift_right(lattice[i][j], tempmz0, bit_len-bit_diff);
1942 }
1943
1944 k = 1;
1945 reduction_steps = 0;
1946 swapping_steps = 0;
1947 correction_steps = 0;
1948 vectsize = columns;
1949 bfl_assign_zero_bfl(vect);
1950 bfl_scalprod_bfl(vect[0], lattice[0], lattice[0]);
1951
1952 while (k < rows) {
1953 Fc = true;
1954 repeat_loop = false;
1955 while (Fc) {
1956 //
1957 // begin of orthogonalization
1958 //
1959 debug_handler("bigfloat_lattice_gensys", "Tr_lll_trans_var_bfl(Bi, Factor, bit_prec) [2]");
1960 bfl_scalprod_bfl(vect[k], lattice[k], lattice[k]);
1961 sqrt(vect_sq[k], vect[k]);
1962 for (lidia_size_t j = 0; j < k; ++j) {
1963 bfl_scalprod_bfl(tempmz0, lattice[k], lattice[j]);
1964 LiDIA::multiply(tempmz1, invbound, vect_sq[k]);
1965 LiDIA::multiply(tempmz1, tempmz1, vect_sq[j]);
1966 //
1967 // Correction Step I, only if needed, |lattice[k]*lattice[j]| to big
1968 //
1969 if (abs(tempmz0).compare(tempmz1) < 0) {
1970 //
1971 // Compute the correct scalar product
1972 //
1973 bfl_scalprod_bin(tempbiz, BiAddr[k], BiAddr[j]);
1974 bi_bit_len = tempbiz.bit_length();
1975 if (bi_bit_len > bit_prec) {
1976 bit_diff = bi_bit_len-bit_prec;
1977 shift_right(tempbiz, tempbiz, bit_diff);
1978 tempmz0.assign(tempbiz);
1979 }
1980 else {
1981 bit_diff = 0;
1982 tempmz0.assign(tempbiz);
1983 }
1984 shift_right(my[j][k], tempmz0, 2*bit_len-bit_diff);
1985 }
1986 else
1987 my[j][k].assign(tempmz0);
1988
1989 for (lidia_size_t i = 0; i < j; ++i) {
1990 LiDIA::multiply(tempmz1, my[i][j], my[i][k]);
1991 LiDIA::multiply(tempmz1, tempmz1, vect[i]);
1992 LiDIA::subtract(my[j][k], my[j][k], tempmz1);
1993 }
1994
1995 tempmz0.assign(my[j][k]);
1996 LiDIA::divide(my[j][k], my[j][k], vect[j]);
1997 LiDIA::multiply(tempmz1, my[j][k], tempmz0);
1998 LiDIA::subtract(vect[k], vect[k], tempmz1);
1999 }
2000 //
2001 // end of orthogonalization
2002 //
2003 Fc = false;
2004 korr = false;
2005 //
2006 // begin of reduction
2007 //
2008 debug_handler("bigfloat_lattice_gensys", "Tr_lll_trans_var_bfl(Bi, Tr, Factor, bit_prec) [3]");
2009 if (!repeat_loop) {
2010 for (lidia_size_t j = k-1; j >= 0; j--)
2011 if (abs(my[j][k]).compare(halb) > 0) {
2012 ++reduction_steps;
2013 korr = true;
2014 round (tempmz0, my[j][k]);
2015 if (abs(tempmz0).compare(bound) > 0)
2016 Fc = true;
2017
2018 tempmz0.bigintify(tempbiz);
2019 bfl_scalmul_bin(tempvect1, tempbiz, BiAddr[j]);
2020 bfl_subtract_bin(BiAddr[k], BiAddr[k], tempvect1);
2021 vectsize = rows;
2022 bfl_scalmul_bin(tempvect1, tempbiz, T[j]);
2023 bfl_subtract_bin(T[k], T[k], tempvect1);
2024 vectsize = columns;
2025
2026 for (lidia_size_t i = 0; i < j; ++i) {
2027 LiDIA::multiply(tempmz1, tempmz0, my[i][j]);
2028 LiDIA::subtract(my[i][k], my[i][k], tempmz1);
2029 }
2030 LiDIA::subtract(my[j][k], my[j][k], tempmz0);
2031
2032 }
2033 }
2034
2035 //
2036 // Correction Step II
2037 // lattice=(bigfloat )value
2038 //
2039 if (korr) {
2040 ++correction_steps;
2041 for (lidia_size_t j = 0; j < columns; j++) {
2042 bi_bit_len = BiAddr[k][j].bit_length();
2043 if (bi_bit_len > bit_prec) {
2044 bit_diff = bi_bit_len-bit_prec;
2045 shift_right(tempbiz, BiAddr[k][j], bit_diff);
2046 tempmz0.assign(tempbiz);
2047 }
2048 else {
2049 bit_diff = 0;
2050 tempmz0.assign(BiAddr[k][j]);
2051 }
2052 shift_right(lattice[k][j], tempmz0, bit_len-bit_diff);
2053 }
2054 }
2055
2056 //
2057 // end of reduction
2058 //
2059 if (Fc) {
2060 if (k > 1)
2061 --k;
2062 }
2063 else {
2064 Fc = (!repeat_loop) && korr;
2065 repeat_loop = !repeat_loop;
2066 }
2067 }
2068
2069 debug_handler("bigfloat_lattice_gensys", "Tr_lll_trans_var_bfl(Bi, Tr, Factor, bit_prec) [4]");
2070 //
2071 // Check lll - condition for parameter y_par
2072 //
2073 tempmz0.assign(y_par);
2074 LiDIA::multiply(tempmz0, tempmz0, vect[k-1]);
2075 LiDIA::square(tempmz1, my[k-1][k]);
2076 LiDIA::multiply(tempmz1, tempmz1, vect[k-1]);
2077 LiDIA::add(tempmz1, tempmz1, vect[k]);
2078
2079 if (tempmz0.compare(tempmz1) > 0) {
2080 //
2081 // ::swap vectors k and k-1
2082 //
2083 ++swapping_steps;
2084 bfl_swap_bin(BiAddr[k], BiAddr[k-1]);
2085 bfl_swap_bfl(lattice[k], lattice[k-1]);
2086 bfl_swap_bin(T[k], T[k-1]);
2087
2088 if (k > 1)
2089 --k;
2090 if (k == 1) {
2091 bfl_scalprod_bfl(vect[0], lattice[0], lattice[0]);
2092 sqrt(vect_sq[0], vect[0]);
2093 }
2094 }
2095 else
2096 ++k;
2097 }
2098 //
2099 // restore precision and
2100 // generate bigfloats
2101 //
2102 bigfloat::set_precision(old_prec);
2103 for (lidia_size_t i = 0; i < rows; i++)
2104 for (lidia_size_t j = 0; j < columns; j++) {
2105 bi_bit_len = BiAddr[i][j].bit_length();
2106 if (bi_bit_len > old_bit_prec) {
2107 bit_diff = bi_bit_len-old_bit_prec;
2108 shift_right(tempbiz, BiAddr[i][j], bit_diff);
2109 tempmz0.assign(tempbiz);
2110 }
2111 else {
2112 bit_diff = 0;
2113 tempmz0.assign(BiAddr[i][j]);
2114 }
2115 shift_right(value[i][j], tempmz0, bit_len-bit_diff);
2116 }
2117
2118 //
2119 // Copy into math_matrix< bigint > Tr
2120 //
2121 Tr.set_no_of_rows(rows);
2122 Tr.set_no_of_columns(rows);
2123 TrAddr = Tr.get_data_address();
2124 for (lidia_size_t i = 0; i < rows; i++)
2125 for (lidia_size_t j = 0; j < rows; j++)
2126 if (trans_flag)
2127 LiDIA::swap(TrAddr[j][i], T[i][j]);
2128 else
2129 LiDIA::swap(TrAddr[j][i], T[j][i]);
2130
2131 //
2132 // Free allocated storage
2133 //
2134 delete[] Tdel;
2135 delete[] T;
2136 delete[] mydel;
2137 delete[] my;
2138 }
2139
2140
2141
Tr_lll_var_bfl_gensys(math_matrix<bigint> & Bi,sdigit bit_len,lidia_size_t & rank,sdigit bit_prec)2142 void bigfloat_lattice_gensys::Tr_lll_var_bfl_gensys(math_matrix< bigint > & Bi,
2143 sdigit bit_len,
2144 lidia_size_t& rank, sdigit bit_prec)
2145 {
2146 debug_handler("bigfloat_lattice_gensys", "Tr_lll_var_bfl_gensys(Bi, Factor, rank, bit_len, bit_prec) [1]");
2147 bool Fc;
2148 bool korr;
2149 bool break_flag;
2150 bool vector_is_zero;
2151 bool repeat_loop;
2152 sdigit old_prec;
2153 sdigit old_bit_prec;
2154 sdigit bit_diff;
2155 sdigit bi_bit_len;
2156 sdigit short_prec;
2157 lidia_size_t k;
2158 lidia_size_t mom_rows = rows;
2159 bigint tempbiz;
2160 bigint *tempvect1;
2161 bigint **BiAddr;
2162 bigfloat bound = std::exp(std::log(2.0)*bit_prec/2);
2163 bigfloat invbound;
2164 bigfloat *mydel;
2165 bigfloat **my;
2166 bigfloat **lattice;
2167 bigfloat *vect;
2168 bigfloat *vect_sq;
2169 bigfloat halb(0.5);
2170
2171 LiDIA::divide(invbound, 1, bound);
2172
2173 //
2174 // Compute needed precision for exact arithmetik
2175 //
2176 old_prec = bigfloat::get_precision();
2177 old_bit_prec = static_cast<sdigit>(static_cast<double>(old_prec)/std::log(2.0)*std::log(10.0)+1);
2178 short_prec = static_cast<sdigit>(static_cast<double>(bit_prec)*std::log(2.0)/std::log(10.0));
2179 BiAddr = Bi.get_data_address();
2180
2181 vectsize = columns;
2182
2183 //
2184 // Allocating Memory for matrix
2185 //
2186 my = new bigfloat*[rows*2];
2187 memory_handler(my, "bigfloat_lattice_gensys", "Tr_lll_var_bfl_gensys(Bi, Factor, rank, bit_prec) :: "
2188 "not enough memory !");
2189 lattice = &my[rows];
2190 //
2191 // think again
2192 //
2193 mydel = new bigfloat[rows*(rows+columns+2)+columns];
2194 memory_handler(mydel, "bigfloat_lattice_gensys", "Tr_lll_var_bfl_gensys(Bi, Factor, rank, bit_prec) :: "
2195 "not enough memory !");
2196 for (lidia_size_t i = 0; i < rows; i++) {
2197 my[i] = &mydel[i*rows]; // rows x rows
2198 lattice[i] = &mydel[(rows*rows)+(i*columns)]; // rows x columns
2199 }
2200
2201 //
2202 // Allocating Memory for vector
2203 //
2204
2205 //
2206 // think again Part II
2207 //
2208 tempvect1 = new bigint[columns];
2209 memory_handler(tempvect1, "bigfloat_lattice_gensys", "Tr_lll_var_bfl_gensys(Bi, Factor, rank, bit_prec) :: "
2210 "not enough memory !");
2211
2212 vect = &mydel[(rows+columns)*rows];
2213 vect_sq = &vect[rows];
2214
2215
2216 //
2217 // Start
2218 //
2219 // Copy into floating point
2220 //
2221 bigfloat::set_precision(short_prec);
2222 for (lidia_size_t i = 0; i < rows; i++)
2223 for (lidia_size_t j = 0; j < columns; j++) {
2224 bi_bit_len = BiAddr[i][j].bit_length();
2225 if (bi_bit_len > bit_prec) {
2226 bit_diff = bi_bit_len-bit_prec;
2227 shift_right(tempbiz, BiAddr[i][j], bit_diff);
2228 tempmz0.assign(tempbiz);
2229 }
2230 else {
2231 bit_diff = 0;
2232 tempmz0.assign(BiAddr[i][j]);
2233 }
2234 shift_right(lattice[i][j], tempmz0, bit_len-bit_diff);
2235 }
2236
2237 k = 1;
2238 reduction_steps = 0;
2239 swapping_steps = 0;
2240 correction_steps = 0;
2241 vectsize = columns;
2242 break_flag = false;
2243 bfl_assign_zero_bfl(vect);
2244 bfl_scalprod_bfl(vect[0], lattice[0], lattice[0]);
2245
2246 while (k < mom_rows) {
2247 Fc = true;
2248 repeat_loop = false;
2249 while (Fc) {
2250 //
2251 // begin of orthogonalization
2252 //
2253 debug_handler("bigfloat_lattice_gensys", "Tr_lll_var_bfl_gensys(Bi, Factor, rank, bit_prec) [2]");
2254 bfl_scalprod_bfl(vect[k], lattice[k], lattice[k]);
2255 sqrt(vect_sq[k], vect[k]);
2256 for (lidia_size_t j = 0; j < k; ++j) {
2257 bfl_scalprod_bfl(tempmz0, lattice[k], lattice[j]);
2258 LiDIA::multiply(tempmz1, invbound, vect_sq[k]);
2259 LiDIA::multiply(tempmz1, tempmz1, vect_sq[j]);
2260 //
2261 // Correction Step I, only if needed, |lattice[k]*lattice[j]| to big
2262 //
2263 if (abs(tempmz0).compare(tempmz1) < 0) {
2264 //
2265 // Compute the correct scalar product
2266 //
2267 bfl_scalprod_bin(tempbiz, BiAddr[k], BiAddr[j]);
2268 bi_bit_len = tempbiz.bit_length();
2269 if (bi_bit_len > bit_prec) {
2270 bit_diff = bi_bit_len-bit_prec;
2271 shift_right(tempbiz, tempbiz, bit_diff);
2272 tempmz0.assign(tempbiz);
2273 }
2274 else {
2275 bit_diff = 0;
2276 tempmz0.assign(tempbiz);
2277 }
2278 shift_right(my[j][k], tempmz0, 2*bit_len-bit_diff);
2279 // ::divide(my[j][k], tempbiz, SqrFactor);
2280 }
2281 else
2282 my[j][k].assign(tempmz0);
2283
2284 for (lidia_size_t i = 0; i < j; ++i) {
2285 LiDIA::multiply(tempmz1, my[i][j], my[i][k]);
2286 LiDIA::multiply(tempmz1, tempmz1, vect[i]);
2287 LiDIA::subtract(my[j][k], my[j][k], tempmz1);
2288 }
2289
2290 tempmz0.assign(my[j][k]);
2291 LiDIA::divide(my[j][k], my[j][k], vect[j]);
2292 LiDIA::multiply(tempmz1, my[j][k], tempmz0);
2293 LiDIA::subtract(vect[k], vect[k], tempmz1);
2294 }
2295 //
2296 // end of orthogonalization
2297 //
2298 Fc = false;
2299 korr = false;
2300 //
2301 // begin of reduction
2302 //
2303 debug_handler("bigfloat_lattice_gensys", "Tr_lll_var_bfl_gensys(Bi, Factor, rank, bit_prec) [3]");
2304 if (!repeat_loop) {
2305 for (lidia_size_t j = k-1; j >= 0; j--)
2306 if (abs(my[j][k]).compare(halb) > 0) {
2307 ++reduction_steps;
2308 korr = true;
2309 round (tempmz0, my[j][k]);
2310 if (abs(tempmz0).compare(bound) > 0)
2311 Fc = true;
2312
2313 tempmz0.bigintify(tempbiz);
2314 bfl_scalmul_bin(tempvect1, tempbiz, BiAddr[j]);
2315 bfl_subtract_bin(BiAddr[k], BiAddr[k], tempvect1);
2316
2317 for (lidia_size_t i = 0; i < j; ++i) {
2318 LiDIA::multiply(tempmz1, tempmz0, my[i][j]);
2319 LiDIA::subtract(my[i][k], my[i][k], tempmz1);
2320 }
2321 LiDIA::subtract(my[j][k], my[j][k], tempmz0);
2322
2323 }
2324 }
2325
2326 //
2327 // Correction Step II
2328 // lattice=(bigfloat )value
2329 //
2330 if (korr) {
2331 ++correction_steps;
2332 for (lidia_size_t j = 0; j < columns; j++)
2333 if (!(BiAddr[k][j].is_zero())) {
2334 bi_bit_len = BiAddr[k][j].bit_length();
2335 if (bi_bit_len > bit_prec) {
2336 bit_diff = bi_bit_len-bit_prec;
2337 shift_right(tempbiz, BiAddr[k][j], bit_diff);
2338 tempmz0.assign(tempbiz);
2339 }
2340 else {
2341 bit_diff = 0;
2342 tempmz0.assign(BiAddr[k][j]);
2343 }
2344 shift_right(lattice[k][j], tempmz0, bit_len-bit_diff);
2345 }
2346 else
2347 lattice[k][j].assign_zero();
2348 }
2349
2350 //
2351 // Check if it`s 0 - vector
2352 // delete then
2353 //
2354 vector_is_zero = true;
2355 for (lidia_size_t j = 0; j < columns; j++)
2356 if (!(BiAddr[k][j].is_zero())) {
2357 vector_is_zero = false;
2358 break;
2359 }
2360
2361 if (vector_is_zero) {
2362 for (lidia_size_t j = k; j < mom_rows-1; j++) {
2363 bfl_swap_bfl(lattice[j], lattice[j+1]);
2364 bfl_swap_bin(BiAddr[j], BiAddr[j+1]);
2365 }
2366 bfl_assign_zero_bin(BiAddr[--mom_rows]);
2367 bfl_assign_zero_bfl(lattice[mom_rows]);
2368 k = 1;
2369 bfl_assign_zero_bfl(vect);
2370 bfl_scalprod_bfl(vect[0], lattice[0], lattice[0]);
2371 sqrt(vect_sq[0], vect[0]);
2372 vector_is_zero = false;
2373 break_flag = true;
2374 korr = false;
2375 Fc = false;
2376 break;
2377 }
2378
2379 //
2380 // end of reduction
2381 //
2382 if (Fc) {
2383 if (k > 1)
2384 --k;
2385 }
2386 else {
2387 Fc = (!repeat_loop) && korr;
2388 repeat_loop = !repeat_loop;
2389 }
2390 }
2391
2392 if (!break_flag) {
2393 debug_handler("bigfloat_lattice_gensys", "Tr_lll_var_bfl_gensys(Bi, Factor, rank, bit_prec) [4]");
2394 //
2395 // Check lll - condition for parameter y_par
2396 //
2397 tempmz0.assign(y_par);
2398 LiDIA::multiply(tempmz0, tempmz0, vect[k-1]);
2399 LiDIA::square(tempmz1, my[k-1][k]);
2400 LiDIA::multiply(tempmz1, tempmz1, vect[k-1]);
2401 LiDIA::add(tempmz1, tempmz1, vect[k]);
2402
2403 if (tempmz0.compare(tempmz1) > 0) {
2404 //
2405 // ::swap vectors k and k-1
2406 //
2407 ++swapping_steps;
2408 bfl_swap_bin(BiAddr[k], BiAddr[k-1]);
2409 bfl_swap_bfl(lattice[k], lattice[k-1]);
2410
2411 if (k > 1)
2412 --k;
2413 if (k == 1) {
2414 bfl_scalprod_bfl(vect[0], lattice[0], lattice[0]);
2415 sqrt(vect_sq[0], vect[0]);
2416 }
2417 }
2418 else
2419 ++k;
2420 }
2421 break_flag = false;
2422 }
2423 rank = mom_rows;
2424 //
2425 // restore precision and
2426 // generate bigfloats
2427 //
2428 bigfloat::set_precision(old_prec);
2429 for (lidia_size_t i = 0; i < rows; i++)
2430 for (lidia_size_t j = 0; j < columns; j++) {
2431 bi_bit_len = BiAddr[i][j].bit_length();
2432 if (bi_bit_len > old_bit_prec) {
2433 bit_diff = bi_bit_len-old_bit_prec;
2434 shift_right(tempbiz, BiAddr[i][j], bit_diff);
2435 tempmz0.assign(tempbiz);
2436 }
2437 else {
2438 bit_diff = 0;
2439 tempmz0.assign(BiAddr[i][j]);
2440 }
2441 shift_right(value[i][j], tempmz0, bit_len-bit_diff);
2442 }
2443
2444 //
2445 // Free allocated storage
2446 //
2447 bigfloat::set_precision(old_prec);
2448 delete[] tempvect1;
2449 delete[] mydel;
2450 delete[] my;
2451 }
2452
2453
2454
Tr_lll_trans_var_bfl_gensys(math_matrix<bigint> & Bi,math_matrix<bigint> & Tr,sdigit bit_len,lidia_size_t & rank,sdigit bit_prec)2455 void bigfloat_lattice_gensys::Tr_lll_trans_var_bfl_gensys(math_matrix< bigint > & Bi,
2456 math_matrix< bigint > &Tr,
2457 sdigit bit_len,
2458 lidia_size_t& rank, sdigit bit_prec)
2459 {
2460 debug_handler("bigfloat_lattice_gensys", "Tr_lll_trans_var_bfl_gensys(Bi, Tr, "
2461 " Factor, bit_len, rank, bit_prec) [1]");
2462 bool Fc;
2463 bool korr;
2464 bool break_flag;
2465 bool vector_is_zero;
2466 bool repeat_loop;
2467 sdigit old_prec;
2468 sdigit old_bit_prec;
2469 sdigit bit_diff;
2470 sdigit bi_bit_len;
2471 sdigit short_prec;
2472 lidia_size_t k = 1;
2473 lidia_size_t mom_rows = rows;
2474 bigint tempbiz;
2475 bigint *tempvect1;
2476 bigint **BiAddr;
2477 bigint *Tdel;
2478 bigint **T;
2479 bigint **TrAddr;
2480 bigfloat bound = std::exp(std::log(2.0)*bit_prec/2);
2481 bigfloat invbound;
2482 bigfloat *mydel;
2483 bigfloat **my;
2484 bigfloat **lattice;
2485 bigfloat *vect;
2486 bigfloat *vect_sq;
2487 bigfloat halb(0.5);
2488
2489 LiDIA::divide(invbound, 1, bound);
2490
2491 //
2492 // Compute needed precision for exact arithmetik
2493 //
2494 old_prec = bigfloat::get_precision();
2495 old_bit_prec = static_cast<sdigit>(static_cast<double>(old_prec)/std::log(2.0)*std::log(10.0)+1);
2496 short_prec = static_cast<sdigit>(static_cast<double>(bit_prec)*std::log(2.0)/std::log(10.0));
2497 BiAddr = Bi.get_data_address();
2498
2499 vectsize = columns;
2500
2501 //
2502 // Allocating Memory for matrix
2503 //
2504 my = new bigfloat*[rows*2];
2505 memory_handler(my, "bigfloat_lattice_gensys", "Tr_lll_trans_var_bfl_gensys(Bi, Tr, Factor, bit_prec) :: "
2506 "not enough memory !");
2507 lattice = &my[rows];
2508 T = new bigint*[rows];
2509 memory_handler(T, "bigfloat_lattice_gensys", "Tr_lll_trans_var_bfl_gensys(Bi, Tr, Factor, bit_prec) :: "
2510 "Not enough memory !");
2511
2512
2513 //
2514 // think again
2515 //
2516 mydel = new bigfloat[rows*(rows+columns+2)+columns];
2517 memory_handler(mydel, "bigfloat_lattice_gensys", "Tr_lll_trans_var_bfl_gensys(Bi, Tr, Factor, bit_prec) :: "
2518 "not enough memory !");
2519 Tdel = new bigint[rows*rows+((rows > columns)?rows:columns)];
2520 memory_handler(Tdel, "bigfloat_lattice_gensys", "Tr_lll_trans_var_bfl_gensys(Bi, Tr, Factor, bit_prec) :: "
2521 "Not enough memory !");
2522
2523 for (lidia_size_t i = 0; i < rows; i++) {
2524 my[i] = &mydel[i*rows]; // rows x rows
2525 T[i] = &Tdel[i*rows];
2526 T[i][i].assign_one();
2527 lattice[i] = &mydel[(rows*rows)+(i*columns)]; // rows x columns
2528 }
2529
2530 //
2531 // Allocating Memory for vector
2532 //
2533
2534 //
2535 // think again Part II
2536 //
2537
2538 tempvect1 = &Tdel[rows*rows];
2539 vect = &mydel[(rows+columns)*rows];
2540 vect_sq = &vect[rows];
2541
2542 //
2543 // Start
2544 //
2545 // Copy into floating point
2546 //
2547 bigfloat::set_precision(short_prec);
2548 for (lidia_size_t i = 0; i < rows; i++)
2549 for (lidia_size_t j = 0; j < columns; j++) {
2550 bi_bit_len = BiAddr[i][j].bit_length();
2551 if (bi_bit_len > bit_prec) {
2552 bit_diff = bi_bit_len-bit_prec;
2553 shift_right(tempbiz, BiAddr[i][j], bit_diff);
2554 tempmz0.assign(tempbiz);
2555 }
2556 else {
2557 bit_diff = 0;
2558 tempmz0.assign(BiAddr[i][j]);
2559 }
2560 shift_right(lattice[i][j], tempmz0, bit_len-bit_diff);
2561 }
2562
2563 k = 1;
2564 reduction_steps = 0;
2565 swapping_steps = 0;
2566 correction_steps = 0;
2567 vectsize = columns;
2568 break_flag = false;
2569 bfl_assign_zero_bfl(vect);
2570 bfl_scalprod_bfl(vect[0], lattice[0], lattice[0]);
2571
2572 while (k < mom_rows) {
2573 Fc = true;
2574 repeat_loop = false;
2575 while (Fc) {
2576 //
2577 // begin of orthogonalization
2578 //
2579 debug_handler("bigfloat_lattice_gensys", "Tr_lll_trans_var_bfl_gensys(Bi, Factor, bit_prec) [2]");
2580 bfl_scalprod_bfl(vect[k], lattice[k], lattice[k]);
2581 sqrt(vect_sq[k], vect[k]);
2582 for (lidia_size_t j = 0; j < k; ++j) {
2583 bfl_scalprod_bfl(tempmz0, lattice[k], lattice[j]);
2584 LiDIA::multiply(tempmz1, invbound, vect_sq[k]);
2585 LiDIA::multiply(tempmz1, tempmz1, vect_sq[j]);
2586 //
2587 // Correction Step I, only if needed, |lattice[k]*lattice[j]| to big
2588 //
2589 if (abs(tempmz0).compare(tempmz1) < 0) {
2590 //
2591 // Compute the correct scalar product
2592 //
2593 bfl_scalprod_bin(tempbiz, BiAddr[k], BiAddr[j]);
2594 bi_bit_len = tempbiz.bit_length();
2595 if (bi_bit_len > bit_prec) {
2596 bit_diff = bi_bit_len-bit_prec;
2597 shift_right(tempbiz, tempbiz, bit_diff);
2598 tempmz0.assign(tempbiz);
2599 }
2600 else {
2601 bit_diff = 0;
2602 tempmz0.assign(tempbiz);
2603 }
2604 shift_right(my[j][k], tempmz0, 2*bit_len-bit_diff);
2605 }
2606 else
2607 my[j][k].assign(tempmz0);
2608
2609 for (lidia_size_t i = 0; i < j; ++i) {
2610 LiDIA::multiply(tempmz1, my[i][j], my[i][k]);
2611 LiDIA::multiply(tempmz1, tempmz1, vect[i]);
2612 LiDIA::subtract(my[j][k], my[j][k], tempmz1);
2613 }
2614
2615 tempmz0.assign(my[j][k]);
2616 LiDIA::divide(my[j][k], my[j][k], vect[j]);
2617 LiDIA::multiply(tempmz1, my[j][k], tempmz0);
2618 LiDIA::subtract(vect[k], vect[k], tempmz1);
2619 }
2620 //
2621 // end of orthogonalization
2622 //
2623 Fc = false;
2624 korr = false;
2625 //
2626 // begin of reduction
2627 //
2628 debug_handler("bigfloat_lattice_gensys", "Tr_lll_trans_var_bfl_gensys(Bi, Tr, Factor, bit_prec) [3]");
2629 if (!repeat_loop) {
2630 for (lidia_size_t j = k-1; j >= 0; j--)
2631 if (abs(my[j][k]).compare(halb) > 0) {
2632 ++reduction_steps;
2633 korr = true;
2634 round (tempmz0, my[j][k]);
2635 if (abs(tempmz0).compare(bound) > 0)
2636 Fc = true;
2637
2638 tempmz0.bigintify(tempbiz);
2639 bfl_scalmul_bin(tempvect1, tempbiz, BiAddr[j]);
2640 bfl_subtract_bin(BiAddr[k], BiAddr[k], tempvect1);
2641 vectsize = rows;
2642 bfl_scalmul_bin(tempvect1, tempbiz, T[j]);
2643 bfl_subtract_bin(T[k], T[k], tempvect1);
2644 vectsize = columns;
2645
2646 for (lidia_size_t i = 0; i < j; ++i) {
2647 LiDIA::multiply(tempmz1, tempmz0, my[i][j]);
2648 LiDIA::subtract(my[i][k], my[i][k], tempmz1);
2649 }
2650 LiDIA::subtract(my[j][k], my[j][k], tempmz0);
2651
2652 }
2653 }
2654
2655 //
2656 // Correction Step II
2657 // lattice=(bigfloat )value
2658 //
2659 if (korr) {
2660 ++correction_steps;
2661 for (lidia_size_t j = 0; j < columns; j++)
2662 if (!(BiAddr[k][j].is_zero())) {
2663 bi_bit_len = BiAddr[k][j].bit_length();
2664 if (bi_bit_len > bit_prec) {
2665 bit_diff = bi_bit_len-bit_prec;
2666 shift_right(tempbiz, BiAddr[k][j], bit_diff);
2667 tempmz0.assign(tempbiz);
2668 }
2669 else {
2670 bit_diff = 0;
2671 tempmz0.assign(BiAddr[k][j]);
2672 }
2673 shift_right(lattice[k][j], tempmz0, bit_len-bit_diff);
2674 }
2675 else
2676 lattice[k][j].assign_zero();
2677 bfl_scalquad_bfl(vect_sq[k], lattice[k]);
2678 sqrt(vect_sq[k], vect_sq[k]);
2679 korr = false;
2680 }
2681 //
2682 // Check if it`s 0 - vector
2683 // delete then
2684 //
2685 vector_is_zero = true;
2686 for (lidia_size_t j = 0; j < columns; j++)
2687 if (!(BiAddr[k][j].is_zero())) {
2688 vector_is_zero = false;
2689 break;
2690 }
2691
2692 if (vector_is_zero) {
2693 for (lidia_size_t j = k; j < mom_rows-1; j++) {
2694 bfl_swap_bfl(lattice[j], lattice[j+1]);
2695 bfl_swap_bin(BiAddr[j], BiAddr[j+1]);
2696 bfl_swap_bin(T[j], T[j+1]);
2697 }
2698 bfl_assign_zero_bin(BiAddr[--mom_rows]);
2699 bfl_assign_zero_bfl(lattice[mom_rows]);
2700 vectsize = rows;
2701 bfl_assign_zero_bin(T[mom_rows]);
2702 vectsize = columns;
2703 k = 1;
2704 bfl_assign_zero_bfl(vect);
2705 bfl_scalprod_bfl(vect[0], lattice[0], lattice[0]);
2706 sqrt(vect_sq[0], vect[0]);
2707 vector_is_zero = false;
2708 break_flag = true;
2709 korr = false;
2710 Fc = false;
2711 break;
2712 }
2713
2714 //
2715 // end of reduction
2716 //
2717 if (Fc) {
2718 if (k > 1)
2719 --k;
2720 }
2721 else {
2722 Fc = (!repeat_loop) && korr;
2723 repeat_loop = !repeat_loop;
2724 }
2725 }
2726
2727 if (!break_flag) {
2728 debug_handler("bigfloat_lattice_gensys", "Tr_lll_trans_var_bfl_gensys(Bi, Tr, Factor, bit_prec) [4]");
2729 //
2730 // Check lll - condition for parameter y_par
2731 //
2732 tempmz0.assign(y_par);
2733 LiDIA::multiply(tempmz0, tempmz0, vect[k-1]);
2734 LiDIA::square(tempmz1, my[k-1][k]);
2735 LiDIA::multiply(tempmz1, tempmz1, vect[k-1]);
2736 LiDIA::add(tempmz1, tempmz1, vect[k]);
2737
2738 if (tempmz0.compare(tempmz1) > 0) {
2739 //
2740 // ::swap vectors k and k-1
2741 //
2742 ++swapping_steps;
2743 bfl_swap_bin(BiAddr[k], BiAddr[k-1]);
2744 bfl_swap_bfl(lattice[k], lattice[k-1]);
2745 bfl_swap_bin(T[k], T[k-1]);
2746
2747 if (k > 1)
2748 --k;
2749 if (k == 1) {
2750 bfl_scalprod_bfl(vect[0], lattice[0], lattice[0]);
2751 sqrt(vect_sq[0], vect[0]);
2752 }
2753 }
2754 else
2755 ++k;
2756 }
2757 break_flag = false;
2758 }
2759 rank = mom_rows;
2760 //
2761 // restore precision and
2762 // generate bigfloats
2763 //
2764 bigfloat::set_precision(old_prec);
2765 for (lidia_size_t i = 0; i < rows; i++)
2766 for (lidia_size_t j = 0; j < columns; j++) {
2767 bi_bit_len = BiAddr[i][j].bit_length();
2768 if (bi_bit_len > old_bit_prec) {
2769 bit_diff = bi_bit_len-old_bit_prec;
2770 shift_right(tempbiz, BiAddr[i][j], bit_diff);
2771 tempmz0.assign(tempbiz);
2772 }
2773 else {
2774 bit_diff = 0;
2775 tempmz0.assign(BiAddr[i][j]);
2776 }
2777 shift_right(value[i][j], tempmz0, bit_len-bit_diff);
2778 }
2779
2780 //
2781 // Copy into math_matrix< bigint > Tr
2782 //
2783 Tr.set_no_of_rows(rows);
2784 Tr.set_no_of_columns(rows);
2785 TrAddr = Tr.get_data_address();
2786 for (lidia_size_t i = 0; i < rows; i++)
2787 for (lidia_size_t j = 0; j < rows; j++)
2788 if (trans_flag)
2789 LiDIA::swap(TrAddr[j][i], T[i][j]);
2790 else
2791 LiDIA::swap(TrAddr[j][i], T[j][i]);
2792
2793 //
2794 // Free allocated storage
2795 //
2796 delete[] Tdel;
2797 delete[] T;
2798 delete[] mydel;
2799 delete[] my;
2800 }
2801
2802
2803
2804 #ifdef LIDIA_NAMESPACE
2805 } // end of namespace LiDIA
2806 #endif
2807