source: git/ntl/src/mat_ZZ_p.c @ 2cfffe

spielwiese
Last change on this file since 2cfffe was 2cfffe, checked in by Hans Schönemann <hannes@…>, 21 years ago
This commit was generated by cvs2svn to compensate for changes in r6316, which included commits to RCS files with non-trunk default branches. git-svn-id: file:///usr/local/Singular/svn/trunk@6317 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 15.6 KB
Line 
1
2#include <NTL/mat_ZZ_p.h>
3#include <NTL/vec_ZZVec.h>
4#include <NTL/vec_long.h>
5
6#include <NTL/new.h>
7
8NTL_START_IMPL
9
10NTL_matrix_impl(ZZ_p,vec_ZZ_p,vec_vec_ZZ_p,mat_ZZ_p)
11NTL_io_matrix_impl(ZZ_p,vec_ZZ_p,vec_vec_ZZ_p,mat_ZZ_p)
12NTL_eq_matrix_impl(ZZ_p,vec_ZZ_p,vec_vec_ZZ_p,mat_ZZ_p)
13
14
15 
16void add(mat_ZZ_p& X, const mat_ZZ_p& A, const mat_ZZ_p& B) 
17{ 
18   long n = A.NumRows(); 
19   long m = A.NumCols(); 
20 
21   if (B.NumRows() != n || B.NumCols() != m)   
22      Error("matrix add: dimension mismatch"); 
23 
24   X.SetDims(n, m); 
25 
26   long i, j; 
27   for (i = 1; i <= n; i++)   
28      for (j = 1; j <= m; j++) 
29         add(X(i,j), A(i,j), B(i,j)); 
30} 
31 
32void sub(mat_ZZ_p& X, const mat_ZZ_p& A, const mat_ZZ_p& B) 
33{ 
34   long n = A.NumRows(); 
35   long m = A.NumCols(); 
36 
37   if (B.NumRows() != n || B.NumCols() != m) 
38      Error("matrix sub: dimension mismatch"); 
39 
40   X.SetDims(n, m); 
41 
42   long i, j; 
43   for (i = 1; i <= n; i++) 
44      for (j = 1; j <= m; j++) 
45         sub(X(i,j), A(i,j), B(i,j)); 
46} 
47
48void negate(mat_ZZ_p& X, const mat_ZZ_p& A) 
49{ 
50   long n = A.NumRows(); 
51   long m = A.NumCols(); 
52 
53 
54   X.SetDims(n, m); 
55 
56   long i, j; 
57   for (i = 1; i <= n; i++) 
58      for (j = 1; j <= m; j++) 
59         negate(X(i,j), A(i,j)); 
60} 
61 
62void mul_aux(mat_ZZ_p& X, const mat_ZZ_p& A, const mat_ZZ_p& B) 
63{ 
64   long n = A.NumRows(); 
65   long l = A.NumCols(); 
66   long m = B.NumCols(); 
67 
68   if (l != B.NumRows()) 
69      Error("matrix mul: dimension mismatch"); 
70 
71   X.SetDims(n, m); 
72 
73   long i, j, k; 
74   ZZ acc, tmp; 
75 
76   for (i = 1; i <= n; i++) { 
77      for (j = 1; j <= m; j++) { 
78         clear(acc); 
79         for(k = 1; k <= l; k++) { 
80            mul(tmp, rep(A(i,k)), rep(B(k,j))); 
81            add(acc, acc, tmp); 
82         } 
83         conv(X(i,j), acc); 
84      } 
85   } 
86} 
87 
88 
89void mul(mat_ZZ_p& X, const mat_ZZ_p& A, const mat_ZZ_p& B) 
90{ 
91   if (&X == &A || &X == &B) { 
92      mat_ZZ_p tmp; 
93      mul_aux(tmp, A, B); 
94      X = tmp; 
95   } 
96   else 
97      mul_aux(X, A, B); 
98} 
99 
100 
101static
102void mul_aux(vec_ZZ_p& x, const mat_ZZ_p& A, const vec_ZZ_p& b) 
103{ 
104   long n = A.NumRows(); 
105   long l = A.NumCols(); 
106 
107   if (l != b.length()) 
108      Error("matrix mul: dimension mismatch"); 
109 
110   x.SetLength(n); 
111 
112   long i, k; 
113   ZZ acc, tmp; 
114 
115   for (i = 1; i <= n; i++) { 
116      clear(acc); 
117      for (k = 1; k <= l; k++) { 
118         mul(tmp, rep(A(i,k)), rep(b(k))); 
119         add(acc, acc, tmp); 
120      } 
121      conv(x(i), acc); 
122   } 
123} 
124 
125 
126void mul(vec_ZZ_p& x, const mat_ZZ_p& A, const vec_ZZ_p& b) 
127{ 
128   if (&b == &x || A.position(b) != -1) {
129      vec_ZZ_p tmp;
130      mul_aux(tmp, A, b);
131      x = tmp;
132   }
133   else
134      mul_aux(x, A, b);
135} 
136
137static
138void mul_aux(vec_ZZ_p& x, const vec_ZZ_p& a, const mat_ZZ_p& B) 
139{ 
140   long n = B.NumRows(); 
141   long l = B.NumCols(); 
142 
143   if (n != a.length()) 
144      Error("matrix mul: dimension mismatch"); 
145 
146   x.SetLength(l); 
147 
148   long i, k; 
149   ZZ acc, tmp; 
150 
151   for (i = 1; i <= l; i++) { 
152      clear(acc); 
153      for (k = 1; k <= n; k++) { 
154         mul(tmp, rep(a(k)), rep(B(k,i)));
155         add(acc, acc, tmp); 
156      } 
157      conv(x(i), acc); 
158   } 
159} 
160
161void mul(vec_ZZ_p& x, const vec_ZZ_p& a, const mat_ZZ_p& B)
162{
163   if (&a == &x || B.position(a) != -1) {
164      vec_ZZ_p tmp;
165      mul_aux(tmp, a, B);
166      x = tmp;
167   }
168   else
169      mul_aux(x, a, B);
170}
171
172     
173 
174void ident(mat_ZZ_p& X, long n) 
175{ 
176   X.SetDims(n, n); 
177   long i, j; 
178 
179   for (i = 1; i <= n; i++) 
180      for (j = 1; j <= n; j++) 
181         if (i == j) 
182            set(X(i, j)); 
183         else 
184            clear(X(i, j)); 
185} 
186
187
188void determinant(ZZ_p& d, const mat_ZZ_p& M_in)
189{
190   long k, n;
191   long i, j;
192   long pos;
193   ZZ t1, t2;
194   ZZ *x, *y;
195
196   const ZZ& p = ZZ_p::modulus();
197
198   n = M_in.NumRows();
199
200   if (M_in.NumCols() != n)
201      Error("determinant: nonsquare matrix");
202
203   if (n == 0) {
204      set(d);
205      return;
206   }
207
208   vec_ZZVec M;
209   sqr(t1, p);
210   mul(t1, t1, n);
211
212   M.SetLength(n);
213   for (i = 0; i < n; i++) {
214      M[i].SetSize(n, t1.size());
215      for (j = 0; j < n; j++)
216         M[i][j] = rep(M_in[i][j]);
217   }
218
219   ZZ det;
220   set(det);
221
222   for (k = 0; k < n; k++) {
223      pos = -1;
224      for (i = k; i < n; i++) {
225         rem(t1, M[i][k], p);
226         M[i][k] = t1;
227         if (pos == -1 && !IsZero(t1))
228            pos = i;
229      }
230
231      if (pos != -1) {
232         if (k != pos) {
233            swap(M[pos], M[k]);
234            NegateMod(det, det, p);
235         }
236
237         MulMod(det, det, M[k][k], p);
238
239         // make M[k, k] == -1 mod p, and make row k reduced
240
241         InvMod(t1, M[k][k], p);
242         NegateMod(t1, t1, p);
243         for (j = k+1; j < n; j++) {
244            rem(t2, M[k][j], p);
245            MulMod(M[k][j], t2, t1, p);
246         }
247
248         for (i = k+1; i < n; i++) {
249            // M[i] = M[i] + M[k]*M[i,k]
250
251            t1 = M[i][k];   // this is already reduced
252
253            x = M[i].elts() + (k+1);
254            y = M[k].elts() + (k+1);
255
256            for (j = k+1; j < n; j++, x++, y++) {
257               // *x = *x + (*y)*t1
258
259               mul(t2, *y, t1);
260               add(*x, *x, t2);
261            }
262         }
263      }
264      else {
265         clear(d);
266         return;
267      }
268   }
269
270   conv(d, det);
271}
272
273long IsIdent(const mat_ZZ_p& A, long n)
274{
275   if (A.NumRows() != n || A.NumCols() != n)
276      return 0;
277
278   long i, j;
279
280   for (i = 1; i <= n; i++)
281      for (j = 1; j <= n; j++)
282         if (i != j) {
283            if (!IsZero(A(i, j))) return 0;
284         }
285         else {
286            if (!IsOne(A(i, j))) return 0;
287         }
288
289   return 1;
290}
291           
292
293void transpose(mat_ZZ_p& X, const mat_ZZ_p& A)
294{
295   long n = A.NumRows();
296   long m = A.NumCols();
297
298   long i, j;
299
300   if (&X == & A) {
301      if (n == m)
302         for (i = 1; i <= n; i++)
303            for (j = i+1; j <= n; j++)
304               swap(X(i, j), X(j, i));
305      else {
306         mat_ZZ_p tmp;
307         tmp.SetDims(m, n);
308         for (i = 1; i <= n; i++)
309            for (j = 1; j <= m; j++)
310               tmp(j, i) = A(i, j);
311         X.kill();
312         X = tmp;
313      }
314   }
315   else {
316      X.SetDims(m, n);
317      for (i = 1; i <= n; i++)
318         for (j = 1; j <= m; j++)
319            X(j, i) = A(i, j);
320   }
321}
322   
323
324void solve(ZZ_p& d, vec_ZZ_p& X, 
325           const mat_ZZ_p& A, const vec_ZZ_p& b)
326
327{
328   long n = A.NumRows();
329   if (A.NumCols() != n)
330      Error("solve: nonsquare matrix");
331
332   if (b.length() != n)
333      Error("solve: dimension mismatch");
334
335   if (n == 0) {
336      set(d);
337      X.SetLength(0);
338      return;
339   }
340
341   long i, j, k, pos;
342   ZZ t1, t2;
343   ZZ *x, *y;
344
345   const ZZ& p = ZZ_p::modulus();
346
347   vec_ZZVec M;
348   sqr(t1, p);
349   mul(t1, t1, n);
350
351   M.SetLength(n);
352
353   for (i = 0; i < n; i++) {
354      M[i].SetSize(n+1, t1.size());
355      for (j = 0; j < n; j++) 
356         M[i][j] = rep(A[j][i]);
357      M[i][n] = rep(b[i]);
358   }
359
360   ZZ det;
361   set(det);
362
363   for (k = 0; k < n; k++) {
364      pos = -1;
365      for (i = k; i < n; i++) {
366         rem(t1, M[i][k], p);
367         M[i][k] = t1;
368         if (pos == -1 && !IsZero(t1)) {
369            pos = i;
370         }
371      }
372
373      if (pos != -1) {
374         if (k != pos) {
375            swap(M[pos], M[k]);
376            NegateMod(det, det, p);
377         }
378
379         MulMod(det, det, M[k][k], p);
380
381         // make M[k, k] == -1 mod p, and make row k reduced
382
383         InvMod(t1, M[k][k], p);
384         NegateMod(t1, t1, p);
385         for (j = k+1; j <= n; j++) {
386            rem(t2, M[k][j], p);
387            MulMod(M[k][j], t2, t1, p);
388         }
389
390         for (i = k+1; i < n; i++) {
391            // M[i] = M[i] + M[k]*M[i,k]
392
393            t1 = M[i][k];   // this is already reduced
394
395            x = M[i].elts() + (k+1);
396            y = M[k].elts() + (k+1);
397
398            for (j = k+1; j <= n; j++, x++, y++) {
399               // *x = *x + (*y)*t1
400
401               mul(t2, *y, t1);
402               add(*x, *x, t2);
403            }
404         }
405      }
406      else {
407         clear(d);
408         return;
409      }
410   }
411
412   X.SetLength(n);
413   for (i = n-1; i >= 0; i--) {
414      clear(t1);
415      for (j = i+1; j < n; j++) {
416         mul(t2, rep(X[j]), M[i][j]);
417         add(t1, t1, t2);
418      }
419      sub(t1, t1, M[i][n]);
420      conv(X[i], t1);
421   }
422
423   conv(d, det);
424}
425
426void inv(ZZ_p& d, mat_ZZ_p& X, const mat_ZZ_p& A)
427{
428   long n = A.NumRows();
429   if (A.NumCols() != n)
430      Error("inv: nonsquare matrix");
431
432   if (n == 0) {
433      set(d);
434      X.SetDims(0, 0);
435      return;
436   }
437
438   long i, j, k, pos;
439   ZZ t1, t2;
440   ZZ *x, *y;
441
442   const ZZ& p = ZZ_p::modulus();
443
444   vec_ZZVec M;
445   sqr(t1, p);
446   mul(t1, t1, n);
447
448   M.SetLength(n);
449
450   for (i = 0; i < n; i++) {
451      M[i].SetSize(2*n, t1.size());
452      for (j = 0; j < n; j++) {
453         M[i][j] = rep(A[i][j]);
454         clear(M[i][n+j]);
455      }
456      set(M[i][n+i]);
457   }
458
459   ZZ det;
460   set(det);
461
462   for (k = 0; k < n; k++) {
463      pos = -1;
464      for (i = k; i < n; i++) {
465         rem(t1, M[i][k], p);
466         M[i][k] = t1;
467         if (pos == -1 && !IsZero(t1)) {
468            pos = i;
469         }
470      }
471
472      if (pos != -1) {
473         if (k != pos) {
474            swap(M[pos], M[k]);
475            NegateMod(det, det, p);
476         }
477
478         MulMod(det, det, M[k][k], p);
479
480         // make M[k, k] == -1 mod p, and make row k reduced
481
482         InvMod(t1, M[k][k], p);
483         NegateMod(t1, t1, p);
484         for (j = k+1; j < 2*n; j++) {
485            rem(t2, M[k][j], p);
486            MulMod(M[k][j], t2, t1, p);
487         }
488
489         for (i = k+1; i < n; i++) {
490            // M[i] = M[i] + M[k]*M[i,k]
491
492            t1 = M[i][k];   // this is already reduced
493
494            x = M[i].elts() + (k+1);
495            y = M[k].elts() + (k+1);
496
497            for (j = k+1; j < 2*n; j++, x++, y++) {
498               // *x = *x + (*y)*t1
499
500               mul(t2, *y, t1);
501               add(*x, *x, t2);
502            }
503         }
504      }
505      else {
506         clear(d);
507         return;
508      }
509   }
510
511   X.SetDims(n, n);
512   for (k = 0; k < n; k++) {
513      for (i = n-1; i >= 0; i--) {
514         clear(t1);
515         for (j = i+1; j < n; j++) {
516            mul(t2, rep(X[j][k]), M[i][j]);
517            add(t1, t1, t2);
518         }
519         sub(t1, t1, M[i][n+k]);
520         conv(X[i][k], t1);
521      }
522   }
523
524   conv(d, det);
525}
526
527
528
529long gauss(mat_ZZ_p& M_in, long w)
530{
531   long k, l;
532   long i, j;
533   long pos;
534   ZZ t1, t2, t3;
535   ZZ *x, *y;
536
537   long n = M_in.NumRows();
538   long m = M_in.NumCols();
539
540   if (w < 0 || w > m)
541      Error("gauss: bad args");
542
543   const ZZ& p = ZZ_p::modulus();
544
545   vec_ZZVec M;
546   sqr(t1, p);
547   mul(t1, t1, n);
548
549   M.SetLength(n);
550   for (i = 0; i < n; i++) {
551      M[i].SetSize(m, t1.size());
552      for (j = 0; j < m; j++) {
553         M[i][j] = rep(M_in[i][j]);
554      }
555   }
556
557   l = 0;
558   for (k = 0; k < w && l < n; k++) {
559
560      pos = -1;
561      for (i = l; i < n; i++) {
562         rem(t1, M[i][k], p);
563         M[i][k] = t1;
564         if (pos == -1 && !IsZero(t1)) {
565            pos = i;
566         }
567      }
568
569      if (pos != -1) {
570         swap(M[pos], M[l]);
571
572         InvMod(t3, M[l][k], p);
573         NegateMod(t3, t3, p);
574
575         for (j = k+1; j < m; j++) {
576            rem(M[l][j], M[l][j], p);
577         }
578
579         for (i = l+1; i < n; i++) {
580            // M[i] = M[i] + M[l]*M[i,k]*t3
581
582            MulMod(t1, M[i][k], t3, p);
583
584            clear(M[i][k]);
585
586            x = M[i].elts() + (k+1);
587            y = M[l].elts() + (k+1);
588
589            for (j = k+1; j < m; j++, x++, y++) {
590               // *x = *x + (*y)*t1
591
592               mul(t2, *y, t1);
593               add(t2, t2, *x);
594               *x = t2;
595            }
596         }
597
598         l++;
599      }
600   }
601   
602   for (i = 0; i < n; i++)
603      for (j = 0; j < m; j++)
604         conv(M_in[i][j], M[i][j]);
605
606   return l;
607}
608
609long gauss(mat_ZZ_p& M)
610{
611   return gauss(M, M.NumCols());
612}
613
614void image(mat_ZZ_p& X, const mat_ZZ_p& A)
615{
616   mat_ZZ_p M;
617   M = A;
618   long r = gauss(M);
619   M.SetDims(r, M.NumCols());
620   X = M;
621}
622
623void kernel(mat_ZZ_p& X, const mat_ZZ_p& A)
624{
625   long m = A.NumRows();
626   long n = A.NumCols();
627
628   mat_ZZ_p M;
629   long r;
630
631   transpose(M, A);
632   r = gauss(M);
633
634   X.SetDims(m-r, m);
635
636   long i, j, k, s;
637   ZZ t1, t2;
638
639   ZZ_p T3;
640
641   vec_long D;
642   D.SetLength(m);
643   for (j = 0; j < m; j++) D[j] = -1;
644
645   vec_ZZ_p inverses;
646   inverses.SetLength(m);
647
648   j = -1;
649   for (i = 0; i < r; i++) {
650      do {
651         j++;
652      } while (IsZero(M[i][j]));
653
654      D[j] = i;
655      inv(inverses[j], M[i][j]); 
656   }
657
658   for (k = 0; k < m-r; k++) {
659      vec_ZZ_p& v = X[k];
660      long pos = 0;
661      for (j = m-1; j >= 0; j--) {
662         if (D[j] == -1) {
663            if (pos == k)
664               set(v[j]);
665            else
666               clear(v[j]);
667            pos++;
668         }
669         else {
670            i = D[j];
671
672            clear(t1);
673
674            for (s = j+1; s < m; s++) {
675               mul(t2, rep(v[s]), rep(M[i][s]));
676               add(t1, t1, t2);
677            }
678
679            conv(T3, t1);
680            mul(T3, T3, inverses[j]);
681            negate(v[j], T3);
682         }
683      }
684   }
685}
686   
687void mul(mat_ZZ_p& X, const mat_ZZ_p& A, const ZZ_p& b_in)
688{
689   NTL_ZZ_pRegister(b);
690   b = b_in;
691   long n = A.NumRows();
692   long m = A.NumCols();
693
694   X.SetDims(n, m);
695
696   long i, j;
697   for (i = 0; i < n; i++)
698      for (j = 0; j < m; j++)
699         mul(X[i][j], A[i][j], b);
700}
701   
702void mul(mat_ZZ_p& X, const mat_ZZ_p& A, long b_in)
703{
704   NTL_ZZ_pRegister(b);
705   b = b_in;
706   long n = A.NumRows();
707   long m = A.NumCols();
708
709   X.SetDims(n, m);
710
711   long i, j;
712   for (i = 0; i < n; i++)
713      for (j = 0; j < m; j++)
714         mul(X[i][j], A[i][j], b);
715}
716
717void diag(mat_ZZ_p& X, long n, const ZZ_p& d_in) 
718{ 
719   ZZ_p d = d_in;
720   X.SetDims(n, n); 
721   long i, j; 
722 
723   for (i = 1; i <= n; i++) 
724      for (j = 1; j <= n; j++) 
725         if (i == j) 
726            X(i, j) = d; 
727         else 
728            clear(X(i, j)); 
729} 
730
731long IsDiag(const mat_ZZ_p& A, long n, const ZZ_p& d)
732{
733   if (A.NumRows() != n || A.NumCols() != n)
734      return 0;
735
736   long i, j;
737
738   for (i = 1; i <= n; i++)
739      for (j = 1; j <= n; j++)
740         if (i != j) {
741            if (!IsZero(A(i, j))) return 0;
742         }
743         else {
744            if (A(i, j) != d) return 0;
745         }
746
747   return 1;
748}
749
750
751long IsZero(const mat_ZZ_p& a)
752{
753   long n = a.NumRows();
754   long i;
755
756   for (i = 0; i < n; i++)
757      if (!IsZero(a[i]))
758         return 0;
759
760   return 1;
761}
762
763void clear(mat_ZZ_p& x)
764{
765   long n = x.NumRows();
766   long i;
767   for (i = 0; i < n; i++)
768      clear(x[i]);
769}
770
771
772mat_ZZ_p operator+(const mat_ZZ_p& a, const mat_ZZ_p& b)
773{
774   mat_ZZ_p res;
775   add(res, a, b);
776   NTL_OPT_RETURN(mat_ZZ_p, res);
777}
778
779mat_ZZ_p operator*(const mat_ZZ_p& a, const mat_ZZ_p& b)
780{
781   mat_ZZ_p res;
782   mul_aux(res, a, b);
783   NTL_OPT_RETURN(mat_ZZ_p, res);
784}
785
786mat_ZZ_p operator-(const mat_ZZ_p& a, const mat_ZZ_p& b)
787{
788   mat_ZZ_p res;
789   sub(res, a, b);
790   NTL_OPT_RETURN(mat_ZZ_p, res);
791}
792
793
794mat_ZZ_p operator-(const mat_ZZ_p& a)
795{
796   mat_ZZ_p res;
797   negate(res, a);
798   NTL_OPT_RETURN(mat_ZZ_p, res);
799}
800
801
802vec_ZZ_p operator*(const mat_ZZ_p& a, const vec_ZZ_p& b)
803{
804   vec_ZZ_p res;
805   mul_aux(res, a, b);
806   NTL_OPT_RETURN(vec_ZZ_p, res);
807}
808
809vec_ZZ_p operator*(const vec_ZZ_p& a, const mat_ZZ_p& b)
810{
811   vec_ZZ_p res;
812   mul_aux(res, a, b);
813   NTL_OPT_RETURN(vec_ZZ_p, res);
814}
815
816void inv(mat_ZZ_p& X, const mat_ZZ_p& A)
817{
818   ZZ_p d;
819   inv(d, X, A);
820   if (d == 0) Error("inv: non-invertible matrix");
821}
822
823void power(mat_ZZ_p& X, const mat_ZZ_p& A, const ZZ& e)
824{
825   if (A.NumRows() != A.NumCols()) Error("power: non-square matrix");
826
827   if (e == 0) {
828      ident(X, A.NumRows());
829      return;
830   }
831
832   mat_ZZ_p T1, T2;
833   long i, k;
834
835   k = NumBits(e);
836   T1 = A;
837
838   for (i = k-2; i >= 0; i--) {
839      sqr(T2, T1);
840      if (bit(e, i))
841         mul(T1, T2, A);
842      else
843         T1 = T2;
844   }
845
846   if (e < 0)
847      inv(X, T1);
848   else
849      X = T1;
850}
851
852NTL_END_IMPL
Note: See TracBrowser for help on using the repository browser.