Data Structures | Typedefs | Functions
linalg_from_matpol.cc File Reference

Go to the source code of this file.

Data Structures

class  row_col_weight
 
class  mp_permmatrix
 

Typedefs

typedef int perm[100]
 

Functions

static void mpReplace (int j, int n, int &sign, int *perm)
 perform replacement for pivot strategy in Bareiss algorithm change sign of determinant More...
 
static int mpNextperm (perm *z, int max)
 
static poly mp_Leibnitz (matrix a, const ring)
 
static poly minuscopy (poly p, const ring)
 
static poly p_Insert (poly p1, poly p2, const ring)
 
static void mp_PartClean (matrix, int, int, const ring)
 
static int mp_PrepareRow (matrix, int, int, const ring)
 
static int mp_PreparePiv (matrix, int, int, const ring)
 
static int mp_PivBar (matrix, int, int, const ring)
 
static int mp_PivRow (matrix, int, int, const ring)
 
static float mp_PolyWeight (poly, const ring)
 
static void mp_SwapRow (matrix, int, int, int, const ring)
 
static void mp_SwapCol (matrix, int, int, int, const ring)
 
static void mp_ElimBar (matrix, matrix, poly, int, int, const ring)
 
static void mpSwapRow (matrix a, int pos, int lr, int lc)
 
static void mpSwapCol (matrix a, int pos, int lr, int lc)
 
poly mp_Det (matrix m, const ring R)
 

Typedef Documentation

§ perm

typedef int perm[100]

Definition at line 1 of file linalg_from_matpol.cc.

Function Documentation

§ minuscopy()

static poly minuscopy ( poly  p,
const ring  R 
)
static

Definition at line 232 of file linalg_from_matpol.cc.

233 {
234  poly w;
235  number e;
236  e = n_Init(-1, R);
237  w = p_Copy(p, R);
238  p_Mult_nn(w, e, R);
239  n_Delete(&e, R);
240  return w;
241 }
return P p
Definition: myNF.cc:203
static FORCE_INLINE number n_Init(long i, const coeffs r)
a number representing i in the given coeff field/ring r
Definition: coeffs.h:542
static poly p_Copy(poly p, const ring r)
returns a copy of p
Definition: p_polys.h:804
const ring R
Definition: DebugPrint.cc:36
static poly p_Mult_nn(poly p, number n, const ring r)
Definition: p_polys.h:895
const CanonicalForm & w
Definition: facAbsFact.cc:55
static FORCE_INLINE void n_Delete(number *p, const coeffs r)
delete 'p'
Definition: coeffs.h:459
polyrec * poly
Definition: hilb.h:10

§ mp_Det()

poly mp_Det ( matrix  m,
const ring  R 
)

Definition at line 248 of file linalg_from_matpol.cc.

