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