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_FP_MATRIX_ALGORITHMS_CC_GUARD_
20 #define LIDIA_FP_MATRIX_ALGORITHMS_CC_GUARD_
21
22
23 #include "LiDIA/modular_operations.inl"
24 #ifndef LIDIA_FP_MATRIX_ALGORITHMS_H_GUARD_
25 # include "LiDIA/matrix/fp_matrix_algorithms.h"
26 #endif
27
28
29
30 #ifdef LIDIA_NAMESPACE
31 # ifndef IN_NAMESPACE_LIDIA
32 namespace LiDIA {
33 # endif
34 #endif
35
36
37
38 //
39 // column step form
40 //
41
42 template< class T, class MATRIX_TYPE >
43 int
STF(MATRIX_TYPE & A,const T & mod) const44 fp_matrix_algorithms< T, MATRIX_TYPE >::STF (MATRIX_TYPE &A, const T & mod) const
45 {
46 //
47 // Input: **value = values of matrix
48 // r = number of rows
49 // c = number of columns
50 // mod = modulus for Fp - class
51 // Task: ex = STF(value, r, c, mod);
52 // = > matrix (value, r, c) in column step form
53 // = > ex = (-1)^Number_of_column_exchanges
54 // Version: bigint 1.9
55 //
56
57 register lidia_size_t index = 0, i = 0, j = 0, z;
58
59 lidia_size_t startr = A.rows;
60 lidia_size_t startc = A.columns;
61
62 T TMP, TMP1, TMP2;
63 T *tmp;
64
65 // Step 1 - 4
66 register int exchange = 1;
67
68 // Step 5 - 8
69 for (--startc, --startr; startr >= 0 && startc >= 0; startc--, startr--) {
70 tmp = A.value[startr];
71
72 // Step 9 - 13
73 for (index = startc; index >= 0 && tmp[index] == A.Zero; index--);
74
75 // Step 14 - 26
76 if (index != -1) {
77 if (index != startc) {
78 exchange = -exchange;
79
80 // exchange column j0 with column j
81 for (z = 0; z <= startr; z++)
82 swap(A.value[z][startc], A.value[z][index]);
83 }
84
85 inv_mod(TMP, tmp[startc], mod);
86
87 // Step 27 - 29
88 for (j = startc - 1; j >= 0; j--) {
89 if (tmp[j] != A.Zero) {
90 // Step 30 - 40
91 mult_mod(TMP1, tmp[j], TMP, mod);
92
93 for (i = 0; i <= startr; i++) {
94 mult_mod(TMP2, A.value[i][startc], TMP1, mod);
95 sub_mod(A.value[i][j], A.value[i][j], TMP2, mod);
96 }
97 }
98 }
99 }
100 else {
101 startc++;
102
103 }
104 }
105 return exchange;
106 }
107
108
109
110 //
111 // extended column step form
112 //
113
114 template< class T, class MATRIX_TYPE >
115 const T *
STF_extended(MATRIX_TYPE & A,const T & mod) const116 fp_matrix_algorithms< T, MATRIX_TYPE >::STF_extended (MATRIX_TYPE &A, const T & mod) const
117 {
118 //
119 // Input: **value = values of matrix
120 // r = number of rows
121 // c = number of columns
122 // mod = modulus for Fp - class
123 // Task: ex = STF_extended(A, mod);
124 // = > matrix A in column step form
125 // = > ex = (-1)^Number_of_column_exchanges
126 // Version: bigint 1.9
127 //
128
129 register lidia_size_t index = 0, i = 0, j = 0, z;
130 lidia_size_t startr = A.rows;
131 lidia_size_t startc = A.columns;
132
133 lidia_size_t line = A.columns - A.rows;
134
135 T *RES = new T[line + 1];
136 memory_handler(RES, "STF_extended", "Error in memory allocation (RES)");
137
138 for (i = 0; i <= line; i++)
139 RES[i] = 1;
140
141 T TMP, TMP1, TMP2;
142 T *tmp;
143 bool SW = true;
144
145 // Step 1 - 4
146 register int exchange = 1;
147
148 // Step 5 - 8
149 for (--startc, --startr; startr >= 0 && startc >= 0; startc--, startr--) {
150 tmp = A.value[startr];
151
152 // Step 9 - 13
153 for (index = startc; index > line && tmp[index] == A.Zero; index--);
154
155 // Step 14 - 26
156 if (index > line) {
157 if (index != startc) {
158 exchange = -exchange;
159
160 // exchange column j0 with column j
161 for (z = 0; z <= startr; z++)
162 swap(A.value[z][startc], A.value[z][index]);
163 }
164
165 for (i = 0; i <= line; i++)
166 ::mult_mod(RES[i], RES[i], tmp[startc], mod);
167
168 inv_mod(TMP, tmp[startc], mod);
169
170 // Step 27 - 29
171 for (j = startc - 1; j >= 0; j--) {
172 if (tmp[j] != A.Zero) {
173 // Step 30 - 40
174 mult_mod(TMP1, tmp[j], TMP, mod);
175
176 for (i = 0; i < startr; i++) {
177 mult_mod(TMP2, A.value[i][startc], TMP1, mod);
178 sub_mod(A.value[i][j], A.value[i][j], TMP2, mod);
179 }
180 }
181 }
182 }
183 else {
184 if (SW) {
185 for (i = 0; i <= line; i++)
186 ::mult_mod(RES[i], RES[i], tmp[i], mod);
187 exchange = -exchange;
188 SW = false;
189 }
190 else {
191 for (i = 0; i <= line; i++)
192 RES[i] = 0;
193 }
194 startc++;
195 }
196 }
197
198 if (exchange == -1)
199 for (i = 0; i <= line; i++)
200 RES[i] = -RES[i];
201 return RES;
202 }
203
204
205
206 //
207 // rank
208 //
209
210 template< class T, class MATRIX_TYPE >
211 lidia_size_t
rank(MATRIX_TYPE & A,const T & mod) const212 fp_matrix_algorithms< T, MATRIX_TYPE >::rank (MATRIX_TYPE &A, const T &mod) const
213 {
214 //
215 // Input: **Avalue = values of matrix
216 // r = number of rows
217 // c = number of columns
218 // Gmod = modulus for Fp - class
219 // Task: rank_intern(Avalue, r, c, Gmod) = rank of matrix (Avalue, r, c)
220 // Version: 1.9
221 //
222
223 debug_handler_l("bigint_matrix",
224 "rank(MATRIX_TYPE &A, const T &)", LDBL_MATRIX);
225
226 register lidia_size_t i, j, No = 0;
227
228 // Step 1, 2
229 STF(A, mod);
230
231 // Step 3 - 12
232 for (j = 0; j < A.columns; j++) {
233 for (i = A.rows - 1; i >= 0 && A.value[i][j] == A.Zero; i--);
234 if (i == -1)
235 No++;
236 }
237
238 // Step 13 - 24
239 return (A.columns - No);
240 }
241
242
243
244 //
245 // rank and linearly independent rows or columns
246 //
247
248 template< class T, class MATRIX_TYPE >
249 lidia_size_t *
lininr(MATRIX_TYPE & A,const T & mod) const250 fp_matrix_algorithms< T, MATRIX_TYPE >::lininr (MATRIX_TYPE &A, const T &mod) const
251 {
252 //
253 // Input: **Avalue = values of matrix
254 // r = number of rows
255 // c = number of columns
256 // Gmod = modulus for Fp - class
257 // Task: IND = lininr_intern(Avalue, r, c, Gmod);
258 // = > IND[0] = Rank of matrix (Avalue, r, c)
259 // = > IND[1], ..., IND[IND[0]], such that
260 // row(IND[1]), ..., row(IND[IND[0]])
261 // of matrix (Avalue, r, c) are linearly independent.
262 // Version: 1.9
263 //
264
265 debug_handler("bigint_matrix",
266 "lininr(MATRIX_TYPE &A, const T &)");
267
268 register lidia_size_t i, j;
269 lidia_size_t *l = new lidia_size_t[A.columns];
270 for (i = 0; i < A.columns; l[i] = 0, i++);
271
272 STF(A, mod);
273
274 // Step 3 - 12
275 for (j = 0; j < A.columns; j++) {
276 for (i = A.rows - 1; i >= 0 && A.value[i][j] == A.Zero; i--);
277 l[j] = i;
278 }
279
280 // Step 13 - 24
281 for (i = 0; i < A.columns && l[i] == -1; i++); // i = number of zero-columns
282
283 lidia_size_t TMP = A.columns - i; // rank
284 lidia_size_t *IND = new lidia_size_t[TMP + 1];
285 memory_handler(IND, "bigint_matrix", "lininr_intern - "
286 "Error in memory allocation (IND)");
287 IND[0] = TMP; // rank
288 for (j = 0; j < TMP; j++)
289 IND[TMP - j] = l[j + i];
290 delete[] l;
291 return IND;
292 }
293
294
295
296 template< class T, class MATRIX_TYPE >
297 lidia_size_t *
lininc(MATRIX_TYPE & A,const T & mod) const298 fp_matrix_algorithms< T, MATRIX_TYPE >::lininc (MATRIX_TYPE &A, const T &mod) const
299 {
300 //
301 // Input: **Avalue = values of matrix
302 // r = number of rows
303 // c = number of columns
304 // Gmod = modulus for Fp - class
305 // Task: IND = lininc_intern(Avalue, r, c, Gmod);
306 // = > IND[0] = Rank of matrix (Avalue, r, c)
307 // = > IND[1], ..., IND[IND[0]], such that
308 // column(IND[1]), ..., column(IND[IND[0]])
309 // of matrix (Avalue, r, c) are linearly independent.
310 // Version: 1.9
311 //
312
313 register lidia_size_t i, j;
314 lidia_size_t *l = new lidia_size_t[A.columns + 1];
315 memory_handler(l, "bigint_matrix", "lininc_intern :: "
316 "Error in memory allocation (l)");
317
318 // Step 1, 2
319 STF(A, mod);
320
321 // Step 3 - 12
322 for (j = 0; j < A.columns; j++) {
323 for (i = A.rows - 1; i >= 0 && A.value[i][j] == A.Zero; i--);
324 l[j] = i;
325 }
326
327 // Step 13 - 24
328 for (i = 0; i < A.columns && l[i] == -1; i++); // i = number of zero-columns
329
330 lidia_size_t TMP = A.columns - i;
331 lidia_size_t *IND = new lidia_size_t[TMP + 1];
332 memory_handler(IND, "bigint_matrix", "lininc_intern - "
333 "Error in memory allocation (IND)");
334 IND[0] = TMP; // rank
335 for (j = 0; j < TMP; j++)
336 IND[TMP - j] = l[j + i];
337 delete[] l;
338 return IND;
339 }
340
341
342
343 //
344 // adjoint matrix
345 //
346
347 template< class T, class MATRIX_TYPE >
348 void
adj(MATRIX_TYPE & A,const T & mod) const349 fp_matrix_algorithms< T, MATRIX_TYPE >::adj (MATRIX_TYPE &A, const T &mod) const
350 {
351 //
352 // Input: **value = values of matrix
353 // r = number of rows = number of columns
354 // mod = modulus for Fp - class
355 // Task: adj_intern(value, r, mod);
356 // = > adjoint matrix (value, r, r)
357 // Version: 1.9
358 //
359
360 debug_handler("bigint_matrix", "in inline - function "
361 "adj(MATRIX_TYPE &, cont T &)");
362
363 register lidia_size_t i, j, z;
364 T TMP, TMP1, TMP2;
365 T *tmp, *tmp1, *Btmp, *Btmp1;
366 lidia_size_t exchange = 1;
367 T DET = 1;
368
369 // Step 1, 2
370 T **Bvalue = new T *[A.rows];
371 memory_handler(Bvalue, "bigint_matrix", "adj_intern - T part :: "
372 "Error in memory allocation (Bvalue)");
373 for (i = 0; i < A.rows; i++) {
374 Btmp = new T[A.rows];
375 tmp = A.value[i];
376 memory_handler(Btmp, "bigint_matrix", "adj_intern - T part :: "
377 "Error in memory allocation (Btmp)");
378 for (j = 0; j < A.rows; j++) {
379 Btmp[j] = tmp[j];
380 if (i == j)
381 tmp[j] = 1;
382 else
383 tmp[j] = 0;
384 }
385 Bvalue[i] = Btmp;
386 }
387
388 // Step 3 - 5
389 for (i = A.rows - 1; i >= 0; i--) {
390 Btmp = Bvalue[i];
391
392 // Step 6 - 11
393 for (j = i; j >= 0 && Btmp[j] == 0; j--);
394
395 // Step 12 - 19
396 if (j != i) {
397 exchange = -exchange;
398 for (z = 0; z < A.rows; z++) {
399 Btmp1 = Bvalue[z];
400 tmp1 = A.value[z];
401
402 // A.swap_columns(i, j);
403 TMP = Btmp1[j];
404 Btmp1[j] = Btmp1[i];
405 Btmp1[i] = TMP;
406
407 // B.swap_columns(i, j);
408 TMP = tmp1[i];
409 tmp1[i] = tmp1[j];
410 tmp1[j] = TMP;
411 }
412 }
413
414 inv_mod(TMP1, Btmp[i], mod);
415
416 // Step 20 - 32
417 for (j = 0; j < A.rows; j++) {
418 if (j != i) {
419 mult_mod(TMP2, Btmp[j], TMP1, mod);
420 for (z = 0; z < i; z++) {
421 Btmp1 = Bvalue[z];
422
423 mult_mod(TMP, Btmp1[i], TMP2, mod);
424 sub_mod(Btmp1[j], Btmp1[j], TMP, mod);
425 }
426
427 for (z = 0; z < A.rows; z++) {
428 tmp1 = A.value[z];
429
430 mult_mod(TMP, tmp1[i], TMP2, mod);
431 sub_mod(tmp1[j], tmp1[j], TMP, mod);
432 }
433 }
434 }
435 mult_mod(DET, DET, Btmp[i], mod);
436 for (z = 0; z < i; z++) {
437 Btmp1 = Bvalue[z];
438 mult_mod(Btmp1[i], Btmp1[i], TMP1, mod);
439 }
440 for (z = 0; z < A.rows; z++) {
441 tmp1 = A.value[z];
442 mult_mod(tmp1[i], tmp1[i], TMP1, mod);
443 }
444 }
445
446 // Step 33 - 43
447 if (exchange < 0)
448 DET = -DET;
449 for (j = 0; j < A.rows; j++) {
450 tmp = A.value[j];
451 for (i = 0; i < A.rows; i++)
452 mult_mod(tmp[i], tmp[i], DET, mod);
453 }
454
455 for (j = 0; j < A.rows; j++)
456 delete[] Bvalue[j];
457 delete[] Bvalue;
458 }
459
460
461
462 //
463 // determinant
464 //
465
466 template< class T, class MATRIX_TYPE >
467 const T
det(MATRIX_TYPE & A,const T & mod) const468 fp_matrix_algorithms< T, MATRIX_TYPE >::det (MATRIX_TYPE &A, const T &mod) const
469 {
470 //
471 // Input: **value = values of matrix
472 // r = number of rows = number of columns
473 // mod = modulus for Fp - class
474 // Task: det_intern(value, r, mod) = determinant of matrix (value, r, r);
475 // Version: T 1.9
476 //
477
478 debug_handler("bigint_matrix",
479 "det(MATRIX_TYPE &, T)");
480
481 register lidia_size_t i, j, z;
482 T TMP, TMP1, TMP2;
483 T *tmp, *tmp1;
484
485 // Step 1 - 4
486 lidia_size_t ex = 1;
487 T ret = 1;
488
489 // Step 5 - 8
490 for (i = 0; i < A.rows; i++) {
491
492 // Step 9 - 13
493 for (j = i; j < A.rows && A.value[j][i] == 0; j++);
494
495 // Step 14 - 26
496 if (j == A.rows)
497 return 0;
498 if (j != i) {
499 ex = -ex;
500 tmp1 = A.value[j];
501 A.value[j] = A.value[i];
502 A.value[i] = tmp1;
503 }
504 tmp = A.value[i];
505
506 // Step 27 - 29
507 inv_mod(TMP1, tmp[i], mod);
508
509 for (j = i + 1; j < A.rows; j++) {
510 tmp1 = A.value[j];
511 mult_mod(TMP2, tmp1[i], TMP1, mod);
512 for (z = i + 1; z < A.rows; z++) {
513 mult_mod(TMP, tmp[z], TMP2, mod);
514 sub_mod(tmp1[z], tmp1[z], TMP, mod);
515 }
516 }
517 mult_mod(ret, ret, tmp[i], mod);
518 }
519 if (ex < 0)
520 ret = -ret;
521 return ret;
522 }
523
524
525
526 //
527 // Hessenberg form
528 //
529
530 template< class T, class MATRIX_TYPE >
531 void
HBF(MATRIX_TYPE & A,const T & mod) const532 fp_matrix_algorithms< T, MATRIX_TYPE >::HBF (MATRIX_TYPE &A, const T &mod) const
533 {
534 //
535 // Input: //value = values of matrix
536 // r = number of rows = number of columns
537 // mod = modulus for Fp - class
538 // Task: HBF_intern(value, r, mod);
539 // = > matrix (value, r, r) in Hessenberg form
540 // Version: 1.9
541 //
542
543 debug_handler("bigint_matrix", "in inline - function "
544 "HBF_intern(T **, lidia_size_t, T)");
545
546 // Step 1, 2
547 lidia_size_t i, j, z;
548 T TMP, TMP1, TMP2;
549 T *tmp;
550
551 // Step 3 - 11
552 for (i = A.rows - 1; i >= 1; i--) {
553 for (j = i - 1; j >= 0 && A.value[i][j] == 0; j--);
554 if (j != -1) {
555
556 // Step 12, 13
557 if (j != i - 1) {
558
559 // Step 14 - 18
560 // exchange columns i-1 and j
561 for (z = 0; z < A.rows; z++) {
562 TMP = A.value[z][i - 1];
563 A.value[z][i - 1] = A.value[z][j];
564 A.value[z][j] = TMP;
565 }
566
567 // Step 19 - 24
568 // exchange rows i-1 and j
569 tmp = A.value[i - 1];
570 A.value[i - 1] = A.value[j];
571 A.value[j] = tmp;
572 }
573 tmp = A.value[i];
574
575 // Step 25 - 41
576 inv_mod(TMP2, tmp[i - 1], mod);
577 for (j = i - 2; j >= 0; j--) {
578 mult_mod(TMP1, tmp[j], TMP2, mod);
579 for (z = 0; z < A.rows; z++) {
580 mult_mod(TMP, A.value[z][i - 1], TMP1, mod);
581 sub_mod(A.value[z][j], A.value[z][j], TMP, mod);
582 }
583 for (z = 0; z < A.rows; z++) {
584 mult_mod(TMP, A.value[j][z], TMP1, mod);
585 add_mod(A.value[i - 1][z], A.value[i - 1][z], TMP, mod);
586 }
587 }
588 }
589 }
590 }
591
592
593
594 //
595 // characteristic polynomial
596 //
597
598 template< class T, class MATRIX_TYPE >
599 T *
600 fp_matrix_algorithms< T, MATRIX_TYPE >::charpoly (MATRIX_TYPE &A, const T &mod) const
601 {
602 //
603 // Input: **value = values of matrix
604 // r = number of rows = number of columns
605 // mod = modulus for Fp - class
606 // Task: RES = charpoly_intern(value, r, mod);
607 // = > RES[0], ..., RES[r] are the coefficients of
608 // the characteristic polynomial of matrix (value, r, r)
609 // Version: 1.9
610 //
611
612 debug_handler("bigint_matrix", "in inline - function "
613 "charpoly_intern(T **, lidia_size_t, T)");
614
615 lidia_size_t i, j, z;
616 T TMP;
617 T *tmp;
618 lidia_size_t sign;
619
620 // Step 1 - 5
621 HBF(A, mod);
622
623 T *K = new T[A.rows]; // size = r
624 memory_handler(K, "bigint_matrix", "charpoly_intern - Version T ::"
625 "Error in memory allocation (K)");
626 for (i = 0; i < A.rows; i++)
627 K[i] = 1;
628
629 // Step 6 - 8
630 T **P = new T *[A.rows + 1];
631 memory_handler(P, "bigint_matrix", "charpoly_intern - Version T :: "
632 "Error in memory allocation (P)");
633 for (i = 0; i < A.rows + 1; i++) {
634 tmp = new T[A.rows + 1];
635 memory_handler(tmp, "bigint_matrix", "charpoly_intern - Version T :: "
636 "Error in memory allocation (tmp)");
637 for (j = 0; j < A.rows + 1; j++)
638 tmp[j] = 0;
639 P[i] = tmp;
640 }
641 P[0][0] = 1;
642
643 // Step 9 - 11
644 for (z = 1; z <= A.rows; z++) {
645
646 // Step 12 - 16
647 for (j = 1; j < z; j++)
648 mult_mod(K[j - 1], K[j - 1], A.value[z - 1][z - 2], mod);
649
650 // Step 17 - 23
651 P[z][z] = mod - P[z - 1][z - 1];
652 for (i = 1; i < z; i++) {
653 mult_mod(TMP, A.value[z - 1][z - 1], P[i][z - 1], mod);
654 sub_mod(P[i][z], TMP, P[i - 1][z - 1], mod);
655 }
656 mult_mod(P[0][z], A.value[z - 1][z - 1], P[0][z - 1], mod);
657
658 // Step 24 - 34
659 sign = 1;
660 for (j = z - 1; j >= 1; j--) {
661 sign = -sign;
662 for (i = 0; i <= j - 1; i++) {
663 mult_mod(TMP, sign, P[i][j - 1], mod);
664 mult_mod(TMP, TMP, A.value[j - 1][z - 1], mod);
665 mult_mod(TMP, TMP, K[j - 1], mod);
666 add_mod(P[i][z], P[i][z], TMP, mod);
667 }
668 }
669 }
670
671 // Step 35 - 40
672 T *RES = new T[A.rows + 1];
673 memory_handler(RES, "bigint_matrix", "charpoly_intern - Version T :: "
674 "Error in memory allocation (RES)");
675 for (i = 0; i <= A.rows; i++) {
676 RES[i] = P[i][A.rows];
677 delete[] P[i];
678 }
679 delete[] P;
680 delete[] K;
681 return RES;
682 }
683
684
685
686 #ifdef LIDIA_NAMESPACE
687 # ifndef IN_NAMESPACE_LIDIA
688 } // end of namespace LiDIA
689 # endif
690 #endif
691
692
693
694 #endif // LIDIA_FP_MATRIX_ALGORITHMS_CC_GUARD_
695