249 {
250  int i,j,k,n;
251  poly p,q;
252  matrix a, s;
253  matrix ma[100];
254  number c=NULL, d=NULL, ONE=NULL;
255 
256  n = MATROWS(m);
257  if (n != MATCOLS(m))
258  {
259  Werror("det of %d x %d matrix",n,MATCOLS(m));
260  return NULL;
261  }
262  k=rChar(R);
263  if ((k > 0) && (k <= n))
264  return mp_Leibnitz(m, R);
265  ONE = n_Init(1, R);
266  ma[1]=mp_Copy(m, R);
267  k = (n+1) / 2;
268  s = mpNew(1, n);
269  MATELEM(s,1,1) = mp_Trace(m, R);
270  for (i=2; i<=k; i++)
271  {
272  //ma[i] = mpNew(n,n);
273  ma[i]=mp_Mult(ma[i-1], ma[1], R);
274  MATELEM(s,1,i) = mp_Trace(ma[i], R);
275  p_Test(MATELEM(s,1,i), R);
276  }
277  for (i=k+1; i<=n; i++)
278  {
279  MATELEM(s,1,i) = TraceOfProd(ma[i / 2], ma[(i+1) / 2], n, R);
280  p_Test(MATELEM(s,1,i), R);
281  }
282  for (i=1; i<=k; i++)
283  id_Delete((ideal *)&(ma[i]), R);
284 /* the array s contains the traces of the powers of the matrix m,
285 * these are the power sums of the eigenvalues of m */
286  a = mpNew(1,n);
287  MATELEM(a,1,1) = minuscopy(MATELEM(s,1,1), R);
288  for (i=2; i<=n; i++)
289  {
290  p = p_Copy(MATELEM(s,1,i), R);
291  for (j=i-1; j>=1; j--)
292  {
293  q = pp_Mult_qq(MATELEM(s,1,j), MATELEM(a,1,i-j), R);
294  p_Test(q, R);
295  p = p_Add_q(p,q, R);
296  }
297  // c= -1/i
298  d = n_Init(-(int)i, R);
299  c = n_Div(ONE, d, R);
300  n_Delete(&d, R);
301 
302  p_Mult_nn(p, c, R);
303  p_Test(p, R);
304  MATELEM(a,1,i) = p;
305  n_Delete(&c, R);
306  }
307 /* the array a contains the elementary symmetric functions of the
308 * eigenvalues of m */
309  for (i=1; i<=n-1; i++)
310  {
311  //p_Delete(&(MATELEM(a,1,i)), R);
312  p_Delete(&(MATELEM(s,1,i)), R);
313  }
314  p_Delete(&(MATELEM(s,1,n)), R);
315 /* up to a sign, the determinant is the n-th elementary symmetric function */
316  if ((n/2)*2 < n)
317  {
318  d = n_Init(-1, R);
319  p_Mult_nn(MATELEM(a,1,n), d, R);
320  n_Delete(&d, R);
321  }
322  n_Delete(&ONE, R);
323  id_Delete((ideal *)&s, R);
324  poly result=MATELEM(a,1,n);
325  MATELEM(a,1,n)=NULL;
326  id_Delete((ideal *)&a, R);
327  return result;
328 }
const CanonicalForm int s
Definition: facAbsFact.cc:55
const poly a
Definition: syzextra.cc:212
return P p
Definition: myNF.cc:203
poly mp_Trace(matrix a, const ring R)
Definition: matpol.cc:289
static FORCE_INLINE number n_Init(long i, const coeffs r)
a number representing i in the given coeff field/ring r
Definition: coeffs.h:542
int rChar(ring r)
Definition: ring.cc:684
void id_Delete(ideal *h, ring r)
deletes an ideal/module/matrix
int k
Definition: cfEzgcd.cc:93
static poly p_Copy(poly p, const ring r)
returns a copy of p
Definition: p_polys.h:804
int j
Definition: myNF.cc:70
static poly pp_Mult_qq(poly p, poly q, const ring r)
Definition: p_polys.h:1070
const ring R
Definition: DebugPrint.cc:36
static poly minuscopy(poly p, const ring)
int i
Definition: cfEzgcd.cc:123
static poly p_Mult_nn(poly p, number n, const ring r)
Definition: p_polys.h:895
#define p_Test(p, r)
Definition: p_polys.h:160
matrix mpNew(int r, int c)
create a r x c zero-matrix
Definition: matpol.cc:48
static void p_Delete(poly *p, const ring r)
Definition: p_polys.h:843
matrix mp_Mult(matrix a, matrix b, const ring R)
Definition: matpol.cc:224
#define MATCOLS(i)
Definition: matpol.h:28
#define NULL
Definition: omList.c:10
static FORCE_INLINE number n_Div(number a, number b, const coeffs r)
return the quotient of &#39;a&#39; and &#39;b&#39;, i.e., a/b; raises an error if &#39;b&#39; is not invertible in r exceptio...
Definition: coeffs.h:619
matrix mp_Copy(matrix a, const ring r)
copies matrix a (from ring r to r)
Definition: matpol.cc:75
static FORCE_INLINE void n_Delete(number *p, const coeffs r)
delete &#39;p&#39;
Definition: coeffs.h:459
poly TraceOfProd(matrix a, matrix b, int n, const ring R)
Definition: matpol.cc:303
#define MATROWS(i)
Definition: matpol.h:27
polyrec * poly
Definition: hilb.h:10
static poly p_Add_q(poly p, poly q, const ring r)
Definition: p_polys.h:877
static poly mp_Leibnitz(matrix a, const ring)
void Werror(const char *fmt,...)
Definition: reporter.cc:189
return result
Definition: facAbsBiFact.cc:76
#define MATELEM(mat, i, j)
Definition: matpol.h:29

