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