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 : Patrick Theobald (PT)
14 // Changes : See CVS log
15 //
16 //==============================================================================================
17
18
19 #ifndef LIDIA_HNF_KERNEL_CC_GUARD_
20 #define LIDIA_HNF_KERNEL_CC_GUARD_
21
22
23 #ifndef LIDIA_INFO_H_GUARD_
24 # include "LiDIA/info.h"
25 #endif
26 #ifndef LIDIA_HNF_CONF_H_GUARD_
27 # include "LiDIA/matrix/hnf_conf.h"
28 #endif
29 #ifndef LIDIA_HNF_KERNEL_H_GUARD_
30 # include "LiDIA/matrix/hnf_kernel.h"
31 #endif
32
33 #if defined LIDIA_COMPUTE_TIME_INFO || defined LIDIA_COMPUTE_RUNTIME_INFO
34 # include <fstream>
35 #endif
36 #ifdef LIDIA_COMPUTE_TIME_INFO
37 # ifndef LIDIA_TIMER_H_GUARD_
38 # include "LiDIA/timer.h"
39 # endif
40 # ifndef LIDIA_BIGFLOAT_H_GUARD_
41 # include "LiDIA/bigfloat.h"
42 # endif
43 #endif
44
45
46
47 #ifdef LIDIA_NAMESPACE
48 # ifndef IN_NAMESPACE_LIDIA
49 namespace LiDIA {
50 # endif
51 #endif
52
53
54
55 #define DMESSAGE "hnf_kernel"
56
57
58 ////////////////////////////////////////////////////////////////////////////////
59 // HNF computation
60 ////////////////////////////////////////////////////////////////////////////////
61
62 //
63 // mgcd
64 //
65
66 template< class T, class REP, class REP1, class CONF >
67 inline bool havas_kernel< T, REP, REP1, CONF >::
mgcd(MR<T> & A,lidia_size_t & startr,lidia_size_t & startc,const T & BOUND)68 mgcd(MR< T > &A, lidia_size_t &startr, lidia_size_t &startc, const T &BOUND)
69 {
70 // Init
71 conf.init(A, BOUND);
72 --startc;
73 --startr;
74
75 conf.mgcd(A, startr, startc);
76 return true;
77 }
78
79
80
81 template< class T, class REP, class REP1, class CONF >
82 inline bool havas_kernel< T, REP, REP1, CONF >::
mgcd(MR<T> & A,matrix<bigint> & TR,lidia_size_t & startr,lidia_size_t & startc,const T & BOUND)83 mgcd(MR< T > &A, matrix< bigint > &TR, lidia_size_t &startr,
84 lidia_size_t &startc, const T &BOUND)
85 {
86 // Init
87 conf.init(A, BOUND);
88 --startr;
89 --startc;
90
91 conf.mgcd(A, TR, startr, startc);
92 return true;
93 }
94
95
96
97 //
98 // row elimination
99 //
100
101 template< class T, class REP, class REP1, class CONF >
102 bool havas_kernel< T, REP, REP1, CONF >::
hnf_Z1(MR<T> & A,lidia_size_t & startr,lidia_size_t & startc,const T & BOUND)103 hnf_Z1(MR< T > &A, lidia_size_t &startr, lidia_size_t &startc, const T &BOUND)
104 {
105
106 // Init
107 conf.init(A, BOUND);
108
109 #ifdef LIDIA_COMPUTE_RUNTIME_INFO
110 std::ofstream log1_ofs("maximaler_Eintrag.dat");
111 std::ofstream log2_ofs("Eintragsdichte.dat");
112 bigint Bound = 10;
113 power(Bound, Bound, 100);
114 #endif
115
116 #ifdef LIDIA_COMPUTE_TIME_INFO
117 std::ofstream log3_ofs("Time_pro_Iteration.dat");
118 timer t;
119 #endif
120
121 // row elimination
122 for (--startc, --startr; startr >= 0 && startc >= 0; startc--, startr--) {
123 lidia_xinfo_handler("stf", startr, A.rows);
124
125 #ifdef LIDIA_COMPUTE_RUNTIME_INFO
126 T MAX = 0, Durch = 0;
127 lidia_size_t no_of_elements;
128 modul.kennwerte(A, MAX, no_of_elements, Durch);
129
130 no_of_elements = 0;
131 for (lidia_size_t i = 0; i <= startc; i++)
132 no_of_elements += A.value_counter[i];
133
134 log1_ofs << A.rows - startr << " " << MAX << std::endl;
135 if (startr >= 0 && startc >= 0)
136 log2_ofs << A.rows - startr << " " << bigfloat(100 * no_of_elements)/bigfloat((startr+1) * (startc+1)) << std::endl;
137
138 #if 0
139 if (MAX > Bound) {
140 for (lidia_size_t i = A.rows - startr + 1; i < A.rows; i++) {
141 log1_ofs << i << " " << 0 << std::endl;
142 log2_ofs << i << " " << 0 << std::endl;
143 }
144 log1_ofs.close();
145 log2_ofs.close();
146 return true;
147 }
148 #endif
149 #endif
150
151 #ifdef LIDIA_COMPUTE_TIME_INFO
152 t.start_timer();
153 #endif
154
155 lidia_size_t st = conf.mgcd(A, startr, startc);
156
157 if (st == -1) {
158 startr++;
159 startc++;
160 return false;
161 }
162 else if (st == 0) // All elements are zero
163 startc++;
164 else {
165 if (modul.member(A, startr, startc) < A.Zero)
166 modul.negate_column(A, startc, startr);
167
168 if (!conf.normalize_row(A, startr, A.columns - 1, startc, A.rows - 1)) {
169 startr++;
170 startc++;
171 return false;
172 }
173 }
174 #ifdef LIDIA_COMPUTE_TIME_INFO
175 t.stop_timer();
176 log3_ofs << A.rows - 1 - startr << " " << " " << t.user_time() << std::endl;
177 #endif
178 }
179 startr++;
180 startc++;
181
182 #ifdef LIDIA_COMPUTE_TIME_INFO
183 log3_ofs.close();
184 #endif
185
186 #ifdef LIDIA_COMPUTE_RUNTIME_INFO
187 log1_ofs.close();
188 log2_ofs.close();
189 #endif
190
191 return true;
192 }
193
194
195
196 template< class T, class REP, class REP1, class CONF >
197 bool havas_kernel< T, REP, REP1, CONF >::
hnf_Z1(MR<T> & A,matrix<bigint> & TR,lidia_size_t & startr,lidia_size_t & startc,const T & BOUND)198 hnf_Z1(MR< T > &A, matrix< bigint > &TR, lidia_size_t &startr,
199 lidia_size_t &startc, const T &BOUND)
200 {
201 // Init
202 conf.init(A, BOUND);
203
204
205 // row elimination
206 for (--startc, --startr; startr >= 0 && startc >= 0; startc--, startr--) {
207 lidia_xinfo_handler("stf", startr, A.rows);
208
209 lidia_size_t st = conf.mgcd(A, TR, startr, startc);
210
211 if (st == -1) {
212 startr++;
213 startc++;
214 return false;
215 }
216 else if (st == 0) // All elements are zero
217 startc++;
218 else {
219 if (modul.member(A, startr, startc) < A.Zero) {
220 modul.negate_column(A, startc, startr);
221
222 for (lidia_size_t j = 0; j < A.columns; j++)
223 TR.sto(j, startc, -TR.member(j, startc));
224 }
225
226 if (!conf.normalize_row(A, TR, startr, A.columns - 1, startc, A.rows - 1)) {
227 startr++;
228 startc++;
229 return false;
230 }
231 }
232 }
233 startr++;
234 startc++;
235 return true;
236 }
237
238
239
240 template< class T, class REP, class REP1, class CONF >
241 bool havas_kernel< T, REP, REP1, CONF >::
hnf_Z2(MR<T> & A,lidia_size_t & startr,lidia_size_t & startc,const T & BOUND)242 hnf_Z2(MR< T > &A, lidia_size_t &startr, lidia_size_t &startc, const T &BOUND)
243 {
244
245 #ifdef LIDIA_COMPUTE_TIME_INFO
246 std::ofstream log3_ofs("Time_pro_Iteration.dat");
247 timer t;
248 #endif
249
250 // Init
251 conf.init(A, BOUND);
252 for (--startc, --startr; startr >= 0 && startc >= 0; startc--, startr--) {
253 lidia_xinfo_handler("stf", startr, A.rows);
254
255 #ifdef LIDIA_COMPUTE_TIME_INFO
256 t.start_timer();
257 #endif
258
259 // row elimination
260 lidia_size_t st = conf.mgcd(A, startr, startc);
261
262 // stop computation
263 if (st == -1) {
264 startr++;
265 startc++;
266 return false;
267 }
268 else if (st == 1) {
269 if (modul.member(A, startr, startc) < A.Zero)
270 modul.negate_column(A, startc, startr);
271 }
272
273 #ifdef LIDIA_COMPUTE_TIME_INFO
274 t.stop_timer();
275 log3_ofs << A.rows - 1 - startr << " " << " " << t.user_time() << std::endl;
276 #endif
277
278 }
279
280 #ifdef LIDIA_COMPUTE_TIME_INFO
281 log3_ofs.close();
282 #endif
283
284 startr++;
285 startc++;
286 return true;
287 }
288
289
290
291 template< class T, class REP, class REP1, class CONF >
292 bool havas_kernel< T, REP, REP1, CONF >::
hnf_Z2(MR<T> & A,matrix<bigint> & TR,lidia_size_t & startr,lidia_size_t & startc,const T & BOUND)293 hnf_Z2(MR< T > &A, matrix< bigint > &TR, lidia_size_t &startr, lidia_size_t &startc, const T &BOUND)
294 {
295 // Init
296 conf.init(A, BOUND);
297
298 // row elimination
299 for (--startc, --startr; startr >= 0 && startc >= 0; startc--, startr--) {
300 lidia_xinfo_handler("stf", startr, A.rows);
301
302 lidia_size_t st = conf.mgcd(A, TR, startr, startc);
303
304 if (st == -1) {
305 startr++;
306 startc++;
307 return false;
308 }
309 else if (st == 1) {
310 if (modul.member(A, startr, startc) < A.Zero) {
311 modul.negate_column(A, startc, startr);
312
313 for (lidia_size_t j = 0; j < A.columns; j++)
314 TR.sto(j, startc, -TR.member(j, startc));
315 }
316 }
317 }
318 startr++;
319 startc++;
320 return true;
321 }
322
323
324
325 template< class T, class REP, class REP1, class CONF >
326 inline bool havas_kernel< T, REP, REP1, CONF >::
stf(MR<T> & A,lidia_size_t & startr,lidia_size_t & startc,const T & BOUND)327 stf(MR< T > &A, lidia_size_t &startr, lidia_size_t &startc, const T &BOUND)
328 {
329 // Init
330 conf.init(A, BOUND);
331
332 // row elimination
333 for (--startc, --startr; startr >= 0 && startc >= 0; startc--, startr--) {
334 lidia_xinfo_handler("stf", startr, A.rows);
335
336 lidia_size_t st = conf.mgcd(A, startr, startc);
337 if (st == -1) {
338 startr++;
339 startc++;
340 return false;
341 }
342 else if (st == 0) // All elements are zero
343 startc++;
344 else {
345 if (modul.member(A, startr, startc) < A.Zero)
346 modul.negate_column(A, startc, startr);
347 }
348 }
349 startr++;
350 startc++;
351 return true;
352 }
353
354
355
356 template< class T, class REP, class REP1, class CONF >
357 bool havas_kernel< T, REP, REP1, CONF >::
stf(MR<T> & A,matrix<bigint> & TR,lidia_size_t & startr,lidia_size_t & startc,const T & BOUND)358 stf(MR< T > &A, matrix< bigint > &TR, lidia_size_t &startr,
359 lidia_size_t &startc, const T &BOUND)
360 {
361 // Init
362 conf.init(A, BOUND);
363
364 // row elimination
365 for (--startc, --startr; startr >= 0 && startc >= 0; startc--, startr--) {
366 lidia_xinfo_handler("stf", startr, A.rows);
367 lidia_size_t st = conf.mgcd(A, TR, startr, startc);
368
369 if (st == -1) {
370 startr++;
371 startc++;
372 return false;
373 }
374 else if (st == 0) // All elements are zero
375 startc++;
376 else {
377 if (modul.member(A, startr, startc) < A.Zero) {
378 modul.negate_column(A, startc, startr);
379 modul_TR.negate_column(TR, startc, TR.get_no_of_rows() - 1);
380 }
381 }
382 }
383 startr++;
384 startc++;
385 return true;
386 }
387
388
389
390 #if 0
391 template< class T, class REP, class REP1, class CONF >
392 bool havas_kernel< T, REP, REP1, CONF >::
393 hnf_Z2_mod(MR< T > &A, lidia_size_t &startr, lidia_size_t &startc, const T &DET)
394 {
395
396 #ifdef LIDIA_COMPUTE_TIME_INFO
397 std::ofstream log3_ofs("Time_pro_Iteration.dat");
398 timer t;
399 #endif
400
401 // Init
402 //conf.init(A, BOUND);
403 for (--startc, --startr; startr >= 0 && startc >= 0; startc--, startr--) {
404 lidia_xinfo_handler("stf", startr, A.rows);
405
406 #ifdef LIDIA_COMPUTE_TIME_INFO
407 t.start_timer();
408 #endif
409
410 // row elimination
411 lidia_size_t st = conf.mgcd(A, startr, startc);
412 modul.remainder(A, DET);
413
414 // stop computation
415 if (st == -1) {
416 startr++;
417 startc++;
418 return false;
419 }
420 else if (st == 1) {
421 if (modul.member(A, startr, startc) < A.Zero)
422 modul.negate_column_mod(A, startc, startr, DET);
423 }
424
425 #ifdef LIDIA_COMPUTE_TIME_INFO
426 t.stop_timer();
427 log3_ofs << A.rows - 1 - startr << " " << " " << t.user_time() << std::endl;
428 #endif
429
430 }
431
432 #ifdef LIDIA_COMPUTE_TIME_INFO
433 log3_ofs.close();
434 #endif
435
436 startr++;
437 startc++;
438 return true;
439 }
440 #endif
441
442
443
444 // ADDED BY MJJ
445
446 template< class T, class REP, class REP1, class CONF >
447 int havas_kernel< T, REP, REP1, CONF >::
hnf_Z2(MR<T> & A,trans_matrix & TR,matrix<bigint> & tran,lidia_size_t & startr,lidia_size_t & startc,lidia_size_t rowsorig,const T & BOUND)448 hnf_Z2(MR< T > &A, trans_matrix & TR, matrix< bigint > & tran,
449 lidia_size_t &startr, lidia_size_t &startc, lidia_size_t rowsorig,
450 const T &BOUND)
451 {
452 lidia_size_t n = 2;
453 lidia_size_t dense_stop;
454
455 if ((startr <= 500) &&
456 (tran.get_representation() == matrix_flags::sparse_representation))
457 return -1;
458
459 if (rowsorig >= 4000)
460 dense_stop = 1400;
461 else if (rowsorig >= 3000)
462 dense_stop = 800;
463 else if (rowsorig >= 2000)
464 dense_stop = 600;
465 else if (rowsorig >= 800)
466 dense_stop = 500;
467 else
468 dense_stop = -1;
469
470
471 tran.reset();
472 tran.resize(startc, startc);
473 tran.diag(1, 0);
474
475 conf.init(A, BOUND);
476
477 // Elimination
478 for (--startc, --startr; startr >= 0 && startc >= 0; startc--, startr--) {
479 lidia_xinfo_handler("stf", startr, A.rows);
480
481 lidia_size_t st = conf.mgcd(A, tran, startr, startc);
482 if (st == -1) {
483 startr++;
484 startc++;
485 return 0;
486 }
487 else if (st == 0) {
488 if (modul.member(A, startr, startc) < A.Zero) {
489 modul.negate_column(A, startc, startr);
490 modul_TR.negate_column(tran, startc, tran.get_no_of_rows() - 1);
491 }
492 }
493
494 if (tran.get_representation() == matrix_flags::sparse_representation) {
495 if (startr <= dense_stop) {
496 startr++;
497 startc++;
498 return -1;
499 }
500 }
501 else if ((A.bitfield.get_representation() == matrix_flags::sparse_representation) &&
502 ((startr <= (dense_stop-200)))) {
503 startr++;
504 startc++;
505 return 0;
506 }
507
508 if (startr == (rowsorig/n)) {
509 TR.store_matrix(tran);
510
511 tran.resize(startc+1, startc+1);
512 tran.diag(1, 0);
513
514 ++n;
515 }
516 }
517
518 startr++;
519 startc++;
520 return 1;
521 }
522
523
524
525 ////////////////////////////////////////////////////////////////////////////////
526 // kannan_kernel
527 ////////////////////////////////////////////////////////////////////////////////
528
529 //
530 // column step form
531 //
532
533 template< class T, class REP, class REP1, class CONF >
534 bool kannan_kernel< T, REP, REP1, CONF >::
hnf(MR<T> & A,lidia_size_t & startr,lidia_size_t & startc,const T & BOUND)535 hnf(MR< T > &A, lidia_size_t &startr, lidia_size_t &startc, const T &BOUND)
536 {
537 normalization_kernel< T, REP, REP > normalize_modul;
538
539 #ifdef LIDIA_COMPUTE_RUNTIME_INFO
540 std::ofstream log1_ofs("maximaler_Eintrag.dat");
541 std::ofstream log2_ofs("Eintragsdichte.dat");
542 bigint Bound = 10;
543 power(Bound, Bound, 100);
544 #endif
545
546 #ifdef LIDIA_COMPUTE_TIME_INFO
547 std::ofstream log3_ofs("Time_pro_Iteration.dat");
548 timer t;
549 #endif
550
551 lidia_size_t pos1 = 0, pos2 = 0;
552
553 T TMP;
554 T RES0, RES1, RES2, TMP1, x, y;
555 lidia_size_t i, j, l, k;
556
557 for (i = startc - 2, k = startr - 2; i >= 0; i--, k--) {
558
559 #ifdef LIDIA_COMPUTE_RUNTIME_INFO
560 T MAX = 0, Durch = 0;
561 lidia_size_t no_of_elements;
562 modul.kennwerte(A, MAX, no_of_elements, Durch);
563
564 no_of_elements = 0;
565 for (lidia_size_t f = 0; f <= startc; f++)
566 no_of_elements += A.value_counter[f];
567
568 log1_ofs << A.rows - k << " " << MAX << std::endl;
569 if (i >= 0 && k >= 0)
570 log2_ofs << A.rows - k << " " << bigfloat(100 * no_of_elements)/bigfloat((k+1) * (i+1)) << std::endl;
571 #endif
572
573 #ifdef LIDIA_COMPUTE_TIME_INFO
574 t.start_timer();
575 #endif
576
577 lidia_xinfo_handler("stf", k, A.rows);
578 for (j = startr - 1, l = startc - 1; j >= 0 && l > i; j--, l--) {
579 if (modul.member(A, j, i) != A.Zero) {
580 if (modul.member(A, j, l) == A.Zero)
581 modul.swap_columns(A, i, l);
582 else {
583 RES2 = xgcd(RES0, RES1, modul.member(A, j, i), modul.member(A, j, l));
584 x = modul.member(A, j, l) / RES2;
585 y = modul.member(A, j , i) / RES2;
586
587 for (lidia_size_t len = 0; len <= j; len++) {
588 TMP = modul.member(A, len, i);
589 TMP1 = modul.member(A, len, l);
590
591 modul.sto(A, len, i, TMP*x - TMP1*y);
592 modul.sto(A, len, l, TMP*RES0 + TMP1*RES1);
593 }
594
595 }
596 }
597 }
598 if (k >= 0 && i >= 0) {
599 pos1 = k;
600 pos2 = i;
601 if (modul.member(A, pos1, pos2) < A.Zero) {
602 for (lidia_size_t len = 0; len <= pos1; len++)
603 modul.sto(A, len, pos2, -modul.member(A, len, pos2));
604 }
605 }
606 conf_modul.normalize(A, pos1, pos2);
607
608 #ifdef LIDIA_COMPUTE_TIME_INFO
609 t.stop_timer();
610 log3_ofs << A.rows - 1 - k << " " << " " << t.user_time() << std::endl;
611 #endif
612
613 }
614
615 #ifdef LIDIA_COMPUTE_TIME_INFO
616 log3_ofs.close();
617 #endif
618
619 #ifdef LIDIA_COMPUTE_RUNTIME_INFO
620 log1_ofs.close();
621 log2_ofs.close();
622 #endif
623
624 return true;
625 }
626
627
628
629 #undef DMESSAGE
630
631
632
633 #ifdef LIDIA_NAMESPACE
634 # ifndef IN_NAMESPACE_LIDIA
635 } // end of namespace LiDIA
636 # endif
637 #endif
638
639
640
641 #endif // LIDIA_HNF_KERNEL_CC_GUARD_
642