§ mp_ElimBar()

static void mp_ElimBar ( matrix  a0,
matrix  re,
poly  div,
int  lr,
int  lc,
const ring  R 
)
static

Definition at line 918 of file linalg_from_matpol.cc.

919 {
920  int r=lr-1, c=lc-1;
921  poly *b = a0->m, *x = re->m;
922  poly piv, elim, q1, q2, *ap, *a, *q;
923  int i, j;
924 
925  ap = &b[r*a0->ncols];
926  piv = ap[c];
927  for(j=c-1; j>=0; j--)
928  if (ap[j] != NULL) ap[j] = p_Neg(ap[j], R);
929  for(i=r-1; i>=0; i--)
930  {
931  a = &b[i*a0->ncols];
932  q = &x[i*re->ncols];
933  if (a[c] != NULL)
934  {
935  elim = a[c];
936  for (j=c-1; j>=0; j--)
937  {
938  q1 = NULL;
939  if (a[j] != NULL)
940  {
941  q1 = SM_MULT(a[j], piv, div, R);
942  if (ap[j] != NULL)
943  {
944  q2 = SM_MULT(ap[j], elim, div, R);
945  q1 = p_Add_q(q1,q2, R);
946  }
947  }
948  else if (ap[j] != NULL)
949  q1 = SM_MULT(ap[j], elim, div, R);
950  if (q1 != NULL)
951  {
952  if (div)
953  SM_DIV(q1, div, R);
954  q[j] = q1;
955  }
956  }
957  }
958  else
959  {
960  for (j=c-1; j>=0; j--)
961  {
962  if (a[j] != NULL)
963  {
964  q1 = SM_MULT(a[j], piv, div, R);
965  if (div)
966  SM_DIV(q1, div, R);
967  q[j] = q1;
968  }
969  }
970  }
971  }
972 }
const poly a
Definition: syzextra.cc:212
int ncols
Definition: matpol.h:22
CanonicalForm lc(const CanonicalForm &f)
poly * m
Definition: matpol.h:19
const ring r
Definition: syzextra.cc:208
int j
Definition: myNF.cc:70
const ring R
Definition: DebugPrint.cc:36
int i
Definition: cfEzgcd.cc:123
CF_NO_INLINE CanonicalForm div(const CanonicalForm &, const CanonicalForm &)
CF_INLINE CanonicalForm div, mod ( const CanonicalForm & lhs, const CanonicalForm & rhs ) ...
Definition: cf_inline.cc:553
#define NULL
Definition: omList.c:10
#define SM_MULT
Definition: sparsmat.h:23
Variable x
Definition: cfModGcd.cc:4023
#define SM_DIV
Definition: sparsmat.h:24
static poly p_Neg(poly p, const ring r)
Definition: p_polys.h:1013
polyrec * poly
Definition: hilb.h:10
static poly p_Add_q(poly p, poly q, const ring r)
Definition: p_polys.h:877
const poly b
Definition: syzextra.cc:213

§ mp_Leibnitz()

static poly mp_Leibnitz ( matrix  a,
const ring  R 
)
static

