source: git/ntl/src/mat_ZZ_p.c @ 287cc8

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