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_ALGORITHMS_CC_GUARD_
20 #define LIDIA_SPARSE_FP_MATRIX_ALGORITHMS_CC_GUARD_
21
22
23 #ifndef LIDIA_RANDOM_GENERATOR_H_GUARD_
24 # include "LiDIA/random_generator.h"
25 #endif
26 #include "LiDIA/modular_operations.inl"
27 #ifndef LIDIA_SPARSE_FP_MATRIX_ALGORITHMS_H_GUARD_
28 # include "LiDIA/matrix/sparse_fp_matrix_algorithms.h"
29 #endif
30
31 #ifndef LIDIA_BIGINT_MATRIX_H_GUARD_
32 # include "LiDIA/bigint_matrix.h"
33 #endif
34
35
36
37 #ifdef LIDIA_NAMESPACE
38 # ifndef IN_NAMESPACE_LIDIA
39 namespace LiDIA {
40 # endif
41 #endif
42
43
44
45 //
46 // multiply
47 //
48
49 math_vector< bigint >
operator *(const base_power_product<ring_matrix<bigint>,lidia_size_t> & A,const math_vector<bigint> & v)50 operator * (const base_power_product< ring_matrix < bigint >, lidia_size_t > &A,
51 const math_vector< bigint > &v)
52 {
53 math_vector< bigint > w = v;
54 lidia_size_t len = A.get_no_of_components();
55 for (lidia_size_t i = 0; i < len; i++)
56 w = A.get_base(i) * w;
57 return w;
58 }
59
60
61
62 math_vector< long >
operator *(const base_power_product<ring_matrix<long>,lidia_size_t> & A,const math_vector<long> & v)63 operator * (const base_power_product< ring_matrix < long >, lidia_size_t > &A,
64 const math_vector< long > &v)
65 {
66 math_vector< long > w = v;
67 lidia_size_t len = A.get_no_of_components();
68 for (lidia_size_t i = 0; i < len; i++)
69 w = A.get_base(i) * w;
70 return w;
71 }
72
73
74
75 template< class T, class MATRIX_TYPE >
76 inline const T
multiply_special(const math_vector<T> & w,const math_vector<T> & v,const T & mod) const77 sparse_fp_matrix_algorithms< T, MATRIX_TYPE >::multiply_special(const math_vector< T > &w,
78 const math_vector< T > &v,
79 const T &mod) const
80 {
81 T TMP, RES = 0;
82 for (register lidia_size_t i = 0; i < v.size(); i++) {
83 mult_mod(TMP, v[i], w[i], mod);
84 add_mod(RES, RES, TMP, mod);
85 }
86 return RES;
87 }
88
89
90
91 template< class T, class MATRIX_TYPE >
92 void
multiply_special(math_vector<T> & w,const base_power_product<ring_matrix<T>,lidia_size_t> & B,const T & mod) const93 sparse_fp_matrix_algorithms< T, MATRIX_TYPE >::multiply_special (math_vector< T > &w,
94 const base_power_product< ring_matrix < T >, lidia_size_t > &B,
95 const T &mod) const
96 {
97 lidia_size_t *itmp;
98 T *vtmp;
99
100 T TMP, RES;
101 for (lidia_size_t l = B.get_no_of_components() - 1; l >= 0; l--) {
102 math_vector< T > v((B.get_base(l)).rows, (B.get_base(l)).rows);
103 lidia_size_t i = 0;
104 for (i = 0; i < (B.get_base(l)).rows; i++) {
105 itmp = (B.get_base(l)).index[i];
106 vtmp = (B.get_base(l)).value[i];
107
108 RES = 0;
109 for (lidia_size_t j = 0; j < (B.get_base(l)).value_counter[i]; j++) {
110 if (vtmp[i] == 1)
111 add_mod(RES, RES, w[itmp[j]], mod);
112 else {
113 mult_mod(TMP, vtmp[j], w[itmp[j]], mod);
114 add_mod(RES, RES, TMP, mod);
115 }
116 }
117 v[i] = RES;
118 }
119 w = v;
120 }
121 }
122
123
124
125 template< class T, class MATRIX_TYPE >
126 void
multiply_special(math_vector<T> & w,const MATRIX_TYPE & B,const T & mod) const127 sparse_fp_matrix_algorithms< T, MATRIX_TYPE >::multiply_special (math_vector< T > &w,
128 const MATRIX_TYPE &B,
129 const T &mod) const
130 {
131 lidia_size_t *itmp;
132 T *vtmp;
133 T TMP, RES;
134
135 math_vector< T > v(B.rows, B.rows);
136 lidia_size_t i = 0;
137 for (i = 0; i < B.rows; i++) {
138 itmp = B.index[i];
139 vtmp = B.value[i];
140
141 RES = 0;
142 for (lidia_size_t j = 0; j < B.value_counter[i]; j++) {
143 if (vtmp[i] == 1)
144 add_mod(RES, RES, w[itmp[j]], mod);
145 else if (vtmp[i] == -1)
146 sub_mod(RES, RES, w[itmp[j]], mod);
147 else {
148 mult_mod(TMP, vtmp[j], w[itmp[j]], mod);
149 add_mod(RES, RES, TMP, mod);
150 }
151 }
152 v[i] = RES;
153 }
154 w = v;
155 }
156
157
158
159 template< class T, class MATRIX_TYPE >
160 void
multiply_special(math_vector<T> & v,const math_vector<T> & w,const MATRIX_TYPE & B,const T & mod) const161 sparse_fp_matrix_algorithms< T, MATRIX_TYPE >::multiply_special(math_vector< T > &v,
162 const math_vector< T > &w,
163 const MATRIX_TYPE &B,
164 const T &mod) const
165 {
166 if (v.get_size() < w.get_size()) {
167 v.set_capacity(w.get_size());
168 v.set_size(w.get_size());
169 }
170
171 lidia_size_t *itmp;
172 T *vtmp;
173 T TMP, RES;
174 lidia_size_t i = 0;
175 for (i = 0; i < B.rows; i++) {
176 itmp = B.index[i];
177 vtmp = B.value[i];
178
179 RES = 0;
180 for (lidia_size_t j = 0; j < B.value_counter[i]; j++) {
181 if (vtmp[j] == 1)
182 add_mod(RES, RES, w[itmp[j]], mod);
183 else {
184 mult_mod(TMP, vtmp[j], w[itmp[j]], mod);
185 add_mod(RES, RES, TMP, mod);
186 }
187 }
188 v[i] = RES;
189 }
190 }
191
192
193
194 //
195 // berlekamp massay algorithm
196 //
197
198 template< class T, class MATRIX_TYPE >
199 T *
200 sparse_fp_matrix_algorithms< T, MATRIX_TYPE >::berlekamp_massay(T *s, lidia_size_t n, const T &mod) const
201 {
202 // Memory allocation
203 T D, Dold, TMP2;
204 T *C = new T[n + 1];
205 T *B = new T[n + 1];
206 T *tmp = new T[n + 1];
207 lidia_size_t m, N, i, L;
208
209 for (i = 0; i <= n; i++) {
210 C[i] = 0;
211 B[i] = 0;
212 tmp[i] = 0;
213 }
214
215 // Initialization
216 C[0] = 1;
217 L = 0;
218 m = -1;
219 B[0] = 1;
220 N = 0;
221 Dold = 1;
222
223 // Iteration
224 while (N < n) {
225 //Compute the next discrepancy
226 best_remainder(D, s[N], mod);
227 for (i = 1; i <= L; i++) {
228 //D += (C[i]*s[N-i]); D %= mod;
229 mult_mod(TMP2, C[i], s[N-i], mod);
230 add_mod(D, D, TMP2, mod);
231 }
232
233 if (D != 0) {
234 for (i = 0; i <= L; i++)
235 tmp[i] = C[i];
236
237 for (i = 0; i <= L; i++) {
238 inv_mod(TMP2, Dold, mod);
239 mult_mod(TMP2, D, TMP2, mod);
240 mult_mod(TMP2, TMP2, B[i], mod);
241 sub_mod(C[i + N - m], C[i + N - m], TMP2, mod); // -= (D*TMP2)*B[i]; C[i + N - m] %= mod;
242 }
243 if (L*2 <= N+1) {
244 L = N + 1 - L;
245 m = N;
246 for (i = 0; i <= L; i++)
247 B[i] = tmp[i];
248 Dold = D;
249 }
250 }
251 N++;
252 }
253 delete[] B;
254 delete[] tmp;
255
256 T *RES = new T[L+2];
257 RES[0] = L;
258
259 for (i = 0; i <= L; i++)
260 RES[i+1] = C[L-i];
261 delete[] C;
262 return RES;
263 }
264
265
266
267 //
268 // column step form
269 //
270
271 template< class T, class MATRIX_TYPE >
272 inline int
STF(MATRIX_TYPE &,const T &) const273 sparse_fp_matrix_algorithms< T, MATRIX_TYPE >::STF(MATRIX_TYPE &, const T &) const
274 {
275 return 0;
276 }
277
278
279
280 //
281 // rank
282 //
283
284 template< class T, class MATRIX_TYPE >
285 inline lidia_size_t
rank(MATRIX_TYPE & A,const T &) const286 sparse_fp_matrix_algorithms< T, MATRIX_TYPE >::rank(MATRIX_TYPE &A, const T &) const
287 {
288 return A.columns;
289 }
290
291
292
293 //
294 // ran kand linearly independent rows or columns
295 //
296
297 template< class T, class MATRIX_TYPE >
298 inline lidia_size_t *
lininr(MATRIX_TYPE &,const T &) const299 sparse_fp_matrix_algorithms< T, MATRIX_TYPE >::lininr(MATRIX_TYPE &, const T &) const
300 {
301 return NULL;
302 }
303
304
305
306 template< class T, class MATRIX_TYPE >
307 inline lidia_size_t *
lininc(MATRIX_TYPE &,const T &) const308 sparse_fp_matrix_algorithms< T, MATRIX_TYPE >::lininc(MATRIX_TYPE &, const T &) const
309 {
310 return NULL;
311 }
312
313
314
315 //
316 // adjoint matrix
317 //
318
319 template< class T, class MATRIX_TYPE >
320 inline void
adj(MATRIX_TYPE &,const T &) const321 sparse_fp_matrix_algorithms< T, MATRIX_TYPE >::adj(MATRIX_TYPE &, const T &) const
322 {
323
324 }
325
326
327
328 //
329 // determinant
330 //
331
332 template< class T, class MATRIX_TYPE >
333 inline const T
det(MATRIX_TYPE & A,const T & mod) const334 sparse_fp_matrix_algorithms< T, MATRIX_TYPE >::det(MATRIX_TYPE &A, const T &mod) const
335 {
336 T *data = charpoly(A, mod);
337 T ret;
338 if (A.rows % 2 == 1)
339 ret = -data[0];
340 else
341 ret = data[0];
342 delete[] data;
343 best_remainder(ret, ret, mod);
344 return ret;
345 }
346
347
348
349 template< class T, class MATRIX_TYPE >
350 inline const T
det(const base_power_product<ring_matrix<T>,lidia_size_t> & A,const T & mod) const351 sparse_fp_matrix_algorithms< T, MATRIX_TYPE >::det (const base_power_product< ring_matrix < T >, lidia_size_t > &A,
352 const T &mod) const
353 {
354 T *data = charpoly(A, mod);
355 T ret;
356 if (A.get_base(0).rows % 2 == 1)
357 ret = -data[0];
358 else
359 ret = data[0];
360 delete[] data;
361 return ret;
362 }
363
364
365
366 //
367 // Hessenberg form
368 //
369
370 template< class T, class MATRIX_TYPE >
371 inline void
HBF(MATRIX_TYPE &,const T &) const372 sparse_fp_matrix_algorithms< T, MATRIX_TYPE >::HBF (MATRIX_TYPE &,
373 const T &) const
374 {
375 }
376
377
378
379 //
380 // characteristic polynomial
381 //
382
383 template< class T, class MATRIX_TYPE >
384 T *
385 sparse_fp_matrix_algorithms< T, MATRIX_TYPE >::charpoly (MATRIX_TYPE &A,
386 const T &mod) const
387 {
388 // A symmetrische n x n Matrix
389 lidia_size_t n = A.rows;
390
391 random_generator gen;
392 long TMPlong;
393 Fp_polynomial res;
394 res.set_modulus(mod);
395 math_vector< T > u(n, n), b(n, n);
396 bool SW = false;
397 math_vector< T > w(n, n);
398 T *s = new T[n+n];
399
400 lidia_size_t i;
401 do {
402 for (i = 0; i < n; i++) {
403 gen >> TMPlong;
404 best_remainder(u[i], TMPlong, mod);
405 gen >> TMPlong;
406 best_remainder(b[i], TMPlong, mod);
407 }
408
409 w = b;
410 //s[0] = (w*u) % mod;
411 s[0] = multiply_special(w, u, mod);
412
413
414 for (i = 1; i < n+n; i++) {
415 //w = (A * w) % mod;
416 multiply_special(w, A, mod);
417 //s[i] = (u * w) % mod;
418 s[i] = multiply_special(u, w, mod);
419 }
420
421 T *c = berlekamp_massay(s, n+n, mod);
422
423
424 bigint LEN = c[0];
425 lidia_size_t len;
426 LEN.sizetify(len);
427
428 polynomial< bigint > va(&(c[1]), len);
429
430 Fp_polynomial vb(va, mod);
431 delete[] c;
432
433 if (res.degree() <= 0)
434 res = vb;
435 else {
436 res = (res * vb)/gcd(res, vb);
437
438 if (res == vb)
439 SW = true;
440 }
441 //std::cout << "grad = " << res.degree() << std::endl;
442 }
443 while (res.degree() != n && SW == false);
444 delete[] s;
445
446 return Fp_polynomial_convert(res, mod);
447 }
448
449
450
451 template< class T, class MATRIX_TYPE >
452 T *
453 sparse_fp_matrix_algorithms< T, MATRIX_TYPE >::charpoly (const base_power_product< ring_matrix < T >, lidia_size_t > &A,
454 const T &mod) const
455 {
456 // A symmetrische n x n Matrix
457 lidia_size_t n = A.get_base(0).rows;
458 random_generator gen;
459 long TMPlong;
460 Fp_polynomial res;
461 res.set_modulus(mod);
462 bool SW = false;
463 math_vector< T > u(n, n), b(n, n);
464
465 math_vector< T > w(n, n);
466 T *s = new T[n+n];
467
468 lidia_size_t i;
469 do {
470 for (i = 0; i < n; i++) {
471 gen >> TMPlong;
472 best_remainder(u[i], TMPlong, mod);
473 gen >> TMPlong;
474 best_remainder(b[i], TMPlong, mod);
475 }
476
477 w = b;
478 //s[0] = (w * u);
479 s[0] = multiply_special(w, u, mod);
480
481 for (i = 1; i < n+n; i++) {
482 //w = A * w;
483 multiply_special(w, A, mod);
484 //s[i] = (u * w) % mod;
485 s[i] = multiply_special(u, w, mod);
486 }
487
488 T *c = berlekamp_massay(s, n+n, mod);
489 bigint LEN = c[0];
490 lidia_size_t len;
491 LEN.sizetify(len);
492
493 polynomial< bigint > va(&(c[1]), len);
494 Fp_polynomial vb(va, mod);
495
496 delete[] c;
497
498 if (res.degree() <= 0)
499 res = vb;
500 else {
501 res = (res * vb)/gcd(res, vb);
502 if (res == vb)
503 SW = true;
504 }
505 }
506 while (res.degree() != n && SW == false);
507
508 delete[] s;
509
510 return Fp_polynomial_convert(res, mod);
511 }
512
513
514
515 template< class T, class MATRIX_TYPE >
516 bool
lanczos(const MATRIX_TYPE & A,math_vector<T> & x,const math_vector<T> & b,T & mod) const517 sparse_fp_matrix_algorithms< T, MATRIX_TYPE >::lanczos (const MATRIX_TYPE &A,
518 math_vector< T > &x,
519 const math_vector< T > &b,
520 T &mod) const
521 {
522 // A symmetrische n x n Matrix
523 lidia_size_t i, n = A.columns;
524
525 math_vector< T > *w = new math_vector< T > [n+n];
526 math_vector< T > *v = new math_vector< T > [n+n];
527
528 for (i = 0; i < n+n; i++) {
529 w[i].set_capacity(n);
530 w[i].set_size(n);
531 v[i].set_capacity(n);
532 v[i].set_size(n);
533 }
534
535 T *wv = new T[n+n];
536
537 // Init
538 w[0] = b % mod;
539 multiply_special(v[1], w[0], A, mod);
540
541 T TMP1, TMP2, TMP3, TMP4, TMP = multiply_special(w[0], v[1], mod);
542 inv_mod(wv[0], TMP, mod);
543
544 TMP1 = multiply_special(v[1], v[1], mod);
545
546 //w[1] = (v[1] - (TMP1*wv[0])*w[0]) % mod;
547 mult_mod(TMP3, TMP1, wv[0], mod);
548 for (i = 0; i < n; i++) {
549 mult_mod(TMP2, TMP3, w[0][i], mod);
550 sub_mod(w[1][i], v[1][i], TMP2, mod);
551 }
552
553 //x = ((w[0]*b)*wv[0]*w[0]) % mod;
554 TMP2 = multiply_special(w[0], b, mod);
555 mult_mod(TMP2, TMP2, wv[0], mod);
556 for (i = 0; i < n; i++)
557 mult_mod(x[i], TMP2, w[0][i], mod);
558
559 i = 1;
560 lidia_size_t j;
561
562 while (TMP != 0) {
563 //std::cout << "Lanczos Iteration: " << i << std::endl;
564 multiply_special(v[i+1], w[i], A, mod);
565
566 TMP = multiply_special(w[i], v[i+1], mod);
567 if (TMP != 0) {
568 inv_mod(wv[i], TMP, mod);
569
570 TMP1 = multiply_special(v[i+1], v[i+1], mod);
571 mult_mod(TMP1, TMP1, wv[i], mod);
572 TMP2 = multiply_special(w[i], v[i+1], mod);
573 mult_mod(TMP2, TMP2, wv[i-1], mod);
574
575 //w[i+1] = (v[i+1] - TMP1*w[i] - TMP2*w[i-1]) % mod;
576 for (j = 0; j < n; j++) {
577 mult_mod(TMP3, TMP1, w[i][j], mod);
578 mult_mod(TMP4, TMP2, w[i-1][j], mod);
579 add_mod(TMP3, TMP3, TMP4, mod);
580 sub_mod(w[i+1][j], v[i+1][j], TMP3, mod);
581 }
582 // Update of solution
583 //x = (x+((w[i]*b)*wv[i])*w[i]) % mod;
584 TMP2 = multiply_special(w[i], b, mod);
585 mult_mod(TMP2, TMP2, wv[i], mod);
586 for (j = 0; j < n; j++) {
587 mult_mod(TMP3, TMP2, w[i][j], mod);
588 add_mod(x[j], x[j], TMP3, mod);
589 }
590 i++;
591 }
592 }
593 delete[] w;
594 delete[] v;
595 delete[] wv;
596 return true;
597 }
598
599
600
601 template< class T, class MATRIX_TYPE >
602 T
603 sparse_fp_matrix_algorithms< T, MATRIX_TYPE >::lanczos_ZmZ (const MATRIX_TYPE &A,
604 math_vector< T > &x,
605 const math_vector< T > &b,
606 T &mod) const
607 {
608 // A symmetrische n x n Matrix
609 lidia_size_t i, n = A.columns;
610 T factor_1, factor_2, factor_3, factor_1_2;
611
612 math_vector< T > *w = new math_vector< T > [n+n];
613 math_vector< T > *v = new math_vector< T > [n+n];
614
615 for (i = 0; i < n+n; i++) {
616 w[i].set_capacity(n);
617 w[i].set_size(n);
618 v[i].set_capacity(n);
619 v[i].set_size(n);
620 }
621
622 // Init
623 w[0] = b % mod;
624
625 multiply_special(v[1], w[0], A, mod);
626
627 T TMP1, TMP2, TMP3, TMP4;
628 T TMP = multiply_special(w[0], v[1], mod);
629
630 // FAKTOR_1, FACTOR_3
631 factor_1 = TMP;
632 factor_3 = factor_1;
633
634 TMP1 = multiply_special(v[1], v[1], mod);
635
636 //w[1]
637 for (i = 0; i < n; i++) {
638 mult_mod(TMP2, TMP1, w[0][i], mod);
639 mult_mod(TMP3, TMP, v[1][i], mod);
640 sub_mod(w[1][i], TMP3, TMP2, mod);
641 }
642
643 //FACTOR_2
644 mult_mod(factor_2, factor_1, factor_1, mod);
645
646 //x
647 TMP2 = multiply_special(w[0], b, mod);
648 for (i = 0; i < n; i++)
649 mult_mod(x[i], TMP2, w[0][i], mod);
650
651 i = 1;
652 lidia_size_t j;
653
654 while (TMP != 0) {
655 multiply_special(v[i+1], w[i], A, mod);
656
657 TMP = multiply_special(w[i], v[i+1], mod);
658
659 if (TMP != 0) {
660
661 TMP1 = multiply_special(v[i+1], v[i+1], mod);
662 TMP2 = multiply_special(w[i], v[i+1], mod);
663
664 //FACTOR_1
665 factor_1 = TMP2;
666 mult_mod(factor_3, factor_3, factor_2, mod);
667 mult_mod(factor_1_2, factor_2, factor_1, mod);
668
669 mult_mod(TMP1, TMP1, factor_2, mod);
670 mult_mod(TMP2, TMP2, factor_1, mod);
671
672 //w[i+1]
673 for (j = 0; j < n; j++) {
674 mult_mod(TMP3, TMP1, w[i][j], mod);
675 mult_mod(TMP4, TMP2, w[i-1][j], mod);
676 add_mod(TMP3, TMP3, TMP4, mod);
677 mult_mod(TMP4, factor_1_2, v[i+1][j], mod);
678 sub_mod(w[i+1][j], TMP4, TMP3, mod);
679 }
680
681 // Update of solution
682 TMP2 = multiply_special(w[i], b, mod);
683 for (j = 0; j < n; j++) {
684 mult_mod(TMP3, TMP2*factor_3, w[i][j], mod);
685 mult_mod(TMP4, factor_1_2, x[j], mod);
686 add_mod(x[j], TMP4, TMP3, mod);
687 }
688
689 //FACTOR_2
690 mult_mod(factor_3, factor_3, factor_1, mod);
691 mult_mod(factor_2, factor_2, factor_1, mod);
692 mult_mod(factor_2, factor_2, factor_1, mod);
693
694 i++;
695 }
696 }
697
698 delete[] w;
699 delete[] v;
700 return factor_3;
701 }
702
703
704
705 //
706 // wiedemann algorithm
707 //
708
709 template< class T, class MATRIX_TYPE >
710 bool
wiedemann(const ring_matrix<T> & A,math_vector<T> & x,const math_vector<T> & b,T & mod) const711 sparse_fp_matrix_algorithms< T, MATRIX_TYPE >::wiedemann (const ring_matrix< T > &A,
712 math_vector< T > &x,
713 const math_vector< T > &b,
714 T &mod) const
715 {
716 // A symmetrische n x n Matrix
717 lidia_size_t n = A.columns;
718
719 // vector u
720 math_vector< T > u(n, n);
721 long TMPlong;
722 random_generator gen;
723 lidia_size_t i;
724
725 for (i = 0; i < n; i++) {
726 gen >> TMPlong;
727 u[i] = TMPlong % mod;
728 }
729
730 math_vector< T > *w = new math_vector< T > [n+n];
731 T *s = new T[n+n];
732
733 w[0] = b;
734 s[0] = (w[0]*u);
735
736 for (i = 1; i < n+n; i++) {
737 w[i] = A * w[i-1] % mod;
738 s[i] = (u* w[i]) % mod;
739 }
740
741 T *c = berlekamp_massay(s, n+n, mod);
742
743 for (i = 0; i < x.size(); i++)
744 x[i] = 0;
745 for (i = 2; i <= c[0]+1; i++) {
746 x += c[i]*w[i-2];
747 }
748
749 T TMP;
750 inv_mod(TMP, c[1], mod);
751 x = x * (-TMP) % mod;
752 return true;
753 }
754
755
756
757 template< class T, class MATRIX_TYPE >
758 bool
wiedemann(const base_power_product<ring_matrix<T>,lidia_size_t> & A,math_vector<T> & x,const math_vector<T> & b,T & mod) const759 sparse_fp_matrix_algorithms< T, MATRIX_TYPE >::wiedemann (const base_power_product< ring_matrix < T >, lidia_size_t > &A,
760 math_vector< T > &x,
761 const math_vector< T > &b,
762 T &mod) const
763 {
764 // A symmetrische n x n Matrix
765 lidia_size_t n = A.get_base(0).columns;
766
767 // vector u
768 math_vector< T > u(n, n);
769 long TMPlong;
770 random_generator gen;
771 lidia_size_t i;
772
773 for (i = 0; i < n; i++) {
774 gen >> TMPlong;
775 u[i] = TMPlong % mod;
776 }
777
778 math_vector< T > *w = new math_vector< T > [n+n];
779 T *s = new T[n+n];
780
781 w[0] = b;
782 s[0] = (w[0]*u);
783
784 for (i = 1; i < n+n; i++) {
785 w[i] = A * w[i-1] % mod;
786 s[i] = (u* w[i]) % mod;
787 }
788
789 T *c = berlekamp_massay(s, n+n, mod);
790
791 for (i = 0; i < x.size(); i++)
792 x[i] = 0;
793
794 for (i = 2; i <= c[0]+1; i++) {
795 x += c[i]*w[i-2];
796 }
797
798 T TMP;
799 inv_mod(TMP, c[1], mod);
800 x = x * (-TMP) % mod;
801 return true;
802 }
803
804
805
806 template< class T, class MATRIX_TYPE >
807 bool
conjugate_gradient(const ring_matrix<T> & A,math_vector<T> & x,const math_vector<T> & b,T & mod) const808 sparse_fp_matrix_algorithms< T, MATRIX_TYPE >::conjugate_gradient (const ring_matrix< T > &A,
809 math_vector< T > &x,
810 const math_vector< T > &b,
811 T &mod) const
812 {
813 // A symmetrische n x n Matrix
814 lidia_size_t n = A.columns;
815
816 // vector u
817 math_vector< T > u(n, n), p(n, n), r(n, n);
818 long TMPlong;
819 random_generator gen;
820 for (register lidia_size_t i = 0; i < n; i++) {
821 gen >> TMPlong;
822 u[i] = TMPlong % mod;
823 }
824
825 T TMP, TMP1, a, beta;
826 // Init
827 p = (b - A*u) % mod;
828 r = p;
829
830 TMP = (r*r);
831 while (TMP != 0) {
832 // PART 1
833 TMP1 = (p * (A * p)) % mod;
834
835 inv_mod(TMP1, TMP1, mod);
836 a = TMP*TMP1;
837
838 u = (u + a * p) % mod;
839 r = (r - a*A*p) % mod;
840
841 // PART 2
842 inv_mod(TMP1, TMP, mod);
843 TMP = (r*r);
844 beta = (TMP*TMP1) % mod;
845
846 p = (r + beta*p) % mod;
847
848 }
849 x = u;
850 return true;
851 }
852
853
854
855 #ifdef LIDIA_NAMESPACE
856 # ifndef IN_NAMESPACE_LIDIA
857 } // end of namespace LiDIA
858 # endif
859 #endif
860
861
862
863 #endif // LIDIA_SPARSE_FP_MATRIX_ALGORITHMS_CC_GUARD_
864