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