1 //==============================================================================================
2 //
3 // This file is part of LiDIA --- a library for computational number theory
4 //
5 // Copyright (c) 1994--2001 the LiDIA Group. All rights reserved.
6 //
7 // See http://www.informatik.tu-darmstadt.de/TI/LiDIA/
8 //
9 //----------------------------------------------------------------------------------------------
10 //
11 // $Id$
12 //
13 // Author : Patrick Theobald (PT)
14 // Changes : See CVS log
15 //
16 //==============================================================================================
17
18
19 #ifndef LIDIA_DENSE_BIGINT_MATRIX_MODULES_CC_GUARD_
20 #define LIDIA_DENSE_BIGINT_MATRIX_MODULES_CC_GUARD_
21
22
23 #include "LiDIA/modular_operations.inl"
24
25
26
27 #ifdef LIDIA_NAMESPACE
28 # ifndef IN_NAMESPACE_LIDIA
29 namespace LiDIA {
30 # endif
31 #endif
32
33
34
35 #define row_oriented_dense_matrix_modules RODMM
36 #define column_oriented_dense_matrix_modules CODMM
37
38
39
40 //
41 // max
42 //
43
44 template< class T >
45 inline void
max_abs(const MR<T> & RES,T & MAX) const46 row_oriented_dense_matrix_modules< T >::max_abs (const MR< T > &RES, T &MAX) const
47 {
48 MAX.assign(abs(RES.value[0][0]));
49
50 register lidia_size_t i, j;
51 register bigint *tmp;
52
53 for (i = 0; i < RES.rows; i++) {
54 tmp = RES.value[i];
55 for (j = 0; j < RES.columns; j++)
56 if (MAX < abs(tmp[j]))
57 MAX.assign(abs(tmp[j]));
58 }
59 }
60
61
62
63 //
64 // special functions
65 //
66
67 template< class T >
68 inline void
subtract_multiple_of_column(MR<T> & A,lidia_size_t l1,const T & q,lidia_size_t l2,lidia_size_t l) const69 row_oriented_dense_matrix_modules< T >::subtract_multiple_of_column (MR< T > &A,
70 lidia_size_t l1,
71 const T &q,
72 lidia_size_t l2,
73 lidia_size_t l) const
74 {
75 for (register lidia_size_t i = 0; i <= l; i++)
76 LiDIA::subtract(A.value[i][l1], A.value[i][l1], q*A.value[i][l2]);
77 }
78
79
80
81 template< class T >
82 inline void
subtract_multiple_of_column_mod(MR<T> & A,lidia_size_t l1,const T & q,lidia_size_t l2,lidia_size_t l,const T & mod) const83 row_oriented_dense_matrix_modules< T >::subtract_multiple_of_column_mod (MR< T > &A,
84 lidia_size_t l1,
85 const T &q,
86 lidia_size_t l2,
87 lidia_size_t l,
88 const T &mod) const
89 {
90 T TMP;
91 for (register lidia_size_t i = 0; i <= l; i++) {
92 LiDIA::mult_mod(TMP, q, A.value[i][l2], mod);
93 LiDIA::sub_mod(A.value[i][l1], A.value[i][l1], TMP, mod);
94 }
95 }
96
97
98
99 template< class T >
100 inline void
normalize_column_mod(MR<T> & A,lidia_size_t l1,const T & q,lidia_size_t l2,lidia_size_t l,const T & mod) const101 row_oriented_dense_matrix_modules< T >::normalize_column_mod (MR< T > &A,
102 lidia_size_t l1,
103 const T &q,
104 lidia_size_t l2,
105 lidia_size_t l,
106 const T &mod) const
107 {
108 T TMP;
109 for (register lidia_size_t i = 0; i <= l; i++) {
110 LiDIA::mult_mod(TMP, q, A.value[i][l2], mod);
111 LiDIA::sub_mod(A.value[i][l1], A.value[i][l1], TMP, mod);
112 if (A.value[i][l1] < A.Zero)
113 A.value[i][l1] += mod;
114 }
115 }
116
117
118
119 template< class T >
120 inline void
negate_column(MR<T> & A,lidia_size_t index,lidia_size_t l) const121 row_oriented_dense_matrix_modules< T >::negate_column (MR< T > &A,
122 lidia_size_t index,
123 lidia_size_t l) const
124 {
125 for (register lidia_size_t i = 0; i <= l; i++)
126 A.value[i][index] = -A.value[i][index];
127 }
128
129
130
131 template< class T >
132 inline void
negate_column_mod(MR<T> & A,lidia_size_t index,lidia_size_t l,const T & mod) const133 row_oriented_dense_matrix_modules< T >::negate_column_mod (MR< T > &A,
134 lidia_size_t index,
135 lidia_size_t l,
136 const T &mod) const
137 {
138 for (register lidia_size_t i = 0; i <= l; i++) {
139 A.value[i][index] = -A.value[i][index];
140 LiDIA::best_remainder(A.value[i][index], A.value[i][index], mod);
141 }
142 }
143
144
145
146 template< class T >
147 inline T *
init_max_array(const MR<T> & A) const148 row_oriented_dense_matrix_modules< T >::init_max_array (const MR< T > &A) const
149 {
150 T *max_array = new T[A.columns];
151 memory_handler(max_array, DMESSAGE, "init :: "
152 "Error in memory allocation (max_array)");
153 T *tmp;
154 for (register lidia_size_t i = 0; i < A.rows; i++) {
155 tmp = A.value[i];
156 for (register lidia_size_t j = 0; j < A.columns; j++)
157 if (max_array[j] < abs(tmp[j]))
158 max_array[j] = abs(tmp[j]);
159 }
160 return max_array;
161 }
162
163
164
165 template< class T >
166 inline void
update_max_array(const MR<T> & A,lidia_size_t i,T * max_array) const167 row_oriented_dense_matrix_modules< T >::update_max_array (const MR< T > &A,
168 lidia_size_t i,
169 T *max_array) const
170 {
171 for (register lidia_size_t j = 0; j < A.rows; j++)
172 if (max_array[i] < abs(A.value[j][i]))
173 max_array[i] = abs(A.value[j][i]);
174 }
175
176
177
178 //
179 // Hermite normal form: pivot search
180 //
181
182 template< class T >
183 inline lidia_size_t
normal_pivot(const MR<T> & A,lidia_size_t startr,lidia_size_t startc,lidia_size_t & index) const184 row_oriented_dense_matrix_modules< T >::normal_pivot (const MR< T > &A,
185 lidia_size_t startr,
186 lidia_size_t startc,
187 lidia_size_t &index) const
188 {
189 lidia_size_t FOUND = 0;
190
191 T *tmp = A.value[startr];
192
193 for (lidia_size_t i = 0; i <= startc; i++)
194 if (tmp[i] != A.Zero) {
195 FOUND++;
196 index = i;
197 for (lidia_size_t j = i + 1; j <= startc; j++)
198 if (tmp[j] != A.Zero && (j != index)) {
199 FOUND++;
200 if (j != index) {
201 if (abs(tmp[j]) < abs(tmp[i]))
202 index = j;
203 break;
204 }
205 }
206 break;
207 }
208 return FOUND;
209 }
210
211
212
213 template< class T >
214 lidia_size_t
min_abs_of_row(const MR<T> & A,lidia_size_t startr,lidia_size_t startc,lidia_size_t & pos) const215 row_oriented_dense_matrix_modules< T >::min_abs_of_row (const MR< T > &A,
216 lidia_size_t startr,
217 lidia_size_t startc,
218 lidia_size_t &pos) const
219 {
220 lidia_size_t i, FOUND = 0;
221 T *tmp = A.value[startr];
222 T val, minval;
223 bool SW = false;
224
225 for (i = 0; i <= startc; i++) {
226 val = abs(tmp[i]);
227 if (val != A.Zero) {
228 FOUND++;
229
230 if (SW == false) {
231 minval = val;
232 pos = i;
233 SW = true;
234 }
235 else
236 if (minval > val) {
237 minval = val;
238 pos = i;
239 }
240 }
241 }
242 return FOUND;
243 }
244
245
246
247 template< class T >
248 lidia_size_t
pivot_sorting_gcd(const MR<T> & A,lidia_size_t startr,lidia_size_t startc,lidia_size_t & index) const249 row_oriented_dense_matrix_modules< T >::pivot_sorting_gcd (const MR< T > &A,
250 lidia_size_t startr,
251 lidia_size_t startc,
252 lidia_size_t &index) const
253 {
254 lidia_size_t FOUND = 0;
255 lidia_size_t index_largest = 0;
256 T *tmp = A.value[startr];
257 T maxval, val;
258
259 for (lidia_size_t i = 0; i <= startc; i++) {
260 val = abs(tmp[i]);
261 if (val != A.Zero) {
262 if (FOUND == 0) {
263 maxval = val;
264 index_largest = i;
265 }
266 else
267 if (maxval < val) {
268 maxval = val;
269 index_largest = i;
270 }
271 FOUND++;
272 }
273 }
274
275 bool SW = false;
276
277 if (FOUND > 1) {
278 for (lidia_size_t i = 0; i <= startc; i++) {
279 val = abs(tmp[i]);
280 if (val != A.Zero && index_largest != i) {
281 if (SW == false) {
282 maxval = val;
283 index = i;
284 SW = true;
285 }
286 else
287 if (maxval <= val) {
288 maxval = val;
289 index = i;
290 }
291 }
292 }
293 }
294 else
295 index = index_largest;
296
297 return FOUND;
298 }
299
300
301
302 template< class T >
303 lidia_size_t
minimal_norm(MR<T> & A,lidia_size_t startr,lidia_size_t startc,lidia_size_t & index,lidia_size_t norm) const304 row_oriented_dense_matrix_modules< T >::minimal_norm (MR< T > &A,
305 lidia_size_t startr,
306 lidia_size_t startc,
307 lidia_size_t &index,
308 lidia_size_t norm) const
309 {
310 lidia_size_t i, FOUND = 0;
311 bigint min_norm, nnorm, TMP;
312 T val, maxval;
313
314 T *tmp = A.value[startr];
315 bool Pair = false;
316 bool SW = false;
317
318 for (i = 0; i <= startc; i++) {
319 val = abs(tmp[i]);
320 if (val != A.Zero) {
321 FOUND++;
322 if (SW == false) {
323 maxval = val;
324 SW = true;
325 index = i;
326 }
327 else {
328 if (val > maxval) {
329 maxval = val;
330 Pair = false;
331 }
332 else if (val == maxval)
333 Pair = true;
334 }
335 }
336 }
337
338 SW = false;
339
340 for (i = 0; i <= startc; i++) {
341 val = abs(tmp[i]);
342 if (val != A.Zero) {
343 if (SW == false) {
344 if (Pair || val != maxval) {
345 min_norm = column_norm(A, i, norm);
346 index = i;
347 SW = true;
348 }
349 }
350 else {
351 if (Pair) {
352 nnorm = column_norm(A, i, norm);
353 if (min_norm > nnorm) {
354 min_norm = nnorm;
355 index = i;
356 }
357 }
358 else {
359 if (val != maxval) {
360 nnorm = column_norm(A, i, norm);
361 if (min_norm > nnorm) {
362 min_norm = nnorm;
363 index = i;
364 }
365 }
366 }
367 }
368 }
369 }
370 return FOUND;
371 }
372
373
374
375 template< class T >
376 lidia_size_t
min_abs_of_row_plus_minimal_norm(MR<T> & A,lidia_size_t startr,lidia_size_t startc,lidia_size_t & index,lidia_size_t norm) const377 row_oriented_dense_matrix_modules< T >::min_abs_of_row_plus_minimal_norm (MR< T > &A,
378 lidia_size_t startr,
379 lidia_size_t startc,
380 lidia_size_t &index,
381 lidia_size_t norm) const
382 {
383 lidia_size_t i, FOUND = 0;
384 bigint min_norm, nnorm, TMP;
385 T minval = 0, val;
386 T *tmp = A.value[startr];
387 bool SW = false;
388
389 for (i = startc; i >= 0; i--) {
390 val = abs(tmp[i]);
391 if (val != A.Zero) {
392 FOUND++;
393
394 if (SW == false) {
395 minval = val;
396 index = i;
397 SW = true;
398 }
399 else {
400 if (minval > val) {
401 index = i;
402 minval = val;
403 }
404 }
405 }
406 }
407
408 if (FOUND > 0) {
409 min_norm = column_norm(A, index, norm);
410
411 // norm computation
412 for (i = index - 1; i >= 0; i--) {
413 if (minval == abs(tmp[i])) {
414 nnorm = column_norm(A, i, norm);
415 if (nnorm < min_norm) {
416 index = i;
417 min_norm = nnorm;
418 }
419 }
420 }
421 }
422 return FOUND;
423 }
424
425
426 template< class T >
427 lidia_size_t
minimal_norm_plus_min_abs_of_row(MR<T> & A,lidia_size_t startr,lidia_size_t startc,lidia_size_t & index,lidia_size_t norm) const428 row_oriented_dense_matrix_modules< T >::minimal_norm_plus_min_abs_of_row (MR< T > &A,
429 lidia_size_t startr,
430 lidia_size_t startc,
431 lidia_size_t &index,
432 lidia_size_t norm) const
433 {
434 lidia_size_t i, FOUND = 0;
435 bigint min_norm, nnorm, TMP;
436 T val, maxval, min_val;
437
438 T *tmp = A.value[startr];
439 bool Pair = false;
440 bool SW = false;
441
442 for (i = 0; i <= startc; i++) {
443 val = abs(tmp[i]);
444 if (val != A.Zero) {
445 FOUND++;
446 if (SW == false) {
447 maxval = val;
448 SW = true;
449 index = i;
450 }
451 else {
452 if (val > maxval) {
453 maxval = val;
454 Pair = false;
455 }
456 else if (val == maxval)
457 Pair = true;
458 }
459 }
460 }
461
462 SW = false;
463
464 for (i = 0; i <= startc; i++) {
465 val = abs(tmp[i]);
466 if (val != A.Zero) {
467 if (SW == false) {
468 if (Pair || val != maxval) {
469 min_val = val;
470 min_norm = column_norm(A, i, norm);
471 index = i;
472 SW = true;
473 }
474 }
475 else {
476 if (Pair) {
477 nnorm = column_norm(A, i, norm);
478 if (min_norm > nnorm) {
479 min_val = val;
480 min_norm = nnorm;
481 index = i;
482 }
483 else if (min_norm == nnorm)
484 if (min_val > val) {
485 min_val = val;
486 min_norm = nnorm;
487 index = i;
488 }
489 }
490 else {
491 if (val != maxval) {
492 nnorm = column_norm(A, i, norm);
493 if (min_norm > nnorm) {
494 min_val = val;
495 min_norm = nnorm;
496 index = i;
497 }
498 else if (min_norm == nnorm)
499 if (min_val > val) {
500 min_val = val;
501 min_norm = nnorm;
502 index = i;
503 }
504 }
505 }
506 }
507 }
508 }
509 return FOUND;
510 }
511
512
513
514 template< class T >
515 lidia_size_t
minimal_norm_plus_sorting_gcd(MR<T> & A,lidia_size_t startr,lidia_size_t startc,lidia_size_t & index,lidia_size_t norm) const516 row_oriented_dense_matrix_modules< T >::minimal_norm_plus_sorting_gcd (MR< T > &A,
517 lidia_size_t startr,
518 lidia_size_t startc,
519 lidia_size_t &index,
520 lidia_size_t norm) const
521 {
522 lidia_size_t i, FOUND = 0;
523 bigint min_norm, nnorm, TMP;
524 T val, maxval, min_val;
525
526 T *tmp = A.value[startr];
527 bool Pair = false;
528 bool SW = false;
529
530 for (i = 0; i <= startc; i++) {
531 val = abs(tmp[i]);
532 if (val != A.Zero) {
533 FOUND++;
534 if (SW == false) {
535 maxval = val;
536 SW = true;
537 index = i;
538 }
539 else {
540 if (val > maxval) {
541 maxval = val;
542 Pair = false;
543 }
544 else if (val == maxval)
545 Pair = true;
546 }
547 }
548 }
549
550 SW = false;
551
552 for (i = 0; i <= startc; i++) {
553 val = abs(tmp[i]);
554 if (val != A.Zero) {
555 if (SW == false) {
556 if (Pair || val != maxval) {
557 min_val = val;
558 min_norm = column_norm(A, i, norm);
559 index = i;
560 SW = true;
561 }
562 }
563 else {
564 if (Pair) {
565 nnorm = column_norm(A, i, norm);
566 if (min_norm > nnorm) {
567 min_val = val;
568 min_norm = nnorm;
569 index = i;
570 }
571 else if (min_norm == nnorm)
572 if (min_val < val) {
573 min_val = val;
574 min_norm = nnorm;
575 index = i;
576 }
577 }
578 else {
579 if (val != maxval) {
580 nnorm = column_norm(A, i, norm);
581 if (min_norm > nnorm) {
582 min_val = val;
583 min_norm = nnorm;
584 index = i;
585 }
586 else if (min_norm == nnorm)
587 if (min_val < val) {
588 min_val = val;
589 min_norm = nnorm;
590 index = i;
591 }
592 }
593 }
594 }
595 }
596 }
597 return FOUND;
598 }
599
600
601
602 template< class T >
603 lidia_size_t
minimal_norm_plus_min_no_of_elements(MR<T> & A,lidia_size_t startr,lidia_size_t startc,lidia_size_t & index,lidia_size_t norm) const604 row_oriented_dense_matrix_modules< T >::minimal_norm_plus_min_no_of_elements (MR< T > &A,
605 lidia_size_t startr,
606 lidia_size_t startc,
607 lidia_size_t &index,
608 lidia_size_t norm) const
609 {
610 lidia_size_t i, FOUND = 0;
611 bigint min_no, no;
612 bigint min_norm, nnorm, TMP;
613 T val, maxval;
614
615 T *tmp = A.value[startr];
616 bool Pair = false;
617 bool SW = false;
618
619 for (i = 0; i <= startc; i++) {
620 val = abs(tmp[i]);
621 if (val != A.Zero) {
622 FOUND++;
623 if (SW == false) {
624 maxval = val;
625 SW = true;
626 index = i;
627 }
628 else {
629 if (val > maxval) {
630 maxval = val;
631 Pair = false;
632 }
633 else if (val == maxval)
634 Pair = true;
635 }
636 }
637 }
638
639 SW = false;
640
641 for (i = 0; i <= startc; i++) {
642 val = abs(tmp[i]);
643 if (val != A.Zero) {
644 if (SW == false) {
645 if (Pair || val != maxval) {
646 min_no = column_norm(A, i, 0);
647 min_norm = column_norm(A, i, norm);
648 index = i;
649 SW = true;
650 }
651 }
652 else {
653 if (Pair) {
654 no = column_norm(A, i, 0);
655 nnorm = column_norm(A, i, norm);
656 if (min_norm > nnorm) {
657 min_no = no;
658 min_norm = nnorm;
659 index = i;
660 }
661 else if (min_norm == nnorm)
662 if (min_no > no) {
663 min_no = no;
664 min_norm = nnorm;
665 index = i;
666 }
667 }
668 else {
669 if (val != maxval) {
670 no = column_norm(A, i, 0);
671 nnorm = column_norm(A, i, norm);
672 if (min_norm > nnorm) {
673 min_no = no;
674 min_norm = nnorm;
675 index = i;
676 }
677 else if (min_norm == nnorm)
678 if (min_no > no) {
679 min_no = no;
680 min_norm = nnorm;
681 index = i;
682 }
683 }
684 }
685 }
686 }
687 }
688 return FOUND;
689 }
690
691
692
693 template< class T >
694 lidia_size_t
sorting_gcd_plus_minimal_norm(MR<T> & A,lidia_size_t startr,lidia_size_t startc,lidia_size_t & index,lidia_size_t norm) const695 row_oriented_dense_matrix_modules< T >::sorting_gcd_plus_minimal_norm (MR< T > &A,
696 lidia_size_t startr,
697 lidia_size_t startc,
698 lidia_size_t &index,
699 lidia_size_t norm) const
700 {
701 lidia_size_t FOUND = 0;
702 lidia_size_t index_largest = 0;
703 T *tmp = A.value[startr];
704 T maxval, val;
705
706 for (lidia_size_t i = 0; i <= startc; i++) {
707 val = abs(tmp[i]);
708 if (val != A.Zero) {
709 if (FOUND == 0) {
710 maxval = val;
711 index_largest = i;
712 }
713 else
714 if (maxval < val) {
715 maxval = val;
716 index_largest = i;
717 }
718 FOUND++;
719 }
720 }
721
722 bool SW = false;
723
724 if (FOUND > 1) {
725 for (lidia_size_t i = 0; i <= startc; i++) {
726 val = abs(tmp[i]);
727 if (val != A.Zero && index_largest != i) {
728 if (SW == false) {
729 maxval = val;
730 index = i;
731 SW = true;
732 }
733 else
734 if (maxval <= val) {
735 maxval = val;
736 index = i;
737 }
738 }
739 }
740 }
741 else {
742 index = index_largest;
743 maxval = abs(tmp[index]);
744 }
745
746 bigint nnorm, min_norm = column_norm(A, index, norm);
747
748 for (lidia_size_t i = 0; i <= startc; i++)
749 if (abs(tmp[i]) == maxval) {
750 nnorm = column_norm(A, index, norm);
751 if (nnorm < min_norm) {
752 min_norm = nnorm;
753 index = i;
754 }
755 }
756
757 return FOUND;
758 }
759
760
761
762 template< class T >
763 lidia_size_t
min_no_of_elements_plus_minimal_norm(MR<T> & A,lidia_size_t startr,lidia_size_t startc,lidia_size_t & index,lidia_size_t norm) const764 row_oriented_dense_matrix_modules< T >::min_no_of_elements_plus_minimal_norm (MR< T > &A,
765 lidia_size_t startr,
766 lidia_size_t startc,
767 lidia_size_t &index,
768 lidia_size_t norm) const
769 {
770 lidia_size_t i, FOUND = 0;
771 bigint min_no, no;
772 bigint min_norm, nnorm, TMP;
773 T val, maxval;
774
775 T *tmp = A.value[startr];
776 bool Pair = false;
777 bool SW = false;
778
779 for (i = 0; i <= startc; i++) {
780 val = abs(tmp[i]);
781 if (val != A.Zero) {
782 FOUND++;
783 if (SW == false) {
784 maxval = val;
785 SW = true;
786 index = i;
787 }
788 else {
789 if (val > maxval) {
790 maxval = val;
791 Pair = false;
792 }
793 else if (val == maxval)
794 Pair = true;
795 }
796 }
797 }
798
799 SW = false;
800
801 for (i = 0; i <= startc; i++) {
802 val = abs(tmp[i]);
803 if (val != A.Zero) {
804 if (SW == false) {
805 if (Pair || val != maxval) {
806 min_no = column_norm(A, i, norm);
807 min_norm = column_norm(A, i, 0);
808 index = i;
809 SW = true;
810 }
811 }
812 else {
813 if (Pair) {
814 no = column_norm(A, i, norm);
815 nnorm = column_norm(A, i, 0);
816 if (min_norm > nnorm) {
817 min_no = no;
818 min_norm = nnorm;
819 index = i;
820 }
821 else if (min_norm == nnorm)
822 if (min_no > no) {
823 min_no = no;
824 min_norm = nnorm;
825 index = i;
826 }
827 }
828 else {
829 if (val != maxval) {
830 no = column_norm(A, i, norm);
831 nnorm = column_norm(A, i, 0);
832 if (min_norm > nnorm) {
833 min_no = no;
834 min_norm = nnorm;
835 index = i;
836 }
837 else if (min_norm == nnorm)
838 if (min_no > no) {
839 min_no = no;
840 min_norm = nnorm;
841 index = i;
842 }
843 }
844 }
845 }
846 }
847 }
848 return FOUND;
849 }
850
851
852
853 //
854 // Norm
855 //
856
857 template< class T >
858 inline bigint
column_norm(MR<T> & A,lidia_size_t pos,lidia_size_t norm) const859 row_oriented_dense_matrix_modules< T >::column_norm (MR< T > &A,
860 lidia_size_t pos,
861 lidia_size_t norm) const
862 {
863 bigint Norm = A.Zero;
864 bigint TMP;
865 lidia_size_t j;
866
867
868 if (norm == 0) {
869 for (j = 0; j < A.rows; j++)
870 if (A.value[j][pos] != A.Zero)
871 Norm++;
872 }
873 else if (norm == 1) {
874 for (j = 0; j < A.rows; j++)
875 LiDIA::add(Norm, Norm, A.value[j][pos]);
876 }
877 else {
878 for (j = 0; j < A.rows; j++) {
879 power(TMP, A.value[j][pos], norm);
880 LiDIA::add(Norm, Norm, abs(TMP));
881 }
882 }
883 return Norm;
884 }
885
886
887
888 template< class T >
889 inline void
kennwerte(MR<T> & RES,T & MAX,lidia_size_t & no_of_elements,T & Durch) const890 row_oriented_dense_matrix_modules< T >::kennwerte (MR< T > &RES,
891 T &MAX,
892 lidia_size_t &no_of_elements,
893 T &Durch) const
894 {
895 register lidia_size_t i, j;
896 bool SW = false;
897 no_of_elements = 0;
898 Durch = 0;
899
900 for (i = 0; i < RES.rows; i++) {
901 for (j = 0; j < RES.columns; j++) {
902 if (RES.value[i][j] != RES.Zero) {
903 no_of_elements++;
904 Durch += abs(RES.value[i][j]);
905 if (SW == false) {
906 SW = true;
907 MAX = abs(RES.value[i][j]);
908 }
909 else
910 if (MAX < abs(RES.value[i][j]))
911 MAX = abs(RES.value[i][j]);
912 }
913 }
914 }
915 Durch /= no_of_elements;
916 }
917
918
919
920 template< class T >
921 inline void
max(MR<T> & RES,T & MAX) const922 row_oriented_dense_matrix_modules< T >::max (MR< T > &RES, T &MAX) const
923 {
924 register lidia_size_t i, j;
925 bool SW = false;
926
927 for (i = 0; i < RES.rows; i++)
928 for (j = 0; j < RES.columns; j++)
929 if (SW == false) {
930 SW = true;
931 MAX = abs(RES.value[i][j]);
932 }
933 else
934 if (MAX < abs(RES.value[i][j]))
935 MAX = abs(RES.value[i][j]);
936 }
937
938
939
940 //
941 // Hermite normal form: elimination
942 //
943
944 template< class T >
945 inline bool
normalize_one_element_of_row(MR<T> & A,lidia_size_t startr,lidia_size_t startc,lidia_size_t index,lidia_size_t len) const946 row_oriented_dense_matrix_modules< T >::normalize_one_element_of_row (MR< T > &A,
947 lidia_size_t startr,
948 lidia_size_t startc,
949 lidia_size_t index,
950 lidia_size_t len) const
951 {
952 T q, TMP;
953 for (lidia_size_t i = 0; i <= startc; i++)
954 if ((i != index) && this->member(A, startr, i) != A.Zero) {
955 pos_div_rem(q, TMP, this->member(A, startr, i), this->member(A, startr, index));
956 if (q != 0)
957 subtract_multiple_of_column(A, i, q, index, len);
958 break;
959 }
960 return true;
961 }
962
963
964
965 template< class T >
966 inline bool
normalize_one_element_of_row(MR<T> & A,matrix<bigint> & TR,lidia_size_t startr,lidia_size_t startc,lidia_size_t index,lidia_size_t len) const967 row_oriented_dense_matrix_modules< T >::normalize_one_element_of_row (MR< T > &A,
968 matrix< bigint > &TR,
969 lidia_size_t startr,
970 lidia_size_t startc,
971 lidia_size_t index,
972 lidia_size_t len) const
973 {
974 T q, TMP;
975 for (lidia_size_t i = 0; i <= startc; i++)
976 if ((i != index) && this->member(A, startr, i) != A.Zero) {
977 pos_div_rem(q, TMP, this->member(A, startr, i), this->member(A, startr, index));
978 if (q != 0) {
979 subtract_multiple_of_column(A, i, q, index, len);
980 subtract_multiple_of_column(TR, i, q, index, TR.get_no_of_rows() - 1);
981 }
982 break;
983 }
984 return true;
985 }
986
987
988
989 template< class T >
990 inline bool
normalize_one_element_of_row(MR<T> & A,math_vector<bigfloat> & vec,lidia_size_t startr,lidia_size_t startc,lidia_size_t index,lidia_size_t len) const991 row_oriented_dense_matrix_modules< T >::normalize_one_element_of_row (MR< T > &A,
992 math_vector< bigfloat > &vec,
993 lidia_size_t startr,
994 lidia_size_t startc,
995 lidia_size_t index,
996 lidia_size_t len) const
997 {
998 T q, TMP;
999 bigfloat rtemp;
1000 for (lidia_size_t i = 0; i <= startc; i++)
1001 if ((i != index) && this->member(A, startr, i) != A.Zero) {
1002 pos_div_rem(q, TMP, this->member(A, startr, i), this->member(A, startr, index));
1003 if (q != 0) {
1004 subtract_multiple_of_column(A, i, q, index, len);
1005 LiDIA::multiply(rtemp, bigfloat(q), vec.member(index));
1006 LiDIA::subtract(vec[i], vec[i], rtemp);
1007 }
1008 break;
1009 }
1010 return true;
1011 }
1012
1013
1014
1015 template< class T >
1016 bool
mgcd_linear(MR<T> & A,lidia_size_t startr,lidia_size_t startc,lidia_size_t len) const1017 row_oriented_dense_matrix_modules< T >::mgcd_linear (MR< T > &A,
1018 lidia_size_t startr,
1019 lidia_size_t startc,
1020 lidia_size_t len) const
1021 {
1022 T TMP, RES0, RES1, RES2, TMP1, x, y;
1023
1024 T *tmp = A.value[startr];
1025
1026 // Init
1027 lidia_size_t index;
1028 for (index = startc; tmp[index] == A.Zero && index >= 0; index--);
1029 if (index < 0)
1030 return false;
1031 if (index != startc)
1032 this->swap_columns(A, startc, index);
1033
1034 for (lidia_size_t i = index - 1; i >= 0; i--)
1035 if (tmp[i] != A.Zero) {
1036 RES2 = xgcd(RES0, RES1, tmp[i], tmp[startc]);
1037 x = tmp[startc] / RES2;
1038 y = tmp[i] / RES2;
1039
1040 for (lidia_size_t l = 0; l <= len; l++) {
1041 TMP = A.value[l][i];
1042 TMP1 = A.value[l][startc];
1043
1044 A.value[l][i] = (TMP*x - TMP1*y);
1045 A.value[l][startc] = (TMP*RES0 + TMP1*RES1);
1046 }
1047 }
1048 return true;
1049 }
1050
1051
1052
1053 template< class T >
1054 bool
mgcd_linear(MR<T> & A,matrix<bigint> & TR,lidia_size_t startr,lidia_size_t startc,lidia_size_t len) const1055 row_oriented_dense_matrix_modules< T >::mgcd_linear (MR< T > &A,
1056 matrix< bigint > &TR,
1057 lidia_size_t startr,
1058 lidia_size_t startc,
1059 lidia_size_t len) const
1060 {
1061 T TMP, RES0, RES1, RES2, TMP1, x, y;
1062
1063 bigint TMP_1, TMP_2;
1064 T *tmp = A.value[startr];
1065
1066 // Init
1067 lidia_size_t index;
1068 for (index = startc; tmp[index] == A.Zero && index >= 0; index--);
1069 if (index < 0)
1070 return false;
1071 if (index != startc) {
1072 this->swap_columns(A, startc, index);
1073 TR.swap_columns(startc, index);
1074 }
1075
1076 for (lidia_size_t i = index - 1; i >= 0; i--)
1077 if (tmp[i] != A.Zero) {
1078 RES2 = xgcd(RES0, RES1, tmp[i], tmp[startc]);
1079 x = tmp[startc] / RES2;
1080 y = tmp[i] / RES2;
1081
1082 for (lidia_size_t l = 0; l <= len; l++) {
1083 TMP = A.value[l][i];
1084 TMP1 = A.value[l][startc];
1085
1086 A.value[l][i] = TMP*x - TMP1*y;
1087 A.value[l][startc] = TMP*RES0 + TMP1*RES1;
1088 }
1089
1090 for (lidia_size_t l = 0; l < TR.get_no_of_rows(); l++) {
1091 TMP_1 = TR.member(l, i);
1092 TMP_2 = TR.member(l, startc);
1093
1094 TR.sto(l, i, (TMP_1*x - TMP_2*y));
1095 TR.sto(l, startc, (TMP_1*RES0 + TMP_2*RES1));
1096 }
1097 }
1098 return true;
1099 }
1100
1101
1102
1103 template< class T >
1104 bool
mgcd_bradley(MR<T> & A,lidia_size_t startr,lidia_size_t startc,lidia_size_t len) const1105 row_oriented_dense_matrix_modules< T >::mgcd_bradley (MR< T > &A,
1106 lidia_size_t startr,
1107 lidia_size_t startc,
1108 lidia_size_t len) const
1109 {
1110 lidia_size_t n = startc + 1;
1111 register lidia_size_t i, j;
1112
1113
1114 matrix< bigint > Tr(A.columns, A.columns);
1115 Tr.diag(1, 0);
1116
1117 T *a = A.value[startr];
1118
1119 if (n <= 0)
1120 lidia_error_handler("multiple_gcd", "mgcd :: Error in parameter !!");
1121
1122 if (n > 1) {
1123 // Step 1 - 8
1124 bigint y, z, TMP, TMP1, t1, t2;
1125 bigint g_old = a[0];
1126 bigint g = g_old;
1127 for (i = 1; i < n; i++) {
1128 if (a[i] != A.Zero) {
1129 g = xgcd(y, z, g_old, a[i]);
1130 TMP = -a[i]/g;
1131 TMP1 = g_old/g;
1132 for (j = 0; j < Tr.rows; j++) {
1133 t1 = Tr.value[j][i-1];
1134 t2 = Tr.value[j][i];
1135 Tr.value[j][i-1] = TMP*t1+TMP1*t2;
1136 Tr.value[j][i] = y*t1+z*t2;
1137 }
1138 g_old = g;
1139 }
1140 else {
1141 Tr.swap_columns(i-1, i);
1142 }
1143 }
1144 if (g == 0)
1145 return false;
1146
1147 // Balaciertes Rueckwaertseinsetzen
1148 bigint q, r;
1149 lidia_size_t l, k;
1150 for (i = n - 1, j = n - 2; i >= 0 && j >= 0; i--, j--)
1151 for (k = j + 1; k < n; k++) {
1152 if (Tr.value[i][j] != Tr.Zero) {
1153 div_rem(q, r, Tr.value[i][k], Tr.value[i][j]);
1154
1155 Tr.value[i][k] = Tr.value[i][k] - q*Tr.value[i][j];
1156 TMP = Tr.value[i][k] + Tr.value[i][j];
1157 TMP1 = Tr.value[i][k] - Tr.value[i][j];
1158 if (abs(TMP) < abs(Tr.value[i][k])) {
1159 q--;
1160 Tr.value[i][k] = TMP;
1161 }
1162 else if (abs(TMP1) < abs(Tr.value[i][k])) {
1163 q++;
1164 Tr.value[i][k] = TMP1;
1165 }
1166
1167 if (abs(q) > 0)
1168 for (l = 0; l < i; l++)
1169 Tr.value[l][k] = Tr.value[l][k] - q*Tr.value[l][j];
1170 }
1171 }
1172 }
1173 //multiply(A, A, Tr);
1174 T TMP;
1175 T *tmp2 = new T[A.columns];
1176 for (i = 0; i < A.rows; i++) {
1177
1178 for (lidia_size_t p = 0; p < A.columns; p++)
1179 tmp2[p] = this->member(A, i, p);
1180
1181 for (j = 0; j < A.columns; j++) {
1182 TMP = A.Zero;
1183 for (lidia_size_t l = 0; l < A.columns; l++)
1184 if (tmp2[l] != A.Zero && Tr.member(l, j) != A.Zero)
1185 TMP += T(tmp2[l] * Tr.member(l, j));
1186
1187 this->sto(A, i, j, TMP);
1188 }
1189 }
1190 delete[] tmp2;
1191
1192 return true;
1193 }
1194
1195
1196
1197 template< class T >
1198 bool
mgcd_bradley(MR<T> & A,matrix<bigint> & TR,lidia_size_t startr,lidia_size_t startc,lidia_size_t len) const1199 row_oriented_dense_matrix_modules< T >::mgcd_bradley (MR< T > &A,
1200 matrix< bigint > &TR,
1201 lidia_size_t startr,
1202 lidia_size_t startc,
1203 lidia_size_t len) const
1204 {
1205 lidia_size_t n = startc + 1;
1206 register lidia_size_t i, j;
1207
1208
1209 matrix< bigint > Tr(A.columns, A.columns);
1210 Tr.diag(1, 0);
1211
1212 T *a = A.value[startr];
1213
1214 if (n <= 0)
1215 lidia_error_handler("multiple_gcd", "mgcd :: Error in parameter !!");
1216
1217 if (n > 1) {
1218 // Step 1 - 8
1219 bigint y, z, TMP, TMP1, t1, t2;
1220 bigint g_old = a[0];
1221 bigint g = g_old;
1222 for (i = 1; i < n; i++) {
1223 if (a[i] != A.Zero) {
1224 g = xgcd(y, z, g_old, a[i]);
1225 TMP = -a[i]/g;
1226 TMP1 = g_old/g;
1227 for (j = 0; j < n; j++) {
1228 t1 = Tr.value[j][i-1];
1229 t2 = Tr.value[j][i];
1230 Tr.value[j][i-1] = TMP*t1+TMP1*t2;
1231 Tr.value[j][i] = y*t1+z*t2;
1232 }
1233 g_old = g;
1234 }
1235 else {
1236 Tr.swap_columns(i-1, i);
1237 }
1238 }
1239 if (g == 0)
1240 return false;
1241
1242 // Balaciertes Rueckwaertseinsetzen
1243 bigint q, r;
1244 lidia_size_t l, k;
1245 for (i = n - 1, j = n - 2; i >= 0 && j >= 0; i--, j--)
1246 for (k = j + 1; k < n; k++) {
1247 if (Tr.value[i][j] != Tr.Zero) {
1248 div_rem(q, r, Tr.value[i][k], Tr.value[i][j]);
1249
1250 Tr.value[i][k] = Tr.value[i][k] - q*Tr.value[i][j];
1251 TMP = Tr.value[i][k] + Tr.value[i][j];
1252 TMP1 = Tr.value[i][k] - Tr.value[i][j];
1253 if (abs(TMP) < abs(Tr.value[i][k])) {
1254 q--;
1255 Tr.value[i][k] = TMP;
1256 }
1257 else if (abs(TMP1) < abs(Tr.value[i][k])) {
1258 q++;
1259 Tr.value[i][k] = TMP1;
1260 }
1261
1262 if (abs(q) > 0)
1263 for (l = 0; l < i; l++)
1264 Tr.value[l][k] = Tr.value[l][k] - q*Tr.value[l][j];
1265 }
1266 }
1267 }
1268 //multiply(A, A, Tr);
1269 T TMP;
1270 T *tmp2 = new T[A.columns];
1271 for (i = 0; i < A.rows; i++) {
1272
1273 for (lidia_size_t p = 0; p < A.columns; p++)
1274 tmp2[p] = this->member(A, i, p);
1275
1276 for (j = 0; j < A.columns; j++) {
1277 TMP = A.Zero;
1278 for (lidia_size_t l = 0; l < A.columns; l++)
1279 if (tmp2[l] != A.Zero && Tr.member(l, j) != A.Zero)
1280 TMP += T(tmp2[l] * Tr.member(l, j));
1281
1282 this->sto(A, i, j, TMP);
1283 }
1284 }
1285 delete[] tmp2;
1286
1287 LiDIA::multiply(TR, TR, Tr);
1288 return true;
1289 }
1290
1291
1292
1293 template< class T >
1294 bool
mgcd_ilio(MR<T> & A,lidia_size_t startr,lidia_size_t startc,lidia_size_t len) const1295 row_oriented_dense_matrix_modules< T >::mgcd_ilio (MR< T > &A,
1296 lidia_size_t startr,
1297 lidia_size_t startc,
1298 lidia_size_t len) const
1299 {
1300 lidia_size_t n = startc + 1;
1301
1302 bigint *a = A.value[startr];
1303
1304 register lidia_size_t i, j;
1305
1306 if (n <= 0)
1307 lidia_error_handler("multiple_gcd", "mgcd :: Error in parameter !!");
1308
1309 T g = A.Zero;
1310 if (n > 1) {
1311 bigint y, z, TMP, TMP1, t1, t2;
1312
1313 lidia_size_t step = 1, step2 = 2;
1314 while (step < n) {
1315 for (i = 0; i < n; i += step2) {
1316 if (i + step < n) {
1317 g = xgcd(y, z, a[i], a[i+step]);
1318 if (g != A.Zero) {
1319 TMP = -a[i+step]/g;
1320 TMP1 = a[i]/g;
1321
1322 for (j = 0; j < A.rows; j++) {
1323 t1 = A.value[j][i];
1324 t2 = A.value[j][i+step];
1325 A.value[j][i] = y*t1+z*t2;
1326 A.value[j][i+step] = TMP*t1+TMP1*t2;
1327 }
1328 }
1329 }
1330 }
1331 step = step2;
1332 step2 = 2*step2;
1333 }
1334 for (j = 0; j < A.rows; j++)
1335 LiDIA::swap(A.value[j][0], A.value[j][n-1]);
1336 if (A.value[startr][startc] < 0)
1337 for (j = 0; j < A.rows; j++)
1338 A.value[j][startc].negate();
1339 }
1340 if (g == A.Zero)
1341 return false;
1342 return true;
1343 }
1344
1345
1346
1347 template< class T >
1348 bool
mgcd_ilio(MR<T> & A,matrix<bigint> & Tr,lidia_size_t startr,lidia_size_t startc,lidia_size_t len) const1349 row_oriented_dense_matrix_modules< T >::mgcd_ilio (MR< T > &A,
1350 matrix< bigint > &Tr,
1351 lidia_size_t startr,
1352 lidia_size_t startc,
1353 lidia_size_t len) const
1354 {
1355 lidia_size_t n = startc + 1;
1356
1357 bigint *a = A.value[startr];
1358
1359 register lidia_size_t i, j;
1360
1361 if (n <= 0)
1362 lidia_error_handler("multiple_gcd", "mgcd :: Error in parameter !!");
1363
1364 T g = A.Zero;
1365 if (n > 1) {
1366 bigint y, z, TMP, TMP1, t1, t2;
1367
1368 lidia_size_t step = 1, step2 = 2;
1369 while (step < n) {
1370 for (i = 0; i < n; i += step2) {
1371 if (i + step < n) {
1372 g = xgcd(y, z, a[i], a[i+step]);
1373 if (g != A.Zero) {
1374 TMP = -a[i+step]/g;
1375 TMP1 = a[i]/g;
1376 for (j = 0; j < A.columns; j++) {
1377 t1 = Tr.member(j, i);
1378 t2 = Tr.member(j, i + step);
1379 Tr.sto(j, i, y*t1+z*t2);
1380 Tr.sto(j, i + step, TMP*t1+TMP1*t2);
1381 }
1382
1383 for (j = 0; j < A.rows; j++) {
1384 t1 = A.value[j][i];
1385 t2 = A.value[j][i+step];
1386 A.value[j][i] = y*t1+z*t2;
1387 A.value[j][i+step] = TMP*t1+TMP1*t2;
1388 }
1389 }
1390 }
1391 }
1392 step = step2;
1393 step2 = 2*step2;
1394 }
1395 for (j = 0; j < A.rows; j++)
1396 LiDIA::swap(A.value[j][0], A.value[j][n-1]);
1397 Tr.swap_columns(0, n-1);
1398 if (A.value[startr][startc] < 0)
1399 {
1400 for (j = 0; j < A.rows; j++)
1401 A.value[j][startc].negate();
1402 for (j = 0; j < A.columns; j++)
1403 Tr.sto(j, startc, -Tr.member(j, startc));
1404 }
1405 }
1406 if (g == A.Zero)
1407 return false;
1408 return true;
1409 }
1410
1411
1412
1413 template< class T >
1414 bool
mgcd_opt(MR<T> & A,lidia_size_t startr,lidia_size_t startc,lidia_size_t len) const1415 row_oriented_dense_matrix_modules< T >::mgcd_opt (MR< T > &A,
1416 lidia_size_t startr,
1417 lidia_size_t startc,
1418 lidia_size_t len) const
1419 {
1420 T TMP, RES0, RES1, RES2, TMP1, x, y;
1421
1422 T *tmp = A.value[startr];
1423
1424 // Init
1425 lidia_size_t index;
1426 for (index = startc; tmp[index] == A.Zero && index >= 0; index--);
1427 if (index < 0)
1428 return false;
1429 if (index != startc)
1430 this->swap_columns(A, startc, index);
1431
1432 // sort
1433 if (startc >= 2) {
1434 lidia_size_t gcdindex = startc - 1;
1435 T gcd_min = gcd(tmp[startc], tmp[startc - 1]), gcd_tmp;
1436 for (lidia_size_t i = startc - 2; i >= 0; i--) {
1437 gcd_tmp = gcd(tmp[startc], tmp[i]);
1438 if (gcd_tmp < gcd_min)
1439 gcdindex = i;
1440 }
1441 this->swap_columns(A, gcdindex, startc - 1);
1442 }
1443
1444 for (lidia_size_t i = startc - 1; i >= 0; i--)
1445 if (tmp[i] != A.Zero) {
1446 RES2 = xgcd(RES0, RES1, tmp[i], tmp[startc]);
1447 x = tmp[startc] / RES2;
1448 y = tmp[i] / RES2;
1449
1450 for (lidia_size_t l = 0; l <= len; l++) {
1451 TMP = A.value[l][i];
1452 TMP1 = A.value[l][startc];
1453
1454 A.value[l][i] = (TMP*x - TMP1*y);
1455 A.value[l][startc] = (TMP*RES0 + TMP1*RES1);
1456 }
1457 }
1458 return true;
1459 }
1460
1461
1462
1463 template< class T >
1464 bool
mgcd_opt(MR<T> & A,matrix<bigint> & TR,lidia_size_t startr,lidia_size_t startc,lidia_size_t len) const1465 row_oriented_dense_matrix_modules< T >::mgcd_opt (MR< T > &A,
1466 matrix< bigint > &TR,
1467 lidia_size_t startr,
1468 lidia_size_t startc,
1469 lidia_size_t len) const
1470 {
1471 T TMP, RES0, RES1, RES2, TMP1, x, y;
1472
1473 bigint TMP_1, TMP_2;
1474 T *tmp = A.value[startr];
1475
1476 // Init
1477 lidia_size_t index;
1478 for (index = startc; tmp[index] == A.Zero && index >= 0; index--);
1479 if (index < 0)
1480 return false;
1481 if (index != startc) {
1482 this->swap_columns(A, startc, index);
1483 TR.swap_columns(startc, index);
1484 }
1485
1486 // sort
1487 if (startc >= 2) {
1488 lidia_size_t gcdindex = startc - 1;
1489 T gcd_min = gcd(tmp[startc], tmp[startc - 1]), gcd_tmp;
1490 for (lidia_size_t i = startc - 2; i >= 0; i--) {
1491 gcd_tmp = gcd(tmp[index], tmp[i]);
1492 if (gcd_tmp < gcd_min)
1493 gcdindex = i;
1494 }
1495 this->swap_columns(A, gcdindex, startc - 1);
1496 TR.swap_columns(gcdindex, startc - 1);
1497 }
1498
1499
1500 for (lidia_size_t i = startc - 1; i >= 0; i--)
1501 if (tmp[i] != A.Zero) {
1502 RES2 = xgcd(RES0, RES1, tmp[i], tmp[startc]);
1503 x = tmp[startc] / RES2;
1504 y = tmp[i] / RES2;
1505
1506 for (lidia_size_t l = 0; l <= len; l++) {
1507 TMP = A.value[l][i];
1508 TMP1 = A.value[l][startc];
1509
1510 A.value[l][i] = TMP*x - TMP1*y;
1511 A.value[l][startc] = TMP*RES0 + TMP1*RES1;
1512 }
1513
1514 for (lidia_size_t l = 0; l < TR.get_no_of_rows(); l++) {
1515 TMP_1 = TR.member(l, i);
1516 TMP_2 = TR.member(l, startc);
1517
1518 TR.sto(l, i, (TMP_1*x - TMP_2*y));
1519 TR.sto(l, startc, (TMP_1*RES0 + TMP_2*RES1));
1520 }
1521 }
1522 return true;
1523 }
1524
1525
1526
1527 template< class T >
1528 bool
mgcd_storjohann(MR<T> & A,lidia_size_t startr,lidia_size_t startc,lidia_size_t len) const1529 row_oriented_dense_matrix_modules< T >::mgcd_storjohann (MR< T > &A,
1530 lidia_size_t startr,
1531 lidia_size_t startc,
1532 lidia_size_t len) const
1533 {
1534 T q, res, TMP;
1535 T *tmp;
1536 T g, N, a, qa, ra, qb, rb, t;
1537 T A3, A2, A1, A0, x, y;
1538 T TMP1, TMP2, TMP3;
1539
1540 register lidia_size_t i, j, l;
1541 lidia_size_t index;
1542
1543 //
1544 // init
1545 //
1546
1547 tmp = A.value[startr];
1548 for (index = startc; index >= 0 && tmp[index].is_zero(); index--);
1549
1550 if (index < 0)
1551 return false;
1552 else {
1553 //
1554 // conditioning
1555 //
1556
1557 N = tmp[index];
1558 for (j = index - 1; j >= 0 && tmp[j].is_zero(); j--);
1559 if (j >= 0) {
1560 a = tmp[j];
1561 t = 0;
1562
1563 for (i = j - 1; i >= 0; i--) {
1564 if (tmp[i] != 0) {
1565 g = gcd(tmp[i], a);
1566
1567 div_rem(qa, ra, a/g, N);
1568 div_rem(qb, rb, tmp[i]/g, N);
1569
1570 for (t = 0; gcd((ra + t*rb), N) != 1; t++);
1571
1572 a = a + t * tmp[i]; // NEW
1573
1574 //M.sto(1, i, t);
1575 for (l = 0; l <= startr; l++)
1576 LiDIA::add(A.value[l][j], A.value[l][j], t*A.value[l][i]);
1577 }
1578 }
1579
1580 //
1581 // elimination
1582 //
1583
1584 A2 = xgcd(A0, A1, tmp[j], tmp[index]);
1585 div_rem(x, A3, tmp[index], A2);
1586 div_rem(y, A3, tmp[j], A2);
1587
1588 for (l = 0; l <= startr; l++) {
1589 TMP = A.value[l][j];
1590 TMP1 = A.value[l][index];
1591
1592
1593 // Atmp1[j] = ((TMP * x) + (TMP1 * y))
1594 LiDIA::multiply(TMP2, TMP, x);
1595 LiDIA::multiply(TMP3, TMP1, y);
1596 LiDIA::subtract(A.value[l][j], TMP2, TMP3);
1597
1598 // Atmp1[n-m+i] = ((TMP * A0) + (TMP1 * A1))
1599 LiDIA::multiply(TMP2, TMP, A0);
1600 LiDIA::multiply(TMP3, TMP1, A1);
1601 LiDIA::add(A.value[l][index], TMP2, TMP3);
1602 }
1603
1604 for (i = j - 1; i >= 0; i--)
1605 if (!tmp[i].is_zero()) {
1606 div_rem(TMP, TMP1, tmp[i], tmp[index]);
1607 for (l = 0; l <= startr; l++)
1608 LiDIA::subtract(A.value[l][i], A.value[l][i], TMP*A.value[l][index]);
1609 }
1610
1611 }
1612
1613 if (tmp[index].is_lt_zero())
1614 for (i = 0; i <= startr; i++)
1615 A.value[i][index].negate();
1616 if (index != startc)
1617 this->swap_columns(A, startc, index);
1618 }
1619 return true;
1620 }
1621
1622
1623
1624 template< class T >
1625 bool
mgcd_storjohann(MR<T> & A,matrix<bigint> & TR,lidia_size_t startr,lidia_size_t startc,lidia_size_t len) const1626 row_oriented_dense_matrix_modules< T >::mgcd_storjohann (MR< T > &A,
1627 matrix< bigint > &TR,
1628 lidia_size_t startr,
1629 lidia_size_t startc,
1630 lidia_size_t len) const
1631 {
1632 T q, res, TMP;
1633 T *tmp;
1634 T g, N, a, qa, ra, qb, rb, t;
1635 T A3, A2, A1, A0, x, y;
1636 T TMP1, TMP2, TMP3;
1637
1638 register lidia_size_t i, j, l;
1639 lidia_size_t index;
1640
1641 //
1642 // init
1643 //
1644
1645 tmp = A.value[startr];
1646 for (index = startc; index >= 0 && tmp[index].is_zero(); index--);
1647
1648 if (index < 0)
1649 return false;
1650 else {
1651 //
1652 // conditioning
1653 //
1654
1655 N = tmp[index];
1656 for (j = index - 1; j >= 0 && tmp[j].is_zero(); j--);
1657 if (j >= 0) {
1658 a = tmp[j];
1659 t = 0;
1660
1661 for (i = j - 1; i >= 0; i--) {
1662 if (tmp[i] != 0) {
1663 g = gcd(tmp[i], a);
1664
1665 div_rem(qa, ra, a/g, N);
1666 div_rem(qb, rb, tmp[i]/g, N);
1667
1668 for (t = 0; gcd((ra + t*rb), N) != 1; t++);
1669
1670 a = a + t * tmp[i]; // NEW
1671
1672 //M.sto(1, i, t);
1673 for (l = 0; l <= startr; l++)
1674 LiDIA::add(A.value[l][j], A.value[l][j], t*A.value[l][i]);
1675 for (l = 0; l < A.columns; l++)
1676 TR.sto(l, j, TR.member(l, j) + t * TR.member(l, i));
1677 }
1678 }
1679
1680 //
1681 // elimination
1682 //
1683
1684 A2 = xgcd(A0, A1, tmp[j], tmp[index]);
1685 div_rem(x, A3, tmp[index], A2);
1686 div_rem(y, A3, tmp[j], A2);
1687
1688 for (l = 0; l <= startr; l++) {
1689 TMP = A.value[l][j];
1690 TMP1 = A.value[l][index];
1691
1692
1693 // Atmp1[j] = ((TMP * x) + (TMP1 * y))
1694 LiDIA::multiply(TMP2, TMP, x);
1695 LiDIA::multiply(TMP3, TMP1, y);
1696 LiDIA::subtract(A.value[l][j], TMP2, TMP3);
1697
1698 // Atmp1[n-m+i] = ((TMP * A0) + (TMP1 * A1))
1699 LiDIA::multiply(TMP2, TMP, A0);
1700 LiDIA::multiply(TMP3, TMP1, A1);
1701 LiDIA::add(A.value[l][index], TMP2, TMP3);
1702 }
1703
1704 bigint TTMP, TTMP1, TTMP2, TTMP3;
1705 for (l = 0; l < A.columns; l++) {
1706 TTMP = TR.member(l, j);
1707 TTMP1 = TR.member(l, index);
1708
1709
1710 // Atmp1[j] = ((TMP * x) + (TMP1 * y))
1711 LiDIA::multiply(TTMP2, TTMP, x);
1712 LiDIA::multiply(TTMP3, TTMP1, y);
1713 TR.sto(l, j, TTMP2 - TTMP3);
1714
1715 // Atmp1[n-m+i] = ((TMP * A0) + (TMP1 * A1))
1716 LiDIA::multiply(TTMP2, TTMP, A0);
1717 LiDIA::multiply(TTMP3, TTMP1, A1);
1718 TR.sto(l, index, TTMP2 + TTMP3);
1719 }
1720
1721 for (i = j - 1; i >= 0; i--)
1722 if (!tmp[i].is_zero()) {
1723 div_rem(TMP, TMP1, tmp[i], tmp[index]);
1724 for (l = 0; l <= startr; l++)
1725 LiDIA::subtract(A.value[l][i], A.value[l][i], TMP*A.value[l][index]);
1726 for (l = 0; l < A.columns; l++)
1727 TR.sto(l, i, TR.member(l, i) - TMP * TR.member(l, index));
1728 }
1729
1730 }
1731
1732 if (tmp[index].is_lt_zero()) {
1733 for (i = 0; i <= startr; i++)
1734 A.value[i][index].negate();
1735 for (i = 0; i < A.columns; i++)
1736 TR.sto(i, index, -TR.member(i, index));
1737 }
1738
1739
1740 if (index != startc) {
1741 this->swap_columns(A, startc, index);
1742 TR.swap_columns(startc, index);
1743 }
1744 }
1745 return true;
1746 }
1747
1748
1749
1750 template< class T >
1751 bool
xgcd_elimination(MR<T> & A,lidia_size_t startr,lidia_size_t startc,lidia_size_t index,lidia_size_t len) const1752 row_oriented_dense_matrix_modules< T >::xgcd_elimination (MR< T > &A,
1753 lidia_size_t startr,
1754 lidia_size_t startc,
1755 lidia_size_t index,
1756 lidia_size_t len) const
1757 {
1758 T TMP;
1759 T RES0, RES1, RES2, TMP1, x, y;
1760 T TMP2, TMP3;
1761
1762 T *tmp = A.value[startr];
1763
1764 for (lidia_size_t i = 0; i <= startc; i++)
1765 if ((i != index) && tmp[i] != A.Zero) {
1766 RES2 = xgcd(RES0, RES1, tmp[i], tmp[index]);
1767 x = tmp[index] / RES2;
1768 y = tmp[i] / RES2;
1769
1770 for (lidia_size_t l = 0; l <= len; l++) {
1771 TMP = A.value[l][i];
1772 TMP1 = A.value[l][index];
1773
1774 LiDIA::multiply(TMP2, TMP, x);
1775 LiDIA::multiply(TMP3, TMP1, y);
1776 LiDIA::subtract(A.value[l][i], TMP2, TMP3);
1777
1778 LiDIA::multiply(TMP2, TMP, RES0);
1779 LiDIA::multiply(TMP3, TMP1, RES1);
1780 LiDIA::add(A.value[l][index], TMP2, TMP3);
1781 }
1782 }
1783 return true;
1784 }
1785
1786
1787
1788 template< class T >
1789 bool
xgcd_elimination(MR<T> & A,matrix<bigint> & TR,lidia_size_t startr,lidia_size_t startc,lidia_size_t index,lidia_size_t len) const1790 row_oriented_dense_matrix_modules< T >::xgcd_elimination (MR< T > &A,
1791 matrix< bigint > &TR,
1792 lidia_size_t startr,
1793 lidia_size_t startc,
1794 lidia_size_t index,
1795 lidia_size_t len) const
1796 {
1797 T TMP;
1798 T RES0, RES1, RES2, TMP1, x, y;
1799 T TMP2, TMP3;
1800
1801 T *tmp = A.value[startr];
1802
1803 for (lidia_size_t i = 0; i <= startc; i++)
1804 if ((i != index) && tmp[i] != A.Zero) {
1805 RES2 = xgcd(RES0, RES1, tmp[i], tmp[index]);
1806 x = tmp[index] / RES2;
1807 y = tmp[i] / RES2;
1808
1809 for (lidia_size_t l = 0; l <= len; l++) {
1810 TMP = A.value[l][i];
1811 TMP1 = A.value[l][index];
1812
1813 LiDIA::multiply(TMP2, TMP, x);
1814 LiDIA::multiply(TMP3, TMP1, y);
1815 LiDIA::subtract(A.value[l][i], TMP2, TMP3);
1816
1817 LiDIA::multiply(TMP2, TMP, RES0);
1818 LiDIA::multiply(TMP3, TMP1, RES1);
1819 LiDIA::add(A.value[l][index], TMP2, TMP3);
1820 }
1821
1822 bigint TTMP, TTMP1, TTMP2, TTMP3, TTMP4;
1823 for (lidia_size_t l = 0; l < A.columns; l++) {
1824 TTMP = TR.member(l, i);
1825 TTMP1 = TR.member(l, index);
1826
1827 LiDIA::multiply(TTMP2, TTMP, x);
1828 LiDIA::multiply(TTMP3, TTMP1, y);
1829 TR.sto(l, i, TTMP2 - TTMP3);
1830
1831 LiDIA::multiply(TTMP2, TTMP, RES0);
1832 LiDIA::multiply(TTMP3, TTMP1, RES1);
1833 TR.sto(l, index, TTMP2 + TTMP3);
1834 }
1835 }
1836 return true;
1837 }
1838
1839
1840
1841 template< class T >
1842 bool
xgcd_elimination_mod(MR<T> & A,lidia_size_t startr,lidia_size_t startc,lidia_size_t index,lidia_size_t len,const T & mod) const1843 row_oriented_dense_matrix_modules< T >::xgcd_elimination_mod (MR< T > &A,
1844 lidia_size_t startr,
1845 lidia_size_t startc,
1846 lidia_size_t index,
1847 lidia_size_t len,
1848 const T &mod) const
1849 {
1850 T TMP;
1851 T RES0, RES1, RES2, TMP1, x, y;
1852 T TMP2, TMP3;
1853
1854 T *tmp = A.value[startr];
1855
1856 for (lidia_size_t i = 0; i <= startc; i++)
1857 if ((i != index) && tmp[i] != A.Zero) {
1858 RES2 = xgcd(RES0, RES1, tmp[i], tmp[index]);
1859 x = tmp[index] / RES2;
1860 y = tmp[i] / RES2;
1861
1862 for (lidia_size_t l = 0; l <= len; l++) {
1863 TMP = A.value[l][i];
1864 TMP1 = A.value[l][index];
1865
1866 LiDIA::mult_mod(TMP2, TMP, x, mod);
1867 LiDIA::mult_mod(TMP3, TMP1, y, mod);
1868 LiDIA::sub_mod(A.value[l][i], TMP2, TMP3, mod);
1869
1870 LiDIA::mult_mod(TMP2, TMP, RES0, mod);
1871 LiDIA::mult_mod(TMP3, TMP1, RES1, mod);
1872 LiDIA::add_mod(A.value[l][index], TMP2, TMP3, mod);
1873 }
1874 }
1875 return true;
1876 }
1877
1878
1879
1880 template< class T >
1881 inline bool
normalize_row(MR<T> & A,lidia_size_t startr,lidia_size_t startc,lidia_size_t index,lidia_size_t len) const1882 row_oriented_dense_matrix_modules< T >::normalize_row (MR< T > &A,
1883 lidia_size_t startr,
1884 lidia_size_t startc,
1885 lidia_size_t index,
1886 lidia_size_t len) const
1887 {
1888 T q, TMP;
1889 for (lidia_size_t i = 0; i <= startc; i++)
1890 if ((i != index) && this->member(A, startr, i) != A.Zero) {
1891 pos_div_rem(q, TMP, this->member(A, startr, i), this->member(A, startr, index));
1892 if (q != 0)
1893 subtract_multiple_of_column(A, i, q, index, len);
1894 }
1895 return true;
1896 }
1897
1898
1899
1900 template< class T >
1901 inline bool
normalize_row_mod(MR<T> & A,lidia_size_t startr,lidia_size_t startc,lidia_size_t index,lidia_size_t len,const T & DET) const1902 row_oriented_dense_matrix_modules< T >::normalize_row_mod (MR< T > &A,
1903 lidia_size_t startr,
1904 lidia_size_t startc,
1905 lidia_size_t index,
1906 lidia_size_t len,
1907 const T &DET) const
1908 {
1909 T q, TMP;
1910 for (lidia_size_t i = 0; i <= startc; i++)
1911 if ((i != index) && this->member(A, startr, i) != A.Zero) {
1912 pos_div_rem(q, TMP, this->member(A, startr, i), this->member(A, startr, index));
1913 if (q != 0)
1914 subtract_multiple_of_column_mod(A, i, q, index, len, DET);
1915 }
1916 return true;
1917 }
1918
1919
1920
1921 template< class T >
1922 inline bool
normalize_row(MR<T> & A,matrix<bigint> & TR,lidia_size_t startr,lidia_size_t startc,lidia_size_t index,lidia_size_t len) const1923 row_oriented_dense_matrix_modules< T >::normalize_row (MR< T > &A,
1924 matrix< bigint > &TR,
1925 lidia_size_t startr,
1926 lidia_size_t startc,
1927 lidia_size_t index,
1928 lidia_size_t len) const
1929 {
1930 T q, TMP;
1931 for (lidia_size_t i = 0; i <= startc; i++)
1932 if ((i != index) && this->member(A, startr, i) != A.Zero) {
1933 pos_div_rem(q, TMP, this->member(A, startr, i), this->member(A, startr, index));
1934
1935 // Check
1936 if (q != 0) {
1937 subtract_multiple_of_column(A, i, q, index, len);
1938
1939 for (lidia_size_t j = 0; j < A.columns; j++)
1940 TR.sto(j, i, TR.member(j, i) - q * TR.member(j, index));
1941 }
1942 }
1943 return true;
1944 }
1945
1946
1947
1948 template< class T >
1949 inline bool
normalize_row(MR<T> & A,lidia_size_t startr,lidia_size_t startc,lidia_size_t index,lidia_size_t len,T * max_array,const T & BOUND) const1950 row_oriented_dense_matrix_modules< T >::normalize_row (MR< T > &A,
1951 lidia_size_t startr,
1952 lidia_size_t startc,
1953 lidia_size_t index,
1954 lidia_size_t len,
1955 T *max_array,
1956 const T &BOUND) const
1957 {
1958 T q, TMP;
1959 for (lidia_size_t i = 0; i <= startc; i++)
1960 if ((i != index) && this->member(A, startr, i) != A.Zero) {
1961 pos_div_rem(q, TMP, this->member(A, startr, i), this->member(A, startr, index));
1962
1963 // Check
1964 if (q != 0) {
1965 if (abs(bigint(bigint(max_array[index]) * bigint(q)) + bigint(max_array[i])) > bigint(BOUND))
1966 return false;
1967
1968 subtract_multiple_of_column(A, i, q, index, len);
1969 update_max_array(A, i, max_array);
1970 }
1971 }
1972 return true;
1973 }
1974
1975
1976
1977 template< class T >
1978 bool
normalize_row(MR<T> & A,matrix<bigint> & TR,lidia_size_t startr,lidia_size_t startc,lidia_size_t index,lidia_size_t len,T * max_array,const T & BOUND) const1979 row_oriented_dense_matrix_modules< T >::normalize_row (MR< T > &A,
1980 matrix< bigint > &TR,
1981 lidia_size_t startr,
1982 lidia_size_t startc,
1983 lidia_size_t index,
1984 lidia_size_t len,
1985 T *max_array,
1986 const T &BOUND) const
1987 {
1988 T q, TMP;
1989 for (lidia_size_t i = 0; i <= startc; i++)
1990 if ((i != index) && this->member(A, startr, i) != A.Zero) {
1991 pos_div_rem(q, TMP, this->member(A, startr, i), this->member(A, startr, index));
1992
1993 // Check
1994 if (q != 0) {
1995 if (abs(bigint(bigint(max_array[index]) * q) + max_array[i]) > bigint(BOUND))
1996 return false;
1997 subtract_multiple_of_column(A, i, q, index, len);
1998
1999 for (lidia_size_t j = 0; j < A.columns; j++)
2000 TR.sto(j, i, TR.member(j, i) - q* TR.member(j, index));
2001
2002 update_max_array(A, i, max_array);
2003 }
2004 }
2005 return true;
2006 }
2007
2008
2009
2010 //
2011 // special functions
2012 //
2013
2014 template< class T >
2015 inline const T &
member(const MR<T> & A,lidia_size_t x,lidia_size_t y) const2016 column_oriented_dense_matrix_modules< T >::member (const MR< T > &A,
2017 lidia_size_t x,
2018 lidia_size_t y) const
2019 {
2020 return A.value[y][x];
2021 }
2022
2023
2024
2025 template< class T >
2026 inline void
swap_columns(const MR<T> & A,lidia_size_t l1,lidia_size_t l2) const2027 column_oriented_dense_matrix_modules< T >::swap_columns (const MR< T > &A,
2028 lidia_size_t l1,
2029 lidia_size_t l2) const
2030 {
2031 T *tmp = A.value[l1];
2032 A.value[l1] = A.value[l2];
2033 A.value[l2] = tmp;
2034 }
2035
2036
2037
2038 template< class T >
2039 inline void
subtract_multiple_of_column(MR<T> & A,lidia_size_t l1,const T & q,lidia_size_t l2,lidia_size_t l) const2040 column_oriented_dense_matrix_modules< T >::subtract_multiple_of_column (MR< T > &A,
2041 lidia_size_t l1,
2042 const T &q,
2043 lidia_size_t l2,
2044 lidia_size_t l) const
2045 {
2046 for (register lidia_size_t i = 0; i <= l; i++)
2047 LiDIA::subtract(A.value[l1][i], A.value[l1][i], q*A.value[l2][i]);
2048 }
2049
2050
2051
2052 template< class T >
2053 inline void
subtract_multiple_of_column_mod(MR<T> & A,lidia_size_t l1,const T & q,lidia_size_t l2,lidia_size_t l,const T & mod) const2054 column_oriented_dense_matrix_modules< T >::subtract_multiple_of_column_mod (MR< T > &A,
2055 lidia_size_t l1,
2056 const T &q,
2057 lidia_size_t l2,
2058 lidia_size_t l,
2059 const T &mod) const
2060 {
2061 T TMP;
2062 for (register lidia_size_t i = 0; i <= l; i++) {
2063 LiDIA::mult_mod(TMP, q, A.value[l2][i], mod);
2064 LiDIA::sub_mod(A.value[l1][i], A.value[l1][i], TMP, mod);
2065 }
2066 }
2067
2068
2069
2070 template< class T >
2071 inline void
normalize_column_mod(MR<T> & A,lidia_size_t l1,const T & q,lidia_size_t l2,lidia_size_t l,const T & mod) const2072 column_oriented_dense_matrix_modules< T >::normalize_column_mod (MR< T > &A,
2073 lidia_size_t l1,
2074 const T &q,
2075 lidia_size_t l2,
2076 lidia_size_t l,
2077 const T &mod) const
2078 {
2079 T TMP;
2080 for (register lidia_size_t i = 0; i <= l; i++) {
2081 LiDIA::mult_mod(TMP, q, A.value[l2][i], mod);
2082 LiDIA::sub_mod(A.value[l1][i], A.value[l1][i], TMP, mod);
2083 if (A.value[l1][i] < A.Zero)
2084 A.value[l1][i] += mod;
2085 }
2086 }
2087
2088
2089
2090 template< class T >
2091 inline void
negate_column(MR<T> & A,lidia_size_t index,lidia_size_t l) const2092 column_oriented_dense_matrix_modules< T >::negate_column (MR< T > &A,
2093 lidia_size_t index,
2094 lidia_size_t l) const
2095 {
2096 for (register lidia_size_t i = 0; i <= l; i++)
2097 A.value[index][i] = -A.value[index][i];
2098 }
2099
2100
2101
2102 template< class T >
2103 inline void
negate_column_mod(MR<T> & A,lidia_size_t index,lidia_size_t l,const T & mod) const2104 column_oriented_dense_matrix_modules< T >::negate_column_mod (MR< T > &A,
2105 lidia_size_t index,
2106 lidia_size_t l,
2107 const T &mod) const
2108 {
2109 for (register lidia_size_t i = 0; i <= l; i++) {
2110 A.value[index][i] = -A.value[index][i];
2111 LiDIA::best_remainder(A.value[index][i], A.value[index][i], mod);
2112 }
2113 }
2114
2115
2116
2117 template< class T >
2118 inline T *
init_max_array(const MR<T> & A) const2119 column_oriented_dense_matrix_modules< T >::init_max_array (const MR< T > &A) const
2120 {
2121 T *max_array = new T[A.columns];
2122 memory_handler(max_array, DMESSAGE, "init :: "
2123 "Error in memory allocation (max_array)");
2124 T *tmp;
2125 for (register lidia_size_t i = 0; i < A.columns; i++) {
2126 tmp = A.value[i];
2127 for (register lidia_size_t j = 0; j < A.rows; j++)
2128 if (max_array[j] < abs(tmp[j]))
2129 max_array[j] = abs(tmp[j]);
2130 }
2131 return max_array;
2132 }
2133
2134
2135
2136 template< class T >
2137 inline void
update_max_array(const MR<T> & A,lidia_size_t i,T * max_array) const2138 column_oriented_dense_matrix_modules< T >::update_max_array (const MR< T > &A,
2139 lidia_size_t i,
2140 T *max_array) const
2141 {
2142 for (register lidia_size_t j = 0; j < A.rows; j++)
2143 if (max_array[i] < abs(A.value[i][j]))
2144 max_array[i] = abs(A.value[i][j]);
2145 }
2146
2147
2148
2149 //
2150 // Pivot search
2151 //
2152
2153 template< class T >
2154 inline lidia_size_t
normal_pivot(const MR<T> & A,lidia_size_t startr,lidia_size_t startc,lidia_size_t & index) const2155 column_oriented_dense_matrix_modules< T >::normal_pivot (const MR< T > &A,
2156 lidia_size_t startr,
2157 lidia_size_t startc,
2158 lidia_size_t &index) const
2159 {
2160 lidia_size_t FOUND = 0;
2161
2162 for (lidia_size_t i = 0; i <= startc; i++)
2163 if (A.value[i][startr] != A.Zero) {
2164 FOUND++;
2165 index = i;
2166 for (lidia_size_t j = 0; j <= startc; j++)
2167 if (A.value[j][startr] != A.Zero && (j != index)) {
2168 FOUND++;
2169 if (j != index) {
2170 if (abs(A.value[j][startr]) < abs(A.value[i][startr]))
2171 index = j;
2172 break;
2173 }
2174 }
2175 break;
2176 }
2177 return FOUND;
2178 }
2179
2180
2181
2182 template< class T >
2183 inline lidia_size_t
min_abs_of_row(const MR<T> & A,lidia_size_t startr,lidia_size_t startc,lidia_size_t & pos) const2184 column_oriented_dense_matrix_modules< T >::min_abs_of_row (const MR< T > &A,
2185 lidia_size_t startr,
2186 lidia_size_t startc,
2187 lidia_size_t &pos) const
2188 {
2189 lidia_size_t FOUND = 0;
2190 for (; startc >= 0 && !FOUND; startc--)
2191 if (A.value[startc][startr] != A.Zero) {
2192 pos = startc;
2193 FOUND++;
2194 }
2195
2196 for (; startc >= 0; startc--)
2197 if (A.value[startc][startr] != A.Zero) {
2198 FOUND++;
2199 if (abs(A.value[pos][startr]) >= abs(A.value[startc][startr]))
2200 pos = startc;
2201 }
2202
2203 return FOUND;
2204 }
2205
2206
2207
2208 template< class T >
2209 lidia_size_t
pivot_sorting_gcd(const MR<T> & A,lidia_size_t startr,lidia_size_t startc,lidia_size_t & index) const2210 column_oriented_dense_matrix_modules< T >::pivot_sorting_gcd (const MR< T > &A,
2211 lidia_size_t startr,
2212 lidia_size_t startc,
2213 lidia_size_t &index) const
2214 {
2215 lidia_size_t FOUND = 0;
2216 lidia_size_t index_largest = 0;
2217 T val, maxval;
2218
2219 for (lidia_size_t i = 0; i <= startc; i++) {
2220 val = abs(A.value[i][startr]);
2221 if (val != A.Zero) {
2222 if (FOUND == 0) {
2223 maxval = val;
2224 index_largest = i;
2225 }
2226 else
2227 if (maxval < val) {
2228 maxval = val;
2229 index_largest = i;
2230 }
2231
2232 FOUND++;
2233 }
2234 }
2235
2236 bool SW = false;
2237 if (FOUND > 1) {
2238 for (lidia_size_t i = 0; i <= startc; i++) {
2239 val = abs(A.value[i][startr]);
2240 if (val != A.Zero && index_largest != i) {
2241 if (SW == false) {
2242 maxval = val;
2243 index = i;
2244 }
2245 else
2246 if (maxval < val) {
2247 maxval = val;
2248 index = i;
2249 }
2250
2251 SW = true;
2252 }
2253 }
2254 }
2255 else
2256 index = index_largest;
2257
2258 return FOUND;
2259 }
2260
2261
2262
2263 template< class T >
2264 lidia_size_t
minimal_norm(const MR<T> & A,lidia_size_t startr,lidia_size_t startc,lidia_size_t & index,lidia_size_t norm) const2265 column_oriented_dense_matrix_modules< T >::minimal_norm (const MR< T > &A,
2266 lidia_size_t startr,
2267 lidia_size_t startc,
2268 lidia_size_t &index,
2269 lidia_size_t norm) const
2270 {
2271 lidia_size_t i, j, FOUND = 0;
2272 bigint min_norm, nnorm, TMP;
2273 T val, maxval, minval;
2274
2275 bool Pair = false;
2276 bool SW = false;
2277
2278 for (i = 0; i <= startc; i++) {
2279 val = abs(A.value[i][startr]);
2280 if (val != A.Zero) {
2281 FOUND++;
2282 if (SW == false) {
2283 maxval = val;
2284 SW = true;
2285 index = i;
2286 }
2287 else {
2288 if (val > maxval) {
2289 maxval = val;
2290 Pair = false;
2291 }
2292 else if (maxval == val)
2293 Pair = true;
2294 }
2295 }
2296 }
2297
2298 SW = false;
2299
2300 for (i = 0; i <= startc; i++) {
2301 val = abs(A.value[i][startr]);
2302 if (val != A.Zero) {
2303 if (SW == false) {
2304 if (Pair || maxval != val) {
2305 min_norm = column_norm(A, i, norm);
2306 index = i;
2307 SW = true;
2308 }
2309 }
2310 else {
2311 if (Pair) {
2312 nnorm = column_norm(A, i, norm);
2313 if (min_norm > nnorm) {
2314 min_norm = nnorm;
2315 index = i;
2316 }
2317 }
2318 else {
2319 if (val != maxval) {
2320 nnorm = column_norm(A, i, norm);
2321 if (min_norm > nnorm) {
2322 min_norm = nnorm;
2323 index = i;
2324 }
2325 }
2326 }
2327 }
2328 }
2329 }
2330 return FOUND;
2331 }
2332
2333
2334
2335 template< class T >
2336 lidia_size_t
min_abs_of_row_plus_minimal_norm(const MR<T> & A,lidia_size_t startr,lidia_size_t startc,lidia_size_t & index,lidia_size_t norm) const2337 column_oriented_dense_matrix_modules< T >::min_abs_of_row_plus_minimal_norm (const MR< T > &A,
2338 lidia_size_t startr,
2339 lidia_size_t startc,
2340 lidia_size_t &index,
2341 lidia_size_t norm) const
2342 {
2343 lidia_size_t i, j, FOUND = 0;
2344 bigint min_norm, nnorm, TMP;
2345 T min_val = 0, val;
2346 bool SW = false;
2347 bool Pair = false;
2348
2349 for (i = startc; i >= 0; i--) {
2350 val = abs(A.value[i][startr]);
2351
2352 if (val != A.Zero) {
2353 FOUND++;
2354 if (SW == false) {
2355 index = i;
2356 min_val = val;
2357 SW = true;
2358 }
2359 else {
2360 if (min_val > val) {
2361 index = i;
2362 min_val = val;
2363 }
2364 }
2365 }
2366 }
2367
2368 min_norm = column_norm(A, index, norm);
2369
2370 // norm computation
2371 for (i = index - 1; i >= 0; i--) {
2372 if (min_val == abs(A.value[i][startr])) {
2373 nnorm = column_norm(A, i, norm);
2374 if (nnorm < min_norm) {
2375 index = i;
2376 min_norm = nnorm;
2377 }
2378 }
2379 }
2380 return FOUND;
2381 }
2382
2383
2384
2385 template< class T >
2386 lidia_size_t
minimal_norm_plus_min_abs_of_row(const MR<T> & A,lidia_size_t startr,lidia_size_t startc,lidia_size_t & index,lidia_size_t norm) const2387 column_oriented_dense_matrix_modules< T >::minimal_norm_plus_min_abs_of_row (const MR< T > &A,
2388 lidia_size_t startr,
2389 lidia_size_t startc,
2390 lidia_size_t &index,
2391 lidia_size_t norm) const
2392 {
2393 lidia_size_t i, j, FOUND = 0;
2394 bigint min_norm, nnorm, TMP;
2395 T val, maxval, minval;
2396
2397 bool SW = false;
2398 bool Pair = false;
2399
2400 for (i = 0; i <= startc; i++) {
2401 val = abs(A.value[i][startr]);
2402 if (val != A.Zero) {
2403 FOUND++;
2404 if (SW == false) {
2405 maxval = val;
2406 SW = true;
2407 index = i;
2408 }
2409 else {
2410 if (val > maxval) {
2411 maxval = val;
2412 Pair = false;
2413 }
2414 else if (maxval == val)
2415 Pair = true;
2416 }
2417 }
2418 }
2419
2420 SW = false;
2421
2422 for (i = 0; i <= startc; i++) {
2423 val = abs(A.value[i][startr]);
2424 if (val != A.Zero) {
2425 if (SW == false) {
2426 if (Pair || val != maxval) {
2427 minval = val;
2428 min_norm = column_norm(A, i, norm);
2429 index = i;
2430 SW = true;
2431 }
2432 }
2433 else {
2434 if (Pair) {
2435 nnorm = column_norm(A, i, norm);
2436 if (min_norm > nnorm) {
2437 min_norm = nnorm;
2438 index = i;
2439 }
2440 else if (min_norm == nnorm)
2441 if (minval > val) {
2442 minval = val;
2443 min_norm = nnorm;
2444 index = i;
2445 }
2446 }
2447 else {
2448 if (val != maxval) {
2449 nnorm = column_norm(A, i, norm);
2450 if (min_norm > nnorm) {
2451 min_norm = nnorm;
2452 index = i;
2453 }
2454 else if (min_norm == nnorm)
2455 if (minval > val) {
2456 minval = val;
2457 min_norm = nnorm;
2458 index = i;
2459 }
2460 }
2461 }
2462 }
2463 }
2464 }
2465 return FOUND;
2466 }
2467
2468
2469
2470 template< class T >
2471 lidia_size_t
sorting_gcd_plus_minimal_norm(const MR<T> & A,lidia_size_t startr,lidia_size_t startc,lidia_size_t & index,lidia_size_t norm) const2472 column_oriented_dense_matrix_modules< T >::sorting_gcd_plus_minimal_norm (const MR< T > &A,
2473 lidia_size_t startr,
2474 lidia_size_t startc,
2475 lidia_size_t &index,
2476 lidia_size_t norm) const
2477 {
2478 lidia_size_t FOUND = 0;
2479 lidia_size_t index_largest = 0;
2480 T val, maxval;
2481
2482 for (lidia_size_t i = 0; i <= startc; i++) {
2483 val = abs(A.value[i][startr]);
2484 if (val != A.Zero) {
2485 if (FOUND == 0) {
2486 maxval = val;
2487 index_largest = i;
2488 }
2489 else
2490 if (maxval < val) {
2491 maxval = val;
2492 index_largest = i;
2493 }
2494
2495 FOUND++;
2496 }
2497 }
2498
2499 bool SW = false;
2500 if (FOUND > 1) {
2501 for (lidia_size_t i = 0; i <= startc; i++) {
2502 val = abs(A.value[i][startr]);
2503 if (val != A.Zero && index_largest != i) {
2504 if (SW == false) {
2505 maxval = val;
2506 index = i;
2507 }
2508 else
2509 if (maxval < val) {
2510 maxval = val;
2511 index = i;
2512 }
2513
2514 SW = true;
2515 }
2516 }
2517 }
2518 else {
2519 index = index_largest;
2520 maxval = abs(A.value[index][startr]);
2521 }
2522
2523 bigint nnorm, min_norm = column_norm(A, index, norm);
2524
2525 for (lidia_size_t i = 0; i <= startc; i++)
2526 if (abs(A.value[i][startr]) == maxval) {
2527 nnorm = column_norm(A, i, norm);
2528 if (nnorm < min_norm) {
2529 min_norm = nnorm;
2530 index = i;
2531 }
2532 }
2533 return FOUND;
2534 }
2535
2536
2537
2538 //
2539 // Norm
2540 //
2541
2542 template< class T >
2543 inline bigint
column_norm(const MR<T> & A,lidia_size_t pos,lidia_size_t norm) const2544 column_oriented_dense_matrix_modules< T >::column_norm (const MR< T > &A,
2545 lidia_size_t pos,
2546 lidia_size_t norm) const
2547 {
2548 bigint Norm = A.Zero;
2549 T * tmp = A.value[pos];
2550 bigint TMP;
2551
2552 if (norm == 0) {
2553 for (lidia_size_t j = 0; j < A.rows; j++)
2554 if (tmp[j] != A.Zero)
2555 Norm++;
2556 }
2557 else if (norm == 1) {
2558 for (lidia_size_t j = 0; j < A.rows; j++)
2559 LiDIA::add(Norm, Norm, tmp[j]);
2560 }
2561 else {
2562 for (lidia_size_t j = 0; j < A.rows; j++) {
2563 power(TMP, tmp[j], norm);
2564 LiDIA::add(Norm, Norm, abs(TMP));
2565 }
2566 }
2567 return Norm;
2568 }
2569
2570
2571
2572 template< class T >
2573 inline void
kennwerte(MR<T> & RES,T & MAX,lidia_size_t & no_of_elements,T & Durch) const2574 column_oriented_dense_matrix_modules< T >::kennwerte (MR< T > &RES,
2575 T &MAX,
2576 lidia_size_t &no_of_elements,
2577 T &Durch) const
2578 {
2579 register lidia_size_t i, j;
2580 bool SW = false;
2581 no_of_elements = 0;
2582 Durch = 0;
2583
2584 for (i = 0; i < RES.columns; i++) {
2585 for (j = 0; j < RES.rows; j++) {
2586 if (RES.value[i][j] != RES.Zero) {
2587 no_of_elements++;
2588 Durch += abs(RES.value[i][j]);
2589 if (SW == false) {
2590 SW = true;
2591 MAX = abs(RES.value[i][j]);
2592 }
2593 else
2594 if (MAX < abs(RES.value[i][j]))
2595 MAX = abs(RES.value[i][j]);
2596 }
2597 }
2598 }
2599 Durch /= no_of_elements;
2600 }
2601
2602
2603
2604 template< class T >
2605 inline void
max(MR<T> & RES,T & MAX) const2606 column_oriented_dense_matrix_modules< T >::max (MR< T > &RES, T &MAX) const
2607 {
2608 register lidia_size_t i, j;
2609 bool SW = false;
2610
2611 for (i = 0; i < RES.columns; i++)
2612 for (j = 0; j < RES.rows; j++)
2613 if (SW == false) {
2614 SW = true;
2615 MAX = abs(RES.value[i][j]);
2616 }
2617 else
2618 if (MAX < abs(RES.value[i][j]))
2619 MAX = abs(RES.value[i][j]);
2620 }
2621
2622
2623
2624 //
2625 // Hermite normal form: elimination
2626 //
2627
2628 template< class T >
2629 inline bool
normalize_one_element_of_row(MR<T> & A,lidia_size_t startr,lidia_size_t startc,lidia_size_t index,lidia_size_t len) const2630 column_oriented_dense_matrix_modules< T >::normalize_one_element_of_row (MR< T > &A,
2631 lidia_size_t startr,
2632 lidia_size_t startc,
2633 lidia_size_t index,
2634 lidia_size_t len) const
2635 {
2636 T q, TMP;
2637 for (lidia_size_t i = 0; i <= startc; i++)
2638 if ((i != index) && member(A, startr, i) != A.Zero) {
2639 pos_div_rem(q, TMP, member(A, startr, i), member(A, startr, index));
2640 if (q != 0)
2641 subtract_multiple_of_column(A, i, q, index, len);
2642 break;
2643 }
2644 return true;
2645 }
2646
2647
2648
2649 template< class T >
2650 inline bool
normalize_one_element_of_row(MR<T> & A,matrix<bigint> & TR,lidia_size_t startr,lidia_size_t startc,lidia_size_t index,lidia_size_t len) const2651 column_oriented_dense_matrix_modules< T >::normalize_one_element_of_row (MR< T > &A,
2652 matrix< bigint > &TR,
2653 lidia_size_t startr,
2654 lidia_size_t startc,
2655 lidia_size_t index,
2656 lidia_size_t len) const
2657 {
2658 T q, TMP;
2659 for (lidia_size_t i = 0; i <= startc; i++)
2660 if ((i != index) && member(A, startr, i) != A.Zero) {
2661 pos_div_rem(q, TMP, member(A, startr, i), member(A, startr, index));
2662 if (q != 0) {
2663 subtract_multiple_of_column(A, i, q, index, len);
2664 subtract_multiple_of_column(TR, i, q, index, TR.get_no_of_rows() - 1);
2665 }
2666 break;
2667 }
2668 return true;
2669 }
2670
2671
2672
2673 template< class T >
2674 inline bool
normalize_one_element_of_row(MR<T> & A,math_vector<bigfloat> & vec,lidia_size_t startr,lidia_size_t startc,lidia_size_t index,lidia_size_t len) const2675 column_oriented_dense_matrix_modules< T >::normalize_one_element_of_row (MR< T > &A,
2676 math_vector< bigfloat > &vec,
2677 lidia_size_t startr,
2678 lidia_size_t startc,
2679 lidia_size_t index,
2680 lidia_size_t len) const
2681 {
2682 T q, TMP;
2683 bigfloat rtemp;
2684 for (lidia_size_t i = 0; i <= startc; i++)
2685 if ((i != index) && member(A, startr, i) != A.Zero) {
2686 pos_div_rem(q, TMP, member(A, startr, i), member(A, startr, index));
2687 if (q != 0) {
2688 subtract_multiple_of_column(A, i, q, index, len);
2689 multiply(rtemp, bigfloat(q), vec.member(index));
2690 subtract(vec[i], vec[i], rtemp);
2691 }
2692 break;
2693 }
2694 return true;
2695 }
2696
2697
2698
2699 template< class T >
2700 bool
mgcd_linear(MR<T> & A,lidia_size_t startr,lidia_size_t startc,lidia_size_t len) const2701 column_oriented_dense_matrix_modules< T >::mgcd_linear (MR< T > &A,
2702 lidia_size_t startr,
2703 lidia_size_t startc,
2704 lidia_size_t len) const
2705 {
2706 T TMP;
2707 T RES0, RES1, RES2, TMP1, x, y;
2708
2709 // Init
2710 lidia_size_t index;
2711 for (index = startc; A.value[index][startr] == A.Zero && index >= 0; index--);
2712 if (index < 0)
2713 return false;
2714 if (index != startc)
2715 swap_columns(A, startc, index);
2716
2717 for (lidia_size_t i = index - 1; i >= 0; i--)
2718 if (A.value[i][startr] != A.Zero) {
2719 RES2 = xgcd(RES0, RES1, A.value[i][startr], A.value[startc][startr]);
2720 x = A.value[startc][startr] / RES2;
2721 y = A.value[i][startr] / RES2;
2722
2723 for (lidia_size_t l = 0; l <= len; l++) {
2724 TMP = A.value[i][l];
2725 TMP1 = A.value[startc][l];
2726
2727 A.value[i][l] = TMP*x - TMP1*y;
2728 A.value[startc][l] = TMP*RES0 + TMP1*RES1;
2729 }
2730 }
2731 return true;
2732 }
2733
2734
2735
2736 template< class T >
2737 bool
mgcd_linear(MR<T> & A,matrix<bigint> & TR,lidia_size_t startr,lidia_size_t startc,lidia_size_t len) const2738 column_oriented_dense_matrix_modules< T >::mgcd_linear (MR< T > &A,
2739 matrix< bigint > &TR,
2740 lidia_size_t startr,
2741 lidia_size_t startc,
2742 lidia_size_t len) const
2743 {
2744 T TMP;
2745 T RES0, RES1, RES2, TMP1, TMP2, x, y;
2746
2747 // Init
2748 lidia_size_t index;
2749 for (index = startc; A.value[index][startr] == A.Zero && index >= 0; index--);
2750 if (index < 0)
2751 return false;
2752 if (index != startc) {
2753 swap_columns(A, startc, index);
2754 TR.swap_columns(startc, index);
2755 }
2756
2757 for (lidia_size_t i = index - 1; i >= 0; i--)
2758 if (A.value[i][startr] != A.Zero) {
2759 RES2 = xgcd(RES0, RES1, A.value[i][startr], A.value[startc][startr]);
2760 x = A.value[startc][startr] / RES2;
2761 y = A.value[i][startr] / RES2;
2762
2763 for (lidia_size_t l = 0; l <= len; l++) {
2764 TMP = A.value[i][l];
2765 TMP1 = A.value[startc][l];
2766
2767 A.value[i][l] = TMP*x - TMP1*y;
2768 A.value[startc][l] = TMP*RES0 + TMP1*RES1;
2769 }
2770
2771 for (lidia_size_t l = 0; l < TR.get_no_of_rows(); l++) {
2772 TMP1 = TR(l, i);
2773 TMP2 = TR(l, startc);
2774
2775 TR.sto(l, i, (TMP1*x - TMP2*y));
2776 TR.sto(l, startc, (TMP1*RES0 + TMP2*RES1));
2777 }
2778 }
2779 return true;
2780 }
2781
2782
2783
2784 template< class T >
2785 bool
xgcd_elimination_mod(MR<T> & A,lidia_size_t startr,lidia_size_t startc,lidia_size_t index,lidia_size_t len,const T & mod) const2786 column_oriented_dense_matrix_modules< T >::xgcd_elimination_mod (MR< T > &A,
2787 lidia_size_t startr,
2788 lidia_size_t startc,
2789 lidia_size_t index,
2790 lidia_size_t len,
2791 const T &mod) const
2792 {
2793 T TMP;
2794 T RES0, RES1, RES2, TMP1, x, y;
2795 T TMP2, TMP3;
2796
2797 for (lidia_size_t i = 0; i <= startc; i++)
2798 if ((i != index) && member(A, startr, i) != A.Zero) {
2799 RES2 = xgcd(RES0, RES1, member(A, startr, i), member(A, startr, index));
2800 x = member(A, startr, index) / RES2;
2801 y = member(A, startr, i) / RES2;
2802
2803 for (lidia_size_t l = 0; l <= len; l++) {
2804 TMP = member(A, l, i);
2805 TMP1 = member(A, l, index);
2806
2807 LiDIA::mult_mod(TMP2, TMP, x, mod);
2808 LiDIA::mult_mod(TMP3, TMP1, y, mod);
2809 LiDIA::sub_mod(A.value[i][l], TMP2, TMP3, mod);
2810
2811 LiDIA::mult_mod(TMP2, TMP, RES0, mod);
2812 LiDIA::mult_mod(TMP3, TMP1, RES1, mod);
2813 LiDIA::add_mod(A.value[index][l], TMP2, TMP3, mod);
2814 }
2815 }
2816 return true;
2817 }
2818
2819
2820
2821 template< class T >
2822 inline bool
normalize_row(MR<T> & A,lidia_size_t startr,lidia_size_t startc,lidia_size_t index,lidia_size_t len) const2823 column_oriented_dense_matrix_modules< T >::normalize_row (MR< T > &A,
2824 lidia_size_t startr,
2825 lidia_size_t startc,
2826 lidia_size_t index,
2827 lidia_size_t len) const
2828 {
2829 T q, TMP;
2830 for (lidia_size_t i = 0; i <= startc; i++)
2831 if ((i != index) && member(A, startr, i) != A.Zero) {
2832 pos_div_rem(q, TMP, member(A, startr, i), member(A, startr, index));
2833 if (q != 0)
2834 subtract_multiple_of_column(A, i, q, index, len);
2835 }
2836 return true;
2837 }
2838
2839
2840
2841 template< class T >
2842 inline bool
normalize_row_mod(MR<T> & A,lidia_size_t startr,lidia_size_t startc,lidia_size_t index,lidia_size_t len,const T & DET) const2843 column_oriented_dense_matrix_modules< T >::normalize_row_mod (MR< T > &A,
2844 lidia_size_t startr,
2845 lidia_size_t startc,
2846 lidia_size_t index,
2847 lidia_size_t len,
2848 const T &DET) const
2849 {
2850 T q, TMP;
2851 for (lidia_size_t i = 0; i <= startc; i++)
2852 if ((i != index) && member(A, startr, i) != A.Zero) {
2853 pos_div_rem(q, TMP, member(A, startr, i), member(A, startr, index));
2854 if (q != 0)
2855 subtract_multiple_of_column_mod(A, i, q, index, len, DET);
2856 }
2857 return true;
2858 }
2859
2860
2861
2862 template< class T >
2863 inline bool
normalize_row(MR<T> & A,matrix<bigint> & TR,lidia_size_t startr,lidia_size_t startc,lidia_size_t index,lidia_size_t len) const2864 column_oriented_dense_matrix_modules< T >::normalize_row (MR< T > &A,
2865 matrix< bigint > &TR,
2866 lidia_size_t startr,
2867 lidia_size_t startc,
2868 lidia_size_t index,
2869 lidia_size_t len) const
2870 {
2871 T q, TMP;
2872 for (lidia_size_t i = 0; i <= startc; i++)
2873 if ((i != index) && member(A, startr, i) != A.Zero) {
2874 pos_div_rem(q, TMP, member(A, startr, i), member(A, startr, index));
2875
2876 // Check
2877 if (q != 0) {
2878 subtract_multiple_of_column(A, i, q, index, len);
2879 subtract_multiple_of_column(TR, i, q, index, TR.get_no_of_rows() - 1);
2880 }
2881 }
2882 return true;
2883 }
2884
2885
2886
2887 template< class T >
2888 inline bool
normalize_row(MR<T> & A,lidia_size_t startr,lidia_size_t startc,lidia_size_t index,lidia_size_t len,T * max_array,const T & BOUND) const2889 column_oriented_dense_matrix_modules< T >::normalize_row (MR< T > &A,
2890 lidia_size_t startr,
2891 lidia_size_t startc,
2892 lidia_size_t index,
2893 lidia_size_t len,
2894 T *max_array,
2895 const T &BOUND) const
2896 {
2897 T q, TMP;
2898 for (lidia_size_t i = 0; i <= startc; i++)
2899 if ((i != index) && member(A, startr, i) != A.Zero) {
2900 pos_div_rem(q, TMP, member(A, startr, i), member(A, startr, index));
2901
2902 // Check
2903 if (q != 0) {
2904 if (abs(bigint(bigint(max_array[index]) * bigint(q)) + bigint(max_array[i])) > bigint(BOUND))
2905 return false;
2906
2907 subtract_multiple_of_column(A, i, q, index, len);
2908 update_max_array(A, i, max_array);
2909 }
2910 }
2911 return true;
2912 }
2913
2914
2915
2916 template< class T >
2917 inline bool
normalize_row(MR<T> & A,matrix<bigint> & TR,lidia_size_t startr,lidia_size_t startc,lidia_size_t index,lidia_size_t len,T * max_array,const T & BOUND) const2918 column_oriented_dense_matrix_modules< T >::normalize_row (MR< T > &A,
2919 matrix< bigint > &TR,
2920 lidia_size_t startr,
2921 lidia_size_t startc,
2922 lidia_size_t index,
2923 lidia_size_t len,
2924 T *max_array,
2925 const T &BOUND) const
2926 {
2927 T q, TMP;
2928 for (lidia_size_t i = 0; i <= startc; i++)
2929 if ((i != index) && member(A, startr, i) != A.Zero) {
2930 pos_div_rem(q, TMP, member(A, startr, i), member(A, startr, index));
2931
2932 // Check
2933 if (q != 0) {
2934 if (abs(bigint(bigint(max_array[index]) * q) + max_array[i]) > bigint(BOUND))
2935 return false;
2936 subtract_multiple_of_column(A, i, q, index, len);
2937 subtract_multiple_of_column(TR, i, q, index, TR.get_no_of_rows() - 1);
2938 update_max_array(A, i, max_array);
2939 }
2940 }
2941 return true;
2942 }
2943
2944
2945
2946 #undef row_oriented_dense_matrix_modules
2947 #undef column_oriented_dense_matrix_modules
2948
2949
2950
2951 #ifdef LIDIA_NAMESPACE
2952 # ifndef IN_NAMESPACE_LIDIA
2953 } // end of namespace LiDIA
2954 # endif
2955 #endif
2956
2957
2958
2959 #endif // LIDIA_DENSE_BIGINT_MATRIX_MODULES_CC_GUARD_
2960