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 #ifdef HAVE_CONFIG_H
20 # include "config.h"
21 #endif
22 #include "LiDIA/timer.h"
23 #include "LiDIA/error.h"
24 #include "LiDIA/info.h"
25 #include "LiDIA/bigint_matrix.h"
26 #include "LiDIA/base_power_product.h"
27 #include "LiDIA/matrix/dense_bigint_matrix_modules.h"
28 #include "LiDIA/matrix/sparse_bigint_matrix_modules.h"
29 #include "LiDIA/matrix/hnf_conf.h"
30 #include "LiDIA/matrix/hnf_kernel.h"
31 #include "LiDIA/matrix/crt_and_prime_handling.h"
32 #include "LiDIA/matrix/dense_fp_matrix_kernel.h"
33 #include "LiDIA/matrix/sparse_fp_matrix_kernel.h"
34 #include "LiDIA/matrix/sparse_fp_matrix_algorithms.h"
35 #include "LiDIA/matrix/bigint_matrix_algorithms.h"
36 #include "LiDIA/matrix/modular_arithmetic.h"
37
38
39
40 #ifdef LIDIA_NAMESPACE
41 namespace LiDIA {
42 #endif
43
44
45
46 //
47 // debug defines / error defines
48 //
49
50 extern const char *PRT;
51 extern const char *matrix_error_msg[];
52
53 #define DVALUE LDBL_MATRIX // Debug value
54 #define DMESSAGE "bigint_matrix" // Debug message
55 #define EMESSAGE matrix_error_msg // Error message
56
57 //
58 // debug level
59 //
60 // 0 : remainder
61 // 1 : divide
62 // 2 : norms and bounds
63 // 3 : randomize
64 // 4 : Linear algebra PART 1
65 // 5 : Linear algebra PART 2
66 //
67
68 #define DRMKex DRMK< bigint >
69 #define SRMKex SRMK< bigint >
70
71 //
72 // constructor
73 //
74
75 matrix< bigint >::
matrix(const base_matrix<long> & M)76 matrix(const base_matrix< long > &M)
77 {
78 if (M.bitfield.get_representation() == matrix_flags::dense_representation) {
79 rows = M.rows;
80 sparse_rows = M.sparse_rows;
81
82 columns = M.columns;
83 sparse_columns = M.sparse_columns;
84
85 allocated = NULL;
86 index = NULL;
87 value_counter = NULL;
88
89 if (M.rows == 0)
90 value = NULL;
91 else {
92 value = new bigint *[rows];
93 memory_handler(value, DMESSAGE,
94 "constructor(MR< T > &, "
95 "const MR< T > &) :: "
96 "Error in memory allocation (value)");
97 }
98
99 register bigint *Atmp;
100 register long *Mtmp;
101 for (register lidia_size_t i = 0; i < rows; i++) {
102 Atmp = new bigint[columns];
103 Mtmp = M.value[i];
104 memory_handler(Atmp, DMESSAGE,
105 "constructor(MR< T > &, "
106 "const MR< T > &) :: "
107 "Error in memory allocation (Atmp)");
108
109 for (register lidia_size_t j = 0; j < columns; j++)
110 Atmp[j] = bigint(Mtmp[j]);
111 value[i] = Atmp;
112 }
113 }
114 else {
115 rows = M.rows;
116 columns = M.columns;
117
118 sparse_rows = M.sparse_rows;
119 sparse_columns = M.sparse_columns;
120
121 value = new bigint *[rows];
122 memory_handler(value, DMESSAGE,
123 "constructor((MR< T > &, const MR< T > &) ::"
124 "Error in memory allocation (value)");
125
126 index = new lidia_size_t *[rows];
127 memory_handler(index, DMESSAGE,
128 "constructor((MR< T > &, const MR< T > &) ::"
129 "Error in memory allocation (index)");
130
131 value_counter = new lidia_size_t[rows];
132 memory_handler(value_counter, DMESSAGE,
133 "constructor((MR< T > &, const MR< T > &) ::"
134 "Error in memory allocation (value_counter)");
135
136 allocated = new lidia_size_t[rows];
137 memory_handler(allocated, DMESSAGE,
138 "constructor((MR< T > &, const MR< T > &) ::"
139 "Error in memory allocation (allocated)");
140
141 for (register lidia_size_t i = 0; i < rows; i++) {
142 register lidia_size_t size = M.allocated[i];
143 if (size) {
144 index[i] = new lidia_size_t[size];
145 memory_handler(index[i], DMESSAGE,
146 "constructor((MR< T > &, const MR< T > &) ::"
147 "Error in memory allocation (index[i])");
148 value[i] = new bigint[size];
149 memory_handler(value[i], DMESSAGE,
150 "constructor((MR< T > &, const MR< T > &) ::"
151 "Error in memory allocation (value[i])");
152 value_counter[i] = M.value_counter[i];
153 allocated[i] = size;
154 for (register lidia_size_t p = 0; p < value_counter[i]; p++) {
155 value[i][p] = bigint(M.value[i][p]);
156 index[i][p] = M.index[i][p];
157 }
158 }
159 }
160 }
161
162 bitfield = M.bitfield;
163 }
164
165
166
167 //
168 // remainder
169 //
170
171 void
remainder(matrix<bigint> & RES,const matrix<bigint> & M,const bigint & mod)172 remainder(matrix< bigint > &RES, const matrix< bigint > & M, const bigint & mod)
173 {
174 //
175 // Task: RES.remainder(M,mod);
176 // => RES = M % mod,
177 // Version: 2.0
178 //
179
180 debug_handler_l(DMESSAGE, "remainder(matrix< bigint >, const matrix< bigint > &, "
181 "const bigint &)", DVALUE);
182
183 const bigint_matrix_algorithms< DRMKex, DRMKex, DRMKex > DDD_bigint_modul;
184 const bigint_matrix_algorithms< SRMKex, DRMKex, DRMKex > SDD_bigint_modul;
185 const bigint_matrix_algorithms< DRMKex, SRMKex, DRMKex > DSD_bigint_modul;
186 const bigint_matrix_algorithms< SRMKex, SRMKex, DRMKex > SSD_bigint_modul;
187
188 if (RES.rows != M.rows)
189 RES.set_no_of_rows(M.rows);
190 if (RES.columns != M.columns)
191 RES.set_no_of_columns(M.columns);
192
193 if (M.bitfield.get_representation() == matrix_flags::dense_representation) {
194 if (RES.bitfield.get_representation() == matrix_flags::dense_representation)
195 DDD_bigint_modul.remainder(RES, M, mod);
196 else
197 SDD_bigint_modul.remainder(RES, M, mod);
198 }
199 else {
200 if (RES.bitfield.get_representation() == matrix_flags::dense_representation)
201 DSD_bigint_modul.remainder(RES, M, mod);
202 else
203 SSD_bigint_modul.remainder(RES, M, mod);
204 }
205 }
206
207
208
209 void
remainder(matrix<long> & RES,const matrix<bigint> & M,long mod)210 remainder(matrix< long > &RES, const matrix< bigint > & M, long mod)
211 {
212 //
213 // Task: RES.remainder(M,mod);
214 // => RES = M % mod,
215 // Version: 2.0
216 //
217
218 debug_handler_l(DMESSAGE, "remainder(matrix< long > &, const matrix< bigint > &, "
219 "const bigint &)", DVALUE);
220
221 const bigint_matrix_algorithms< DRMKex, DRMKex, DRMKex > DDD_bigint_modul;
222 const bigint_matrix_algorithms< SRMKex, DRMKex, DRMKex > SDD_bigint_modul;
223 const bigint_matrix_algorithms< DRMKex, SRMKex, DRMKex > DSD_bigint_modul;
224 const bigint_matrix_algorithms< SRMKex, SRMKex, DRMKex > SSD_bigint_modul;
225
226 if (RES.get_no_of_rows() != M.rows)
227 RES.set_no_of_rows(M.rows);
228 if (RES.get_no_of_columns() != M.columns)
229 RES.set_no_of_columns(M.columns);
230
231 if (RES.bitfield.get_representation() == matrix_flags::dense_representation) {
232 if (RES.bitfield.get_representation() == matrix_flags::dense_representation)
233 DDD_bigint_modul.remainder(RES, M, mod);
234 else
235 SDD_bigint_modul.remainder(RES, M, mod);
236 }
237 else {
238 if (RES.bitfield.get_representation() == matrix_flags::dense_representation)
239 DSD_bigint_modul.remainder(RES, M, mod);
240 else
241 SSD_bigint_modul.remainder(RES, M, mod);
242 }
243 }
244
245
246
247 //
248 // pseudo-division
249 //
250
251 void matrix< bigint >::
divide(const matrix<bigint> & A,const bigint & k)252 divide(const matrix< bigint > &A, const bigint &k)
253 {
254 //
255 // Task: RES.divide(A, k)
256 // => RES.value[x][y] = A.value[x][y] / k,
257 // x=0,...,A.rows-1, y=0,...,A.columns-1
258 // Version: 2.0
259 //
260
261 debug_handler_l(DMESSAGE, "divide(const matrix< bigint > &, const T &)", DVALUE + 1);
262
263 const bigint_matrix_algorithms< DRMKex, SRMKex, DRMKex > DSD_bigint_modul;
264 const bigint_matrix_algorithms< SRMKex, DRMKex, SRMKex > SDD_bigint_modul;
265 const bigint_matrix_algorithms< SRMKex, SRMKex, DRMKex > SSD_bigint_modul;
266
267 if (rows != A.rows)
268 set_no_of_rows(A.rows);
269 if (columns != A.columns)
270 set_no_of_columns(A.columns);
271
272 if (bitfield.get_representation() == matrix_flags::dense_representation) {
273 if (A.bitfield.get_representation() == matrix_flags::dense_representation)
274 D_bigint_modul.divide(*this, A, k);
275 else
276 DSD_bigint_modul.divide(*this, A, k);
277 }
278 else {
279 if (A.bitfield.get_representation() == matrix_flags::dense_representation)
280 SDD_bigint_modul.divide(*this, A, k);
281 else
282 SSD_bigint_modul.divide(*this, A, k);
283 }
284 }
285
286
287
288 void matrix< bigint >::
compwise_divide(const matrix<bigint> & A,const matrix<bigint> & B)289 compwise_divide(const matrix< bigint > &A, const matrix< bigint > &B)
290 {
291 //
292 // Task: RES.compwise_divide(A, B)
293 // => RES.value[x][y] = A.value[x][y] / B.value[x][y],
294 // x=0,...,A.rows-1, y=0,...,A.columns-1
295 // Conditions: A.rows == A.rows and
296 // A.columns == B.columns
297 // Version: 2.0
298 //
299
300 debug_handler_l(DMESSAGE, "compwise_divide(const matrix< bigint > &, "
301 "const matrix< bigint > &)", DVALUE + 1);
302
303 if (A.rows != B.rows || A.columns != B.columns)
304 precondition_error_handler(A.rows, "A.rows", "A.rows == B.rows",
305 B.rows, "B.rows", "A.rows == B.rows",
306 A.columns, "A.columns", "A.columns == B.columns",
307 B.columns, "B.columns", "B.columns == B.columns",
308 "void matrix< bigint >::"
309 "compwise_divide(const matrix< bigint > &A, "
310 "const matrix< bigint > &B)",
311 DMESSAGE, EMESSAGE[4]);
312
313 const bigint_matrix_algorithms< DRMKex, DRMKex, SRMKex > DDS_bigint_modul;
314 const bigint_matrix_algorithms< DRMKex, SRMKex, DRMKex > DSD_bigint_modul;
315 const bigint_matrix_algorithms< SRMKex, DRMKex, SRMKex > SDD_bigint_modul;
316 const bigint_matrix_algorithms< DRMKex, SRMKex, SRMKex > DSS_bigint_modul;
317 const bigint_matrix_algorithms< SRMKex, DRMKex, SRMKex > SDS_bigint_modul;
318 const bigint_matrix_algorithms< SRMKex, SRMKex, DRMKex > SSD_bigint_modul;
319 const bigint_matrix_algorithms< SRMKex, SRMKex, SRMKex > SSS_bigint_modul;
320
321 if (rows != A.rows)
322 set_no_of_rows(A.rows);
323 if (columns != A.columns)
324 set_no_of_columns(A.columns);
325
326 if (bitfield.get_representation() == matrix_flags::dense_representation) {
327 if (A.bitfield.get_representation() ==
328 matrix_flags::dense_representation)
329 if (B.bitfield.get_representation() ==
330 matrix_flags::dense_representation)
331 D_bigint_modul.compwise_divide(*this, A, B);
332 else
333 DDS_bigint_modul.compwise_divide(*this, A, B);
334 else
335 if (B.bitfield.get_representation() ==
336 matrix_flags::dense_representation)
337 DSD_bigint_modul.compwise_divide(*this, A, B);
338 else
339 DSS_bigint_modul.compwise_divide(*this, A, B);
340 }
341 else {
342 if (A.bitfield.get_representation() ==
343 matrix_flags::dense_representation)
344 if (B.bitfield.get_representation() ==
345 matrix_flags::dense_representation)
346 SDD_bigint_modul.compwise_divide(*this, A, B);
347 else
348 SDS_bigint_modul.compwise_divide(*this, A, B);
349 else
350 if (B.bitfield.get_representation() ==
351 matrix_flags::dense_representation)
352 SSD_bigint_modul.compwise_divide(*this, A, B);
353 else
354 SSS_bigint_modul.compwise_divide(*this, A, B);
355 }
356 }
357
358
359
360 //
361 // norms and bounds
362 //
363
364 void matrix< bigint >::
max(bigint & MAX) const365 max(bigint &MAX) const
366 {
367 //
368 // Task: A.max(res);
369 // => res = maximum of all members of matrix A
370 // Version: 2.0
371 //
372
373 debug_handler_l(DMESSAGE, "max(bigint &)", DVALUE + 2);
374
375 if (bitfield.get_representation() == matrix_flags::dense_representation)
376 D_bigint_modul.max(*this, MAX);
377 else
378 S_bigint_modul.max(*this, MAX);
379 }
380
381
382
383 void matrix< bigint >::
max_abs(bigint & MAX) const384 max_abs(bigint &MAX) const
385 {
386 //
387 // Task: A.max_abs(res);
388 // => res = absolute maximum of all members of matrix A
389 // Version: 2.0
390 //
391
392 debug_handler_l(DMESSAGE, "max_abs(bigint &)", DVALUE + 2);
393
394 if (bitfield.get_representation() == matrix_flags::dense_representation)
395 D_bigint_modul.max_abs(*this, MAX);
396 else
397 S_bigint_modul.max_abs(*this, MAX);
398 }
399
400
401
402 void matrix< bigint >::
max_pos(bigint & MAX,lidia_size_t & x,lidia_size_t & y) const403 max_pos(bigint & MAX, lidia_size_t & x, lidia_size_t & y) const
404 {
405 //
406 // Task: A.max_pos(res, x, y);
407 // => res = maximum of all members of matrix A
408 // => res = A.value[x][y]
409 // Version: 2.0
410 //
411
412 debug_handler_l(DMESSAGE, "max_pos(bigint &, lidia_size_t &, lidia_size_t &)",
413 DVALUE + 2);
414
415 if (bitfield.get_representation() == matrix_flags::dense_representation)
416 D_bigint_modul.max_pos(*this, MAX, x, y);
417 else
418 S_bigint_modul.max_pos(*this, MAX, x, y);
419 }
420
421
422
423 void matrix< bigint >::
max_abs_pos(bigint & MAX,lidia_size_t & x,lidia_size_t & y) const424 max_abs_pos(bigint & MAX, lidia_size_t & x, lidia_size_t & y) const
425 {
426 //
427 // Task: A.max_abs_pos(res, x, y);
428 // => res = absolute maximum of all members of matrix A
429 // => res = A.value[x][y]
430 // Version: 2.0
431 //
432
433 debug_handler_l(DMESSAGE, "max_abs_pos(bigint &, lidia_size_t &, lidia_size_t &)",
434 DVALUE + 2);
435
436 if (bitfield.get_representation() == matrix_flags::dense_representation)
437 D_bigint_modul.max_abs_pos(*this, MAX, x, y);
438 else
439 S_bigint_modul.max_abs_pos(*this, MAX, x, y);
440 }
441
442
443
444 void matrix< bigint >::
min(bigint & MIN) const445 min(bigint &MIN) const
446 {
447 //
448 // Task: A.min(MIN);
449 // => MIN = minimum of all members of matrix A
450 // Version: 2.0
451 //
452
453 debug_handler_l(DMESSAGE, "min(bigint &)", DVALUE + 2);
454
455 if (bitfield.get_representation() == matrix_flags::dense_representation)
456 D_bigint_modul.min(*this, MIN);
457 else
458 S_bigint_modul.min(*this, MIN);
459 }
460
461
462
463 void matrix< bigint >::
min_abs(bigint & MIN) const464 min_abs(bigint &MIN) const
465 {
466 //
467 // Task: A.min_abs(MIN);
468 // => MIN = absolute minimum of all members of matrix A
469 // Version: 2.0
470 //
471
472 debug_handler_l(DMESSAGE, "min_abs(bigint &)", DVALUE + 2);
473
474 if (bitfield.get_representation() == matrix_flags::dense_representation)
475 D_bigint_modul.min_abs(*this, MIN);
476 else
477 S_bigint_modul.min_abs(*this, MIN);
478 }
479
480
481
482 void matrix< bigint >::
min_pos(bigint & MIN,lidia_size_t & x,lidia_size_t & y) const483 min_pos(bigint & MIN, lidia_size_t & x, lidia_size_t & y) const
484 {
485 //
486 // Task: A.min_pos(MIN, x, y);
487 // => MIN = minimum of all members of matrix A
488 // => MIN = A.value[x][y]
489 // Version: 2.0
490 //
491
492 debug_handler_l(DMESSAGE, "min_pos(bigint &, lidia_size_t &, lidia_size_t &)",
493 DVALUE + 2);
494
495 if (bitfield.get_representation() == matrix_flags::dense_representation)
496 D_bigint_modul.min_pos(*this, MIN, x, y);
497 else
498 S_bigint_modul.min_pos(*this, MIN, x, y);
499 }
500
501
502
503 void matrix< bigint >::
min_abs_pos(bigint & MIN,lidia_size_t & x,lidia_size_t & y) const504 min_abs_pos(bigint & MIN, lidia_size_t & x, lidia_size_t & y) const
505 {
506 //
507 // Task: A.min_abs_pos(MIN, x, y);
508 // => MIN = absolute minimum of all members of matrix A
509 // => MIN = A.value[x][y]
510 // Version: 2.0
511 //
512
513 debug_handler_l(DMESSAGE, "min_abs_pos(bigint &, lidia_size_t &, lidia_size_t &)",
514 DVALUE + 2);
515
516 if (bitfield.get_representation() == matrix_flags::dense_representation)
517 D_bigint_modul.min_abs_pos(*this, MIN, x, y);
518 else
519 S_bigint_modul.min_abs_pos(*this, MIN, x, y);
520
521 }
522
523
524
525 void matrix< bigint >::
hadamard(bigint & H) const526 hadamard(bigint & H) const
527 {
528 //
529 // Task: A.hadamard(H);
530 // => H = hadamard bound of matrix A
531 // Version: 2.0
532 //
533
534 debug_handler_l(DMESSAGE, "hadamard(bigint &)", DVALUE + 2);
535
536 if (bitfield.get_representation() == matrix_flags::dense_representation)
537 D_bigint_modul.hadamard(*this, H);
538 else
539 S_bigint_modul.hadamard(*this, H);
540 }
541
542
543
544 void matrix< bigint >::
binary_hadamard(lidia_size_t & H) const545 binary_hadamard(lidia_size_t &H) const
546 {
547 //
548 // Task: A.hadamard2(H);
549 // => H = hadamard bound of matrix A
550 // Version: 2.0
551 //
552
553 debug_handler_l(DMESSAGE, "binary_hadamard(lidia_size_t &)", DVALUE + 2);
554
555 if (bitfield.get_representation() == matrix_flags::dense_representation)
556 D_bigint_modul.binary_hadamard(*this, H);
557 else
558 S_bigint_modul.binary_hadamard(*this, H);
559 }
560
561
562
563 void matrix< bigint >::
row_norm(bigint & RES,lidia_size_t pos,long art) const564 row_norm(bigint & RES, lidia_size_t pos, long art) const
565 {
566 //
567 // Task: A.row_norm(RES, pos, art);
568 // => RES = A(pos, 0)^art + ..... + A(pos, A.columns - 1)^art
569 // Version: 2.0
570 //
571
572 debug_handler_l(DMESSAGE, "row_norm(bigint &, lidia_size_t, long)",
573 DVALUE + 2);
574
575 if (pos< 0 || pos >= rows)
576 precondition_error_handler(pos, "pos", "0 <= pos < rows",
577 rows, "rows", "",
578 "void matrix< bigint >::"
579 "row_norm(bigint & RES, lidia_size_t pos, long art) const",
580 DMESSAGE, EMESSAGE[1]);
581
582 if (bitfield.get_representation() == matrix_flags::dense_representation)
583 D_bigint_modul.row_norm(*this, RES, pos, art);
584 else
585 S_bigint_modul.row_norm(*this, RES, pos, art);
586 }
587
588
589
590 void matrix< bigint >::
column_norm(bigint & RES,lidia_size_t pos,long art) const591 column_norm(bigint & RES, lidia_size_t pos, long art) const
592 {
593 //
594 // Task: A.column_norm(RES, pos , art);
595 // => RES = A(0, pos)^art + ..... + A(A.rows - 1, pos)^art
596 // Version: 2.0
597 //
598
599 debug_handler_l(DMESSAGE, "column_norm(bigint &, lidia_size_t, long)",
600 DVALUE + 2);
601
602 if (pos< 0 || pos >= columns)
603 precondition_error_handler(pos, "pos", "0 <= pos < columns",
604 columns, "columns", "",
605 "void matrix< bigint >::"
606 "column_norm(bigint & RES, lidia_size_t pos, long art) const",
607 DMESSAGE, EMESSAGE[1]);
608
609 if (bitfield.get_representation() == matrix_flags::dense_representation)
610 D_bigint_modul.column_norm(*this, RES, pos, art);
611 else
612 S_bigint_modul.column_norm(*this, RES, pos, art);
613 }
614
615
616
617 //
618 // randomize
619 //
620
621 void matrix< bigint >::
randomize(const bigint & S)622 randomize(const bigint & S)
623 {
624 //
625 // Task: RES.randomize(S);
626 // => 0 <= RES.value[i][j] <= S, i=0,...,RES.rows-1,
627 // j=0,...,RES.columns-1; random values
628 // Version: 2.0
629 //
630
631 debug_handler_l(DMESSAGE, "randomize(const bigint &)", DVALUE + 3);
632
633 if (bitfield.get_representation() == matrix_flags::dense_representation)
634 D_bigint_modul.randomize(*this, S);
635 else
636 S_bigint_modul.randomize(*this, S);
637 }
638
639
640
641 void matrix< bigint >::
randomize_with_det(const bigint & S,const bigint & DET)642 randomize_with_det(const bigint & S, const bigint & DET)
643 {
644 //
645 // Task: RES.randomize_with_det(S, DET);
646 // => 0 <= RES.value[i][j] <= S, i=0,...,RES.rows-1,
647 // j=0,...,RES.columns-1; random values
648 // => det(RES) == DET
649 // Version: 2.0
650 //
651
652 debug_handler_l(DMESSAGE, "randomize_with_det(const bigint &, const bigint &)",
653 DVALUE + 3);
654
655 if (rows != columns)
656 precondition_error_handler(rows, "rows", "rows == columns",
657 columns, "columns", "rows == columns",
658 "void matrix< bigint >::"
659 "randomize_with_det(const bigint & S, const bigint & DET)",
660 DMESSAGE, EMESSAGE[7]);
661
662 if (bitfield.get_representation() == matrix_flags::dense_representation)
663 D_bigint_modul.randomize_with_det(*this, S, DET);
664 else
665 S_bigint_modul.randomize_with_det(*this, S, DET);
666 }
667
668
669
670 void matrix< bigint >::
randomize(const bigint & S,const long d)671 randomize(const bigint & S, const long d)
672 {
673 //
674 // Task: RES.randomize(S, d);
675 // => 0 <= RES.value[i][j] <= S, i=0,...,RES.rows-1,
676 // j=0,...,RES.columns-1; random values (d %)
677 // Version: 2.0
678 //
679
680 debug_handler_l(DMESSAGE, "randomize(const bigint &, const long)", DVALUE + 3);
681
682 if (bitfield.get_representation() == matrix_flags::dense_representation)
683 D_bigint_modul.randomize(*this, S, d);
684 else
685 S_bigint_modul.randomize(*this, S, d);
686 }
687
688
689
690 ///////////////////////////
691 // BEGIN: Linear algebra //
692 // PART 1 //
693 ///////////////////////////
694
695 //
696 // rank
697 //
698
699 lidia_size_t matrix< bigint >::
rank(const bigint & H) const700 rank(const bigint &H) const
701 {
702 //
703 // Task: A.rank(H) = Rank of matrix A with
704 // H = hadamard(A).
705 // Version: 2.0
706 //
707
708 debug_handler_l(DMESSAGE, "rank()", DVALUE + 4);
709
710 const modular_arithmetic< DRMK< bigint >,
711 dense_fp_matrix_kernel< long, MR< long > >,
712 dense_fp_matrix_kernel< bigint, MR< bigint > > > Dm_bigint_modul;
713
714 const modular_arithmetic< SRMK< bigint >,
715 sparse_fp_matrix_kernel< long, MR< long > >,
716 sparse_fp_matrix_kernel< bigint, MR< bigint > > > Sm_bigint_modul;
717
718 if (bitfield.get_representation() == matrix_flags::dense_representation)
719 return Dm_bigint_modul.rank(*this, H);
720 else
721 return Sm_bigint_modul.rank(*this, H);
722 }
723
724
725
726 lidia_size_t matrix< bigint >::
rank() const727 rank() const
728 {
729 //
730 // Task: A.rank() = Rank of matrix A.
731 // Version: 2.0
732 //
733
734 debug_handler_l(DMESSAGE, "rank()", DVALUE + 4);
735
736 bigint H;
737
738 const modular_arithmetic< DRMK< bigint >,
739 dense_fp_matrix_kernel< long, MR< long > >,
740 dense_fp_matrix_kernel< bigint, MR< bigint > > > Dm_bigint_modul;
741
742 const modular_arithmetic< SRMK< bigint >,
743 sparse_fp_matrix_kernel< long, MR< long > >,
744 sparse_fp_matrix_kernel< bigint, MR< bigint > > > Sm_bigint_modul;
745
746 if (bitfield.get_representation() == matrix_flags::dense_representation) {
747 D_bigint_modul.hadamard(*this, H);
748 return Dm_bigint_modul.rank(*this, H);
749 }
750 else {
751 S_bigint_modul.hadamard(*this, H);
752 return Sm_bigint_modul.rank(*this, H);
753 }
754 }
755
756
757
758 //
759 // rank and linearly independent rows
760 //
761
762 lidia_size_t *matrix< bigint >::
lininr1(const bigint & H) const763 lininr1(const bigint &H) const
764 {
765 //
766 // Task: RES[0] = Rank of matrix (Avalue,r,c).
767 // RES[1],...,RES[RES[0]], such that
768 // row(RES[1]),...,row(RES[RES[0]])
769 // of matrix *this are linearly independent.
770 // Version: 2.0
771 //
772
773 debug_handler_l(DMESSAGE, "lininr1(const bigint &)", DVALUE + 4);
774
775 const modular_arithmetic< DRMK< bigint >,
776 dense_fp_matrix_kernel< long, MR< long > >,
777 dense_fp_matrix_kernel< bigint, MR< bigint > > > Dm_bigint_modul;
778
779 const modular_arithmetic< SRMK< bigint >,
780 sparse_fp_matrix_kernel< long, MR< long > >,
781 sparse_fp_matrix_kernel< bigint, MR< bigint > > > Sm_bigint_modul;
782
783 if (bitfield.get_representation() == matrix_flags::dense_representation)
784 return Dm_bigint_modul.lininr1(*this, H);
785 else
786 return Sm_bigint_modul.lininr1(*this, H);
787 }
788
789
790
791 lidia_size_t *matrix< bigint >::
lininr1() const792 lininr1() const
793 {
794 //
795 // Task: RES[0] = Rank of matrix (Avalue,r,c).
796 // RES[1],...,RES[RES[0]], such that
797 // row(RES[1]),...,row(RES[RES[0]])
798 // of matrix *this are linearly independent.
799 // Version: 2.0
800 //
801
802 debug_handler_l(DMESSAGE, "lininr1()", DVALUE + 4);
803
804 bigint H;
805
806 const modular_arithmetic< DRMK< bigint >,
807 dense_fp_matrix_kernel< long, MR< long > >,
808 dense_fp_matrix_kernel< bigint, MR< bigint > > > Dm_bigint_modul;
809
810 const modular_arithmetic< SRMK< bigint >,
811 sparse_fp_matrix_kernel< long, MR< long > >,
812 sparse_fp_matrix_kernel< bigint, MR< bigint > > > Sm_bigint_modul;
813
814 if (bitfield.get_representation() == matrix_flags::dense_representation) {
815 D_bigint_modul.hadamard(*this, H);
816 return Dm_bigint_modul.lininr1(*this, H);
817 }
818 else {
819 S_bigint_modul.hadamard(*this, H);
820 return Sm_bigint_modul.lininr1(*this, H);
821 }
822 }
823
824
825
826 void matrix< bigint >::
lininr1(base_vector<lidia_size_t> & RES) const827 lininr1(base_vector< lidia_size_t > &RES) const
828 {
829 //
830 // Task: RES[0] = Rank of matrix (Avalue,r,c).
831 // RES[1],...,RES[RES[0]], such that
832 // row(RES[1]),...,row(RES[RES[0]])
833 // of matrix *this are linearly independent.
834 // Version: 2.0
835 //
836
837 debug_handler_l(DMESSAGE, "lininr1(base_vector< lidia_size_t > &)", DVALUE + 4);
838
839 bigint H;
840 lidia_size_t *tmp;
841
842 const modular_arithmetic< DRMK< bigint >,
843 dense_fp_matrix_kernel< long, MR< long > >,
844 dense_fp_matrix_kernel< bigint, MR< bigint > > > Dm_bigint_modul;
845
846 const modular_arithmetic< SRMK< bigint >,
847 sparse_fp_matrix_kernel< long, MR< long > >,
848 sparse_fp_matrix_kernel< bigint, MR< bigint > > > Sm_bigint_modul;
849
850 if (bitfield.get_representation() == matrix_flags::dense_representation) {
851 D_bigint_modul.hadamard(*this, H);
852 tmp = Dm_bigint_modul.lininr1(*this, H);
853 }
854 else {
855 S_bigint_modul.hadamard(*this, H);
856 tmp = Sm_bigint_modul.lininr1(*this, H);
857 }
858
859 if (RES.capacity() < tmp[0])
860 RES.set_capacity(tmp[0]);
861 if (RES.size() != tmp[0])
862 RES.set_size(tmp[0]);
863
864 for (register lidia_size_t i = 0; i <= tmp[0]; i++)
865 RES[i] = tmp[i];
866
867 delete[] tmp;
868 }
869
870
871
872 lidia_size_t *matrix< bigint >::
lininr2(const bigint & H) const873 lininr2(const bigint &H) const
874 {
875 //
876 // Task: RES[0] = Rank of matrix (Avalue,r,c).
877 // RES[1],...,RES[RES[0]], such that
878 // row(RES[1]),...,row(RES[RES[0]])
879 // of matrix *this are linearly independent.
880 // Version: 2.0
881 //
882
883 debug_handler_l(DMESSAGE, "lininr2(const bigint &)", DVALUE + 4);
884
885 const modular_arithmetic< DRMK< bigint >,
886 dense_fp_matrix_kernel< long, MR< long > >,
887 dense_fp_matrix_kernel< bigint, MR< bigint > > > Dm_bigint_modul;
888
889 const modular_arithmetic< SRMK< bigint >,
890 sparse_fp_matrix_kernel< long, MR< long > >,
891 sparse_fp_matrix_kernel< bigint, MR< bigint > > > Sm_bigint_modul;
892
893 if (bitfield.get_representation() == matrix_flags::dense_representation)
894 return Dm_bigint_modul.lininr2(*this, H);
895 else
896 return Sm_bigint_modul.lininr2(*this, H);
897 }
898
899
900
901 lidia_size_t *matrix< bigint >::
lininr2() const902 lininr2() const
903 {
904 //
905 // Task: RES[0] = Rank of matrix (Avalue,r,c).
906 // RES[1],...,RES[RES[0]], such that
907 // row(RES[1]),...,row(RES[RES[0]])
908 // of matrix *this are linearly independent.
909 // Version: 2.0
910 //
911
912 debug_handler_l(DMESSAGE, "lininr2()", DVALUE + 4);
913
914 bigint H;
915
916 const modular_arithmetic< DRMK< bigint >,
917 dense_fp_matrix_kernel< long, MR< long > >,
918 dense_fp_matrix_kernel< bigint, MR< bigint > > > Dm_bigint_modul;
919
920 const modular_arithmetic< SRMK< bigint >,
921 sparse_fp_matrix_kernel< long, MR< long > >,
922 sparse_fp_matrix_kernel< bigint, MR< bigint > > > Sm_bigint_modul;
923
924 if (bitfield.get_representation() == matrix_flags::dense_representation) {
925 D_bigint_modul.hadamard(*this, H);
926 return Dm_bigint_modul.lininr2(*this, H);
927 }
928 else {
929 S_bigint_modul.hadamard(*this, H);
930 return Sm_bigint_modul.lininr2(*this, H);
931 }
932 }
933
934
935
936 void matrix< bigint >::
lininr2(base_vector<lidia_size_t> & RES) const937 lininr2(base_vector< lidia_size_t > &RES) const
938 {
939 //
940 // Task: RES[0] = Rank of matrix (Avalue,r,c).
941 // RES[1],...,RES[RES[0]], such that
942 // row(RES[1]),...,row(RES[RES[0]])
943 // of matrix *this are linearly independent.
944 // Version: 2.0
945 //
946
947 debug_handler_l(DMESSAGE, "lininr2(base_vector< lidia_size_t > &)", DVALUE + 4);
948
949 bigint H;
950 lidia_size_t *tmp;
951
952 const modular_arithmetic< DRMK< bigint >,
953 dense_fp_matrix_kernel< long, MR< long > >,
954 dense_fp_matrix_kernel< bigint, MR< bigint > > > Dm_bigint_modul;
955
956 const modular_arithmetic< SRMK< bigint >,
957 sparse_fp_matrix_kernel< long, MR< long > >,
958 sparse_fp_matrix_kernel< bigint, MR< bigint > > > Sm_bigint_modul;
959
960 if (bitfield.get_representation() == matrix_flags::dense_representation) {
961 D_bigint_modul.hadamard(*this, H);
962 tmp = Dm_bigint_modul.lininr2(*this, H);
963 }
964 else {
965 S_bigint_modul.hadamard(*this, H);
966 tmp = Sm_bigint_modul.lininr2(*this, H);
967 }
968
969 if (RES.capacity() < tmp[0])
970 RES.set_capacity(tmp[0]);
971 if (RES.size() != tmp[0])
972 RES.set_size(tmp[0]);
973
974 for (register lidia_size_t i = 0; i <= tmp[0]; i++)
975 RES[i] = tmp[i];
976
977 delete[] tmp;
978 }
979
980
981
982 //
983 // rank linearly independent columns
984 //
985
986 lidia_size_t *matrix< bigint >::
lininc1(const bigint & H) const987 lininc1(const bigint &H) const
988 {
989 //
990 // Task: RES[0] = Rank of matrix (Avalue,r,c).
991 // RES[1],...,RES[RES[0]], such that
992 // column(RES[1]),...,column(RES[RES[0]])
993 // of matrix *this are linearly independent.
994 // Version: 2.0
995 //
996
997 debug_handler_l(DMESSAGE, "lininc1(const bigint &)", DVALUE + 4);
998
999 const modular_arithmetic< DRMK< bigint >,
1000 dense_fp_matrix_kernel< long, MR< long > >,
1001 dense_fp_matrix_kernel< bigint, MR< bigint > > > Dm_bigint_modul;
1002
1003 const modular_arithmetic< SRMK< bigint >,
1004 sparse_fp_matrix_kernel< long, MR< long > >,
1005 sparse_fp_matrix_kernel< bigint, MR< bigint > > > Sm_bigint_modul;
1006
1007 if (bitfield.get_representation() == matrix_flags::dense_representation)
1008 return Dm_bigint_modul.lininc1(*this, H);
1009 else
1010 return Sm_bigint_modul.lininc1(*this, H);
1011 }
1012
1013
1014
1015 lidia_size_t *matrix< bigint >::
lininc1() const1016 lininc1() const
1017 {
1018 //
1019 // Task: RES[0] = Rank of matrix (Avalue,r,c).
1020 // RES[1],...,RES[RES[0]], such that
1021 // column(RES[1]),...,column(RES[RES[0]])
1022 // of matrix *this are linearly independent.
1023 // Version: 2.0
1024 //
1025
1026 debug_handler_l(DMESSAGE, "lininc1()", DVALUE + 4);
1027
1028 bigint H;
1029
1030 const modular_arithmetic< DRMK< bigint >,
1031 dense_fp_matrix_kernel< long, MR< long > >,
1032 dense_fp_matrix_kernel< bigint, MR< bigint > > > Dm_bigint_modul;
1033
1034 const modular_arithmetic< SRMK< bigint >,
1035 sparse_fp_matrix_kernel< long, MR< long > >,
1036 sparse_fp_matrix_kernel< bigint, MR< bigint > > > Sm_bigint_modul;
1037
1038 if (bitfield.get_representation() == matrix_flags::dense_representation) {
1039 D_bigint_modul.hadamard(*this, H);
1040 return Dm_bigint_modul.lininc1(*this, H);
1041 }
1042 else {
1043 S_bigint_modul.hadamard(*this, H);
1044 return Sm_bigint_modul.lininc1(*this, H);
1045 }
1046 }
1047
1048
1049
1050 void matrix< bigint >::
lininc1(base_vector<lidia_size_t> & RES) const1051 lininc1(base_vector< lidia_size_t > &RES) const
1052 {
1053 //
1054 // Task: RES[0] = Rank of matrix (Avalue,r,c).
1055 // RES[1],...,RES[RES[0]], such that
1056 // column(RES[1]),...,column(RES[RES[0]])
1057 // of matrix *this are linearly independent.
1058 // Version: 2.0
1059 //
1060
1061 debug_handler_l(DMESSAGE, "lininc1(base_vector< lidia_size_t > &)",
1062 DVALUE + 4);
1063
1064 bigint H;
1065 lidia_size_t *tmp;
1066
1067 const modular_arithmetic< DRMK< bigint >,
1068 dense_fp_matrix_kernel< long, MR< long > >,
1069 dense_fp_matrix_kernel< bigint, MR< bigint > > > Dm_bigint_modul;
1070
1071 const modular_arithmetic< SRMK< bigint >,
1072 sparse_fp_matrix_kernel< long, MR< long > >,
1073 sparse_fp_matrix_kernel< bigint, MR< bigint > > > Sm_bigint_modul;
1074
1075 if (bitfield.get_representation() == matrix_flags::dense_representation) {
1076 D_bigint_modul.hadamard(*this, H);
1077 tmp = Dm_bigint_modul.lininc1(*this, H);
1078 }
1079 else {
1080 S_bigint_modul.hadamard(*this, H);
1081 tmp = Sm_bigint_modul.lininc1(*this, H);
1082 }
1083
1084 if (RES.capacity() < tmp[0])
1085 RES.set_capacity(tmp[0]);
1086 if (RES.size() != tmp[0])
1087 RES.set_size(tmp[0]);
1088
1089 for (register lidia_size_t i = 0; i <= tmp[0]; i++)
1090 RES[i] = tmp[i];
1091
1092 delete[] tmp;
1093 }
1094
1095
1096
1097 lidia_size_t *matrix< bigint >::
lininc2(const bigint & H) const1098 lininc2(const bigint &H) const
1099 {
1100 //
1101 // Task: RES[0] = Rank of matrix (Avalue,r,c).
1102 // RES[1],...,RES[RES[0]], such that
1103 // column(RES[1]),...,column(RES[RES[0]])
1104 // of matrix *this are linearly independent.
1105 // Version: 2.0
1106 //
1107
1108 debug_handler_l(DMESSAGE, "lininc2(const bigint &)", DVALUE + 4);
1109
1110 const modular_arithmetic< DRMK< bigint >,
1111 dense_fp_matrix_kernel< long, MR< long > >,
1112 dense_fp_matrix_kernel< bigint, MR< bigint > > > Dm_bigint_modul;
1113
1114 const modular_arithmetic< SRMK< bigint >,
1115 sparse_fp_matrix_kernel< long, MR< long > >,
1116 sparse_fp_matrix_kernel< bigint, MR< bigint > > > Sm_bigint_modul;
1117
1118 if (bitfield.get_representation() == matrix_flags::dense_representation)
1119 return Dm_bigint_modul.lininc2(*this, H);
1120 else
1121 return Sm_bigint_modul.lininc2(*this, H);
1122 }
1123
1124
1125
1126 lidia_size_t *matrix< bigint >::
lininc2() const1127 lininc2() const
1128 {
1129 //
1130 // Task: RES[0] = Rank of matrix (Avalue,r,c).
1131 // RES[1],...,RES[RES[0]], such that
1132 // column(RES[1]),...,column(RES[RES[0]])
1133 // of matrix *this are linearly independent.
1134 // Version: 2.0
1135 //
1136
1137 debug_handler_l(DMESSAGE, "lininc2()", DVALUE + 4);
1138
1139 bigint H;
1140
1141 const modular_arithmetic< DRMK< bigint >,
1142 dense_fp_matrix_kernel< long, MR< long > >,
1143 dense_fp_matrix_kernel< bigint, MR< bigint > > > Dm_bigint_modul;
1144
1145 const modular_arithmetic< SRMK< bigint >,
1146 sparse_fp_matrix_kernel< long, MR< long > >,
1147 sparse_fp_matrix_kernel< bigint, MR< bigint > > > Sm_bigint_modul;
1148
1149 if (bitfield.get_representation() == matrix_flags::dense_representation) {
1150 D_bigint_modul.hadamard(*this, H);
1151 return Dm_bigint_modul.lininc2(*this, H);
1152 }
1153 else {
1154 S_bigint_modul.hadamard(*this, H);
1155 return Sm_bigint_modul.lininc2(*this, H);
1156 }
1157 }
1158
1159
1160
1161 void matrix< bigint >::
lininc2(base_vector<lidia_size_t> & RES) const1162 lininc2(base_vector< lidia_size_t > &RES) const
1163 {
1164 //
1165 // Task: RES[0] = Rank of matrix (Avalue,r,c).
1166 // RES[1],...,RES[RES[0]], such that
1167 // column(RES[1]),...,column(RES[RES[0]])
1168 // of matrix *this are linearly independent.
1169 // Version: 2.0
1170 //
1171
1172 debug_handler_l(DMESSAGE, "lininc2(base_vector< lidia_size_t > &)",
1173 DVALUE + 4);
1174
1175 bigint H;
1176 lidia_size_t *tmp;
1177
1178 const modular_arithmetic< DRMK< bigint >,
1179 dense_fp_matrix_kernel< long, MR< long > >,
1180 dense_fp_matrix_kernel< bigint, MR< bigint > > > Dm_bigint_modul;
1181
1182 const modular_arithmetic< SRMK< bigint >,
1183 sparse_fp_matrix_kernel< long, MR< long > >,
1184 sparse_fp_matrix_kernel< bigint, MR< bigint > > > Sm_bigint_modul;
1185
1186 if (bitfield.get_representation() == matrix_flags::dense_representation) {
1187 D_bigint_modul.hadamard(*this, H);
1188 tmp = Dm_bigint_modul.lininc2(*this, H);
1189 }
1190 else {
1191 S_bigint_modul.hadamard(*this, H);
1192 tmp = Sm_bigint_modul.lininc2(*this, H);
1193 }
1194
1195 if (RES.capacity() < tmp[0])
1196 RES.set_capacity(tmp[0]);
1197 if (RES.size() != tmp[0])
1198 RES.set_size(tmp[0]);
1199
1200 for (register lidia_size_t i = 0; i <= tmp[0]; i++)
1201 RES[i] = tmp[i];
1202
1203 delete[] tmp;
1204 }
1205
1206
1207
1208 //
1209 // regular expansion
1210 //
1211
1212 void matrix< bigint >::
regexpansion(const lidia_size_t * v)1213 regexpansion(const lidia_size_t * v)
1214 {
1215 //
1216 // Task: A.regexpansion(v);
1217 // => A = Regular Expansion of the old matrix A relative v.
1218 // Version: 2.0
1219 //
1220
1221 debug_handler_l(DMESSAGE, "regexpansion(const lidia_size_t *)", DVALUE + 4);
1222
1223 if (bitfield.get_representation() == matrix_flags::dense_representation)
1224 D_bigint_modul.regexpansion(*this, v);
1225 else
1226 S_bigint_modul.regexpansion(*this, v);
1227 }
1228
1229
1230
1231 //
1232 // adjoint matrix
1233 //
1234
1235 void matrix< bigint >::
adj1(const matrix<bigint> & A,const bigint & H,const bigint & DET)1236 adj1(const matrix< bigint > & A, const bigint &H, const bigint &DET)
1237 {
1238 //
1239 // Task: B.adj(A);
1240 // => B = adjoint matrix of matrix A
1241 // Conditions: A.columns != A.rows
1242 // Version: 2.0
1243 //
1244
1245 debug_handler_l(DMESSAGE, "adj1(const matrix< bigint > &, const bigint &, const bigint &)",
1246 DVALUE + 4);
1247
1248 if (A.columns != A.rows || DET == 0)
1249 precondition_error_handler(A.columns, "A.columns", "A.columns == A.rows",
1250 A.rows, "A.rows", "A.columns == A.rows",
1251 DET, "DET", "DET != 0",
1252 "void matrix< bigint >::adj1(const matrix< bigint > & A, "
1253 "const bigint &H, const bigint &DET)", DMESSAGE, EMESSAGE[7]);
1254
1255 const modular_arithmetic< DRMK< bigint >,
1256 dense_fp_matrix_kernel< long, MR< long > >,
1257 dense_fp_matrix_kernel< bigint, MR< bigint > > > Dm_bigint_modul;
1258
1259 const modular_arithmetic< SRMK< bigint >,
1260 sparse_fp_matrix_kernel< long, MR< long > >,
1261 sparse_fp_matrix_kernel< bigint, MR< bigint > > > Sm_bigint_modul;
1262
1263 if (bitfield.get_representation() == matrix_flags::dense_representation)
1264 Dm_bigint_modul.adj1(*this, A, H, DET);
1265 else
1266 Sm_bigint_modul.adj1(*this, A, H, DET);
1267 }
1268
1269
1270
1271 void matrix< bigint >::
adj1(const matrix<bigint> & A,const bigint & H)1272 adj1(const matrix< bigint > & A, const bigint &H)
1273 {
1274 //
1275 // Task: B.adj(A);
1276 // => B = adjoint matrix of matrix A
1277 // Conditions: A.columns != A.rows
1278 // Version: 2.0
1279 //
1280
1281 debug_handler_l(DMESSAGE, "adj1(const matrix< bigint > &, const bigint &H)", DVALUE + 4);
1282
1283 if (A.columns != A.rows)
1284 precondition_error_handler(A.columns, "A.columns", "A.columns == A.rows",
1285 A.rows, "A.rows", "A.columns == A.rows",
1286 "void matrix< bigint >::"
1287 "adj1(const matrix< bigint > & A, const bigint &H)", DMESSAGE, EMESSAGE[7]);
1288
1289 const modular_arithmetic< DRMK< bigint >,
1290 dense_fp_matrix_kernel< long, MR< long > >,
1291 dense_fp_matrix_kernel< bigint, MR< bigint > > > Dm_bigint_modul;
1292
1293 const modular_arithmetic< SRMK< bigint >,
1294 sparse_fp_matrix_kernel< long, MR< long > >,
1295 sparse_fp_matrix_kernel< bigint, MR< bigint > > > Sm_bigint_modul;
1296
1297 bigint DET;
1298
1299 if (bitfield.get_representation() == matrix_flags::dense_representation) {
1300 Dm_bigint_modul.det(A, DET, H);
1301 if (DET == 0)
1302 precondition_error_handler(DET, "det(A)", "det(A) != 0", "void matrix< bigint >::"
1303 "adj1(const matrix< bigint > & A, const bigint &H)",
1304 DMESSAGE, EMESSAGE[7]);
1305 Dm_bigint_modul.adj1(*this, A, H, DET);
1306 }
1307 else {
1308 Sm_bigint_modul.det(A, DET, H);
1309 if (DET == 0)
1310 precondition_error_handler(DET, "det(A)", "det(A) != 0", "void matrix< bigint >::"
1311 "adj1(const matrix< bigint > & A, const bigint &H)",
1312 DMESSAGE, EMESSAGE[7]);
1313 Sm_bigint_modul.adj1(*this, A, H, DET);
1314 }
1315 }
1316
1317
1318
1319 void matrix< bigint >::
adj1(const matrix<bigint> & A)1320 adj1(const matrix< bigint > & A)
1321 {
1322 //
1323 // Task: B.adj(A);
1324 // => B = adjoint matrix of matrix A
1325 // Conditions: A.columns != A.rows
1326 // Version: 2.0
1327 //
1328
1329 debug_handler_l(DMESSAGE, "adj1(const matrix< bigint > &)", DVALUE + 4);
1330
1331 if (A.columns != A.rows)
1332 precondition_error_handler(A.columns, "A.columns", "A.columns == A.rows",
1333 A.rows, "A.rows", "A.columns == A.rows",
1334 "void matrix< bigint >::"
1335 "adj1(const matrix< bigint > & A)", DMESSAGE, EMESSAGE[7]);
1336
1337 bigint H, DET;
1338
1339 const modular_arithmetic< DRMK< bigint >,
1340 dense_fp_matrix_kernel< long, MR< long > >,
1341 dense_fp_matrix_kernel< bigint, MR< bigint > > > Dm_bigint_modul;
1342
1343 const modular_arithmetic< SRMK< bigint >,
1344 sparse_fp_matrix_kernel< long, MR< long > >,
1345 sparse_fp_matrix_kernel< bigint, MR< bigint > > > Sm_bigint_modul;
1346
1347 if (bitfield.get_representation() == matrix_flags::dense_representation) {
1348 D_bigint_modul.hadamard(A, H);
1349 Dm_bigint_modul.det(A, DET, H);
1350 if (DET == 0)
1351 precondition_error_handler(DET, "det(A)", "det(A) != 0", "void matrix< bigint >::"
1352 "adj1(const matrix< bigint > & A)", DMESSAGE, EMESSAGE[7]);
1353 Dm_bigint_modul.adj1(*this, A, H, DET);
1354 }
1355 else {
1356 S_bigint_modul.hadamard(A, H);
1357 Sm_bigint_modul.det(A, DET, H);
1358 if (DET == 0)
1359 precondition_error_handler(DET, "det(A)", "det(A) != 0", "void matrix< bigint >::"
1360 "adj1(const matrix< bigint > & A)", DMESSAGE, EMESSAGE[7]);
1361 Sm_bigint_modul.adj1(*this, A, H, DET);
1362 }
1363 }
1364
1365
1366
1367 void matrix< bigint >::
adj2(const matrix<bigint> & A,const bigint & H,const bigint & DET)1368 adj2(const matrix< bigint > & A, const bigint &H, const bigint &DET)
1369 {
1370 //
1371 // Task: B.adj(A);
1372 // => B = adjoint matrix of matrix A
1373 // Conditions: A.columns != A.rows
1374 // Version: 2.0
1375 //
1376
1377 debug_handler_l(DMESSAGE, "adj2(const matrix< bigint > &, const bigint &, const bigint &)",
1378 DVALUE + 4);
1379
1380 if (A.columns != A.rows || DET == 0)
1381 precondition_error_handler(A.columns, "A.columns", "A.columns == A.rows",
1382 A.rows, "A.rows", "A.columns == A.rows",
1383 DET, "DET", "DET != 0",
1384 "void matrix< bigint >::adj2(const matrix< bigint > & A, "
1385 "const bigint &H, const bigint &DET)", DMESSAGE, EMESSAGE[7]);
1386
1387 const modular_arithmetic< DRMK< bigint >,
1388 dense_fp_matrix_kernel< long, MR< long > >,
1389 dense_fp_matrix_kernel< bigint, MR< bigint > > > Dm_bigint_modul;
1390
1391 const modular_arithmetic< SRMK< bigint >,
1392 sparse_fp_matrix_kernel< long, MR< long > >,
1393 sparse_fp_matrix_kernel< bigint, MR< bigint > > > Sm_bigint_modul;
1394
1395 if (bitfield.get_representation() == matrix_flags::dense_representation)
1396 Dm_bigint_modul.adj2(*this, A, H, DET);
1397 else
1398 Sm_bigint_modul.adj2(*this, A, H, DET);
1399 }
1400
1401
1402
1403 void matrix< bigint >::
adj2(const matrix<bigint> & A,const bigint & H)1404 adj2(const matrix< bigint > & A, const bigint &H)
1405 {
1406 //
1407 // Task: B.adj(A);
1408 // => B = adjoint matrix of matrix A
1409 // Conditions: A.columns != A.rows
1410 // Version: 2.0
1411 //
1412
1413 debug_handler_l(DMESSAGE, "adj2(const matrix< bigint > &, const bigint &H)", DVALUE + 4);
1414
1415 if (A.columns != A.rows)
1416 precondition_error_handler(A.columns, "A.columns", "A.columns == A.rows",
1417 A.rows, "A.rows", "A.columns == A.rows",
1418 "void matrix< bigint >::"
1419 "adj2(const matrix< bigint > & A, const bigint &H)", DMESSAGE, EMESSAGE[7]);
1420
1421 const modular_arithmetic< DRMK< bigint >,
1422 dense_fp_matrix_kernel< long, MR< long > >,
1423 dense_fp_matrix_kernel< bigint, MR< bigint > > > Dm_bigint_modul;
1424
1425 const modular_arithmetic< SRMK< bigint >,
1426 sparse_fp_matrix_kernel< long, MR< long > >,
1427 sparse_fp_matrix_kernel< bigint, MR< bigint > > > Sm_bigint_modul;
1428
1429 bigint DET;
1430
1431 if (bitfield.get_representation() == matrix_flags::dense_representation) {
1432 Dm_bigint_modul.det(A, DET, H);
1433 if (DET == 0)
1434 precondition_error_handler(DET, "det(A)", "det(A) != 0", "void matrix< bigint >::"
1435 "adj2(const matrix< bigint > & A, const bigint &H)",
1436 DMESSAGE, EMESSAGE[7]);
1437 Dm_bigint_modul.adj2(*this, A, H, DET);
1438 }
1439 else {
1440 Sm_bigint_modul.det(A, DET, H);
1441 if (DET == 0)
1442 precondition_error_handler(DET, "det(A)", "det(A) != 0", "void matrix< bigint >::"
1443 "adj2(const matrix< bigint > & A, const bigint &H)",
1444 DMESSAGE, EMESSAGE[7]);
1445 Sm_bigint_modul.adj2(*this, A, H, DET);
1446 }
1447 }
1448
1449
1450
1451 void matrix< bigint >::
adj2(const matrix<bigint> & A)1452 adj2(const matrix< bigint > & A)
1453 {
1454 //
1455 // Task: B.adj(A);
1456 // => B = adjoint matrix of matrix A
1457 // Conditions: A.columns != A.rows
1458 // Version: 2.0
1459 //
1460
1461 debug_handler_l(DMESSAGE, "adj2(const matrix< bigint > &)", DVALUE + 4);
1462
1463 if (A.columns != A.rows)
1464 precondition_error_handler(A.columns, "A.columns", "A.columns == A.rows",
1465 A.rows, "A.rows", "A.columns == A.rows",
1466 "void matrix< bigint >::"
1467 "adj2(const matrix< bigint > & A)", DMESSAGE, EMESSAGE[7]);
1468
1469 bigint H, DET;
1470
1471 const modular_arithmetic< DRMK< bigint >,
1472 dense_fp_matrix_kernel< long, MR< long > >,
1473 dense_fp_matrix_kernel< bigint, MR< bigint > > > Dm_bigint_modul;
1474
1475 const modular_arithmetic< SRMK< bigint >,
1476 sparse_fp_matrix_kernel< long, MR< long > >,
1477 sparse_fp_matrix_kernel< bigint, MR< bigint > > > Sm_bigint_modul;
1478
1479 if (bitfield.get_representation() == matrix_flags::dense_representation) {
1480 D_bigint_modul.hadamard(A, H);
1481 Dm_bigint_modul.det(A, DET, H);
1482 if (DET == 0)
1483 precondition_error_handler(DET, "det(A)", "det(A) != 0", "void matrix< bigint >::"
1484 "adj2(const matrix< bigint > & A)", DMESSAGE, EMESSAGE[7]);
1485 Dm_bigint_modul.adj2(*this, A, H, DET);
1486 }
1487 else {
1488 S_bigint_modul.hadamard(A, H);
1489 Sm_bigint_modul.det(A, DET, H);
1490 if (DET == 0)
1491 precondition_error_handler(DET, "det(A)", "det(A) != 0", "void matrix< bigint >::"
1492 "adj2(const matrix< bigint > & A)", DMESSAGE, EMESSAGE[7]);
1493 Sm_bigint_modul.adj2(*this, A, H, DET);
1494 }
1495 }
1496
1497
1498
1499 //
1500 // lattice determinant
1501 //
1502
1503 void matrix< bigint >::
latticedet1(bigint & DET,const bigint & H) const1504 latticedet1(bigint & DET, const bigint &H) const
1505 {
1506 //
1507 // Task: A.latticedet1(DET, H);
1508 // => DET = lattice determinant
1509 // of the lattice formed by the columns of matrix A
1510 // => H = hadamard(A)
1511 // Version: 2.0
1512 //
1513
1514 debug_handler_l(DMESSAGE, "latticedet1(bigint &, const bigint &)", DVALUE + 4);
1515
1516 const modular_arithmetic< DRMK< bigint >,
1517 dense_fp_matrix_kernel< long, MR< long > >,
1518 dense_fp_matrix_kernel< bigint, MR< bigint > > > Dm_bigint_modul;
1519
1520 const modular_arithmetic< SRMK< bigint >,
1521 sparse_fp_matrix_kernel< long, MR< long > >,
1522 sparse_fp_matrix_kernel< bigint, MR< bigint > > > Sm_bigint_modul;
1523
1524 if (bitfield.get_representation() == matrix_flags::dense_representation)
1525 Dm_bigint_modul.latticedet1(*this, DET, H);
1526 else
1527 Sm_bigint_modul.latticedet1(*this, DET, H);
1528 }
1529
1530
1531
1532 void matrix< bigint >::
latticedet1(bigint & DET) const1533 latticedet1(bigint & DET) const
1534 {
1535 //
1536 // Task: A.latticedet1(DET);
1537 // => DET = lattice determinant
1538 // of the lattice formed by the columns of matrix A
1539 // Version: 2.0
1540 //
1541
1542 debug_handler_l(DMESSAGE, "latticedet1(bigint &)", DVALUE + 4);
1543
1544 bigint H;
1545
1546 const modular_arithmetic< DRMK< bigint >,
1547 dense_fp_matrix_kernel< long, MR< long > >,
1548 dense_fp_matrix_kernel< bigint, MR< bigint > > > Dm_bigint_modul;
1549
1550 const modular_arithmetic< SRMK< bigint >,
1551 sparse_fp_matrix_kernel< long, MR< long > >,
1552 sparse_fp_matrix_kernel< bigint, MR< bigint > > > Sm_bigint_modul;
1553
1554 if (bitfield.get_representation() == matrix_flags::dense_representation) {
1555 D_bigint_modul.hadamard(*this, H);
1556 Dm_bigint_modul.latticedet1(*this, DET, H);
1557 }
1558 else {
1559 S_bigint_modul.hadamard(*this, H);
1560 Sm_bigint_modul.latticedet1(*this, DET, H);
1561 }
1562 }
1563
1564
1565
1566 void matrix< bigint >::
latticedet2(bigint & DET,const bigint & H) const1567 latticedet2(bigint & DET, const bigint &H) const
1568 {
1569 //
1570 // Task: A.latticedet2(DET, H);
1571 // => DET = lattice determinant
1572 // of the lattice formed by the columns of matrix A
1573 // => H = hadamard(A)
1574 // Version: 2.0
1575 //
1576
1577 debug_handler_l(DMESSAGE, "latticedet2(bigint &, const bigint &)", DVALUE + 4);
1578
1579 const modular_arithmetic< DRMK< bigint >,
1580 dense_fp_matrix_kernel< long, MR< long > >,
1581 dense_fp_matrix_kernel< bigint, MR< bigint > > > Dm_bigint_modul;
1582
1583 const modular_arithmetic< SRMK< bigint >,
1584 sparse_fp_matrix_kernel< long, MR< long > >,
1585 sparse_fp_matrix_kernel< bigint, MR< bigint > > > Sm_bigint_modul;
1586
1587 if (bitfield.get_representation() == matrix_flags::dense_representation)
1588 Dm_bigint_modul.latticedet2(*this, DET, H);
1589 else
1590 Sm_bigint_modul.latticedet2(*this, DET, H);
1591 }
1592
1593
1594
1595 void matrix< bigint >::
latticedet2(bigint & DET) const1596 latticedet2(bigint & DET) const
1597 {
1598 //
1599 // Task: A.latticedet2(DET);
1600 // => DET = lattice determinant
1601 // of the lattice formed by the columns of matrix A
1602 // Version: 2.0
1603 //
1604
1605 debug_handler_l(DMESSAGE, "latticedet2(bigint &)", DVALUE + 4);
1606
1607 bigint H;
1608
1609 const modular_arithmetic< DRMK< bigint >,
1610 dense_fp_matrix_kernel< long, MR< long > >,
1611 dense_fp_matrix_kernel< bigint, MR< bigint > > > Dm_bigint_modul;
1612
1613 const modular_arithmetic< SRMK< bigint >,
1614 sparse_fp_matrix_kernel< long, MR< long > >,
1615 sparse_fp_matrix_kernel< bigint, MR< bigint > > > Sm_bigint_modul;
1616
1617 if (bitfield.get_representation() == matrix_flags::dense_representation) {
1618 D_bigint_modul.hadamard(*this, H);
1619 Dm_bigint_modul.latticedet2(*this, DET, H);
1620 }
1621 else {
1622 S_bigint_modul.hadamard(*this, H);
1623 Sm_bigint_modul.latticedet2(*this, DET, H);
1624 }
1625 }
1626
1627
1628
1629 void matrix< bigint >::
latticedet3(bigint & DET,const bigint & H) const1630 latticedet3(bigint & DET, const bigint &H) const
1631 {
1632 //
1633 // Task: A.latticedet3(DET, H);
1634 // => DET = lattice determinant
1635 // of the lattice formed by the columns of matrix A
1636 // => H = hadamard(A)
1637 // Version: 2.0
1638 //
1639
1640 debug_handler_l(DMESSAGE, "latticedet3(bigint &, const bigint &)", DVALUE + 4);
1641
1642 const modular_arithmetic< DRMK< bigint >,
1643 dense_fp_matrix_kernel< long, MR< long > >,
1644 dense_fp_matrix_kernel< bigint, MR< bigint > > > Dm_bigint_modul;
1645
1646 const modular_arithmetic< SRMK< bigint >,
1647 sparse_fp_matrix_kernel< long, MR< long > >,
1648 sparse_fp_matrix_kernel< bigint, MR< bigint > > > Sm_bigint_modul;
1649
1650 if (bitfield.get_representation() == matrix_flags::dense_representation)
1651 Dm_bigint_modul.latticedet3(*this, DET, H);
1652 else
1653 Sm_bigint_modul.latticedet3(*this, DET, H);
1654 }
1655
1656
1657
1658 void matrix< bigint >::
latticedet3(bigint & DET) const1659 latticedet3(bigint & DET) const
1660 {
1661 //
1662 // Task: A.latticedet3(DET);
1663 // => DET = lattice determinant
1664 // of the lattice formed by the columns of matrix A
1665 // Version: 2.0
1666 //
1667
1668 debug_handler_l(DMESSAGE, "latticedet3(bigint &)", DVALUE + 4);
1669
1670 bigint H;
1671
1672 const modular_arithmetic< DRMK< bigint >,
1673 dense_fp_matrix_kernel< long, MR< long > >,
1674 dense_fp_matrix_kernel< bigint, MR< bigint > > > Dm_bigint_modul;
1675
1676 const modular_arithmetic< SRMK< bigint >,
1677 sparse_fp_matrix_kernel< long, MR< long > >,
1678 sparse_fp_matrix_kernel< bigint, MR< bigint > > > Sm_bigint_modul;
1679
1680 if (bitfield.get_representation() == matrix_flags::dense_representation) {
1681 D_bigint_modul.hadamard(*this, H);
1682 Dm_bigint_modul.latticedet3(*this, DET, H);
1683 }
1684 else {
1685 S_bigint_modul.hadamard(*this, H);
1686 Sm_bigint_modul.latticedet3(*this, DET, H);
1687 }
1688 }
1689
1690
1691
1692 void matrix< bigint >::
latticedet_special(bigint & DET) const1693 latticedet_special(bigint & DET) const
1694 {
1695 //
1696 // Task: A.latticedet3(DET);
1697 // => DET = lattice determinant
1698 // of the lattice formed by the columns of matrix A
1699 // Version: 2.0
1700 //
1701
1702 debug_handler_l(DMESSAGE, "latticedet3(bigint &)", DVALUE + 4);
1703
1704 bigint H;
1705
1706 const modular_arithmetic< DRMK< bigint >,
1707 dense_fp_matrix_kernel< long, MR< long > >,
1708 dense_fp_matrix_kernel< bigint, MR< bigint > > > Dm_bigint_modul;
1709
1710 const modular_arithmetic< SRMK< bigint >,
1711 sparse_fp_matrix_kernel< long, MR< long > >,
1712 sparse_fp_matrix_kernel< bigint, MR< bigint > > > Sm_bigint_modul;
1713
1714 if (bitfield.get_representation() == matrix_flags::dense_representation) {
1715 D_bigint_modul.hadamard(*this, H);
1716 Dm_bigint_modul.latticedet_special(*this, DET, H, 5);
1717 }
1718 else {
1719 S_bigint_modul.hadamard(*this, H);
1720 Sm_bigint_modul.latticedet_special(*this, DET, H, 5);
1721 }
1722 }
1723
1724
1725
1726 void matrix< bigint >::
real_latticedet(bigint & DET,const bigint & H) const1727 real_latticedet(bigint & DET, const bigint &H) const
1728 {
1729 //
1730 // Task: A.real_latticedet(DET, H);
1731 // => DET = lattice determinant
1732 // of the lattice formed by the columns of matrix A
1733 // => H = hadamard(A)
1734 // Version: 2.0
1735 //
1736
1737 debug_handler_l(DMESSAGE, "real_latticedet(bigint &, const bigint &)",
1738 DVALUE + 4);
1739
1740 const modular_arithmetic< DRMK< bigint >,
1741 dense_fp_matrix_kernel< long, MR< long > >,
1742 dense_fp_matrix_kernel< bigint, MR< bigint > > > Dm_bigint_modul;
1743
1744 const modular_arithmetic< SRMK< bigint >,
1745 sparse_fp_matrix_kernel< long, MR< long > >,
1746 sparse_fp_matrix_kernel< bigint, MR< bigint > > > Sm_bigint_modul;
1747
1748 if (bitfield.get_representation() == matrix_flags::dense_representation) {
1749 matrix< bigint > A(columns, rows);
1750 A.trans(*this);
1751 LiDIA::multiply(A, A, *this);
1752 Dm_bigint_modul.det(A, DET, H*H);
1753 }
1754 else {
1755 sparse_fp_matrix_algorithms< long, ring_matrix< long > > long_modul;
1756 sparse_fp_matrix_algorithms< bigint, ring_matrix< bigint > > bigint_modul;
1757
1758 // Step 1,2
1759 bigint *PRIM = get_primes(bigint(H*H), bigint(1));
1760 long n, Modlong;
1761 PRIM[0].longify(n);
1762
1763 bigint Gmod;
1764
1765 // Step 3
1766 bigint *chininput = new bigint[n];
1767 memory_handler(chininput, DMESSAGE, "det :: "
1768 "Error in memory allocation (chininput)");
1769
1770 lidia_size_t i;
1771 for (i = 1; i <= n; i++) {
1772 Gmod.assign(PRIM[i]);
1773 if (Gmod.bit_length() >= bigint::bits_per_digit()) {
1774 // bigint part
1775 base_power_product< ring_matrix< bigint >, lidia_size_t > bpp;
1776 matrix< bigint > A(rows, columns);
1777 matrix< bigint > AT(columns, rows);
1778 remainder(A, *this, Gmod);
1779 AT.trans(A);
1780
1781 bpp.append(AT);
1782 bpp.append(A);
1783
1784 //chininput[i - 1] = bigint_modul.det(bpp, Gmod);
1785 chininput[i - 1] = bigint_modul.det(A, Gmod);
1786 }
1787 else {
1788 // long part
1789 Gmod.longify(Modlong);
1790 base_power_product< ring_matrix< long >, lidia_size_t > bpp;
1791 matrix< long > A(rows, columns);
1792 matrix< long > AT(columns, rows);
1793 A.set_zero_element(0);
1794 AT.set_zero_element(0);
1795 remainder(A, *this, Modlong);
1796 AT.trans(A);
1797
1798 bpp.append(AT);
1799 bpp.append(A);
1800
1801 //std::cout << "Prime " << i << "/" << n << std::endl;
1802 //chininput[i - 1] = long_modul.det(bpp, Modlong);
1803 chininput[i - 1] = long_modul.det(A, Modlong);
1804 }
1805 }
1806
1807 // Step 4
1808 LiDIA::chinrest(DET, static_cast<const bigint *>(chininput), static_cast<const bigint *>(PRIM));
1809 delete[] chininput;
1810 }
1811 }
1812
1813
1814
1815 //
1816 // determinant
1817 //
1818
1819 void matrix< bigint >::
det(bigint & DET,const bigint & H) const1820 det(bigint & DET, const bigint &H) const
1821 {
1822 //
1823 // Task: A.det(DET, H)
1824 // => DET = determinant of matrix A
1825 // => H = hadamard(A)
1826 // Conditions: columns != rows
1827 // Version: 2.0
1828 //
1829
1830 debug_handler_l(DMESSAGE, "det(bigint &, const bigint &)", DVALUE + 4);
1831
1832 if (rows != columns)
1833 precondition_error_handler(rows, "rows", "rows == columns",
1834 columns, "columns", "rows == columns",
1835 "void matrix< bigint >::"
1836 "det(bigint &DET, const bigint &H) const",
1837 DMESSAGE, EMESSAGE[7]);
1838
1839 const modular_arithmetic< DRMK< bigint >,
1840 dense_fp_matrix_kernel< long, MR< long > >,
1841 dense_fp_matrix_kernel< bigint, MR< bigint > > > Dm_bigint_modul;
1842
1843 const modular_arithmetic< SRMK< bigint >,
1844 sparse_fp_matrix_kernel< long, MR< long > >,
1845 sparse_fp_matrix_kernel< bigint, MR< bigint > > > Sm_bigint_modul;
1846
1847 if (bitfield.get_representation() == matrix_flags::dense_representation)
1848 Dm_bigint_modul.det(*this, DET, H);
1849 else
1850 Sm_bigint_modul.det(*this, DET, H);
1851 }
1852
1853
1854
1855 void matrix< bigint >::
det(bigint & DET,const bigint & H,int num) const1856 det(bigint & DET, const bigint &H, int num) const
1857 {
1858 //
1859 // Task: A.det(DET, H, num)
1860 // => DET = determinant of matrix A
1861 // Conditions: columns != rows
1862 // Version: 2.0
1863 //
1864
1865 debug_handler_l(DMESSAGE, "det(bigint &)", DVALUE + 4);
1866
1867 if (rows != columns)
1868 precondition_error_handler(rows, "rows", "rows == columns",
1869 columns, "columns", "rows == columns",
1870 "void matrix< bigint >::"
1871 "det(bigint & DET) const",
1872 DMESSAGE, EMESSAGE[7]);
1873
1874 const modular_arithmetic< DRMK< bigint >,
1875 dense_fp_matrix_kernel< long, MR< long > >,
1876 dense_fp_matrix_kernel< bigint, MR< bigint > > > Dm_bigint_modul;
1877
1878 const modular_arithmetic< SRMK< bigint >,
1879 sparse_fp_matrix_kernel< long, MR< long > >,
1880 sparse_fp_matrix_kernel< bigint, MR< bigint > > > Sm_bigint_modul;
1881
1882 if (bitfield.get_representation() == matrix_flags::dense_representation)
1883 Dm_bigint_modul.det(*this, DET, H, num);
1884 else
1885 Sm_bigint_modul.det(*this, DET, H, num);
1886 }
1887
1888
1889
1890 void matrix< bigint >::
det(bigint & DET) const1891 det(bigint & DET) const
1892 {
1893 //
1894 // Task: A.det(DET)
1895 // => DET = determinant of matrix A
1896 // Conditions: columns != rows
1897 // Version: 2.0
1898 //
1899
1900 debug_handler_l(DMESSAGE, "det(bigint &)", DVALUE + 4);
1901
1902 if (rows != columns)
1903 precondition_error_handler(rows, "rows", "rows == columns",
1904 columns, "columns", "rows == columns",
1905 "void matrix< bigint >::"
1906 "det(bigint & DET) const",
1907 DMESSAGE, EMESSAGE[7]);
1908
1909 bigint H;
1910
1911 const modular_arithmetic< DRMK< bigint >,
1912 dense_fp_matrix_kernel< long, MR< long > >,
1913 dense_fp_matrix_kernel< bigint, MR< bigint > > > Dm_bigint_modul;
1914
1915 const modular_arithmetic< SRMK< bigint >,
1916 sparse_fp_matrix_kernel< long, MR< long > >,
1917 sparse_fp_matrix_kernel< bigint, MR< bigint > > > Sm_bigint_modul;
1918
1919 if (bitfield.get_representation() == matrix_flags::dense_representation) {
1920 D_bigint_modul.hadamard(*this, H);
1921 Dm_bigint_modul.det(*this, DET, H);
1922 }
1923 else {
1924 S_bigint_modul.hadamard(*this, H);
1925 Sm_bigint_modul.det(*this, DET, H);
1926 }
1927 }
1928
1929
1930
1931 //
1932 // characteristic polynomial
1933 //
1934
1935 bigint *matrix< bigint >::
charpoly() const1936 charpoly() const
1937 {
1938 //
1939 // Task: RES = A.charpoly();
1940 // => RES[0],...,RES[r] are the coefficients of
1941 // the characteristic polynomial of matrix A
1942 // Conditions: rows != columns
1943 // Version: 2.0
1944 //
1945
1946 debug_handler_l(DMESSAGE, "charpoly()", DVALUE + 4);
1947
1948 if (columns != rows)
1949 precondition_error_handler(rows, "rows", "rows == columns",
1950 columns, "columns", "rows == columns",
1951 "bigint *matrix< bigint >::"
1952 "charpoly() const",
1953 DMESSAGE, EMESSAGE[7]);
1954
1955 bigint H;
1956 bigint *RES = new bigint[columns + 1];
1957 memory_handler(RES, DMESSAGE, "charpoly :: "
1958 "Error in memory allocation (RES)");
1959
1960 const modular_arithmetic< DRMK< bigint >,
1961 dense_fp_matrix_kernel< long, MR< long > >,
1962 dense_fp_matrix_kernel< bigint, MR< bigint > > > Dm_bigint_modul;
1963
1964 const modular_arithmetic< SRMK< bigint >,
1965 sparse_fp_matrix_kernel< long, MR< long > >,
1966 sparse_fp_matrix_kernel< bigint, MR< bigint > > > Sm_bigint_modul;
1967
1968 if (bitfield.get_representation() == matrix_flags::dense_representation) {
1969 D_bigint_modul.hadamard(*this, H);
1970 Dm_bigint_modul.charpoly(*this, RES, H);
1971 }
1972 else {
1973 S_bigint_modul.hadamard(*this, H);
1974 Sm_bigint_modul.charpoly(*this, RES, H);
1975 }
1976
1977 return RES;
1978 }
1979
1980
1981
1982 void matrix< bigint >::
charpoly(base_vector<bigint> & RES) const1983 charpoly(base_vector< bigint > &RES) const
1984 {
1985 //
1986 // Task: A.charpoly(RES);
1987 // => RES[0],...,RES[r] are the coefficients of
1988 // the characteristic polynomial of matrix A
1989 // Conditions: rows == columns
1990 // Version: 2.0
1991 //
1992
1993 debug_handler_l(DMESSAGE, "charpoly(base_vector< bigint > &)", DVALUE + 4);
1994
1995 if (columns != rows)
1996 precondition_error_handler(rows, "rows", "rows == columns",
1997 columns, "columns", "rows == columns",
1998 "void matrix< bigint >::"
1999 "charpoly(base_vector< bigint > &) const",
2000 DMESSAGE, EMESSAGE[7]);
2001
2002 bigint H;
2003 RES.set_size(columns + 1);
2004
2005 const modular_arithmetic< DRMK< bigint >,
2006 dense_fp_matrix_kernel< long, MR< long > >,
2007 dense_fp_matrix_kernel< bigint, MR< bigint > > > Dm_bigint_modul;
2008
2009 const modular_arithmetic< SRMK< bigint >,
2010 sparse_fp_matrix_kernel< long, MR< long > >,
2011 sparse_fp_matrix_kernel< bigint, MR< bigint > > > Sm_bigint_modul;
2012
2013 if (bitfield.get_representation() == matrix_flags::dense_representation) {
2014 D_bigint_modul.hadamard(*this, H);
2015 Dm_bigint_modul.charpoly(*this, RES.get_data_address(), H);
2016 }
2017 else {
2018 S_bigint_modul.hadamard(*this, H);
2019 Sm_bigint_modul.charpoly(*this, RES.get_data_address(), H);
2020 }
2021 }
2022
2023
2024
2025 /////////////////////////
2026 // END: Linear algebra //
2027 // PART 1 //
2028 /////////////////////////
2029
2030 ///////////////////////////
2031 // BEGIN: Linear algebra //
2032 // PART 2 //
2033 ///////////////////////////
2034
2035 //
2036 // Hermite normal form
2037 //
2038
2039 void matrix< bigint >::
hnfmod_dkt(const bigint & mod)2040 hnfmod_dkt(const bigint &mod)
2041 {
2042 //
2043 // Task: A.hnfmod_dkt(mod);
2044 // => A in Hermite normal form
2045 // => h = lattice determinant of lattice formed
2046 // by the columns of matrix A
2047 // Conditions: rank != rows
2048 // Version: 2.0
2049 //
2050
2051 debug_handler_l(DMESSAGE, "hnfmod_dkt(const bigint &)", DVALUE + 5);
2052
2053 lidia_size_t r = rank();
2054 if (r != rows)
2055 precondition_error_handler(r, "rank", "rank == rows",
2056 rows, "rows", "rank == rows",
2057 "void matrix< bigint >::"
2058 "hnfmod_dkt(const bigint &mod)",
2059 DMESSAGE, EMESSAGE[10]);
2060
2061 if (bitfield.get_representation() == matrix_flags::dense_representation)
2062 D_bigint_modul.hnfmod_dkt(*this, mod);
2063 else
2064 lidia_error_handler("matrix<bigint>",
2065 "hnfmod_dkt: dense matrices supported only");
2066 }
2067
2068
2069
2070 void matrix< bigint >::
hnfmod_dkt(matrix<bigint> & TR,const bigint & mod)2071 hnfmod_dkt(matrix< bigint > &TR, const bigint &mod)
2072 {
2073 //
2074 // Task: A.hnfmod_dkt(TR, mod);
2075 // => A in Hermite normal form
2076 // => h = lattice determinant of lattice formed
2077 // by the columns of matrix A
2078 // Conditions: rank != rows
2079 // Version: 2.0
2080 //
2081
2082 debug_handler_l(DMESSAGE, "hnfmod_dkt(const bigint &)", DVALUE + 5);
2083
2084 if (rank() != rows)
2085 precondition_error_handler(rank(), "rank", "rank == rows",
2086 rows, "rows", "rank == rows",
2087 "void matrix< bigint >::"
2088 "hnfmod_dkt(const bigint &mod)",
2089 DMESSAGE, EMESSAGE[10]);
2090
2091 TR.resize(rows + columns, rows + columns);
2092 TR.diag(1, 0);
2093
2094 if (bitfield.get_representation() == matrix_flags::dense_representation)
2095 D_bigint_modul.hnfmod_dkt(*this, TR, mod);
2096 else
2097 lidia_error_handler("matrix<bigint>",
2098 "hnfmod_dkt: dense matrices supported only");
2099 }
2100
2101
2102
2103 void matrix< bigint >::
hnfmod_cohen(const bigint & D)2104 hnfmod_cohen(const bigint & D)
2105 {
2106 debug_handler_l(DMESSAGE, "hnfmod_cohen(const bigint &)", DVALUE + 5);
2107
2108 if (rank() != rows)
2109 precondition_error_handler(rank(), "rank", "rows == rank",
2110 rows, "rows", "rows == rank",
2111 "void matrix< bigint >::"
2112 "hnfmod_cohen(const bigint & D)",
2113 DMESSAGE, EMESSAGE[10]);
2114
2115 if (bitfield.get_representation() == matrix_flags::dense_representation)
2116 D_bigint_modul.hnfmod_cohen(*this, D);
2117 else
2118 lidia_error_handler("matrix<bigint>",
2119 "hnfmod_cohen: dense matrices supported only");
2120 }
2121
2122
2123
2124 void matrix< bigint >::
hnfmod_mueller(matrix<bigint> & TRANS)2125 hnfmod_mueller(matrix< bigint > & TRANS)
2126 {
2127 //
2128 // Task: A.hnfmod_mueller(TRANS);
2129 // => A in Hermite normal form
2130 // => TRANS = transformtion matrix A*TRANS = HNFmod(A)
2131 // Conditions: TRANS.rows != TRANS.columns or TRANS.rows != columns or
2132 // rank != rows
2133 // Version: 2.0
2134 //
2135
2136 debug_handler_l(DMESSAGE, "hnfmod_mueller(matrix< bigint > &, const bigint &)", DVALUE + 5);
2137
2138 const modular_arithmetic< DRMK< bigint >, dense_fp_matrix_kernel< long, MR< long > >,
2139 dense_fp_matrix_kernel< bigint, MR< bigint > > > Dm_bigint_modul;
2140
2141 lidia_size_t *linuz = lininc();
2142 if (linuz[0] != rows)
2143 precondition_error_handler(linuz[0], "rank", "rank == rows",
2144 rows, "rows", "rank == rows",
2145 "void matrix< bigint >::"
2146 "hnfmod_mueller(matrix< bigint > & TRANS)",
2147 DMESSAGE, EMESSAGE[10]);
2148
2149 if (bitfield.get_representation() == matrix_flags::dense_representation) {
2150 bigint DET, H;
2151 D_bigint_modul.hadamard(*this, H);
2152 Dm_bigint_modul.latticedet1(*this, DET, H);
2153 D_bigint_modul.hnfmod_mueller(*this, TRANS, DET);
2154 }
2155 else
2156 lidia_error_handler("matrix<bigint>",
2157 "hnfmod_mueller: dense matrices supported only");
2158 }
2159
2160
2161
2162 void matrix< bigint >::
hnf_ref()2163 hnf_ref()
2164 {
2165 if (bitfield.get_representation() == matrix_flags::dense_representation)
2166 hnf_ref_modul.hnf_havas(*this);
2167 else
2168 lidia_error_handler("matrix<bigint>",
2169 "hnf_ref: dense matrices supported only");
2170 }
2171
2172
2173
2174
2175 //
2176 // hnf_storjohann
2177 //
2178
2179 void matrix< bigint >::
hnf_storjohann()2180 hnf_storjohann()
2181 {
2182 //
2183 // Task: HNF Computation
2184 // Algorithm: Gauss with reduction
2185 // IMPROVEMENTS: Theory of Havas / best reaminder
2186 // Version: 2.0
2187 //
2188
2189 debug_handler_l(DMESSAGE, "hnf_storjohann()", DVALUE + 5);
2190
2191 if (bitfield.get_representation() == matrix_flags::dense_representation) {
2192 D_bigint_modul.hnf_storjohann(*this);
2193 }
2194 else {
2195 S_bigint_modul.hnf_storjohann(*this);
2196 }
2197 }
2198
2199
2200
2201 void matrix< bigint >::
hnf_storjohann(matrix<bigint> & T,matrix<bigint> & C,matrix<bigint> & Q)2202 hnf_storjohann(matrix< bigint > &T, matrix< bigint > &C, matrix< bigint > &Q)
2203 {
2204 //
2205 // Task: HNF Computation
2206 // Algorithm: Gauss with reduction
2207 // IMPROVEMENTS: Theory of Havas / best reaminder
2208 // Version: 2.0
2209 //
2210
2211 debug_handler_l(DMESSAGE, "hnf_storjohann()", DVALUE + 5);
2212
2213 lidia_size_t i = rows, j = columns;
2214 if (T.columns != columns)
2215 T.set_no_of_columns(columns);
2216 if (T.rows != columns)
2217 T.set_no_of_rows(columns);
2218
2219 if (C.columns != columns)
2220 C.set_no_of_columns(columns);
2221 if (C.rows != columns)
2222 C.set_no_of_rows(columns);
2223
2224 if (Q.columns != columns)
2225 Q.set_no_of_columns(columns);
2226 if (Q.rows != columns)
2227 Q.set_no_of_rows(columns);
2228
2229 T.diag(1, 0);
2230 C.diag(1, 0);
2231 Q.diag(1, 0);
2232
2233 if (bitfield.get_representation() == matrix_flags::dense_representation)
2234 D_bigint_modul.hnf_storjohann(*this, T, C, Q);
2235 else {
2236
2237 set_orientation(matrix_flags::column_oriented);
2238 i = 0;
2239 j = 0;
2240 // hnf_smodul2.hnf(*this, i, j);
2241 set_orientation(matrix_flags::row_oriented);
2242 }
2243 }
2244
2245
2246
2247 //
2248 // hnf_havas
2249 //
2250
2251 #define schema_havas_call_sequenz(name) \
2252 {if (bitfield.get_representation() == matrix_flags::dense_representation)\
2253 { \
2254 havas_kernel< bigint, RODMM< bigint >, RODMM< bigint >, \
2255 name< bigint, RODMM< bigint >, RODMM< bigint > > > modul; \
2256 modul.hnf_Z1(*this, i, j); \
2257 } \
2258 else { \
2259 havas_kernel< bigint, COSMM< bigint >, COSMM< bigint >, \
2260 name< bigint, COSMM< bigint >, COSMM< bigint > > > modul; \
2261 \
2262 set_orientation(matrix_flags::column_oriented); \
2263 modul.hnf_Z1(*this, i, j); \
2264 set_orientation(matrix_flags::row_oriented); \
2265 } \
2266 break; }
2267
2268 #define schema_havas_stf_call_sequenz(name) \
2269 {if (bitfield.get_representation() == matrix_flags::dense_representation)\
2270 { \
2271 havas_kernel< bigint, RODMM< bigint >, RODMM< bigint >, \
2272 name< bigint, RODMM< bigint >, RODMM< bigint > > > modul; \
2273 modul.hnf_Z2(*this, i, j); \
2274 normalization_kernel< bigint, RODMM< bigint >, \
2275 RODMM< bigint > > nmodul; \
2276 \
2277 switch(normalizeModul) { \
2278 case 0: \
2279 { \
2280 nmodul.normalize_Std(*this, 0, columns - rows); \
2281 break; \
2282 } \
2283 case 1: \
2284 { \
2285 nmodul.normalize_ChouCollins(*this, 0, columns - rows); \
2286 break; \
2287 } \
2288 case 2: \
2289 { \
2290 nmodul.normalizeMod_Std(*this, 0, columns - rows); \
2291 break; \
2292 } \
2293 case 3: \
2294 { \
2295 nmodul.normalizeMod_ChouCollins(*this, 0, columns - rows); \
2296 break; \
2297 } \
2298 case 4: \
2299 { \
2300 nmodul.normalizeHybrid_Std(*this, 0, columns - rows); \
2301 break; \
2302 } \
2303 case 5: \
2304 { \
2305 nmodul.normalizeHybrid_ChouCollins(*this, 0, columns - rows); \
2306 break; \
2307 } \
2308 default: \
2309 lidia_error_handler("bigint_matrix", "hnf_havas :: " \
2310 "Mode not supported! (normalize " \
2311 "Modul not defined)"); \
2312 } \
2313 } \
2314 else { \
2315 havas_kernel< bigint, COSMM< bigint >, COSMM< bigint >, \
2316 name< bigint, COSMM< bigint >, COSMM< bigint > > > modul; \
2317 \
2318 set_orientation(matrix_flags::column_oriented); \
2319 modul.hnf_Z2(*this, i, j); \
2320 normalization_kernel< bigint, COSMM< bigint >, \
2321 COSMM< bigint > > nmodul; \
2322 \
2323 switch(normalizeModul) { \
2324 case 0: \
2325 { \
2326 nmodul.normalize_Std(*this, 0, columns - rows); \
2327 break; \
2328 } \
2329 case 1: \
2330 { \
2331 nmodul.normalize_ChouCollins(*this, 0, columns - rows); \
2332 break; \
2333 } \
2334 case 2: \
2335 { \
2336 nmodul.normalize_Std(*this, 0, columns - rows); \
2337 break; \
2338 } \
2339 case 3: \
2340 { \
2341 nmodul.normalize_Std(*this, 0, columns - rows); \
2342 break; \
2343 } \
2344 case 4: \
2345 { \
2346 nmodul.normalizeHybrid_Std(*this, 0, columns - rows); \
2347 break; \
2348 } \
2349 case 5: \
2350 { \
2351 nmodul.normalizeHybrid_ChouCollins(*this, 0, columns - rows); \
2352 break; \
2353 } \
2354 default: \
2355 lidia_error_handler("bigint_matrix", "hnf_havas :: " \
2356 "Mode not supported! (normalize " \
2357 "Modul not defined)"); \
2358 } \
2359 set_orientation(matrix_flags::row_oriented); \
2360 } \
2361 break; }
2362
2363 void matrix< bigint >::
hnf_havas(lidia_size_t KernAlgo,lidia_size_t mgcdModul,lidia_size_t normalizeModul)2364 hnf_havas(lidia_size_t KernAlgo, lidia_size_t mgcdModul, lidia_size_t normalizeModul)
2365 {
2366 //
2367 // Task: HNF Computation
2368 // Algorithm: Gauss with reduction
2369 // IMPROVEMENTS: Theory of Havas / best reaminder
2370 // Version: 2.0
2371 //
2372
2373 debug_handler_l(DMESSAGE, "hnf_havas()", DVALUE + 5);
2374
2375 lidia_size_t i = rows, j = columns;
2376
2377 switch(KernAlgo) {
2378 case 0:
2379 {
2380 switch(mgcdModul) {
2381 case 0:
2382 schema_havas_call_sequenz(hermite)
2383 case 1:
2384 schema_havas_call_sequenz(Bradley)
2385 case 2:
2386 schema_havas_call_sequenz(Iliopoulos)
2387 case 3:
2388 schema_havas_call_sequenz(opt)
2389 case 4:
2390 schema_havas_call_sequenz(Blankinship)
2391 case 5:
2392 schema_havas_call_sequenz(Best_remainder)
2393 case 6:
2394 schema_havas_call_sequenz(havas_best_remainder)
2395 case 7:
2396 schema_havas_call_sequenz(havas_euclidean_norm)
2397 case 8:
2398 schema_havas_call_sequenz(havas_minimal_norm)
2399 case 9:
2400 schema_havas_call_sequenz(havas_sorting_gcd)
2401 case 10:
2402 schema_havas_call_sequenz(havas_min_no_of_elements)
2403 case 11:
2404 schema_havas_call_sequenz(havas_min_abs_of_row_plus_minimal_euclidean_norm)
2405 case 12:
2406 schema_havas_call_sequenz(havas_min_abs_of_row_plus_minimal_norm)
2407 case 13:
2408 schema_havas_call_sequenz(havas_min_abs_of_row_plus_min_no_of_elements)
2409 case 14:
2410 schema_havas_call_sequenz(havas_minimal_norm_plus_min_abs_of_row)
2411 case 15:
2412 schema_havas_call_sequenz(havas_minimal_norm_plus_sorting_gcd)
2413 case 16:
2414 schema_havas_call_sequenz(havas_minimal_norm_plus_min_no_of_elements)
2415 case 17:
2416 schema_havas_call_sequenz(havas_sorting_gcd_plus_minimal_euclidean_norm)
2417 case 18:
2418 schema_havas_call_sequenz(havas_sorting_gcd_plus_minimal_norm)
2419 case 19:
2420 schema_havas_call_sequenz(havas_sorting_gcd_plus_min_no_of_elements)
2421 case 20:
2422 schema_havas_call_sequenz(havas_min_no_of_elements_plus_min_abs_of_row)
2423 case 21:
2424 schema_havas_call_sequenz(havas_min_no_of_elements_plus_minimal_euclidean_norm)
2425 case 22:
2426 schema_havas_call_sequenz(havas_min_no_of_elements_plus_minimal_norm)
2427 case 23:
2428 schema_havas_call_sequenz(havas_min_no_of_elements_plus_sorting_gcd)
2429 case 24:
2430 schema_havas_call_sequenz(Storjohann)
2431 case 25:
2432 schema_havas_call_sequenz(Heuristik)
2433 default:
2434 lidia_error_handler("bigint_matrix", "hnf_havas :: Mode not supported! (mgcd Modul not defined)");
2435 }
2436 break;
2437 }
2438 case 1:
2439 {
2440 switch(mgcdModul) {
2441 case 0:
2442 schema_havas_stf_call_sequenz(hermite)
2443 case 1:
2444 schema_havas_stf_call_sequenz(Bradley)
2445 case 2:
2446 schema_havas_stf_call_sequenz(Iliopoulos)
2447 case 3:
2448 schema_havas_stf_call_sequenz(opt)
2449 case 4:
2450 schema_havas_stf_call_sequenz(Blankinship)
2451 case 5:
2452 schema_havas_stf_call_sequenz(Best_remainder)
2453 case 6:
2454 schema_havas_stf_call_sequenz(havas_best_remainder)
2455 case 7:
2456 schema_havas_stf_call_sequenz(havas_euclidean_norm)
2457 case 8:
2458 schema_havas_stf_call_sequenz(havas_minimal_norm)
2459 case 9:
2460 schema_havas_stf_call_sequenz(havas_sorting_gcd)
2461 case 10:
2462 schema_havas_stf_call_sequenz(havas_min_no_of_elements)
2463 case 11:
2464 schema_havas_stf_call_sequenz(havas_min_abs_of_row_plus_minimal_euclidean_norm)
2465 case 12:
2466 schema_havas_stf_call_sequenz(havas_min_abs_of_row_plus_minimal_norm)
2467 case 13:
2468 schema_havas_stf_call_sequenz(havas_min_abs_of_row_plus_min_no_of_elements)
2469 case 14:
2470 schema_havas_stf_call_sequenz(havas_minimal_norm_plus_min_abs_of_row)
2471 case 15:
2472 schema_havas_stf_call_sequenz(havas_minimal_norm_plus_sorting_gcd)
2473 case 16:
2474 schema_havas_stf_call_sequenz(havas_minimal_norm_plus_min_no_of_elements)
2475 case 17:
2476 schema_havas_stf_call_sequenz(havas_sorting_gcd_plus_minimal_euclidean_norm)
2477 case 18:
2478 schema_havas_stf_call_sequenz(havas_sorting_gcd_plus_minimal_norm)
2479 case 19:
2480 schema_havas_stf_call_sequenz(havas_sorting_gcd_plus_min_no_of_elements)
2481 case 20:
2482 schema_havas_stf_call_sequenz(havas_min_no_of_elements_plus_min_abs_of_row)
2483 case 21:
2484 schema_havas_stf_call_sequenz(havas_min_no_of_elements_plus_minimal_euclidean_norm)
2485 case 22:
2486 schema_havas_stf_call_sequenz(havas_min_no_of_elements_plus_minimal_norm)
2487 case 23:
2488 schema_havas_stf_call_sequenz(havas_min_no_of_elements_plus_sorting_gcd)
2489 case 24:
2490 schema_havas_stf_call_sequenz(Storjohann)
2491 case 25:
2492 schema_havas_stf_call_sequenz(Heuristik)
2493 default:
2494 lidia_error_handler("bigint_matrix", "hnf_havas :: Mode not supported! (mgcd Modul not defined)");
2495 }
2496 }
2497 }
2498 }
2499
2500
2501
2502 #define schema_havas_call_sequenz_ex(name) \
2503 {if (bitfield.get_representation() == matrix_flags::dense_representation)\
2504 { \
2505 havas_kernel< bigint, RODMM< bigint >, RODMM< bigint >, \
2506 name< bigint, RODMM< bigint >, RODMM< bigint > > > modul; \
2507 modul.hnf_Z1(*this, T, i, j); \
2508 } \
2509 else \
2510 { \
2511 havas_kernel< bigint, COSMM< bigint >, COSMM< bigint >, \
2512 name< bigint, COSMM< bigint >, COSMM< bigint > > > modul; \
2513 \
2514 set_orientation(matrix_flags::column_oriented); \
2515 modul.hnf_Z1(*this, T, i, j); \
2516 set_orientation(matrix_flags::row_oriented); \
2517 } \
2518 break; }
2519
2520 #define schema_havas_stf_call_sequenz_ex(name) \
2521 {if (bitfield.get_representation() == matrix_flags::dense_representation)\
2522 { \
2523 havas_kernel< bigint, RODMM< bigint >, RODMM< bigint >, \
2524 name< bigint, RODMM< bigint >, RODMM< bigint > > > modul; \
2525 modul.hnf_Z2(*this, T, i, j); \
2526 normalization_kernel< bigint, RODMM< bigint >, \
2527 RODMM< bigint > > nmodul; \
2528 \
2529 switch(normalizeModul) \
2530 { \
2531 case 0: \
2532 { \
2533 nmodul.normalize_Std(*this, T, 0, columns - rows); \
2534 break; \
2535 } \
2536 case 1: \
2537 { \
2538 nmodul.normalize_ChouCollins(*this, T, 0, columns - rows); \
2539 break; \
2540 } \
2541 case 2: \
2542 { \
2543 nmodul.normalizeHybrid_Std(*this, T, 0, columns - rows); \
2544 break; \
2545 } \
2546 case 3: \
2547 { \
2548 nmodul.normalizeHybrid_ChouCollins(*this, T, 0, columns - rows); \
2549 break; \
2550 } \
2551 default: \
2552 lidia_error_handler("bigint_matrix", "hnf_havas :: " \
2553 "Mode not supported! (normalize " \
2554 "Modul not defined)"); \
2555 } \
2556 } \
2557 else \
2558 { \
2559 havas_kernel< bigint, COSMM< bigint >, COSMM< bigint >, \
2560 name< bigint, COSMM< bigint >, COSMM< bigint > > > modul; \
2561 \
2562 set_orientation(matrix_flags::column_oriented); \
2563 modul.hnf_Z2(*this, T, i, j); \
2564 normalization_kernel< bigint, COSMM< bigint >, \
2565 COSMM< bigint > > nmodul; \
2566 \
2567 switch(normalizeModul) \
2568 { \
2569 case 0: \
2570 { \
2571 nmodul.normalize_Std(*this, T, 0, columns - rows); \
2572 break; \
2573 } \
2574 case 1: \
2575 { \
2576 nmodul.normalize_ChouCollins(*this, T, 0, columns - rows); \
2577 break; \
2578 } \
2579 case 2: \
2580 { \
2581 nmodul.normalizeHybrid_Std(*this, T, 0, columns - rows); \
2582 break; \
2583 } \
2584 case 3: \
2585 { \
2586 nmodul.normalizeHybrid_ChouCollins(*this, T, 0, columns - rows); \
2587 break; \
2588 } \
2589 default: \
2590 lidia_error_handler("bigint_matrix", "hnf_havas :: " \
2591 "Mode not supported! (normalize " \
2592 "Modul not defined)"); \
2593 } \
2594 set_orientation(matrix_flags::row_oriented); \
2595 } \
2596 break; }
2597
2598 void matrix< bigint >::
hnf_havas(matrix<bigint> & T,lidia_size_t KernAlgo,lidia_size_t mgcdModul,lidia_size_t normalizeModul)2599 hnf_havas(matrix< bigint > & T, lidia_size_t KernAlgo, lidia_size_t mgcdModul, lidia_size_t normalizeModul)
2600 {
2601 //
2602 // Task: HNF Computation
2603 // Algorithm: Gauss with reduction
2604 // IMPROVEMENTS: Theory of Havas / best reaminder
2605 // Version: 2.0
2606 //
2607
2608 debug_handler_l(DMESSAGE, "hnf_havas(matrix< bigint > &)", DVALUE + 5);
2609
2610 lidia_size_t i = rows, j = columns;
2611
2612 if (T.columns != columns)
2613 T.set_no_of_columns(columns);
2614 if (T.rows != columns)
2615 T.set_no_of_rows(columns);
2616
2617 T.diag(1, 0);
2618
2619 switch(KernAlgo) {
2620 case 0:
2621 {
2622 switch(mgcdModul) {
2623 case 0:
2624 schema_havas_call_sequenz_ex(hermite)
2625 case 1:
2626 schema_havas_call_sequenz_ex(Bradley)
2627 case 2:
2628 schema_havas_call_sequenz_ex(Iliopoulos)
2629 case 3:
2630 schema_havas_call_sequenz_ex(opt)
2631 case 4:
2632 schema_havas_call_sequenz_ex(Blankinship)
2633 case 5:
2634 schema_havas_call_sequenz_ex(Best_remainder)
2635 case 6:
2636 schema_havas_call_sequenz_ex(havas_best_remainder)
2637 case 7:
2638 schema_havas_call_sequenz_ex(havas_euclidean_norm)
2639 case 8:
2640 schema_havas_call_sequenz_ex(havas_minimal_norm)
2641 case 9:
2642 schema_havas_call_sequenz_ex(havas_sorting_gcd)
2643 case 10:
2644 schema_havas_call_sequenz_ex(havas_min_no_of_elements)
2645 case 11:
2646 schema_havas_call_sequenz_ex(havas_min_abs_of_row_plus_minimal_euclidean_norm)
2647 case 12:
2648 schema_havas_call_sequenz_ex(havas_min_abs_of_row_plus_minimal_norm)
2649 case 13:
2650 schema_havas_call_sequenz_ex(havas_min_abs_of_row_plus_min_no_of_elements)
2651 case 14:
2652 schema_havas_call_sequenz_ex(havas_minimal_norm_plus_min_abs_of_row)
2653 case 15:
2654 schema_havas_call_sequenz_ex(havas_minimal_norm_plus_sorting_gcd)
2655 case 16:
2656 schema_havas_call_sequenz_ex(havas_minimal_norm_plus_min_no_of_elements)
2657 case 17:
2658 schema_havas_call_sequenz_ex(havas_sorting_gcd_plus_minimal_euclidean_norm)
2659 case 18:
2660 schema_havas_call_sequenz_ex(havas_sorting_gcd_plus_minimal_norm)
2661 case 19:
2662 schema_havas_call_sequenz_ex(havas_sorting_gcd_plus_min_no_of_elements)
2663 case 20:
2664 schema_havas_call_sequenz_ex(havas_min_no_of_elements_plus_min_abs_of_row)
2665 case 21:
2666 schema_havas_call_sequenz_ex(havas_min_no_of_elements_plus_minimal_euclidean_norm)
2667 case 22:
2668 schema_havas_call_sequenz_ex(havas_min_no_of_elements_plus_minimal_norm)
2669 case 23:
2670 schema_havas_call_sequenz_ex(havas_min_no_of_elements_plus_sorting_gcd)
2671 case 24:
2672 schema_havas_call_sequenz_ex(Storjohann)
2673 case 25:
2674 schema_havas_call_sequenz_ex(Heuristik)
2675 default:
2676 lidia_error_handler("bigint_matrix", "hnf_havas :: Mode not supported! (mgcd Modul not defined)");
2677 }
2678 break;
2679 }
2680 case 1:
2681 {
2682 switch(mgcdModul) {
2683 case 0:
2684 schema_havas_stf_call_sequenz_ex(hermite)
2685 case 1:
2686 schema_havas_stf_call_sequenz_ex(Bradley)
2687 case 2:
2688 schema_havas_stf_call_sequenz_ex(Iliopoulos)
2689 case 3:
2690 schema_havas_stf_call_sequenz_ex(opt)
2691 case 4:
2692 schema_havas_stf_call_sequenz_ex(Blankinship)
2693 case 5:
2694 schema_havas_stf_call_sequenz_ex(Best_remainder)
2695 case 6:
2696 schema_havas_stf_call_sequenz_ex(havas_best_remainder)
2697 case 7:
2698 schema_havas_stf_call_sequenz_ex(havas_euclidean_norm)
2699 case 8:
2700 schema_havas_stf_call_sequenz_ex(havas_minimal_norm)
2701 case 9:
2702 schema_havas_stf_call_sequenz_ex(havas_sorting_gcd)
2703 case 10:
2704 schema_havas_stf_call_sequenz_ex(havas_min_no_of_elements)
2705 case 11:
2706 schema_havas_stf_call_sequenz_ex(havas_min_abs_of_row_plus_minimal_euclidean_norm)
2707 case 12:
2708 schema_havas_stf_call_sequenz_ex(havas_min_abs_of_row_plus_minimal_norm)
2709 case 13:
2710 schema_havas_stf_call_sequenz_ex(havas_min_abs_of_row_plus_min_no_of_elements)
2711 case 14:
2712 schema_havas_stf_call_sequenz_ex(havas_minimal_norm_plus_min_abs_of_row)
2713 case 15:
2714 schema_havas_stf_call_sequenz_ex(havas_minimal_norm_plus_sorting_gcd)
2715 case 16:
2716 schema_havas_stf_call_sequenz_ex(havas_minimal_norm_plus_min_no_of_elements)
2717 case 17:
2718 schema_havas_stf_call_sequenz_ex(havas_sorting_gcd_plus_minimal_euclidean_norm)
2719 case 18:
2720 schema_havas_stf_call_sequenz_ex(havas_sorting_gcd_plus_minimal_norm)
2721 case 19:
2722 schema_havas_stf_call_sequenz_ex(havas_sorting_gcd_plus_min_no_of_elements)
2723 case 20:
2724 schema_havas_stf_call_sequenz_ex(havas_min_no_of_elements_plus_min_abs_of_row)
2725 case 21:
2726 schema_havas_stf_call_sequenz_ex(havas_min_no_of_elements_plus_minimal_euclidean_norm)
2727 case 22:
2728 schema_havas_stf_call_sequenz_ex(havas_min_no_of_elements_plus_minimal_norm)
2729 case 23:
2730 schema_havas_stf_call_sequenz_ex(havas_min_no_of_elements_plus_sorting_gcd)
2731 case 24:
2732 schema_havas_stf_call_sequenz_ex(Storjohann)
2733 case 25:
2734 schema_havas_stf_call_sequenz_ex(Heuristik)
2735 default:
2736 lidia_error_handler("bigint_matrix", "hnf_havas :: Mode not supported! (mgcd Modul not defined)");
2737 }
2738 }
2739 }
2740 }
2741
2742
2743
2744 //
2745 // mgcd
2746 //
2747
2748 #define schema_mgcd_call_sequenz(name) \
2749 {if (bitfield.get_representation() == matrix_flags::dense_representation)\
2750 { \
2751 havas_kernel< bigint, RODMM< bigint >, RODMM< bigint >, \
2752 name< bigint, RODMM< bigint >, RODMM< bigint > > > modul; \
2753 modul.mgcd(*this, i, j); \
2754 } \
2755 else \
2756 { \
2757 havas_kernel< bigint, COSMM< bigint >, COSMM< bigint >, \
2758 name< bigint, COSMM< bigint >, COSMM< bigint > > > modul; \
2759 \
2760 set_orientation(matrix_flags::column_oriented); \
2761 modul.mgcd(*this, i, j); \
2762 set_orientation(matrix_flags::row_oriented); \
2763 } \
2764 break; }
2765
2766
2767
2768 void matrix< bigint >::
mgcd(lidia_size_t mgcdModul)2769 mgcd(lidia_size_t mgcdModul)
2770 {
2771 lidia_size_t i = rows, j = columns;
2772
2773 switch(mgcdModul) {
2774 case 0:
2775 schema_mgcd_call_sequenz(hermite)
2776 case 1:
2777 schema_mgcd_call_sequenz(Bradley)
2778 case 2:
2779 schema_mgcd_call_sequenz(Iliopoulos)
2780 case 3:
2781 schema_mgcd_call_sequenz(opt)
2782 case 4:
2783 schema_mgcd_call_sequenz(Blankinship)
2784 case 5:
2785 schema_mgcd_call_sequenz(Best_remainder)
2786 case 6:
2787 schema_mgcd_call_sequenz(havas_best_remainder)
2788 case 7:
2789 schema_mgcd_call_sequenz(havas_euclidean_norm)
2790 case 8:
2791 schema_mgcd_call_sequenz(havas_minimal_norm)
2792 case 9:
2793 schema_mgcd_call_sequenz(havas_sorting_gcd)
2794 case 10:
2795 schema_mgcd_call_sequenz(havas_min_no_of_elements)
2796 case 11:
2797 schema_mgcd_call_sequenz(havas_min_abs_of_row_plus_minimal_euclidean_norm)
2798 case 12:
2799 schema_mgcd_call_sequenz(havas_min_abs_of_row_plus_minimal_norm)
2800 case 13:
2801 schema_mgcd_call_sequenz(havas_min_abs_of_row_plus_min_no_of_elements)
2802 case 14:
2803 schema_mgcd_call_sequenz(havas_minimal_norm_plus_min_abs_of_row)
2804 case 15:
2805 schema_mgcd_call_sequenz(havas_minimal_norm_plus_sorting_gcd)
2806 case 16:
2807 schema_mgcd_call_sequenz(havas_minimal_norm_plus_min_no_of_elements)
2808 case 17:
2809 schema_mgcd_call_sequenz(havas_sorting_gcd_plus_minimal_euclidean_norm)
2810 case 18:
2811 schema_mgcd_call_sequenz(havas_sorting_gcd_plus_minimal_norm)
2812 case 19:
2813 schema_mgcd_call_sequenz(havas_sorting_gcd_plus_min_no_of_elements)
2814 case 20:
2815 schema_mgcd_call_sequenz(havas_min_no_of_elements_plus_min_abs_of_row)
2816 case 21:
2817 schema_mgcd_call_sequenz(havas_min_no_of_elements_plus_minimal_euclidean_norm)
2818 case 22:
2819 schema_mgcd_call_sequenz(havas_min_no_of_elements_plus_minimal_norm)
2820 case 23:
2821 schema_mgcd_call_sequenz(havas_min_no_of_elements_plus_sorting_gcd)
2822 case 24:
2823 schema_mgcd_call_sequenz(Storjohann)
2824 case 25:
2825 schema_mgcd_call_sequenz(Heuristik)
2826 default:
2827 lidia_error_handler("bigint_matrix", "mgcd :: Mode not supported! (mgcd Modul not defined)");
2828 }
2829 }
2830
2831
2832
2833 #define schema_mgcd_call_sequenz_ex(name) \
2834 {if (bitfield.get_representation() == matrix_flags::dense_representation)\
2835 { \
2836 havas_kernel< bigint, RODMM< bigint >, RODMM< bigint >, \
2837 name< bigint, RODMM< bigint >, RODMM< bigint > > > modul; \
2838 modul.mgcd(*this, T, i, j); \
2839 } \
2840 else \
2841 { \
2842 havas_kernel< bigint, COSMM< bigint >, COSMM< bigint >, \
2843 name< bigint, COSMM< bigint >, COSMM< bigint > > > modul; \
2844 \
2845 set_orientation(matrix_flags::column_oriented); \
2846 modul.mgcd(*this, T, i, j); \
2847 set_orientation(matrix_flags::row_oriented); \
2848 } \
2849 break; }
2850
2851
2852
2853 void matrix< bigint >::
mgcd(matrix<bigint> & T,lidia_size_t mgcdModul)2854 mgcd(matrix< bigint > & T, lidia_size_t mgcdModul)
2855 {
2856 lidia_size_t i = rows, j = columns;
2857
2858 if (T.columns != columns)
2859 T.set_no_of_columns(columns);
2860 if (T.rows != columns)
2861 T.set_no_of_rows(columns);
2862
2863 T.diag(1, 0);
2864
2865 switch(mgcdModul) {
2866 case 0:
2867 schema_mgcd_call_sequenz_ex(hermite)
2868 case 1:
2869 schema_mgcd_call_sequenz_ex(Bradley)
2870 case 2:
2871 schema_mgcd_call_sequenz_ex(Iliopoulos)
2872 case 3:
2873 schema_mgcd_call_sequenz_ex(opt)
2874 case 4:
2875 schema_mgcd_call_sequenz_ex(Blankinship)
2876 case 5:
2877 schema_mgcd_call_sequenz_ex(Best_remainder)
2878 case 6:
2879 schema_mgcd_call_sequenz_ex(havas_best_remainder)
2880 case 7:
2881 schema_mgcd_call_sequenz_ex(havas_euclidean_norm)
2882 case 8:
2883 schema_mgcd_call_sequenz_ex(havas_minimal_norm)
2884 case 9:
2885 schema_mgcd_call_sequenz_ex(havas_sorting_gcd)
2886 case 10:
2887 schema_mgcd_call_sequenz_ex(havas_min_no_of_elements)
2888 case 11:
2889 schema_mgcd_call_sequenz_ex(havas_min_abs_of_row_plus_minimal_euclidean_norm)
2890 case 12:
2891 schema_mgcd_call_sequenz_ex(havas_min_abs_of_row_plus_minimal_norm)
2892 case 13:
2893 schema_mgcd_call_sequenz_ex(havas_min_abs_of_row_plus_min_no_of_elements)
2894 case 14:
2895 schema_mgcd_call_sequenz_ex(havas_minimal_norm_plus_min_abs_of_row)
2896 case 15:
2897 schema_mgcd_call_sequenz_ex(havas_minimal_norm_plus_sorting_gcd)
2898 case 16:
2899 schema_mgcd_call_sequenz_ex(havas_minimal_norm_plus_min_no_of_elements)
2900 case 17:
2901 schema_mgcd_call_sequenz_ex(havas_sorting_gcd_plus_minimal_euclidean_norm)
2902 case 18:
2903 schema_mgcd_call_sequenz_ex(havas_sorting_gcd_plus_minimal_norm)
2904 case 19:
2905 schema_mgcd_call_sequenz_ex(havas_sorting_gcd_plus_min_no_of_elements)
2906 case 20:
2907 schema_mgcd_call_sequenz_ex(havas_min_no_of_elements_plus_min_abs_of_row)
2908 case 21:
2909 schema_mgcd_call_sequenz_ex(havas_min_no_of_elements_plus_minimal_euclidean_norm)
2910 case 22:
2911 schema_mgcd_call_sequenz_ex(havas_min_no_of_elements_plus_minimal_norm)
2912 case 23:
2913 schema_mgcd_call_sequenz_ex(havas_min_no_of_elements_plus_sorting_gcd)
2914 case 24:
2915 schema_mgcd_call_sequenz_ex(Storjohann)
2916 case 25:
2917 schema_mgcd_call_sequenz_ex(Heuristik)
2918 default:
2919 lidia_error_handler("bigint_matrix", "mgcd :: Mode not supported! (mgcd Modul not defined)");
2920 }
2921 }
2922
2923
2924
2925 void matrix< bigint >::
normalize(lidia_size_t normalizeModul)2926 normalize(lidia_size_t normalizeModul)
2927 {
2928 //
2929 // Task: Normalization
2930 // Version: 2.0
2931 //
2932
2933 debug_handler_l(DMESSAGE, "normalize(lidia_size_t)", DVALUE + 5);
2934
2935 if (bitfield.get_representation() == matrix_flags::dense_representation) {
2936 normalization_kernel< bigint, RODMM< bigint >,
2937 RODMM< bigint > > nmodul;
2938
2939 switch(normalizeModul) {
2940 case 0:
2941 {
2942 nmodul.normalize_Std(*this, 0, columns - rows);
2943 break;
2944 }
2945 case 1:
2946 {
2947 nmodul.normalize_ChouCollins(*this, 0, columns - rows);
2948 break;
2949 }
2950 case 2:
2951 {
2952 nmodul.normalizeMod_Std(*this, 0, columns - rows);
2953 break;
2954 }
2955 case 3:
2956 {
2957 nmodul.normalizeMod_ChouCollins(*this, 0, columns - rows);
2958 break;
2959 }
2960 case 4:
2961 {
2962 nmodul.normalizeHybrid_Std(*this, 0, columns - rows);
2963 break;
2964 }
2965 case 5:
2966 {
2967 nmodul.normalizeHybrid_ChouCollins(*this, 0, columns - rows);
2968 break;
2969 }
2970 default:
2971 lidia_error_handler("bigint_matrix", "normalize :: "
2972 "Mode not supported! (normalize "
2973 "Modul not defined)");
2974 }
2975 }
2976 else {
2977 set_orientation(matrix_flags::column_oriented);
2978 normalization_kernel< bigint, COSMM< bigint >,
2979 COSMM< bigint > > nmodul;
2980
2981 switch(normalizeModul) {
2982 case 0:
2983 {
2984 nmodul.normalize_Std(*this, 0, columns - rows);
2985 break;
2986 }
2987 case 1:
2988 {
2989 nmodul.normalize_ChouCollins(*this, 0, columns - rows);
2990 break;
2991 }
2992 case 2:
2993 {
2994 nmodul.normalizeMod_Std(*this, 0, columns - rows);
2995 break;
2996 }
2997 case 3:
2998 {
2999 nmodul.normalizeMod_ChouCollins(*this, 0, columns - rows);
3000 break;
3001 }
3002 case 4:
3003 {
3004 nmodul.normalizeHybrid_Std(*this, 0, columns - rows);
3005 break;
3006 }
3007 case 5:
3008 {
3009 nmodul.normalizeHybrid_ChouCollins(*this, 0, columns - rows);
3010 break;
3011 }
3012 default:
3013 lidia_error_handler("bigint_matrix", "normalize :: "
3014 "Mode not supported! (normalize "
3015 "Modul not defined)");
3016 }
3017 set_orientation(matrix_flags::row_oriented);
3018 }
3019 }
3020
3021
3022
3023 //
3024 // hnf_kannan
3025 //
3026
3027 void matrix< bigint >::
hnf_kannan(lidia_size_t SW)3028 hnf_kannan(lidia_size_t SW)
3029 {
3030 //
3031 // Task: HNF Computation
3032 // Algorithm: Gauss with reduction
3033 // IMPROVEMENTS: Algorithm of Kannan and Bachem
3034 // Version: 2.0
3035 //
3036
3037 debug_handler_l(DMESSAGE, "hnf_kannan()", DVALUE + 5);
3038
3039 lidia_size_t i = rows, j = columns;
3040
3041 switch(SW) {
3042 case 0:
3043 {
3044 if (bitfield.get_representation() == matrix_flags::dense_representation) {
3045 kannan_kernel< bigint, RODMM< bigint >, RODMM< bigint >,
3046 Standard_normalization< bigint, RODMM< bigint >, RODMM< bigint > > > modul;
3047
3048 modul.hnf(*this, i, j);
3049 }
3050 else {
3051 kannan_kernel< bigint, COSMM< bigint >, COSMM< bigint >,
3052 Standard_normalization< bigint, COSMM< bigint >, COSMM< bigint > > > modul;
3053
3054 set_orientation(matrix_flags::column_oriented);
3055 modul.hnf(*this, i, j);
3056 set_orientation(matrix_flags::row_oriented);
3057 }
3058 break;
3059 }
3060 case 1:
3061 {
3062 if (bitfield.get_representation() == matrix_flags::dense_representation) {
3063 kannan_kernel< bigint, RODMM< bigint >, RODMM< bigint >,
3064 ChouCollins_normalization< bigint, RODMM< bigint >, RODMM< bigint > > > modul;
3065
3066 modul.hnf(*this, i, j);
3067 }
3068 else {
3069 kannan_kernel< bigint, COSMM< bigint >, COSMM< bigint >,
3070 ChouCollins_normalization< bigint, COSMM< bigint >, COSMM< bigint > > > modul;
3071
3072 set_orientation(matrix_flags::column_oriented);
3073 modul.hnf(*this, i, j);
3074 set_orientation(matrix_flags::row_oriented);
3075 }
3076 break;
3077 }
3078 case 2:
3079 {
3080 if (bitfield.get_representation() == matrix_flags::dense_representation) {
3081 kannan_kernel< bigint, RODMM< bigint >, RODMM< bigint >,
3082 Jacobson_normalization< bigint, RODMM< bigint >, RODMM< bigint > > > modul;
3083
3084 modul.hnf(*this, i, j);
3085 }
3086 else {
3087 kannan_kernel< bigint, COSMM< bigint >, COSMM< bigint >,
3088 Jacobson_normalization< bigint, COSMM< bigint >, COSMM< bigint > > > modul;
3089
3090 set_orientation(matrix_flags::column_oriented);
3091 modul.hnf(*this, i, j);
3092 set_orientation(matrix_flags::row_oriented);
3093 }
3094 break;
3095 }
3096 default:
3097 lidia_error_handler("bigint_matrix", "hnf_kannan :: Mode not supported!");
3098 }
3099 }
3100
3101
3102
3103 //
3104 // hnf_cg
3105 //
3106
3107 void matrix< bigint >::
hnf_cg(const matrix<long> & B,long BOUND_1,const bigint & BOUND_2,int no)3108 hnf_cg(const matrix< long > &B, long BOUND_1, const bigint &BOUND_2, int no)
3109 {
3110 //
3111 // modul definitions
3112 //
3113
3114 COSMM< long > modul3;
3115
3116 havas_kernel< long, COSMM< long >, COSMM< bigint >,
3117 nf_conf3e< long, COSMM< long >, COSMM< bigint > > > hnf_smodul3le;
3118
3119 normalization_kernel< long, COSMM< long >, COSMM< bigint > > normalize_smodul3le;
3120
3121 havas_kernel< bigint, RODMM< bigint >, RODMM< bigint >,
3122 havas_best_remainder_ext< bigint, RODMM< bigint >, RODMM< bigint > > > hnf_dmodul2e;
3123
3124 havas_kernel< bigint, COSMM< bigint >, COSMM< bigint >,
3125 havas_best_remainder< bigint, COSMM< bigint >, COSMM< bigint > > > hnf_smodul2;
3126
3127 normalization_kernel< bigint, COSMM< bigint >, COSMM< bigint > > normalize_smodul2;
3128
3129 const modular_arithmetic< DRMK< bigint >,
3130 dense_fp_matrix_kernel< long, MR< long > >,
3131 dense_fp_matrix_kernel< bigint, MR< bigint > > > Dm_bigint_modul;
3132
3133 const modular_arithmetic< SRMK< bigint >,
3134 sparse_fp_matrix_kernel< long, MR< long > >,
3135 sparse_fp_matrix_kernel< bigint, MR< bigint > > > Sm_bigint_modul;
3136
3137 //
3138 // Variables
3139 //
3140
3141 timer t;
3142 bigint H, DET;
3143 lidia_size_t actual_row = B.rows, actual_column = B.columns;
3144 lidia_size_t i, j;
3145
3146 //
3147 // main algorithm
3148 //
3149
3150 lidia_info_handler(t.start_timer(););
3151
3152 // Change representation of matrix B
3153 set_storage_mode(matrix_flags::sparse_representation);
3154
3155 // Change dimensions of member matrix
3156 resize(B.rows, B.columns);
3157
3158 //
3159 // Stage 1: non-modular long
3160 //
3161
3162 {
3163 // Create a copy of the original matrix
3164 matrix< long > A = B;
3165 A.set_zero_element(0);
3166
3167 // change orientation
3168 A.set_orientation(matrix_flags::column_oriented);
3169
3170 // Compute the Hadamard Bound
3171 modul3.hadamard(A, H);
3172 lidia_info_handler(t.stop_timer();
3173 std::cout << std::endl << "Hadamard Bound = " << H << std::endl;
3174 std::cout << "Time: " << t << std::endl;
3175 t.cont_timer(););
3176
3177 // STF computation
3178 if (hnf_smodul3le.hnf_Z2(A, actual_row, actual_column, BOUND_1)) {
3179 normalize_smodul3le.normalize_ChouCollins(A, 0, actual_column);
3180
3181 set_orientation(matrix_flags::row_oriented);
3182
3183 // copy data
3184 for (i = 0; i < A.columns; i++)
3185 for (j = 0; j < A.value_counter[i]; j++)
3186 sto(A.index[i][j], i, bigint(A.value[i][j]));
3187
3188 lidia_info_handler(t.stop_timer();
3189 std::cout << std::endl << "total: " << t << std::endl;);
3190 return;
3191 }
3192
3193 // copy data
3194 set_orientation(matrix_flags::column_oriented);
3195 for (i = 0; i < A.columns; i++) {
3196 lidia_size_t len = A.value_counter[i];
3197 value[i] = new bigint[len];
3198 index[i] = new lidia_size_t[len];
3199 for (j = 0; j < len; j++) {
3200 value[i][j] = bigint(A.value[i][j]);
3201 index[i][j] = A.index[i][j];
3202 }
3203 value_counter[i] = len;
3204 allocated[i] = len;
3205 }
3206 }
3207
3208 //
3209 // Stage 2: non-modular bigint
3210 //
3211
3212 lidia_info_handler(t.stop_timer();
3213 std::cout << std::endl << "hnf_Z2 (long): " << t << std::endl;
3214 t.cont_timer(););
3215
3216 {
3217 matrix< bigint > A3(actual_row, actual_column);
3218
3219 // insert elements
3220 for (i = 0; i < A3.columns; i++) {
3221 for (j = 0; j < value_counter[i]; j++)
3222 A3.value[index[i][j]][i] = value[i][j];
3223 delete[] value[i];
3224 delete[] index[i];
3225 value_counter[i] = 0;
3226 }
3227
3228 if (!hnf_dmodul2e.hnf_Z2(A3, actual_row, actual_column, BOUND_2)) {
3229 lidia_info_handler(t.stop_timer();
3230 std::cout << std::endl << "hnf_Z2 (bigint): " << t << std::endl;
3231 t.cont_timer(););
3232
3233 Dm_bigint_modul.latticedet2(A3, DET, H, no);
3234 D_bigint_modul.hnfmod_dkt_part(A3, DET);
3235 }
3236
3237 // insert elements back
3238 for (i = 0; i < A3.columns; i++) {
3239 value[i] = new bigint[A3.rows];
3240 index[i] = new lidia_size_t[A3.rows];
3241 allocated[i] = A3.rows;
3242
3243 lidia_size_t k = 0;
3244 for (j = 0; j < A3.rows; j++) {
3245 if (A3.value[j][i] != A3.Zero) {
3246 value[i][k] = A3.value[j][i];
3247 index[i][k] = j;
3248 k++;
3249 }
3250 }
3251 value_counter[i] = k;
3252 }
3253 }
3254
3255 // normalize
3256 normalize_smodul2.normalize_ChouCollins(*this, 0, B.columns - B.rows);
3257
3258 set_orientation(matrix_flags::row_oriented);
3259 lidia_info_handler(t.stop_timer();
3260 std::cout << std::endl << "total: " << t << std::endl;);
3261 }
3262
3263
3264
3265 void matrix< bigint >::
hnf_cg(const matrix<long> & B,matrix<bigint> & TR,long BOUND_1,const bigint & BOUND_2,int no)3266 hnf_cg(const matrix< long > &B, matrix< bigint > &TR, long BOUND_1,
3267 const bigint &BOUND_2, int no)
3268 {
3269 //
3270 // modul definitions
3271 //
3272
3273 COSMM< long > modul3;
3274 COSMM< bigint > modul4;
3275
3276 const modular_arithmetic< DRMK< bigint >,
3277 dense_fp_matrix_kernel< long, MR< long > >,
3278 dense_fp_matrix_kernel< bigint, MR< bigint > > > Dm_bigint_modul;
3279
3280 const modular_arithmetic< SRMK< bigint >,
3281 sparse_fp_matrix_kernel< long, MR< long > >,
3282 sparse_fp_matrix_kernel< bigint, MR< bigint > > > Sm_bigint_modul;
3283
3284
3285
3286 TR.set_representation(matrix_flags::sparse_representation);
3287
3288 if (TR.columns != TR.rows || TR.columns != B.columns)
3289 TR.resize(B.columns, B.columns);
3290
3291 TR.diag(1, 0);
3292
3293 timer t;
3294 lidia_info_handler(t.start_timer(););
3295
3296 lidia_size_t actual_row, actual_column;
3297
3298 bigint H;
3299
3300 set_representation(matrix_flags::sparse_representation);
3301 resize(B.rows, B.columns);
3302
3303 lidia_size_t i, j;
3304 //
3305 // Stage 1: non-modular long
3306 //
3307 {
3308 matrix< long > A;
3309 A.set_zero_element(0);
3310
3311 // Create a copy of the original matrix
3312 A = B;
3313
3314 // change orientation
3315 A.set_orientation(matrix_flags::column_oriented);
3316 TR.set_orientation(matrix_flags::column_oriented);
3317
3318 modul3.hadamard(A, H);
3319 lidia_info_handler(t.stop_timer();
3320 std::cout << std::endl << "Hadamard Bound = " << H << std::endl;
3321 std::cout << "Time: " << t << std::endl;
3322 t.cont_timer(););
3323
3324 havas_kernel< long, COSMM< long >, COSMM< bigint >,
3325 nf_conf3e< long, COSMM< long >, COSMM< bigint > > > hnf_smodul3le;
3326
3327 normalization_kernel< long, COSMM< long >, COSMM< bigint > > normalize_smodul3le;
3328
3329 if (hnf_smodul3le.hnf_Z2(A, TR, actual_row, actual_column, BOUND_1)) {
3330 normalize_smodul3le.normalize_Std(A, TR, 0, actual_column);
3331
3332 set_orientation(matrix_flags::row_oriented);
3333 TR.set_orientation(matrix_flags::row_oriented);
3334
3335 // copy data
3336 for (i = 0; i < A.columns; i++)
3337 for (j = 0; j < A.value_counter[i]; j++)
3338 sto(A.index[i][j], i, bigint(A.value[i][j]));
3339
3340 lidia_info_handler(t.stop_timer();
3341 std::cout << std::endl << "total (Phase 1): " << t << std::endl;);
3342 return;
3343 }
3344
3345 // copying matrix
3346 set_orientation(matrix_flags::column_oriented);
3347
3348 for (i = 0; i < A.columns; i++) {
3349 lidia_size_t len = A.value_counter[i];
3350 value[i] = new bigint[len];
3351 index[i] = new lidia_size_t[len];
3352 for (j = 0; j < len; j++) {
3353 value[i][j] = bigint(A.value[i][j]);
3354 index[i][j] = A.index[i][j];
3355 }
3356 value_counter[i] = len;
3357 allocated[i] = len;
3358 }
3359 }
3360
3361 //
3362 // Stage 2: non-modular bigint
3363 //
3364
3365 lidia_info_handler(t.stop_timer();
3366 std::cout << std::endl << "hnf_Z2 (long): " << t << std::endl;
3367 t.cont_timer(););
3368
3369 {
3370 matrix< bigint > A3(actual_row, actual_column);
3371
3372 // insert elements
3373 for (i = 0; i < A3.columns; i++) {
3374 for (j = 0; j < value_counter[i]; j++)
3375 A3.value[index[i][j]][i] = value[i][j];
3376 delete[] value[i];
3377 delete[] index[i];
3378 value_counter[i] = 0;
3379 }
3380
3381 havas_kernel< bigint, RODMM< bigint >, RODMM< bigint >,
3382 havas_best_remainder_ext< bigint, RODMM< bigint >, RODMM< bigint > > > hnf_dmodul2e;
3383
3384 if (!hnf_dmodul2e.hnf_Z2(A3, actual_row, actual_column, BOUND_2)) {
3385 lidia_info_handler(t.stop_timer();
3386 std::cout << std::endl << "hnf_Z2 (bigint): " << t << std::endl;
3387 t.cont_timer(););
3388 bigint DET;
3389 Dm_bigint_modul.latticedet2(A3, DET, H, no);
3390 D_bigint_modul.hnfmod_dkt_part(A3, DET);
3391 }
3392
3393 // insert elements back
3394 for (i = 0; i < A3.columns; i++) {
3395 value[i] = new bigint[A3.rows];
3396 index[i] = new lidia_size_t[A3.rows];
3397 allocated[i] = A3.rows;
3398
3399 lidia_size_t k = 0;
3400 for (j = 0; j < A3.rows; j++) {
3401 if (A3.value[j][i] != A3.Zero) {
3402 value[i][k] = A3.value[j][i];
3403 index[i][k] = j;
3404 k++;
3405 }
3406 }
3407 value_counter[i] = k;
3408 }
3409 }
3410
3411 // normalize
3412 normalization_kernel< bigint, COSMM< bigint >, COSMM< bigint > > normalize_smodul2;
3413
3414 normalize_smodul2.normalize_Std(*this, 0, B.columns - B.rows);
3415
3416 set_orientation(matrix_flags::row_oriented);
3417 lidia_info_handler(t.stop_timer();
3418 std::cout << std::endl << "total: " << t << std::endl;);
3419 }
3420
3421
3422
3423 //
3424 // hnf_gls_solver
3425 //
3426
3427 void matrix< bigint >::
hnf_gls_solver()3428 hnf_gls_solver()
3429 {
3430 matrix_flags flags;
3431 flags.set_representation(matrix_flags::sparse_representation);
3432
3433 //
3434 // temp variables
3435 //
3436
3437 matrix< bigint > C2(columns, columns, flags);
3438 matrix< bigint > T2(columns, columns, flags);
3439
3440 bool SW = false;
3441 bigint DET, DET2;
3442 bigint H, H1, d;
3443
3444 matrix< bigint > B, D, B2;
3445 B.set_representation(matrix_flags::sparse_representation);
3446
3447 matrix< bigint > CT(rows, rows, flags);
3448 matrix< bigint > TT(rows, rows, flags);
3449
3450 lidia_size_t pos = 0;
3451
3452 math_vector< bigint > b(rows - 1, rows - 1), x(rows, rows);
3453
3454 //
3455 // computing the hadamard bound
3456 //
3457
3458 hadamard(H);
3459 std::cout << "hadamard's bound: " << H << "(" << decimal_length(H) << " digits)" << std::endl;
3460
3461 //
3462 // computation of the lattice determinante
3463 //
3464
3465 det(DET, H);
3466 std::cout << "determinante: " << DET << "(" << decimal_length(DET) << " digits)" << std::endl;
3467
3468 DET2 = DET;
3469 lidia_size_t i, j, l;
3470 while (abs(DET) != 1) {
3471 SW = false;
3472 std::cout << "Rest determinante (" << pos << ") :" << DET << std::endl;
3473
3474 C2.diag(1, 0);
3475 T2.diag(1, 0);
3476
3477 //
3478 // creating linear system
3479 //
3480
3481 b.set_size(rows - pos - 1);
3482 for (i = pos + 1; i < rows; i++) {
3483 b[i - pos - 1] = member(i, pos);
3484 if (b[i - pos - 1] != 0)
3485 SW = true;
3486 }
3487
3488 if (SW) {
3489 B.resize(rows - pos - 1, columns - pos - 1);
3490 for (j = pos + 1; j < rows; j++)
3491 for (i = pos + 1; i < columns; i++)
3492 B.sto(i - pos - 1, j - pos - 1, member(i, j));
3493
3494 //
3495 // new system
3496 //
3497
3498 //D.trans(B);
3499
3500 //B.det(DET2, H);
3501 //std::cout << "partial determinante (sparse): " << DET2 << std::endl;
3502
3503 x.set_size(rows - pos);
3504
3505 //
3506 // size reduction 1
3507 //
3508
3509 //d = b[0];
3510 // for (i = 1; i < rows - pos - 1; i++)
3511 //d = gcd(d, b[i]);
3512
3513 //for (i = 0; i < rows - pos - 1; i++)
3514 //b[i] /= d;
3515 // std::cout << "gcd computation (b) " << std::endl;
3516 // std::cout << "gcd = " << d << std::endl;
3517
3518 //B = D * B;
3519 //b = D * b;
3520
3521 //
3522 // system solving
3523 // PART 1: wiedemann
3524 // PART 2: lanczos
3525 //
3526
3527 std::cout << "System solving:" << std::endl;
3528 B.hadamard(H1);
3529
3530 B.reginvimage(x, b, H1*H1*H, DET2);
3531 std::cout << "Loesung: " << x << std::endl;
3532
3533 //
3534 // size reduction 2
3535 //
3536
3537 d = x[0];
3538 for (i = 1; i < rows - pos; i++)
3539 d = gcd(d, x[i]);
3540 std::cout << "gcd computation" << std::endl;
3541 std::cout << "gcd = " << d << std::endl;
3542
3543 bigint *x2 = new bigint[rows - pos];
3544 for (i = 0; i < rows - pos; i++)
3545 x2[i] = x[i] / d;
3546
3547 //
3548 // conditioning matrix
3549 //
3550
3551 CT.resize(rows - pos, rows - pos);
3552 CT.cond_matrix(x2, rows - pos);
3553
3554 x2 = CT * x2;
3555
3556
3557 //
3558 // transformation matrix
3559 //
3560
3561 TT.resize(rows - pos, rows - pos);
3562 TT.diag(1, 0);
3563
3564 bigint TMP1, TMP2, TMP3;
3565 TMP3 = xgcd(TMP1, TMP2, x2[0], x2[1]);
3566 for (i = 0; i < rows - pos; i++)
3567 TT.sto(i, 0, x2[i]);
3568 TT.sto(0, 1, -TMP2);
3569 TT.sto(1, 1, TMP1);
3570 delete[] x2;
3571
3572 //
3573 // update of A and TR
3574 //
3575
3576 for (i = 2; i < rows - pos; i++)
3577 CT.sto(1, i, -CT.member(1, i));
3578 C2.insert_at(pos, pos, CT, 0, 0, CT.get_no_of_rows(), CT.get_no_of_columns());
3579 T2.insert_at(pos, pos, TT, 0, 0, TT.get_no_of_rows(), TT.get_no_of_columns());
3580 multiply(*this, C2);
3581 multiply(*this, T2);
3582 }
3583
3584 std::cout << "OLDdet = " << DET << std::endl;
3585 if (member(pos, pos) < 0) {
3586 for (l = 0; l < rows; l++)
3587 sto(l, pos, -member(l, pos));
3588 }
3589
3590 std::cout << "Element: " << member(pos, pos) << std::endl;
3591 DET = DET / member(pos, pos);
3592
3593 //
3594 // final reduction
3595 //
3596
3597 std::cout << "NORMALIZE: " << std::endl;
3598 bigint TMPQ1, TMPQ2;
3599 for (i = 0; i < pos; i++)
3600 for (j = i + 1; j < pos; j++) {
3601 pos_div_rem(TMPQ1, TMPQ2, member(i, j), member(i, i));
3602 for (l = 0; l < rows; l++)
3603 sto(l, j, member(l, j) - TMPQ1*member(l, i));
3604 }
3605 pos++;
3606 }
3607 }
3608
3609
3610
3611 void matrix< bigint >::
hnf_gls_solver(matrix<bigint> &)3612 hnf_gls_solver(matrix< bigint > &)
3613 {
3614
3615 matrix_flags flags;
3616 flags.set_representation(matrix_flags::sparse_representation);
3617
3618 //
3619 // temp variables
3620 //
3621
3622 matrix< bigint > C(columns, columns, flags);
3623 matrix< bigint > C2(columns, columns, flags);
3624 matrix< bigint > T(columns, columns, flags);
3625 matrix< bigint > T2(columns, columns, flags);
3626
3627 C.diag(1, 0);
3628 T.diag(1, 0);
3629
3630 bool SW = false;
3631 bigint DET, DET2;
3632 bigint H, H1, d;
3633
3634 matrix< bigint > B, D, B2;
3635 B.set_representation(matrix_flags::sparse_representation);
3636
3637 matrix< bigint > CT(rows, rows, flags);
3638 matrix< bigint > TT(rows, rows, flags);
3639
3640 lidia_size_t pos = 0;
3641
3642 math_vector< bigint > b(rows - 1, rows - 1), x(rows, rows);
3643
3644 //
3645 // computing the hadamard bound
3646 //
3647
3648 hadamard(H);
3649 std::cout << "hadamard's bound: " << H << "(" << decimal_length(H) << " digits)" << std::endl;
3650
3651 //
3652 // computation of the lattice determinante
3653 //
3654
3655 det(DET, H);
3656 std::cout << "determinante: " << DET << "(" << decimal_length(DET) << " digits)" << std::endl;
3657
3658 lidia_size_t i, j, l;
3659 while (abs(DET) != 1) {
3660 SW = false;
3661 std::cout << "rest determinante (" << pos << ") :" << DET << std::endl;
3662
3663 C2.diag(1, 0);
3664 T2.diag(1, 0);
3665
3666 //
3667 // creating linear system
3668 //
3669
3670 b.set_size(rows - pos - 1);
3671 for (i = pos + 1; i < rows; i++) {
3672 b[i - pos - 1] = member(i, pos);
3673 if (b[i - pos - 1] != 0)
3674 SW = true;
3675 }
3676
3677 if (SW) {
3678 B.resize(rows - pos - 1, columns - pos - 1);
3679 for (j = pos + 1; j < rows; j++)
3680 for (i = pos + 1; i < columns; i++)
3681 B.sto(i - pos - 1, j - pos - 1, member(i, j));
3682
3683 //
3684 // new system
3685 //
3686
3687 D.trans(B);
3688
3689 B.det(DET2, H);
3690 std::cout << "partial determinante (sparse): " << DET2 << std::endl;
3691
3692 x.set_size(rows - pos);
3693
3694 //
3695 // size reduction 1
3696 //
3697
3698 d = b[0];
3699 for (i = 1; i < rows - pos - 1; i++)
3700 d = gcd(d, b[i]);
3701
3702 for (i = 0; i < rows - pos - 1; i++)
3703 b[i] /= d;
3704 std::cout << "gcd computation (b) " << std::endl;
3705 std::cout << "gcd = " << d << std::endl;
3706
3707 B = D * B;
3708 b = D * b;
3709
3710 //
3711 // system solving
3712 // PART 1: wiedemann
3713 // PART 2: lanczos
3714 //
3715
3716 std::cout << "System solving:" << std::endl;
3717 B.hadamard(H1);
3718
3719 B.reginvimage(x, b, H1*H1*H, DET2);
3720 std::cout << "Loesung: " << x << std::endl;
3721
3722 //
3723 // size reduction 2
3724 //
3725
3726 d = x[0];
3727 for (i = 1; i < rows - pos; i++)
3728 d = gcd(d, x[i]);
3729 std::cout << "gcd computation" << std::endl;
3730 std::cout << "gcd = " << d << std::endl;
3731
3732 bigint *x2 = new bigint[rows - pos];
3733 for (i = 0; i < rows - pos; i++)
3734 x2[i] = x[i] / d;
3735
3736 //
3737 // conditioning matrix
3738 //
3739
3740 CT.resize(rows - pos, rows - pos);
3741 CT.cond_matrix(x2, rows - pos);
3742
3743 x2 = CT * x2;
3744
3745
3746 //
3747 // transformation matrix
3748 //
3749
3750 TT.resize(rows - pos, rows - pos);
3751 TT.diag(1, 0);
3752
3753 bigint TMP1, TMP2, TMP3;
3754 TMP3 = xgcd(TMP1, TMP2, x2[0], x2[1]);
3755 for (i = 0; i < rows - pos; i++)
3756 TT.sto(i, 0, x2[i]);
3757 TT.sto(0, 1, -TMP2);
3758 TT.sto(1, 1, TMP1);
3759 delete[] x2;
3760
3761 //
3762 // update of A and TR
3763 //
3764
3765 for (i = 2; i < rows - pos; i++)
3766 CT.sto(1, i, -CT.member(1, i));
3767 C2.insert_at(pos, pos, CT, 0, 0, CT.get_no_of_rows(), CT.get_no_of_columns());
3768 T2.insert_at(pos, pos, TT, 0, 0, TT.get_no_of_rows(), TT.get_no_of_columns());
3769 multiply(*this, C2);
3770 multiply(*this, T2);
3771 }
3772
3773 std::cout << "OLDdet = " << DET << std::endl;
3774 if (member(pos, pos) < 0) {
3775 for (l = 0; l < rows; l++)
3776 sto(l, pos, -member(l, pos));
3777 }
3778
3779 std::cout << "Element: " << member(pos, pos) << std::endl;
3780 DET = DET / member(pos, pos);
3781
3782 //
3783 // final reduction
3784 //
3785
3786 std::cout << "NORMALIZE: " << std::endl;
3787 bigint TMPQ1, TMPQ2;
3788 for (i = 0; i < pos; i++)
3789 for (j = i + 1; j < pos; j++) {
3790 pos_div_rem(TMPQ1, TMPQ2, member(i, j), member(i, i));
3791 for (l = 0; l < rows; l++)
3792 sto(l, j, member(l, j) - TMPQ1*member(l, i));
3793 }
3794 pos++;
3795 }
3796 }
3797
3798
3799
3800 //
3801 // Kernel
3802 //
3803
3804 void matrix< bigint >::
kernel1(const matrix<bigint> & A)3805 kernel1(const matrix< bigint > & A)
3806 {
3807 //
3808 // Task: B.kernel1(A);
3809 // => The columns of matrix B form a basis
3810 // of the kernel of matrix A.
3811 // Version: 2.0
3812 //
3813
3814 debug_handler_l(DMESSAGE, "kernel1(const matrix< bigint > &)", DVALUE + 5);
3815
3816 const modular_bigint_matrix_algorithms< DRMK< bigint >,
3817 dense_fp_matrix_kernel< long, MR< long > >,
3818 dense_fp_matrix_kernel< bigint, MR< bigint > > > Dm_bigint_modul;
3819
3820 const modular_bigint_matrix_algorithms< SRMK< bigint >,
3821 sparse_fp_matrix_kernel< long, MR< long > >,
3822 sparse_fp_matrix_kernel< bigint, MR< bigint > > > Sm_bigint_modul;
3823
3824 bigint H;
3825
3826 if (A.bitfield.get_representation() == matrix_flags::dense_representation) {
3827 D_bigint_modul.hadamard(A, H);
3828 Dm_bigint_modul.kernel1(*this, A, H);
3829 }
3830 else {
3831 S_bigint_modul.hadamard(A, H);
3832 Sm_bigint_modul.kernel1(*this, A, H);
3833 }
3834
3835 }
3836
3837
3838
3839 void matrix< bigint >::
kernel2(const matrix<bigint> & A)3840 kernel2(const matrix< bigint > & A)
3841 {
3842 //
3843 // Task: B.kernel2(A);
3844 // => The columns of matrix B form a basis
3845 // of the kernel of matrix A.
3846 // Version: 2.0
3847 //
3848
3849 debug_handler_l(DMESSAGE, "kernel2(const matrix< bigint > &)", DVALUE + 5);
3850
3851 D_bigint_modul.kernel2(*this, A);
3852 }
3853
3854
3855
3856 //
3857 // regular InvImage
3858 //
3859
3860 void matrix< bigint >::
reginvimage(math_vector<bigint> & RES,const math_vector<bigint> & b) const3861 reginvimage(math_vector< bigint > &RES, const math_vector< bigint > & b) const
3862 {
3863 bigint TMP, H = hadamard();
3864
3865 sparse_fp_matrix_algorithms< long, MR< long > > long_modul;
3866 sparse_fp_matrix_algorithms< bigint, MR< bigint > > bigint_modul;
3867
3868 // Step 2
3869 long len;
3870 shift_left(TMP, H, 1);
3871 bigint *PRIM = get_primes(TMP, bigint(1));
3872 PRIM[0].longify(len);
3873
3874 // Step 3
3875 matrix< bigint > U(columns + 1, static_cast<int>(len));
3876
3877 bigint MOD;
3878 long Modlong;
3879
3880 // Step 3
3881 lidia_size_t i, j;
3882 for (i = 1; i <= len; i++) {
3883 MOD.assign(PRIM[i]);
3884 if (MOD.bit_length() > bigint::bits_per_digit()) {
3885 matrix< bigint > B(rows, columns);
3886 remainder(B, *this, MOD);
3887 math_vector< bigint > x(columns, columns);
3888
3889 bigint DET = bigint_modul.det(B, MOD);
3890
3891 math_vector< bigint > b1(rows, rows);
3892 for (j = 0; j < rows; j++) {
3893 LiDIA::remainder(b1[j], b[j], MOD);
3894 LiDIA::mult_mod(b1[j], b1[j], DET, MOD);
3895 }
3896
3897 bigint_modul.lanczos(B, x, b1, MOD);
3898
3899 U.value[0][i - 1] = -DET;
3900 for (j = 1; j < columns + 1; j++)
3901 U.value[j][i - 1].assign(x[j-1]);
3902 }
3903 else {
3904 MOD.longify(Modlong);
3905 matrix< long > B(rows, columns);
3906 B.set_zero_element(0);
3907 remainder(B, *this, Modlong);
3908 math_vector< long > x(columns, columns);
3909
3910 long DET = long_modul.det(B, Modlong);
3911 math_vector< long > b1(rows, rows);
3912 for (j = 0; j < rows; j++) {
3913 LiDIA::remainder(b1[j], b[j], Modlong);
3914 LiDIA::mult_mod(b1[j], b1[j], DET, Modlong);
3915 }
3916
3917 long_modul.lanczos(B, x, b1, Modlong);
3918
3919 U.value[0][i - 1] = -DET;
3920 for (j = 1; j < columns + 1; j++)
3921 U.value[j][i - 1].assign(x[j-1]);
3922 }
3923 }
3924
3925 // Step 4,5
3926 RES.set_size(U.rows);
3927 for (i = 0; i < columns + 1; i++)
3928 LiDIA::chinrest(RES[i], U.value[i], PRIM);
3929 }
3930
3931
3932
3933 void matrix< bigint >::
reginvimage(math_vector<bigint> & RES,const math_vector<bigint> & b,const bigint & H,const bigint & DET) const3934 reginvimage(math_vector< bigint > &RES, const math_vector< bigint > & b,
3935 const bigint &H, const bigint &DET) const
3936 {
3937 bigint TMP; //, H = hadamard();
3938
3939 sparse_fp_matrix_algorithms< long, MR< long > > long_modul;
3940 sparse_fp_matrix_algorithms< bigint, MR< bigint > > bigint_modul;
3941
3942 // Step 2
3943 long len;
3944 shift_left(TMP, H, 1);
3945 //bigint *PRIM = get_primes(TMP, bigint(DET));
3946 bigint *PRIM = get_primes(TMP, bigint(H));
3947 PRIM[0].longify(len);
3948
3949 // Step 3
3950 bigint MOD;
3951 long Modlong;
3952
3953 math_vector< bigint > x(columns, columns);
3954 math_vector< bigint > x2(columns, columns);
3955
3956 matrix< bigint > Bbigint(rows, columns);
3957 math_vector< bigint > xbigint(columns, columns);
3958 math_vector< bigint > b1bigint(rows, rows);
3959
3960
3961 matrix< long > Blong(rows, columns);
3962 Blong.set_zero_element(0);
3963 math_vector< long > xlong(columns, columns);
3964 math_vector< long > b1long(rows, rows);
3965
3966 // Step 3
3967 bigint PROD;
3968 lidia_size_t i, j;
3969 for (i = 1; i <= len; i++) {
3970 MOD.assign(PRIM[i]);
3971 if (MOD.bit_length() > bigint::bits_per_digit()) {
3972 remainder(Bbigint, *this, MOD);
3973
3974 for (j = 0; j < rows; j++)
3975 LiDIA::mult_mod(b1bigint[j], b[j], DET, MOD);
3976
3977 lidia_info_handler(std::cout << "- > in lanczos (" << i << ", " << len << ")" << std::endl;);
3978 bigint_modul.lanczos(Bbigint, xbigint, b1bigint, MOD);
3979
3980 x2 = xbigint;
3981 }
3982 else {
3983 MOD.longify(Modlong);
3984 LiDIA::remainder(Blong, *this, Modlong);
3985
3986 long DETlong;
3987 best_remainder(DETlong, DET, Modlong);
3988 for (j = 0; j < rows; j++) {
3989 LiDIA::best_remainder(b1long[j], b[j], Modlong);
3990 LiDIA::mult_mod(b1long[j], b1long[j], DETlong, Modlong);
3991 }
3992
3993 lidia_info_handler(std::cout << "- > in lanczos (" << i << ", " <<
3994 len << ")" << std::endl;);
3995 long_modul.lanczos(Blong, xlong, b1long, Modlong);
3996
3997 for (j = 0; j < columns; j++)
3998 x2[j] = xlong[j];
3999 }
4000
4001 if (i == 1) {
4002 x = x2;
4003 PROD = MOD;
4004 }
4005 else {
4006
4007 for (j = 0; j < columns; j++) {
4008 x[j] = chinese_remainder(x[j], PROD, x2[j], MOD);
4009 best_remainder(x[j], x[j], (PROD*MOD));
4010 }
4011 PROD *= MOD;
4012 }
4013
4014 if ((*this)*x == (DET*b)) {
4015 RES[0] = -DET;
4016 for (j = 1; j <= columns; j++)
4017 RES[j] = x[j - 1];
4018 return;
4019 }
4020 }
4021 }
4022
4023
4024
4025 void matrix< bigint >::
reginvimage(math_vector<bigint> & RES,const math_vector<bigint> & b,const bigint & H,const bigint & DET,const bigint & MODULUS) const4026 reginvimage(math_vector< bigint > &RES, const math_vector< bigint > & b,
4027 const bigint &H, const bigint &DET, const bigint &MODULUS) const
4028 {
4029 bigint TMP; //, H = hadamard();
4030
4031 sparse_fp_matrix_algorithms< long, MR< long > > long_modul;
4032 sparse_fp_matrix_algorithms< bigint, MR< bigint > > bigint_modul;
4033
4034 // Step 2
4035 long len;
4036 shift_left(TMP, H, 1);
4037 //bigint *PRIM = get_primes(TMP, bigint(DET));
4038 bigint *PRIM = get_primes(TMP, bigint(H));
4039 PRIM[0].longify(len);
4040
4041 // Step 3
4042 bigint MOD;
4043 long Modlong;
4044
4045 math_vector< bigint > x(columns, columns);
4046 math_vector< bigint > x2(columns, columns);
4047
4048 matrix< bigint > Bbigint(rows, columns);
4049 math_vector< bigint > xbigint(columns, columns);
4050 math_vector< bigint > b1bigint(rows, rows);
4051
4052
4053 matrix< long > Blong(rows, columns);
4054 Blong.set_zero_element(0);
4055 math_vector< long > xlong(columns, columns);
4056 math_vector< long > b1long(rows, rows);
4057
4058 // Step 3
4059 lidia_size_t i, j;
4060 bigint PROD;
4061 for (i = 1; i <= len; i++) {
4062 MOD.assign(PRIM[i]);
4063 if (MOD.bit_length() > bigint::bits_per_digit()) {
4064 remainder(Bbigint, *this, MOD);
4065
4066 for (j = 0; j < rows; j++)
4067 LiDIA::mult_mod(b1bigint[j], b[j], DET, MOD);
4068
4069 lidia_info_handler(std::cout << "- > in lanczos (" << i << ", "
4070 << len << ")" << std::endl;);
4071 bigint_modul.lanczos(Bbigint, xbigint, b1bigint, MOD);
4072
4073 x2 = xbigint;
4074 }
4075 else {
4076 MOD.longify(Modlong);
4077 LiDIA::remainder(Blong, *this, Modlong);
4078
4079 long DETlong;
4080 best_remainder(DETlong, DET, Modlong);
4081 for (j = 0; j < rows; j++) {
4082 LiDIA::best_remainder(b1long[j], b[j], Modlong);
4083 LiDIA::mult_mod(b1long[j], b1long[j], DETlong, Modlong);
4084 }
4085
4086 lidia_info_handler(std::cout << "- > in lanczos (" << i << ", " <<
4087 len << ")" << std::endl;);
4088 long_modul.lanczos(Blong, xlong, b1long, Modlong);
4089
4090 for (j = 0; j < columns; j++)
4091 x2[j] = xlong[j];
4092 }
4093
4094 if (i == 1) {
4095 x = x2;
4096 PROD = MOD;
4097 }
4098 else {
4099
4100 for (j = 0; j < columns; j++) {
4101 x[j] = chinese_remainder(x[j], PROD, x2[j], MOD);
4102 best_remainder(x[j], x[j], (PROD*MOD));
4103 }
4104 PROD *= MOD;
4105 }
4106 bool RET = true;
4107 math_vector< bigint > c = (*this)*x - DET*b;
4108 for (j = 0; j < columns; j++)
4109 if (c[j] != 0)
4110 RET = false;
4111 if (RET) {
4112 RES[0] = -DET;
4113 for (j = 1; j <= columns; j++)
4114 RES[j] = x[j - 1];
4115 return;
4116 }
4117 }
4118 }
4119
4120
4121
4122 void matrix< bigint >::
reginvimage_ZmZ(math_vector<bigint> & RES,const math_vector<bigint> & b,const bigint & H,bigint & DET) const4123 reginvimage_ZmZ(math_vector< bigint > &RES, const math_vector< bigint > & b,
4124 const bigint &H, bigint &DET) const
4125 {
4126 bigint TMP; //, H = hadamard();
4127
4128 sparse_fp_matrix_algorithms< long, MR< long > > long_modul;
4129 sparse_fp_matrix_algorithms< bigint, MR< bigint > > bigint_modul;
4130
4131 // Step 2
4132 long len;
4133 shift_left(TMP, H, 1);
4134 bigint *PRIM = get_primes(TMP, bigint(DET));
4135 PRIM[0].longify(len);
4136
4137 // Step 3
4138 bigint MOD, DETtmp;
4139 long Modlong;
4140
4141 math_vector< bigint > x(columns, columns);
4142 math_vector< bigint > x2(columns, columns);
4143
4144 matrix< bigint > Bbigint(rows, columns);
4145 math_vector< bigint > xbigint(columns, columns);
4146 math_vector< bigint > b1bigint(rows, rows);
4147
4148
4149 matrix< long > Blong(rows, columns);
4150 Blong.set_zero_element(0);
4151 math_vector< long > xlong(columns, columns);
4152 math_vector< long > b1long(rows, rows);
4153
4154 // Step 3
4155 lidia_size_t i, j;
4156 bigint PROD;
4157 for (i = 1; i <= len; i++) {
4158 MOD.assign(PRIM[i]);
4159 if (MOD.bit_length() > bigint::bits_per_digit()) {
4160 remainder(Bbigint, *this, MOD);
4161
4162 for (j = 0; j < rows; j++)
4163 LiDIA::best_remainder(b1bigint[j], b[j], MOD);
4164
4165 lidia_info_handler(std::cout << "- > in lanczos (" << i << ", " << len << ")" << std::endl;);
4166 DETtmp = bigint_modul.lanczos_ZmZ(Bbigint, xbigint, b1bigint, MOD);
4167
4168 std::cout << "Test : " << (Bbigint * xbigint - DETtmp * b1bigint) % MOD << std::endl;
4169 x2 = xbigint;
4170 }
4171 else {
4172 MOD.longify(Modlong);
4173 LiDIA::remainder(Blong, *this, Modlong);
4174
4175 for (j = 0; j < rows; j++)
4176 LiDIA::best_remainder(b1long[j], b[j], Modlong);
4177
4178 lidia_info_handler(std::cout << "- > in lanczos (" << i << ", " << len << ")" << std::endl;);
4179 DETtmp = long_modul.lanczos_ZmZ(Blong, xlong, b1long, Modlong);
4180
4181 for (j = 0; j < columns; j++)
4182 x2[j] = xlong[j];
4183 }
4184
4185 if (i == 1) {
4186 x = x2;
4187 DET = DETtmp;
4188 PROD = MOD;
4189 }
4190 else {
4191
4192 DET = chinese_remainder(DET, PROD, DETtmp, MOD);
4193 remainder(DET, DET, PROD*MOD);
4194 std::cout << "->DET = " << DET << " mod " << PROD*MOD << std::endl;
4195 for (j = 0; j < columns; j++) {
4196 x[j] = chinese_remainder(x[j], PROD, x2[j], MOD);
4197 best_remainder(x[j], x[j], (PROD*MOD));
4198 }
4199 PROD *= MOD;
4200 }
4201
4202 if ((*this)*x == (DET*b)) {
4203 RES[0] = -DET;
4204 for (j = 1; j < columns + 1; j++)
4205 RES[j] = x[j - 1];
4206 return;
4207 }
4208 }
4209 }
4210
4211
4212
4213 void matrix< bigint >::
reginvimage1(const matrix<bigint> & A,const matrix<bigint> & B)4214 reginvimage1(const matrix< bigint > & A, const matrix< bigint > & B)
4215 {
4216 //
4217 // Task: C.reginvimage1(A,B);
4218 // => A * C.column(j) = g(j)*B.column(j), j=0,...,B.columns
4219 // => g(j) minimal
4220 // Version: 2.0
4221 //
4222
4223 debug_handler_l(DMESSAGE, "reginvimage1(const matrix< bigint > &, const matrix< bigint > &",
4224 DVALUE + 5);
4225
4226 const modular_bigint_matrix_algorithms< DRMK< bigint >,
4227 dense_fp_matrix_kernel< long, MR< long > >,
4228 dense_fp_matrix_kernel< bigint, MR< bigint > > > Dm_bigint_modul;
4229
4230 const modular_bigint_matrix_algorithms< SRMK< bigint >,
4231 sparse_fp_matrix_kernel< long, MR< long > >,
4232 sparse_fp_matrix_kernel< bigint, MR< bigint > > > Sm_bigint_modul;
4233
4234 bigint H;
4235
4236 D_bigint_modul.hadamard(A, H);
4237 Dm_bigint_modul.reginvimage1(*this, A, B, H);
4238 }
4239
4240
4241
4242 void matrix< bigint >::
reginvimage2(const matrix<bigint> & A,const matrix<bigint> & B)4243 reginvimage2(const matrix< bigint > & A, const matrix< bigint > & B)
4244 {
4245 //
4246 // Task: C.reginvimage2(A,B);
4247 // => A * C.column(j) = g(j)*B.column(j), j=0,...,B.columns
4248 // => g(j) minimal
4249 // Version: 2.0
4250 //
4251
4252 debug_handler_l(DMESSAGE, "reginvimage2(const matrix< bigint > &, const matrix< bigint > &",
4253 DVALUE + 5);
4254
4255 bigint H;
4256
4257 D_bigint_modul.hadamard(A, H);
4258 D_bigint_modul.reginvimage2(*this, A, B);
4259 }
4260
4261
4262
4263 //
4264 // Image
4265 //
4266
4267 void matrix< bigint >::
image1(const matrix<bigint> & A)4268 image1(const matrix< bigint > & A)
4269 {
4270 //
4271 // Task: B.image1(A);
4272 // => The columns of matrix B form a basis
4273 // of the image of matrix A.
4274 // Version: 2.0
4275 //
4276
4277 debug_handler_l(DMESSAGE, "image1(const matrix< bigint > &)", DVALUE + 5);
4278
4279 const modular_bigint_matrix_algorithms< DRMK< bigint >,
4280 dense_fp_matrix_kernel< long, MR< long > >,
4281 dense_fp_matrix_kernel< bigint, MR< bigint > > > Dm_bigint_modul;
4282
4283 const modular_bigint_matrix_algorithms< SRMK< bigint >,
4284 sparse_fp_matrix_kernel< long, MR< long > >,
4285 sparse_fp_matrix_kernel< bigint, MR< bigint > > > Sm_bigint_modul;
4286
4287 bigint H;
4288
4289 D_bigint_modul.hadamard(A, H);
4290 Dm_bigint_modul.image1(*this, A, H);
4291 }
4292
4293
4294
4295 void matrix< bigint >::
image2(const matrix<bigint> & A)4296 image2(const matrix< bigint > & A)
4297 {
4298 //
4299 // Task: B.image2(A);
4300 // => The columns of matrix B form a basis
4301 // of the image of matrix A.
4302 // Version: 2.0
4303 //
4304
4305 debug_handler_l(DMESSAGE, "image2(const matrix< bigint > &)", DVALUE + 5);
4306
4307 const modular_bigint_matrix_algorithms< DRMK< bigint >,
4308 dense_fp_matrix_kernel< long, MR< long > >,
4309 dense_fp_matrix_kernel< bigint, MR< bigint > > > Dm_bigint_modul;
4310
4311 const modular_bigint_matrix_algorithms< SRMK< bigint >,
4312 sparse_fp_matrix_kernel< long, MR< long > >,
4313 sparse_fp_matrix_kernel< bigint, MR< bigint > > > Sm_bigint_modul;
4314
4315 Dm_bigint_modul.image2(*this, A);
4316 }
4317
4318
4319
4320 //
4321 // InvImage
4322 //
4323
4324 void matrix< bigint >::
invimage(const matrix<bigint> & B,const bigint * b)4325 invimage(const matrix< bigint > & B, const bigint * b)
4326 {
4327 //
4328 // Task: v = invimage(B,b);
4329 // => A is a basis of the solution of B*x=b
4330 // Version: 2.0
4331 //
4332
4333 debug_handler_l(DMESSAGE, "invimage(const matrix< bigint > &, const bigint *)", DVALUE + 5);
4334
4335 if (b == NULL)
4336 precondition_error_handler(PRT, "b", "b != NULL",
4337 "void matrix< bigint >::"
4338 "invimage(const matrix<bigint>& B,"
4339 " const bigint * b)",
4340 DMESSAGE, EMESSAGE[1]);
4341
4342 D_bigint_modul.invimage(*this, B, b);
4343 }
4344
4345
4346
4347 void matrix< bigint >::
invimage(const matrix<bigint> & B,const math_vector<bigint> & b)4348 invimage(const matrix< bigint > & B, const math_vector< bigint > &b)
4349 {
4350 //
4351 // Task: v = invimage(B,b);
4352 // => A is a basis of the solution of B*x=b
4353 // Version: 2.0
4354 //
4355
4356 debug_handler_l(DMESSAGE, "invimage(const matrix< bigint > &, const math_vector< bigint > &)", DVALUE + 5);
4357
4358 if (b.size() != B.rows)
4359 precondition_error_handler(b.size(), "b.size", "b.size == B.rows",
4360 B.rows, "B.rows", "b.size == B.rows",
4361 "void matrix<bigint>::"
4362 "invimage(const matrix<bigint>& B,"
4363 " const math_vector< bigint > &b)",
4364 DMESSAGE, EMESSAGE[1]);
4365
4366 bigint *tmp = b.get_data_address();
4367 D_bigint_modul.invimage(*this, B, tmp);
4368 }
4369
4370
4371
4372 //
4373 // Smith normal form
4374 //
4375
4376 void matrix< bigint >::
snf_hartley()4377 snf_hartley()
4378 {
4379 //
4380 // Task: snf Computation
4381 // Algorithm: given in Hartley and Hawkes
4382 // IMPROVEMENTS: Havas, Majewski
4383 // PAPER: Recognizing badly represented Z-modules, Havas
4384 // Version: 2.0
4385 //
4386
4387 debug_handler_l(DMESSAGE, "snf_hartley()", DVALUE + 5);
4388
4389 D_bigint_modul.snf_hartley(*this);
4390 }
4391
4392
4393
4394 void matrix< bigint >::
snf_hartley(matrix<bigint> & T1,matrix<bigint> & T2)4395 snf_hartley(matrix< bigint > & T1, matrix< bigint > & T2)
4396 {
4397 //
4398 // Task: snf Computation
4399 // Algorithm: given in Hartley and Hawkes
4400 // PAPER: Recognizing badly represented Z-modules
4401 // Version: 2.0
4402 //
4403
4404 debug_handler_l(DMESSAGE, "snf_hartley(matrix< bigint > &, matrix< bigint > &)", DVALUE + 5);
4405
4406 if (T1.columns != rows)
4407 T1.set_no_of_columns(rows);
4408 if (T1.rows != rows)
4409 T1.set_no_of_rows(rows);
4410
4411 if (T2.columns != columns)
4412 T2.set_no_of_columns(columns);
4413 if (T2.rows != columns)
4414 T2.set_no_of_rows(columns);
4415
4416 T1.diag(1, 0);
4417 T2.diag(1, 0);
4418
4419 D_bigint_modul.snf_hartley(*this, T1, T2);
4420 }
4421
4422
4423
4424 void matrix< bigint >::
snf_simple()4425 snf_simple()
4426 {
4427 //
4428 // Task: SNF Computation
4429 // Algorithm: given in Hartley and Hawkes
4430 // IMPROVEMENT: Havas
4431 // PAPER: Recognizing badly represented Z-modules
4432 // Version: 2.0
4433 //
4434
4435 debug_handler_l(DMESSAGE, "snf_simple()", DVALUE + 5);
4436
4437 D_bigint_modul.snf_simple(*this);
4438 }
4439
4440
4441
4442 void matrix< bigint >::
snf_simple(matrix<bigint> & T1,matrix<bigint> & T2)4443 snf_simple(matrix< bigint > & T1, matrix< bigint > & T2)
4444 {
4445 //
4446 // Task: SNF Computation
4447 // Algorithm: given in Hartley and Hawkes
4448 // IMPROVEMENT: Havas
4449 // PAPER: Recognizing badly represented Z-modules
4450 // Version: 2.0
4451 //
4452
4453 debug_handler_l(DMESSAGE, "snf_simple(matrix< bigint > &, matrix< bigint > &)", DVALUE + 5);
4454 if (T1.columns != rows)
4455 T1.set_no_of_columns(rows);
4456 if (T1.rows != rows)
4457 T1.set_no_of_rows(rows);
4458
4459 if (T2.columns != columns)
4460 T2.set_no_of_columns(columns);
4461 if (T2.rows != columns)
4462 T2.set_no_of_rows(columns);
4463
4464 T1.diag(1, 0);
4465 T2.diag(1, 0);
4466
4467 D_bigint_modul.snf_simple(*this, T1, T2);
4468 }
4469
4470
4471
4472 void matrix< bigint >::
snf_havas()4473 snf_havas()
4474 {
4475 //
4476 // Task: snf Computation
4477 // Algorithm: given in Hartley and Hawkes
4478 // PAPER: Recognizing badly represented Z-modules
4479 // IMPROVEMENTS: Havas, best reaminder include mgcd
4480 // Version: 2.0
4481 //
4482
4483 debug_handler_l(DMESSAGE, "snf_havas()", DVALUE + 5);
4484
4485 D_bigint_modul.snf_havas(*this);
4486 }
4487
4488
4489
4490 void matrix< bigint >::
snf_havas(matrix<bigint> & T1,matrix<bigint> & T2)4491 snf_havas(matrix< bigint > & T1, matrix< bigint > & T2)
4492 {
4493 //
4494 // Task: snf Computation
4495 // Algorithm: given in Hartley and Hawkes
4496 // PAPER: Recognizing badly represented Z-modules
4497 // IMPROVEMENTS: Havas, best reaminder include mgcd
4498 // Version: 2.0
4499 //
4500
4501 debug_handler_l(DMESSAGE, "snf_havas(matrix< bigint > &, matrix< bigint > &)", DVALUE + 5);
4502
4503 if (T1.columns != rows)
4504 T1.set_no_of_columns(rows);
4505 if (T1.rows != rows)
4506 T1.set_no_of_rows(rows);
4507
4508 if (T2.columns != columns)
4509 T2.set_no_of_columns(columns);
4510 if (T2.rows != columns)
4511 T2.set_no_of_rows(columns);
4512
4513 T1.diag(1, 0);
4514 T2.diag(1, 0);
4515
4516 D_bigint_modul.snf_havas(*this, T1, T2);
4517 }
4518
4519
4520
4521 void matrix< bigint >::
snf_mult(long art)4522 snf_mult(long art)
4523 {
4524 //
4525 // Task: snf Computation
4526 // Algorithm: given in Hartley and Hawkes
4527 // PAPER: Recognizing badly represented Z-modules + pivot selection
4528 // Version: 2.0
4529 //
4530
4531 debug_handler_l(DMESSAGE, "snf_mult(matrix< bigint > &, matrix< bigint > &)", DVALUE + 5);
4532
4533 D_bigint_modul.snf_mult(*this, art);
4534 }
4535
4536
4537
4538 void matrix< bigint >::
snf_mult(matrix<bigint> & T1,matrix<bigint> & T2,long art)4539 snf_mult(matrix< bigint > & T1, matrix< bigint > & T2, long art)
4540 {
4541 //
4542 // Task: snf Computation
4543 // Algorithm: given in Hartley and Hawkes
4544 // PAPER: Recognizing badly represented Z-modules + pivot selection
4545 // Version: 2.0
4546 //
4547
4548 debug_handler_l(DMESSAGE, "snf_mult(matrix< bigint > &, matrix< bigint > &)", DVALUE + 5);
4549
4550 if (T1.columns != rows)
4551 T1.set_no_of_columns(rows);
4552 if (T1.rows != rows)
4553 T1.set_no_of_rows(rows);
4554
4555 if (T2.columns != columns)
4556 T2.set_no_of_columns(columns);
4557 if (T2.rows != columns)
4558 T2.set_no_of_rows(columns);
4559
4560 T1.diag(1, 0);
4561 T2.diag(1, 0);
4562
4563 D_bigint_modul.snf_mult(*this, T1, T2, art);
4564 }
4565
4566
4567
4568 void matrix< bigint >::
snf_add(long art)4569 snf_add(long art)
4570 {
4571 //
4572 // Task: snf Computation
4573 // Algorithm: given in Hartley and Hawkes
4574 // PAPER: Recognizing badly represented Z-modules + pivot selection
4575 // Version: 2.0
4576 //
4577
4578 debug_handler_l(DMESSAGE, "snf_add(long)", DVALUE + 5);
4579
4580 D_bigint_modul.snf_add(*this, art);
4581 }
4582
4583
4584
4585 void matrix< bigint >::
snf_add(matrix<bigint> & T1,matrix<bigint> & T2,long art)4586 snf_add(matrix< bigint > & T1, matrix< bigint > & T2, long art)
4587 {
4588 //
4589 // Task: snf Computation
4590 // Algorithm: given in Hartley and Hawkes
4591 // PAPER: Recognizing badly represented Z-modules + pivot selection
4592 // Version: 2.0
4593 //
4594
4595 debug_handler_l(DMESSAGE, "snf_add(matrix< bigint > &, matrix< bigint > &, long)", DVALUE + 5);
4596
4597 if (T1.columns != rows)
4598 T1.set_no_of_columns(rows);
4599 if (T1.rows != rows)
4600 T1.set_no_of_rows(rows);
4601
4602 if (T2.columns != columns)
4603 T2.set_no_of_columns(columns);
4604 if (T2.rows != columns)
4605 T2.set_no_of_rows(columns);
4606
4607 T1.diag(1, 0);
4608 T2.diag(1, 0);
4609
4610 D_bigint_modul.snf_add(*this, T1, T2, art);
4611 }
4612
4613
4614
4615 void matrix< bigint >::
snf_new(long art)4616 snf_new(long art)
4617 {
4618 //
4619 // Task: snf Computation
4620 // Algorithm: given in Hartley and Hawkes
4621 // PAPER: Recognizing badly represented Z-modules + pivot selection
4622 // Version: 2.0
4623 //
4624
4625 debug_handler_l(DMESSAGE, "snf_new(long)", DVALUE + 5);
4626
4627 D_bigint_modul.snf_new(*this, art);
4628 }
4629
4630
4631
4632 void matrix< bigint >::
snf_new(matrix<bigint> & T1,matrix<bigint> & T2,long art)4633 snf_new(matrix< bigint > & T1, matrix< bigint > & T2, long art)
4634 {
4635 //
4636 // Task: snf Computation
4637 // Algorithm: given in Hartley and Hawkes
4638 // PAPER: Recognizing badly represented Z-modules + pivot selection
4639 // Version: 2.0
4640 //
4641
4642 debug_handler_l(DMESSAGE, "snf_new(matrix< bigint > &, matrix< bigint > &, long)", DVALUE + 5);
4643
4644 if (T1.columns != rows)
4645 T1.set_no_of_columns(rows);
4646 if (T1.rows != rows)
4647 T1.set_no_of_rows(rows);
4648
4649 if (T2.columns != columns)
4650 T2.set_no_of_columns(columns);
4651 if (T2.rows != columns)
4652 T2.set_no_of_rows(columns);
4653
4654 T1.diag(1, 0);
4655 T2.diag(1, 0);
4656
4657 D_bigint_modul.snf_new(*this, T1, T2, art);
4658 }
4659
4660
4661
4662 void matrix< bigint >::
snfmod_dkt(const bigint & mod)4663 snfmod_dkt(const bigint &mod)
4664 {
4665 //
4666 // Task: A.snfmod_dkt(mod);
4667 // => A in Smith normal form
4668 // => h = lattice determinant of lattice formed
4669 // by the columns of matrix A
4670 // PAPER: Asymptotically fast triangulazation of matrices over rings
4671 // IMPROVEMENT: Hafner, McCurly
4672 // Conditions: rank != rows
4673 // Version: 2.0
4674 //
4675
4676 debug_handler_l(DMESSAGE, "snfmod_dkt(const bigint &)", DVALUE + 5);
4677
4678 if (rank() != rows)
4679 precondition_error_handler(rank(), "rank", "rows == rank",
4680 rows, "rows", "rows == rank",
4681 "void matrix< bigint >::"
4682 "snfmod_dkt(const bigint &mod)",
4683 DMESSAGE, EMESSAGE[10]);
4684
4685 D_bigint_modul.snfmod_dkt(*this, mod);
4686 }
4687
4688
4689
4690 void matrix< bigint >::
snfmod_cohen(const bigint & mod)4691 snfmod_cohen(const bigint & mod)
4692 {
4693 //
4694 // Task: A.snfmod_cohen(mod);
4695 // => A in Smith normal form
4696 // => mod = multiple of lattice determinant of lattice formed
4697 // by the columns of matrix A
4698 // Conditions: rank != rows
4699 // Version: 2.0
4700 //
4701
4702 debug_handler_l(DMESSAGE, "snfmod_cohen(const bigint &)", DVALUE + 5);
4703
4704 if (rank() != rows)
4705 precondition_error_handler(rank(), "rank", "rank == rows",
4706 rows, "rows", "rank == rows",
4707 "void matrix< bigint >::"
4708 "snfmod_cohen(const bigint & mod)",
4709 DMESSAGE, EMESSAGE[10]);
4710
4711 D_bigint_modul.snfmod_cohen(*this, mod);
4712 }
4713
4714
4715
4716 //
4717 // END: Linear algebra
4718 // PART 2
4719 //
4720
4721 void matrix< bigint >::
gauss()4722 gauss()
4723 {
4724 debug_handler(DMESSAGE, "gauss()");
4725
4726 D_bigint_modul.gauss(*this);
4727 }
4728
4729
4730
4731 bigint *matrix< bigint >::
mgcd1(const bigint * aconst,lidia_size_t n)4732 mgcd1(const bigint * aconst, lidia_size_t n)
4733 {
4734 //
4735 // DESCRIPTION: RES = T.mgcd1(a, n);
4736 // = > RES[0] = RES[1]*a[0] + ... + RES[n]*a[n-1]
4737 // = > RES[0] = gcd(a[0], ..., a[n-1])
4738 // = > T*a = RES
4739 // ALGORITHM: Blankinship, PIVOT: MINIMUM
4740 // VERSION: 1.8
4741 //
4742
4743 debug_handler("multiple_gcd", "mgcd1(const bigint *, lidia_size_t)");
4744
4745 register lidia_size_t i, j, index, bound;
4746
4747 bigint MIN, TMP, q, r, *Ttmp1, *Ttmp2 = NULL;
4748
4749 if (columns != n)
4750 set_no_of_columns(n);
4751 if (rows != n)
4752 set_no_of_rows(n);
4753
4754 diag(1, 0);
4755
4756 bigint *a = new bigint[n + 1];
4757 memory_handler(a, "multiple_gcd", "mgcd1 :: "
4758 "Error in memory allocation (a)");
4759
4760 for (i = 0; i < n; i++)
4761 a[i].assign(aconst[i]);
4762
4763 // Init
4764 for (index = 0; index < n && a[index].is_zero(); index++);
4765
4766 if (index == n) {
4767 delete[] a;
4768 return new bigint[n];
4769 }
4770 else
4771 bound = index;
4772
4773 do {
4774 // Pivot search: MINIMUM
4775 MIN.assign(a[index]);
4776
4777 for (i = bound; i < n; i++)
4778 if ((abs(MIN) > abs(a[i])) && !a[i].is_zero()) {
4779 MIN.assign(a[i]);
4780 index = i;
4781 }
4782
4783 // first element != 0
4784 for (i = bound; i < n && (a[i].is_zero() || i == index); i++);
4785 if (i < n) {
4786 div_rem(q, r, a[i], MIN);
4787 a[i].assign(r);
4788 Ttmp1 = value[i];
4789 Ttmp2 = value[index];
4790 for (j = 0; j < n; j++) {
4791 LiDIA::multiply(TMP, q, Ttmp2[j]);
4792 LiDIA::subtract(Ttmp1[j], Ttmp1[j], TMP);
4793 }
4794 }
4795 } while (i < n);
4796
4797 Ttmp2 = value[index];
4798
4799 // gcd < 0 ?
4800 if (a[index] < 0) {
4801 a[index].negate();
4802 for (i = 0; i < n; i++)
4803 Ttmp2[i].negate();
4804 }
4805
4806 if (index != 0)
4807 a[0].assign(a[index]);
4808 for (i = 1; i <= n; i++)
4809 a[i].assign(Ttmp2[i - 1]);
4810
4811 return a;
4812 }
4813
4814
4815
4816 bigint *matrix< bigint >::
mgcd2(const bigint * aconst,lidia_size_t n)4817 mgcd2(const bigint * aconst, lidia_size_t n)
4818 {
4819 //
4820 // DESCRIPTION: RES = T.mgcd2(a, n);
4821 // = > RES[0] = RES[1]*a[0] + ... + RES[n]*a[n-1]
4822 // = > RES[0] = gcd(a[0], ..., a[n-1])
4823 // = > T*a = RES
4824 // ALGORITHM: Blankinship
4825 // IMPROVEMENTS: Havas, Majewski, reduction of all elements, MIN assignments
4826 // PAPER: Hermite normal form computation for integer matrices, Havas
4827 // VERSION: 1.8
4828 //
4829
4830 debug_handler("multiple_gcd", "mgcd2(const bigint *, lidia_size_t)");
4831
4832 register lidia_size_t i, j, index, bound, SW;
4833 bigint MIN, TMP, q, r, *Ttmp1, *Ttmp2 = NULL;
4834
4835 if (columns != n)
4836 set_no_of_columns(n);
4837 if (rows != n)
4838 set_no_of_rows(n);
4839 diag(1, 0);
4840
4841 bigint *a = new bigint[n + 1];
4842 memory_handler(a, "multiple_gcd", "mgcd2 :: "
4843 "Error in memory allocation (a)");
4844
4845 for (i = 0; i < n; i++)
4846 a[i].assign(aconst[i]);
4847
4848 // init
4849 for (index = 0; index < n && a[index].is_zero(); index++);
4850
4851 if (index == n) {
4852 delete[] a;
4853 return new bigint[n];
4854 }
4855 else
4856 bound = index;
4857
4858 do {
4859 MIN.assign(a[index]);
4860
4861 // Pivot search: MINIMUM
4862 for (i = bound; i < n; i++)
4863 if ((abs(MIN) > abs(a[i])) && !a[i].is_zero()) {
4864 MIN.assign(a[i]);
4865 index = i;
4866 }
4867
4868 // all elements
4869 SW = 0;
4870 Ttmp2 = value[index];
4871 for (i = bound; i < n; i++)
4872 if ((i != index) && !a[i].is_zero()) {
4873 SW = 1;
4874 Ttmp1 = value[i];
4875 div_rem(q, r, a[i], MIN);
4876 a[i].assign(r);
4877 for (j = 0; j < n; j++) {
4878 LiDIA::multiply(TMP, q, Ttmp2[j]);
4879 LiDIA::subtract(Ttmp1[j], Ttmp1[j], TMP);
4880 }
4881 }
4882 } while (SW == 1);
4883
4884 Ttmp2 = value[index];
4885
4886 // gcd < 0 ?
4887 if (a[index] < 0) {
4888 a[index].negate();
4889 for (i = 0; i < n; i++)
4890 Ttmp2[i].negate();
4891 }
4892
4893 if (index != 0)
4894 a[0].assign(a[index]);
4895 for (i = 1; i <= n; i++)
4896 a[i].assign(Ttmp2[i - 1]);
4897
4898 return a;
4899 }
4900
4901
4902
4903 void matrix< bigint >::
basis_completion(bigint * v,lidia_size_t n)4904 basis_completion(bigint *v, lidia_size_t n)
4905 {
4906 //
4907 // Task: A.hadamard(H);
4908 // => H = hadamard bound of matrix A
4909 // Version: 2.0
4910 //
4911
4912 debug_handler_l(DMESSAGE, "hadamard(bigint &)", DVALUE + 3);
4913
4914 if (bitfield.get_representation() == matrix_flags::dense_representation)
4915 D_bigint_modul.basis_completion(*this, v, n);
4916 else
4917 S_bigint_modul.basis_completion(*this, v, n);
4918 }
4919
4920
4921
4922 void matrix< bigint >::
simple_basis_completion(bigint * v,lidia_size_t n)4923 simple_basis_completion(bigint *v, lidia_size_t n)
4924 {
4925 //
4926 // Task: A.hadamard(H);
4927 // => H = hadamard bound of matrix A
4928 // Version: 2.0
4929 //
4930
4931 debug_handler_l(DMESSAGE, "hadamard(bigint &)", DVALUE + 3);
4932
4933 if (bitfield.get_representation() == matrix_flags::dense_representation)
4934 D_bigint_modul.simple_basis_completion(*this, v, n);
4935 else
4936 S_bigint_modul.simple_basis_completion(*this, v, n);
4937 }
4938
4939
4940
4941 lidia_size_t matrix< bigint >::
cond_matrix(bigint * v,lidia_size_t n)4942 cond_matrix(bigint *v, lidia_size_t n)
4943 {
4944 //
4945 // Task: A.cond_matrix(v, n);
4946 //
4947 // Version: 2.0
4948 //
4949
4950 debug_handler_l(DMESSAGE, "cond_matrix(bigint *, lidia_size_t)", DVALUE + 3);
4951
4952 if (bitfield.get_representation() == matrix_flags::dense_representation)
4953 return D_bigint_modul.cond_matrix(*this, v, n);
4954 else
4955 return S_bigint_modul.cond_matrix(*this, v, n);
4956 }
4957
4958
4959
4960 #undef DRMKex
4961 #undef SRMKex
4962
4963 #undef DVALUE
4964 #undef DMESSAGE
4965 #undef EMESSAGE
4966
4967
4968
4969 #ifdef LIDIA_NAMESPACE
4970 } // end of namespace LiDIA
4971 #endif
4972