source: git/ntl/src/mat_poly_lzz_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: 1.5 KB
Line 
1#include <NTL/mat_poly_lzz_p.h>
2
3#include <NTL/new.h>
4
5NTL_START_IMPL
6
7
8void CharPoly(zz_pX& f, const mat_zz_p& M)
9{
10   long n = M.NumRows();
11   if (M.NumCols() != n)
12      Error("CharPoly: nonsquare matrix");
13
14   if (n == 0) {
15      set(f);
16      return;
17   }
18
19   zz_p t;
20
21   if (n == 1) {
22      SetX(f);
23      negate(t, M(1, 1));
24      SetCoeff(f, 0, t);
25      return;
26   }
27
28   mat_zz_p H;
29
30   H = M;
31
32   long i, j, m;
33   zz_p u, t1;
34
35   for (m = 2; m <= n-1; m++) {
36      i = m;
37      while (i <= n && IsZero(H(i, m-1)))
38         i++;
39
40      if (i <= n) {
41         t = H(i, m-1);
42         if (i > m) {
43            swap(H(i), H(m));
44            // swap columns i and m
45            for (j = 1; j <= n; j++) 
46               swap(H(j, i), H(j, m));
47         }
48
49         for (i = m+1; i <= n; i++) {
50            div(u, H(i, m-1), t);
51            for (j = m; j <= n; j++) {
52               mul(t1, u, H(m, j));
53               sub(H(i, j), H(i, j), t1);
54            }
55
56            for (j = 1; j <= n; j++) {
57               mul(t1, u, H(j, i));
58               add(H(j, m), H(j, m), t1);
59            }
60         }
61      }
62   }
63
64   vec_zz_pX F;
65   F.SetLength(n+1);
66   zz_pX T;
67   T.SetMaxLength(n);
68
69   set(F[0]);
70   for (m = 1; m <= n; m++) {
71      LeftShift(F[m], F[m-1], 1);
72      mul(T, F[m-1], H(m, m));
73      sub(F[m], F[m], T);
74
75      set(t);
76      for (i = 1; i <= m-1; i++) {
77         mul(t, t, H(m-i+1, m-i));
78         mul(t1, t, H(m-i, m));
79         mul(T, F[m-i-1], t1);
80         sub(F[m], F[m], T);
81      }
82   }
83
84   f = F[n];
85}
86
87NTL_END_IMPL
Note: See TracBrowser for help on using the repository browser.