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_NORMALIZE_KERNEL_CC_GUARD_
20 #define LIDIA_NORMALIZE_KERNEL_CC_GUARD_
21
22
23 #include <fstream>
24 #ifndef LIDIA_INFO_H_GUARD_
25 # include "LiDIA/info.h"
26 #endif
27 #ifndef LIDIA_NORMALIZE_KERNEL_H_GUARD_
28 # include "LiDIA/matrix/normalize_kernel.h"
29 #endif
30 #include <cstdlib>
31
32
33 #ifdef LIDIA_NAMESPACE
34 # ifndef IN_NAMESPACE_LIDIA
35 namespace LiDIA {
36 # endif
37 #endif
38
39
40
41 #ifdef LIDIA_NAMESPACE
42 using std::abs;
43 #endif
44
45
46
47 #define DMESSAGE "normalize_kernel"
48
49
50 ////////////////////////////////////////////////////////////////////////////////////////
51 // normalization
52 ////////////////////////////////////////////////////////////////////////////////////////
53
54 //
55 // normalize Standard
56 //
57
58 template< class T, class REP, class REP1 >
59 lidia_size_t *
normalize_Std(MR<T> & A,lidia_size_t startr,lidia_size_t startc) const60 normalization_kernel< T, REP, REP1 >::normalize_Std (MR< T > &A,
61 lidia_size_t startr,
62 lidia_size_t startc) const
63 {
64 lidia_size_t i, j, k = 1;
65 T q, r;
66
67 lidia_size_t *RET = new lidia_size_t[A.rows];
68
69 #ifdef LIDIA_COMPUTE_RUNTIME_INFO
70 std::ofstream profile("normalize_max.profile");
71 T MAX = 0, GMAX = 0;
72 #endif
73
74 for (i = A.columns - 1, j = A.rows - 1; i >= startc && j >= startr; i--, j--) {
75 #ifdef LIDIA_COMPUTE_RUNTIME_INFO
76 modul.max(A, MAX);
77 if (GMAX < MAX)
78 GMAX = MAX;
79 #endif
80 if (modul.member(A, j, i) != A.Zero) {
81 for (lidia_size_t l = i + 1; l < A.columns; l++)
82 if (modul.member(A, j, l) != A.Zero) {
83 pos_div_rem(q, r, modul.member(A, j, l), modul.member(A, j, i));
84 modul.subtract_multiple_of_column(A, l, q, i, A.rows - 1);
85 }
86 }
87 else {
88 RET[k] = j;
89 k++;
90 }
91 }
92 if (k == 1) {
93 delete [] RET;
94 RET = NULL;
95 }
96 else
97 RET[0] = k - 1;
98
99 #ifdef LIDIA_COMPUTE_RUNTIME_INFO
100 profile << GMAX << std::endl;
101 profile.close();
102 #endif
103
104 return RET;
105 }
106
107
108
109 template< class T, class REP, class REP1 >
110 lidia_size_t *
normalize_Std(MR<T> & A,matrix<bigint> & TR,lidia_size_t startr,lidia_size_t startc) const111 normalization_kernel< T, REP, REP1 >::normalize_Std (MR< T > &A,
112 matrix< bigint > &TR,
113 lidia_size_t startr,
114 lidia_size_t startc) const
115 {
116 lidia_size_t i, j, k = 1;
117 T q, r;
118
119 lidia_size_t *RET = new lidia_size_t[A.rows];
120
121 for (i = A.columns - 1, j = A.rows - 1; i >= startc && j >= startr; i--, j--)
122 if (modul.member(A, j, i) != A.Zero) {
123 for (lidia_size_t l = i + 1; l < A.columns; l++)
124 if (modul.member(A, j, l) != A.Zero) {
125 pos_div_rem(q, r, modul.member(A, j, l), modul.member(A, j, i));
126 modul.subtract_multiple_of_column(A, l, q, i, A.rows - 1);
127
128 for (lidia_size_t m = 0; m < A.columns; m++)
129 TR.sto(m, l, TR.member(m, l) - q * TR.member(m, i));
130 }
131 }
132 else {
133 RET[k] = j;
134 k++;
135 }
136
137 if (k == 1) {
138 delete [] RET;
139 RET = NULL;
140 }
141 else
142 RET[0] = k - 1;
143 return RET;
144 }
145
146
147
148 //
149 // normalize_ChouCollins
150 //
151
152 template< class T, class REP, class REP1 >
153 lidia_size_t *
normalize_ChouCollins(MR<T> & A,lidia_size_t startr,lidia_size_t startc) const154 normalization_kernel< T, REP, REP1 >::normalize_ChouCollins (MR< T > &A,
155 lidia_size_t startr,
156 lidia_size_t startc) const
157 {
158 lidia_size_t start = startc, i, k, p, l;
159 T q, r;
160
161 lidia_size_t *RET = new lidia_size_t[A.rows];
162
163 #ifdef LIDIA_COMPUTE_RUNTIME_INFO
164 std::ofstream profile("normalize_max.profile");
165 T MAX = 0, GMAX = 0;
166 #endif
167
168 for (i = start + 1, p = startr; i < A.columns; i++, p++) {
169 // std::cout << "i = " << i << " p = " << p << std::endl;
170 for (k = i - 1, l = p; k >= start; k--, l--) {
171 // std::cout << "\t k = " << k << " l = " << l << std::endl;
172 #ifdef LIDIA_COMPUTE_RUNTIME_INFO
173 modul.max(A, MAX);
174 if (GMAX < MAX)
175 GMAX = MAX;
176 #endif
177
178 if (modul.member(A, l, i) != A.Zero && modul.member(A, l, k) != A.Zero) {
179 pos_div_rem(q, r, modul.member(A, l, i), modul.member(A, l, k));
180 modul.subtract_multiple_of_column(A, i, q, k, A.rows - 1);
181 }
182 }
183 }
184
185 k = 1;
186 for (i = startr; i < A.rows; i++) {
187 if (modul.member(A, i, start + i - startr) == A.Zero) {
188 RET[k] = i;
189 k++;
190 }
191 }
192
193 RET[0] = k - 1;
194
195 if (RET[0] == 0) {
196 delete[] RET;
197 RET = NULL;
198 }
199
200 #ifdef LIDIA_COMPUTE_RUNTIME_INFO
201 profile << GMAX << std::endl;
202 profile.close();
203 #endif
204
205 return RET;
206 }
207
208
209
210 template< class T, class REP, class REP1 >
211 lidia_size_t *
normalize_ChouCollins(MR<T> & A,matrix<bigint> & TR,lidia_size_t startr,lidia_size_t startc) const212 normalization_kernel< T, REP, REP1 >::normalize_ChouCollins (MR< T > &A,
213 matrix< bigint > &TR,
214 lidia_size_t startr,
215 lidia_size_t startc) const
216 {
217 lidia_size_t start = startc, i, k, p, l;
218 T q, r;
219
220 lidia_size_t *RET = new lidia_size_t[A.rows];
221
222 for (i = start + 1, p = startr; i < A.columns; i++, p++) {
223 for (k = i - 1, l = p; k >= start; k--, l--) {
224 if (modul.member(A, l, i) != A.Zero && modul.member(A, l, k) != A.Zero) {
225 pos_div_rem(q, r, modul.member(A, l, i), modul.member(A, l, k));
226 modul.subtract_multiple_of_column(A, i, q, k, A.rows - 1);
227
228 for (lidia_size_t m = 0; m < A.columns; m++)
229 TR.sto(m, i, TR.member(m, i) - q * TR.member(m, k));
230 }
231 }
232 }
233
234 k = 1;
235 for (i = startr; i < A.rows; i++)
236 if (modul.member(A, i, start + i - startr) == A.Zero) {
237 RET[k] = i;
238 k++;
239 }
240 RET[0] = k - 1;
241
242 if (RET[0] == 0) {
243 delete[] RET;
244 RET = NULL;
245 }
246
247 return RET;
248 }
249
250
251
252 template< class T, class REP, class REP1 >
253 lidia_size_t *
normalize_ChouCollins(MR<T> & A,trans_matrix & TR,matrix<bigint> & tran,lidia_size_t startr,lidia_size_t startc) const254 normalization_kernel< T, REP, REP1 >::normalize_ChouCollins (MR< T > &A,
255 trans_matrix &TR,
256 matrix< bigint > & tran,
257 lidia_size_t startr,
258 lidia_size_t startc) const
259 {
260 lidia_size_t start = startc, i, k, p, l;
261 T q, r;
262 matrix< bigint > otran;
263
264 tran.reset();
265 otran = tran;
266 tran.resize(A.columns, A.columns);
267 tran.diag(1, 0);
268
269 lidia_size_t *RET = new lidia_size_t[A.rows];
270
271 for (i = start + 1, p = startr; i < A.columns; i++, p++) {
272 for (k = i - 1, l = p; k >= start; k--, l--) {
273 if (modul.member(A, l, i) != A.Zero && modul.member(A, l, k) != A.Zero) {
274 pos_div_rem(q, r, modul.member(A, l, i), modul.member(A, l, k));
275 modul.subtract_multiple_of_column(A, i, q, k, A.rows - 1);
276
277 for (lidia_size_t m = 0; m < A.columns; m++)
278 tran.sto(m, i, tran.member(m, i) - q * tran.member(m, k));
279 }
280 }
281
282 if (i == (A.rows/2)) {
283 TR.store_matrix(tran);
284
285 tran = otran;
286 tran.resize(A.columns, A.columns);
287 tran.diag(1, 0);
288 }
289 }
290
291 k = 1;
292 for (i = startr; i < A.rows; i++)
293 if (modul.member(A, i, start + i - startr) == A.Zero) {
294 RET[k] = i;
295 k++;
296 }
297 RET[0] = k - 1;
298
299 if (RET[0] == 0) {
300 delete[] RET;
301 RET = NULL;
302 }
303
304 return RET;
305 }
306
307
308
309 //
310 // normalize_ChouCollins_extended
311 //
312
313 template< class T, class REP, class REP1 >
314 lidia_size_t *
normalize_ChouCollins_extended(MR<T> & A,lidia_size_t startr,lidia_size_t startc) const315 normalization_kernel< T, REP, REP1 >::normalize_ChouCollins_extended (MR< T > &A,
316 lidia_size_t startr,
317 lidia_size_t startc) const
318 {
319 lidia_size_t start = startc, i, k, p, l;
320 T q, r;
321
322 lidia_size_t *RET = new lidia_size_t[A.rows];
323
324 for (i = start + 1, p = startr; i < A.columns; i++, p++)
325 for (k = i - 1, l = p; k >= start; k--, l--)
326 if (modul.member(A, l, i) != A.Zero && modul.member(A, l, k) != A.Zero) {
327 pos_div_rem(q, r, modul.member(A, l, i), modul.member(A, l, k));
328 modul.subtract_multiple_of_column(A, i, q, k, A.rows - 1);
329 }
330
331 // correct diagonal
332 for (i = A.columns - 1, l = A.rows - 1; i > start; i--, l--)
333 if (modul.member(A, l, i) == A.Zero) {
334 k = l;
335 while (modul.member(A, k, i) == A.Zero && k > startr)
336 k--;
337 if (k > startr) {
338 if (modul.member(A, k, start + k - startr) == A.Zero) {
339 modul.swap_columns(A, i, start + k - startr);
340 i++;
341 l++;
342 }
343 }
344 }
345
346 k = 1;
347 for (i = startr; i < A.rows; i++)
348 if (modul.member(A, i, start + i - startr) == A.Zero) {
349 RET[k] = i;
350 k++;
351 }
352 RET[0] = k - 1;
353 return RET;
354 }
355
356
357
358 //
359 // normalizeHybrid_Std
360 //
361
362 template< class T, class REP, class REP1 >
363 lidia_size_t *
normalizeHybrid_Std(MR<T> & A,lidia_size_t startr,lidia_size_t startc) const364 normalization_kernel< T, REP, REP1 >::normalizeHybrid_Std (MR< T > &A,
365 lidia_size_t startr,
366 lidia_size_t startc) const
367 {
368 lidia_size_t start = startc, i, k;
369 T q, r;
370
371 #ifdef LIDIA_COMPUTE_RUNTIME_INFO
372 std::ofstream profile("normalize_max.profile");
373 T MAX = 0, GMAX = 0;
374 #endif
375
376 // Phase 1
377 for (i = startr; i < A.rows; i++) {
378
379 #ifdef LIDIA_COMPUTE_RUNTIME_INFO
380 modul.max(A, MAX);
381 if (GMAX < MAX)
382 GMAX = MAX;
383 #endif
384
385 if (modul.member(A, i, start) == 1)
386 for (lidia_size_t l = start + 1; l < A.columns; l++)
387 if (modul.member(A, i, l) != A.Zero)
388 modul.subtract_multiple_of_column(A, l, modul.member(A, i, l), start, A.rows - 1);
389 ++start;
390 }
391
392 // Phase 2
393 lidia_size_t *RET = new lidia_size_t[A.rows];
394 k = 1;
395
396 start--;
397 for (--i; i >= startr; i--) {
398
399 #ifdef LIDIA_COMPUTE_RUNTIME_INFO
400 modul.max(A, MAX);
401 if (GMAX < MAX)
402 GMAX = MAX;
403 #endif
404
405 if (modul.member(A, i, start) == A.Zero)
406 RET[k++] = i;
407 else {
408 if (modul.member(A, i, start) != 1)
409 for (lidia_size_t l = start + 1; l < A.columns; ++l) {
410 if (modul.member(A, i, l) != A.Zero) {
411 pos_div_rem(q, r, modul.member(A, i, l), modul.member(A, i, start));
412 modul.subtract_multiple_of_column(A, l, q, start, A.rows - 1);
413 }
414 }
415 }
416 start--;
417 }
418 if (k == 1) {
419 delete [] RET;
420 RET = NULL;
421 }
422 else
423 RET[0] = k - 1;
424
425 #ifdef LIDIA_COMPUTE_RUNTIME_INFO
426 profile << GMAX << std::endl;
427 profile.close();
428 #endif
429
430 return RET;
431 }
432
433
434
435 template< class T, class REP, class REP1 >
436 lidia_size_t *
normalizeHybrid_Std(MR<T> & A,matrix<bigint> & TR,lidia_size_t startr,lidia_size_t startc) const437 normalization_kernel< T, REP, REP1 >::normalizeHybrid_Std (MR< T > &A,
438 matrix< bigint > &TR,
439 lidia_size_t startr,
440 lidia_size_t startc) const
441 {
442 lidia_size_t start = startc, i, k;
443 T q, r;
444
445 // Phase 1
446 for (i = startr; i < A.rows; ++i) {
447 if (modul.member(A, i, start) == 1)
448 for (lidia_size_t l = start + 1; l < A.columns; ++l)
449 if (modul.member(A, i, l) != A.Zero) {
450 q = modul.member(A, i, l);
451 modul.subtract_multiple_of_column(A, l, q, start, A.rows - 1);
452
453 for (lidia_size_t m = 0; m < A.columns; m++)
454 TR.sto(m, l, TR.member(m, l) - q * TR.member(m, start));
455 }
456 ++start;
457 }
458
459 lidia_size_t *RET = new lidia_size_t[A.rows];
460 k = 1;
461
462 start--;
463
464 // Phase 2
465 for (--i; i >= startr; --i) {
466 if (modul.member(A, i, start) == A.Zero)
467 RET[k++] = i;
468 else {
469 if (modul.member(A, i, start) != 1)
470 for (lidia_size_t l = start + 1; l < A.columns; ++l) {
471 if (modul.member(A, i, l) != A.Zero) {
472 pos_div_rem(q, r, modul.member(A, i, l), modul.member(A, i, start));
473 modul.subtract_multiple_of_column(A, l, q, start, A.rows - 1);
474
475 for (lidia_size_t m = 0; m < A.columns; m++)
476 TR.sto(m, l, TR.member(m, l) - q * TR.member(m, start));
477 }
478 }
479 }
480 --start;
481 }
482
483 if (k == 1) {
484 delete [] RET;
485 RET = NULL;
486 }
487 else
488 RET[0] = k-1;
489
490 return RET;
491 }
492
493
494
495 template< class T, class REP, class REP1 >
496 lidia_size_t *
normalizeHybrid_Std(MR<T> & A,trans_matrix & TR,matrix<bigint> & tran,lidia_size_t s)497 normalization_kernel< T, REP, REP1 >::normalizeHybrid_Std (MR< T > &A,
498 trans_matrix &TR,
499 matrix< bigint > & tran,
500 lidia_size_t s)
501 {
502 lidia_size_t start = s, i, k, p, l;
503 T q, r;
504 matrix< bigint > otran;
505
506 tran.reset();
507 otran = tran;
508 tran.resize(A.columns, A.columns);
509 tran.diag(1, 0);
510
511 lidia_size_t *RET = new lidia_size_t[A.rows];
512
513 for (i = start + 1, p = 0; i < A.columns; i++, p++) {
514 for (k = i - 1, l = p; k >= start; k--, l--) {
515 if (modul.member(A, l, i) != A.Zero && modul.member(A, l, k) != A.Zero) {
516 pos_div_rem(q, r, modul.member(A, l, i), modul.member(A, l, k));
517 modul.subtract_multiple_of_column(A, i, q, k, A.rows - 1);
518
519 for (lidia_size_t m = 0; m < A.columns; m++)
520 tran.sto(m, i, tran.member(m, i) - q * tran.member(m, k));
521 }
522 }
523
524 if (i == (A.rows/2)) {
525 TR.store_matrix(tran);
526
527 tran = otran;
528 tran.resize(A.columns, A.columns);
529 tran.diag(1, 0);
530 }
531 }
532
533 k = 1;
534 for (i = 0; i < A.rows; i++)
535 if (modul.member(A, i, start + i) == A.Zero) {
536 RET[k] = i;
537 k++;
538 }
539 RET[0] = k - 1;
540
541 if (RET[0] == 0) {
542 delete[] RET;
543 RET = NULL;
544 }
545
546 return RET;
547 }
548
549
550
551 //
552 // normalizeHybrid_ChouColllins
553 //
554
555 template< class T, class REP, class REP1 >
556 lidia_size_t *
normalizeHybrid_ChouCollins(MR<T> & A,lidia_size_t startr,lidia_size_t startc) const557 normalization_kernel< T, REP, REP1 >::normalizeHybrid_ChouCollins (MR< T > &A,
558 lidia_size_t startr,
559 lidia_size_t startc) const
560 {
561 lidia_size_t start = startc, i, k;
562 T q, r;
563
564 #ifdef LIDIA_COMPUTE_RUNTIME_INFO
565 std::ofstream profile("normalize_max.profile");
566 T MAX = 0, GMAX = 0;
567 #endif
568
569 // Phase 1
570 for (i = startr; i < A.rows; i++) {
571
572 #ifdef LIDIA_COMPUTE_RUNTIME_INFO
573 modul.max(A, MAX);
574 if (GMAX < MAX)
575 GMAX = MAX;
576 #endif
577
578 if (modul.member(A, i, start) == 1)
579 for (lidia_size_t l = start + 1; l < A.columns; l++)
580 if (modul.member(A, i, l) != A.Zero)
581 modul.subtract_multiple_of_column(A, l, modul.member(A, i, l), start, A.rows - 1);
582 ++start;
583 }
584
585 lidia_size_t *RET = new lidia_size_t[A.rows];
586 k = 1;
587
588 start = startc;
589
590 // Phase 2
591 lidia_size_t p, l;
592
593 for (i = start + 1, p = startr; i < A.columns; i++, p++) {
594 std::cout << "i = " << i << " p = " << p << std::endl;
595 for (k = i - 1, l = p; k >= start; k--, l--) {
596 std::cout << "\t k = " << k << " l = " << l << std::endl;
597 #ifdef LIDIA_COMPUTE_RUNTIME_INFO
598 modul.max(A, MAX);
599 if (GMAX < MAX)
600 GMAX = MAX;
601 #endif
602 if (modul.member(A, l, k) != 1)
603 if (modul.member(A, l, i) != A.Zero && modul.member(A, l, k) != A.Zero) {
604 pos_div_rem(q, r, modul.member(A, l, i), modul.member(A, l, k));
605 modul.subtract_multiple_of_column(A, i, q, k, A.rows - 1);
606 }
607 }
608 }
609
610 k = 1;
611 for (i = startr; i < A.rows; i++) {
612 if (modul.member(A, i, start + i - startr) == A.Zero) {
613 RET[k] = i;
614 k++;
615 }
616 }
617
618 RET[0] = k - 1;
619
620 if (RET[0] == 0) {
621 delete[] RET;
622 RET = NULL;
623 }
624
625 #ifdef LIDIA_COMPUTE_RUNTIME_INFO
626 profile << GMAX << std::endl;
627 profile.close();
628 #endif
629
630 return RET;
631 }
632
633
634
635 template< class T, class REP, class REP1 >
636 lidia_size_t *
normalizeHybrid_ChouCollins(MR<T> & A,matrix<bigint> & TR,lidia_size_t startr,lidia_size_t startc) const637 normalization_kernel< T, REP, REP1 >::normalizeHybrid_ChouCollins (MR< T > &A,
638 matrix< bigint > &TR,
639 lidia_size_t startr,
640 lidia_size_t startc) const
641 {
642 lidia_size_t start = startc, i, k;
643 T q, r;
644
645 // Phase 1
646 for (i = startr; i < A.rows; ++i) {
647 if (modul.member(A, i, start) == 1)
648 for (lidia_size_t l = start + 1; l < A.columns; ++l)
649 if (modul.member(A, i, l) != A.Zero) {
650 q = modul.member(A, i, l);
651 modul.subtract_multiple_of_column(A, l, q, start, A.rows - 1);
652
653 for (lidia_size_t m = 0; m < A.columns; m++)
654 TR.sto(m, l, TR.member(m, l) - q * TR.member(m, start));
655 }
656 ++start;
657 }
658
659 lidia_size_t *RET = new lidia_size_t[A.rows];
660 k = 1;
661
662 start = startc;
663
664 // Phase 2
665 lidia_size_t p, l;
666
667 for (i = start + 1, p = startr; i < A.columns; i++, p++) {
668 for (k = i - 1, l = p; k >= start; k--, l--) {
669 if (modul.member(A, l, k) != 1)
670 if (modul.member(A, l, i) != A.Zero && modul.member(A, l, k) != A.Zero) {
671 pos_div_rem(q, r, modul.member(A, l, i), modul.member(A, l, k));
672 modul.subtract_multiple_of_column(A, i, q, k, A.rows - 1);
673
674 for (lidia_size_t m = 0; m < A.columns; m++)
675 TR.sto(m, i, TR.member(m, i) - q * TR.member(m, k));
676 }
677 }
678 }
679
680 k = 1;
681 for (i = startr; i < A.rows; i++)
682 if (modul.member(A, i, start + i - startr) == A.Zero) {
683 RET[k] = i;
684 k++;
685 }
686 RET[0] = k - 1;
687
688 if (RET[0] == 0) {
689 delete[] RET;
690 RET = NULL;
691 }
692
693 return RET;
694 }
695
696
697
698 template< class T, class REP, class REP1 >
699 lidia_size_t *
normalizeHybrid_ChouCollins(MR<T> & A,trans_matrix & TR,matrix<bigint> & tran,lidia_size_t s)700 normalization_kernel< T, REP, REP1 >::normalizeHybrid_ChouCollins (MR< T > &A,
701 trans_matrix &TR,
702 matrix< bigint > & tran,
703 lidia_size_t s)
704 {
705 lidia_size_t start = s, i, k, p, l;
706 T q, r;
707 matrix< bigint > otran;
708
709 tran.reset();
710 otran = tran;
711 tran.resize(A.columns, A.columns);
712 tran.diag(1, 0);
713
714 lidia_size_t *RET = new lidia_size_t[A.rows];
715
716 for (i = start + 1, p = 0; i < A.columns; i++, p++) {
717 for (k = i - 1, l = p; k >= start; k--, l--) {
718 if (modul.member(A, l, i) != A.Zero && modul.member(A, l, k) != A.Zero) {
719 pos_div_rem(q, r, modul.member(A, l, i), modul.member(A, l, k));
720 modul.subtract_multiple_of_column(A, i, q, k, A.rows - 1);
721
722 for (lidia_size_t m = 0; m < A.columns; m++)
723 tran.sto(m, i, tran.member(m, i) - q * tran.member(m, k));
724 }
725 }
726
727 if (i == (A.rows/2)) {
728 TR.store_matrix(tran);
729
730 tran = otran;
731 tran.resize(A.columns, A.columns);
732 tran.diag(1, 0);
733 }
734 }
735
736 k = 1;
737 for (i = 0; i < A.rows; i++)
738 if (modul.member(A, i, start + i) == A.Zero) {
739 RET[k] = i;
740 k++;
741 }
742 RET[0] = k - 1;
743
744 if (RET[0] == 0) {
745 delete[] RET;
746 RET = NULL;
747 }
748
749 return RET;
750 }
751
752
753
754 //
755 // normalize modular Standard
756 //
757
758 template< class T, class REP, class REP1 >
759 lidia_size_t *
normalizeMod_Std(MR<T> & A,lidia_size_t startr,lidia_size_t startc) const760 normalization_kernel< T, REP, REP1 >::normalizeMod_Std (MR< T > &A,
761 lidia_size_t startr,
762 lidia_size_t startc) const
763 {
764 lidia_size_t i, j, k = 1, l;
765 T q, r;
766
767 #ifdef LIDIA_COMPUTE_RUNTIME_INFO
768 std::ofstream profile("normalize_max.profile");
769 T MAX = 0, GMAX = 0;
770 #endif
771
772 T mod = 2;
773 for (i = startr, j = startc; i < A.rows && j < A.columns; i++, j++)
774 LiDIA::multiply(mod, mod, abs(modul.member(A, i, j)));
775
776 lidia_size_t *RET = new lidia_size_t[A.rows];
777
778 for (i = startr, j = startc; i < A.rows && j < A.columns; i++, j++)
779 for (l = j + 1; l < A.columns; l++) {
780 LiDIA::best_remainder(r, modul.member(A, i, l), mod);
781 modul.sto(A, i, l, r);
782 }
783
784 for (i = A.columns - 1, j = A.rows - 1; i >= startc && j >= startr; i--, j--) {
785
786 #ifdef LIDIA_COMPUTE_RUNTIME_INFO
787 modul.max(A, MAX);
788 if (GMAX < MAX)
789 GMAX = MAX;
790 #endif
791
792 if (modul.member(A, j, i) != A.Zero) {
793 for (lidia_size_t l = i + 1; l < A.columns; l++)
794 if (modul.member(A, j, l) != A.Zero) {
795 pos_div_rem(q, r, modul.member(A, j, l), modul.member(A, j, i));
796 modul.normalize_column_mod(A, l, q, i, A.rows - 1, mod);
797 }
798 }
799 else {
800 RET[k] = j;
801 k++;
802 }
803 }
804
805 if (k == 1) {
806 delete [] RET;
807 RET = NULL;
808 }
809 else
810 RET[0] = k - 1;
811
812 #ifdef LIDIA_COMPUTE_RUNTIME_INFO
813 profile << GMAX << std::endl;
814 profile.close();
815 #endif
816
817 return RET;
818 }
819
820
821
822 //
823 // normalize modular ChouCollins
824 //
825
826 template< class T, class REP, class REP1 >
827 lidia_size_t *
normalizeMod_ChouCollins(MR<T> & A,lidia_size_t startr,lidia_size_t startc) const828 normalization_kernel< T, REP, REP1 >::normalizeMod_ChouCollins (MR< T > &A,
829 lidia_size_t startr,
830 lidia_size_t startc) const
831 {
832 lidia_size_t start = startc, i, j, k, p, l;
833 T q, r;
834
835 #ifdef LIDIA_COMPUTE_RUNTIME_INFO
836 std::ofstream profile("normalize_max.profile");
837 T MAX = 0, GMAX = 0;
838 #endif
839
840 T mod = 2;
841 for (i = startr, p = startc; p < A.columns && i < A.rows; i++, p++)
842 LiDIA::multiply(mod, mod, abs(modul.member(A, i, p)));
843
844 lidia_size_t *RET = new lidia_size_t[A.rows];
845
846 for (i = startr, j = startc; i < A.rows && j < A.columns; i++, j++)
847 for (l = j + 1; l < A.columns; l++) {
848 LiDIA::best_remainder(r, modul.member(A, i, l), mod);
849 modul.sto(A, i, l, r);
850 }
851
852 for (i = start + 1, p = startr; i < A.columns; i++, p++) {
853 for (k = i - 1, l = p; k >= start; k--, l--) {
854
855 #ifdef LIDIA_COMPUTE_RUNTIME_INFO
856 modul.max(A, MAX);
857 if (GMAX < MAX)
858 GMAX = MAX;
859 #endif
860
861 if (modul.member(A, l, i) != A.Zero && modul.member(A, l, k) != A.Zero) {
862 pos_div_rem(q, r, modul.member(A, l, i), modul.member(A, l, k));
863 modul.normalize_column_mod(A, i, q, k, A.rows - 1, mod);
864 }
865 }
866 }
867
868 k = 1;
869 for (i = startr; i < A.rows; i++) {
870 if (modul.member(A, i, start + i - startr) == A.Zero) {
871 RET[k] = i;
872 k++;
873 }
874 }
875
876 RET[0] = k - 1;
877
878 if (RET[0] == 0) {
879 delete[] RET;
880 RET = NULL;
881 }
882
883 #ifdef LIDIA_COMPUTE_RUNTIME_INFO
884 profile << GMAX << std::endl;
885 profile.close();
886 #endif
887
888 return RET;
889 }
890
891
892
893 #undef DMESSAGE
894
895
896
897 #ifdef LIDIA_NAMESPACE
898 # ifndef IN_NAMESPACE_LIDIA
899 } // end of namespace LiDIA
900 # endif
901 #endif
902
903
904
905 #endif // LIDIA_NORMALIZE_KERNEL_CC_GUARD_
906