Definition at line 885 of file linalg_from_matpol.cc.

886 {
887  int i, e, n;
888  poly p, d;
889  perm z;
890 
891  n = MATROWS(a);
892  memset(&z,0,(n+2)*sizeof(int));
893  p = p_One(R);
894  for (i=1; i <= n; i++)
895  p = p_Mult_q(p, p_Copy(MATELEM(a, i, i), R), R);
896  d = p;
897  for (i=1; i<= n; i++)
898  z[i] = i;
899  z[0]=1;
900  e = 1;
901  if (n!=1)
902  {
903  while (e)
904  {
905  e = mpNextperm((perm *)&z, n);
906  p = p_One(R);
907  for (i = 1; i <= n; i++)
908  p = p_Mult_q(p, p_Copy(MATELEM(a, i, z[i]), R), R);
909  if (z[0] > 0)
910  d = p_Add_q(d, p, R);
911  else
912  d = p_Sub(d, p, R);
913  }
914  }
915  return d;
916 }
static int mpNextperm(perm *z, int max)
return P p
Definition: myNF.cc:203
poly p_Sub(poly p1, poly p2, const ring r)
Definition: p_polys.cc:1912
static poly p_Copy(poly p, const ring r)
returns a copy of p
Definition: p_polys.h:804
poly p_One(const ring r)
Definition: p_polys.cc:1313
const ring R
Definition: DebugPrint.cc:36
int i
Definition: cfEzgcd.cc:123
#define MATROWS(i)
Definition: matpol.h:27
polyrec * poly
Definition: hilb.h:10
static poly p_Add_q(poly p, poly q, const ring r)
Definition: p_polys.h:877
int perm[100]
static poly p_Mult_q(poly p, poly q, const ring r)
Definition: p_polys.h:1020
#define MATELEM(mat, i, j)
Definition: matpol.h:29

§ mp_PartClean()

static void mp_PartClean ( matrix  ,
int  ,
int  ,
const ring   
)
static

§ mp_PivBar()

static int mp_PivBar ( matrix  a,
int  lr,
int  lc,
const ring  R 
)
static

Definition at line 49 of file linalg_from_matpol.cc.

50 {
51  float f1, f2;
52  poly *q1;
53  int i,j,io;
54 
55  io = -1;
56  f1 = 1.0e30;
57  for (i=lr-1;i>=0;i--)
58  {
59  q1 = &(a->m)[i*a->ncols];
60  f2 = 0.0;
61  for (j=lc-1;j>=0;j--)
62  {
63  if (q1[j]!=NULL)
64  f2 += mp_PolyWeight(q1[j], R);
65  }
66  if ((f2!=0.0) && (f2<f1))
67  {
68  f1 = f2;
69  io = i;
70  }
71  }
72  if (io<0) return 0;
73  else return io+1;
74 }
int ncols
Definition: matpol.h:22
CanonicalForm lc(const CanonicalForm &f)
poly * m
Definition: matpol.h:19
for(int i=0;i< R->ExpL_Size;i++) Print("%09lx "
Definition: cfEzgcd.cc:66
int j
Definition: myNF.cc:70
const ring R
Definition: DebugPrint.cc:36
int i
Definition: cfEzgcd.cc:123
static float mp_PolyWeight(poly, const ring)
#define NULL
Definition: omList.c:10
polyrec * poly
Definition: hilb.h:10

§ mp_PivRow()

static int mp_PivRow ( matrix  a,
int  lr,
int  lc,
const ring  R 
)
static

Definition at line 79 of file linalg_from_matpol.cc.

80 {
81  float f1, f2;
82  poly *q1;
83  int j,jo;
84 
85  jo = -1;
86  f1 = 1.0e30;
87  q1 = &(a->m)[(lr-1)*a->ncols];
88  for (j=lc-1;j>=0;j--)
89  {
90  if (q1[j]!=NULL)
91  {
92  f2 = mp_PolyWeight(q1[j], R);
93  if (f2<f1)
94  {
95  f1 = f2;
96  jo = j;
97  }
98  }
99  }
100  if (jo<0) return 0;
101  else return jo+1;
102 }
int ncols
Definition: matpol.h:22
CanonicalForm lc(const CanonicalForm &f)
poly * m
Definition: matpol.h:19
int j
Definition: myNF.cc:70
const ring R
Definition: DebugPrint.cc:36
static float mp_PolyWeight(poly, const ring)
#define NULL
Definition: omList.c:10
polyrec * poly
Definition: hilb.h:10

§ mp_PolyWeight()

static float mp_PolyWeight ( poly  p,
const ring  R 
)
static

Definition at line 107 of file linalg_from_matpol.cc.

108 {
109  int i;
110  float res;
111 
112  if (pNext(p) == NULL)
113  {
114  res = (float)n_Size(p_GetCoeff(p, R), R);
115  for (i=rVar(R);i>0;i--)
116  {
117  if(p_GetExp(p,i, R)!=0)
118  {
119  res += 2.0;
120  break;
121  }
122  }
123  }
124  else
125  {
126  res = 0.0;
127  do
128  {
129  res += (float)n_Size(p_GetCoeff(p, R), R) + 2.0;
130  pIter(p);
131  }
132  while (p);
133  }
134  return res;
135 }
return P p
Definition: myNF.cc:203
static short rVar(const ring r)
#define rVar(r) (r->N)
Definition: ring.h:580
#define pIter(p)
Definition: monomials.h:44
poly res
Definition: myNF.cc:322
static long p_GetExp(const poly p, const unsigned long iBitmask, const int VarOffset)
get a single variable exponent : the integer VarOffset encodes:
Definition: p_polys.h:464
const ring R
Definition: DebugPrint.cc:36
int i
Definition: cfEzgcd.cc:123
#define NULL
Definition: omList.c:10
#define pNext(p)
Definition: monomials.h:43
#define p_GetCoeff(p, r)
Definition: monomials.h:57
static FORCE_INLINE int n_Size(number n, const coeffs r)
return a non-negative measure for the complexity of n; return 0 only when n represents zero; (used fo...
Definition: coeffs.h:574

§ mp_PreparePiv()

static int mp_PreparePiv ( matrix  a,
int  lr,
int  lc,
const ring  R 
)
static

Definition at line 36 of file linalg_from_matpol.cc.

37 {
38  int c;
39 
40  c = mp_PivRow(a, lr, lc, R);
41  if(c==0) return 0;
42  if(c<lc) mp_SwapCol(a, c, lr, lc, R);
43  return 1;
44 }
static void mp_SwapCol(matrix, int, int, int, const ring)
static int mp_PivRow(matrix, int, int, const ring)
CanonicalForm lc(const CanonicalForm &f)
const ring R
Definition: DebugPrint.cc:36

§ mp_PrepareRow()

static int mp_PrepareRow ( matrix  a,
int  lr,
int  lc,
const ring  R 
)
static

Definition at line 22 of file linalg_from_matpol.cc.

23 {
24  int r;
25 
26  r = mp_PivBar(a,lr,lc, R);
27  if(r==0) return 0;
28  if(r<lr) mp_SwapRow(a, r, lr, lc, R);
29  return 1;
30 }
static void mp_SwapRow(matrix, int, int, int, const ring)
CanonicalForm lc(const CanonicalForm &f)
const ring r
Definition: syzextra.cc:208
const ring R
Definition: DebugPrint.cc:36
static int mp_PivBar(matrix, int, int, const ring)

§ mp_SwapCol()

static void mp_SwapCol ( matrix  ,
int  ,
int  ,
int  ,
const ring   
)
static

§ mp_SwapRow()

static void mp_SwapRow ( matrix  ,
int  ,
int  ,
int  ,
const ring   
)
static

§ mpNextperm()

static int mpNextperm ( perm z,
int  max 
)
static

Definition at line 833 of file linalg_from_matpol.cc.

834 {
835  int s, i, k, t;
836  s = max;
837  do
838  {
839  s--;
840  }
841  while ((s > 0) && ((*z)[s] >= (*z)[s+1]));
842  if (s==0)
843  return 0;
844  do
845  {
846  (*z)[s]++;
847  k = 0;
848  do
849  {
850  k++;
851  }
852  while (((*z)[k] != (*z)[s]) && (k!=s));
853  }
854  while (k < s);
855  for (i=s+1; i <= max; i++)
856  {
857  (*z)[i]=0;
858  do
859  {
860  (*z)[i]++;
861  k=0;
862  do
863  {
864  k++;
865  }
866  while (((*z)[k] != (*z)[i]) && (k != i));
867  }
868  while (k < i);
869  }
870  s = max+1;
871  do
872  {
873  s--;
874  }
875  while ((s > 0) && ((*z)[s] > (*z)[s+1]));
876  t = 1;
877  for (i=1; i<max; i++)
878  for (k=i+1; k<=max; k++)
879  if ((*z)[k] < (*z)[i])
880  t = -t;
881  (*z)[0] = t;
882  return s;
883 }
const CanonicalForm int s
Definition: facAbsFact.cc:55
int k
Definition: cfEzgcd.cc:93
static int max(int a, int b)
Definition: fast_mult.cc:264
int i
Definition: cfEzgcd.cc:123

§ mpReplace()

static void mpReplace ( int  j,
int  n,
int &  sign,
int *  perm 
)
static

perform replacement for pivot strategy in Bareiss algorithm change sign of determinant

Definition at line 820 of file linalg_from_matpol.cc.

821 {
822  int k;
823 
824  if (j != n)
825  {
826  k = perm[n];
827  perm[n] = perm[j];
828  perm[j] = k;
829  sign = -sign;
830  }
831 }
int k
Definition: cfEzgcd.cc:93
int j
Definition: myNF.cc:70
int perm[100]
static int sign(int x)
Definition: ring.cc:3412

§ mpSwapCol()

static void mpSwapCol ( matrix  a,
int  pos,
int  lr,
int  lc 
)
static

Definition at line 153 of file linalg_from_matpol.cc.

154 {
155  poly sw;
156  int j;
157  poly* a2 = a->m;
158  poly* a1 = &a2[pos-1];
159 
160  a2 = &a2[lc-1];
161  for (j=a->ncols*(lr-1); j>=0; j-=a->ncols)
162  {
163  sw = a1[j];
164  a1[j] = a2[j];
165  a2[j] = sw;
166  }
167 }
int ncols
Definition: matpol.h:22
CanonicalForm lc(const CanonicalForm &f)
poly * m
Definition: matpol.h:19
int j
Definition: myNF.cc:70
polyrec * poly
Definition: hilb.h:10

§ mpSwapRow()

static void mpSwapRow ( matrix  a,
int  pos,
int  lr,
int  lc 
)
static

Definition at line 137 of file linalg_from_matpol.cc.

138 {
139  poly sw;
140  int j;
141  poly* a2 = a->m;
142  poly* a1 = &a2[a->ncols*(pos-1)];
143 
144  a2 = &a2[a->ncols*(lr-1)];
145  for (j=lc-1; j>=0; j--)
146  {
147  sw = a1[j];
148  a1[j] = a2[j];
149  a2[j] = sw;
150  }
151 }
int ncols
Definition: matpol.h:22
CanonicalForm lc(const CanonicalForm &f)
poly * m
Definition: matpol.h:19
int j
Definition: myNF.cc:70
polyrec * poly
Definition: hilb.h:10

§ p_Insert()

static poly p_Insert ( poly  p1,
poly  p2,
const ring   
)
static