My Project
eigenval_ip.cc
Go to the documentation of this file.
1 /*****************************************
2 * Computer Algebra System SINGULAR *
3 *****************************************/
4 /*
5 * ABSTRACT: eigenvalues of constant square matrices
6 */
7 
8 
9 
10 
11 #include "kernel/mod2.h"
12 
13 #ifdef HAVE_EIGENVAL
14 
15 #include "Singular/tok.h"
16 #include "Singular/ipid.h"
17 #include "misc/intvec.h"
18 #include "coeffs/numbers.h"
19 #include "kernel/polys.h"
20 #include "kernel/ideals.h"
21 #include "Singular/lists.h"
22 #include "polys/matpol.h"
23 #include "polys/clapsing.h"
25 #include "Singular/ipshell.h"
26 #include "Singular/eigenval_ip.h"
27 
28 
30 {
31  if(currRing)
32  {
33  const short t[]={3,MATRIX_CMD,INT_CMD,INT_CMD};
34  if (iiCheckTypes(h,t,1))
35  {
36  matrix M=(matrix)h->Data();
37  h=h->next;
38  int i=(int)(long)h->Data();
39  h=h->next;
40  int j=(int)(long)h->Data();
41  res->rtyp=MATRIX_CMD;
42  res->data=(void *)evSwap(mp_Copy(M, currRing),i,j);
43  return FALSE;
44  }
45  return TRUE;
46  }
47  WerrorS("no ring active");
48  return TRUE;
49 }
50 
52 {
53  if(currRing)
54  {
55  const short t[]={4,MATRIX_CMD,INT_CMD,INT_CMD,INT_CMD};
56  if (iiCheckTypes(h,t,1))
57  {
58  matrix M=(matrix)h->CopyD();
59  h=h->next;
60  int i=(int)(long)h->Data();
61  h=h->next;
62  int j=(int)(long)h->Data();
63  h=h->next;
64  int k=(int)(long)h->Data();
65  res->rtyp=MATRIX_CMD;
66  res->data=(void *)evRowElim(M,i,j,k);
67  return FALSE;
68  }
69  return TRUE;
70  }
71  WerrorS("no ring active");
72  return TRUE;
73 }
74 
75 #if 0
77 {
78  if(currRing)
79  {
80  const short t[]={4,MATRIX_CMD,INT_CMD,INT_CMD,INT_CMD};
81  if (iiCheckTypes(h,t,1))
82  {
83  matrix M=(matrix)h->Data();
84  h=h->next;
85  int i=(int)(long)h->Data();
86  h=h->next;
87  int j=(int)(long)h->Data();
88  h=h->next;
89  int k=(int)(long)h->Data();
90  res->rtyp=MATRIX_CMD;
91  res->data=(void *)evColElim(mp_Copy(M, currRing),i,j,k);
92  return FALSE;
93  }
94  return TRUE;
95  }
96  WerrorS("no ring active");
97  return TRUE;
98 }
99 #endif
100 
102 {
103  if(currRing)
104  {
105  if(h&&h->Typ()==MATRIX_CMD)
106  {
107  matrix M=(matrix)h->Data();
108  res->rtyp=MATRIX_CMD;
109  res->data=(void *)evHessenberg(mp_Copy(M, currRing));
110  return FALSE;
111  }
112  WerrorS("<matrix> expected");
113  return TRUE;
114  }
115  WerrorS("no ring active");
116  return TRUE;
117 }
118 
119 
121 {
123  if(MATROWS(M)!=MATCOLS(M))
124  {
125  l->Init(0);
126  return(l);
127  }
128 
129  M=evHessenberg(M);
130 
131  int n=MATCOLS(M);
132  ideal e=idInit(n,1);
133  intvec *m=new intvec(n);
134 
135  poly t=pOne();
136  pSetExp(t,1,1);
137  pSetm(t);
138 
139  for(int j0=1,j=2,k=0;j<=n+1;j0=j,j++)
140  {
141  while(j<=n&&MATELEM(M,j,j-1)!=NULL)
142  j++;
143  if(j==j0+1)
144  {
145  e->m[k]=pHead(MATELEM(M,j0,j0));
146  (*m)[k]=1;
147  k++;
148  }
149  else
150  {
151  int n0=j-j0;
152  matrix M0=mpNew(n0,n0);
153 
154  j0--;
155  for(int i=1;i<=n0;i++)
156  for(int j=1;j<=n0;j++)
157  MATELEM(M0,i,j)=pCopy(MATELEM(M,j0+i,j0+j));
158  for(int i=1;i<=n0;i++)
159  MATELEM(M0,i,i)=pSub(MATELEM(M0,i,i),pCopy(t));
160 
161  intvec *m0;
162  ideal e0=singclap_factorize(mp_DetBareiss(M0,currRing),&m0,2, currRing);
163  if (e0==NULL)
164  {
165  l->Init(0);
166  return(l);
167  }
168 
169  for(int i=0;i<IDELEMS(e0);i++)
170  {
171  if(pNext(e0->m[i])==NULL)
172  {
173  (*m)[k]=(*m0)[i];
174  k++;
175  }
176  else
177  if(pGetExp(e0->m[i],1)<2&&pGetExp(pNext(e0->m[i]),1)<2&&
178  pNext(pNext(e0->m[i]))==NULL)
179  {
180  number e1=nCopy(pGetCoeff(e0->m[i]));
181  e1=nInpNeg(e1);
182  if(pGetExp(pNext(e0->m[i]),1)==0)
183  e->m[k]=pNSet(nDiv(pGetCoeff(pNext(e0->m[i])),e1));
184  else
185  e->m[k]=pNSet(nDiv(e1,pGetCoeff(pNext(e0->m[i]))));
186  nDelete(&e1);
187  pNormalize(e->m[k]);
188  (*m)[k]=(*m0)[i];
189  k++;
190  }
191  else
192  {
193  e->m[k]=e0->m[i];
194  pNormalize(e->m[k]);
195  e0->m[i]=NULL;
196  (*m)[k]=(*m0)[i];
197  k++;
198  }
199  }
200 
201  delete(m0);
202  idDelete(&e0);
203  }
204  }
205 
206  pDelete(&t);
207  idDelete((ideal *)&M);
208 
209  for(int i0=0;i0<n-1;i0++)
210  {
211  for(int i1=i0+1;i1<n;i1++)
212  {
213  if(pEqualPolys(e->m[i0],e->m[i1]))
214  {
215  (*m)[i0]+=(*m)[i1];
216  (*m)[i1]=0;
217  }
218  else
219  {
220  if((e->m[i0]==NULL&&!nGreaterZero(pGetCoeff(e->m[i1])))||
221  (e->m[i1]==NULL&&
222  (nGreaterZero(pGetCoeff(e->m[i0]))||pNext(e->m[i0])!=NULL))||
223  e->m[i0]!=NULL&&e->m[i1]!=NULL&&
224  ((pNext(e->m[i0])!=NULL&&pNext(e->m[i1])==NULL)||
225  (pNext(e->m[i0])==NULL&&pNext(e->m[i1])==NULL&&
226  nGreater(pGetCoeff(e->m[i0]),pGetCoeff(e->m[i1])))))
227  {
228  poly e1=e->m[i0];
229  e->m[i0]=e->m[i1];
230  e->m[i1]=e1;
231  int m1=(*m)[i0];
232  (*m)[i0]=(*m)[i1];
233  (*m)[i1]=m1;
234  }
235  }
236  }
237  }
238 
239  int n0=0;
240  for(int i=0;i<n;i++)
241  if((*m)[i]>0)
242  n0++;
243 
244  ideal e0=idInit(n0,1);
245  intvec *m0=new intvec(n0);
246 
247  for(int i=0,i0=0;i<n;i++)
248  if((*m)[i]>0)
249  {
250  e0->m[i0]=e->m[i];
251  e->m[i]=NULL;
252  (*m0)[i0]=(*m)[i];
253  i0++;
254  }
255 
256  idDelete(&e);
257  delete(m);
258 
259  l->Init(2);
260  l->m[0].rtyp=IDEAL_CMD;
261  l->m[0].data=e0;
262  l->m[1].rtyp=INTVEC_CMD;
263  l->m[1].data=m0;
264 
265  return(l);
266 }
267 
268 
270 {
271  if(currRing)
272  {
273  if(h&&h->Typ()==MATRIX_CMD)
274  {
275  matrix M=(matrix)h->CopyD();
276  res->rtyp=LIST_CMD;
277  res->data=(void *)evEigenvals(M);
278  return FALSE;
279  }
280  WerrorS("<matrix> expected");
281  return TRUE;
282  }
283  WerrorS("no ring active");
284  return TRUE;
285 }
286 
287 #endif /* HAVE_EIGENVAL */
int BOOLEAN
Definition: auxiliary.h:87
#define TRUE
Definition: auxiliary.h:100
#define FALSE
Definition: auxiliary.h:96
int l
Definition: cfEzgcd.cc:100
int m
Definition: cfEzgcd.cc:128
int i
Definition: cfEzgcd.cc:132
int k
Definition: cfEzgcd.cc:99
ideal singclap_factorize(poly f, intvec **v, int with_exps, const ring r)
Definition: clapsing.cc:948
Definition: intvec.h:23
Class used for (list of) interpreter objects.
Definition: subexpr.h:83
Definition: lists.h:24
matrix evColElim(matrix M, int i, int j, int k)
Definition: eigenval.cc:76
lists evEigenvals(matrix M)
Definition: eigenval_ip.cc:120
BOOLEAN evSwap(leftv res, leftv h)
Definition: eigenval_ip.cc:29
BOOLEAN evHessenberg(leftv res, leftv h)
Definition: eigenval_ip.cc:101
BOOLEAN evRowElim(leftv res, leftv h)
Definition: eigenval_ip.cc:51
CanonicalForm res
Definition: facAbsFact.cc:60
int j
Definition: facHensel.cc:110
void WerrorS(const char *s)
Definition: feFopen.cc:24
@ IDEAL_CMD
Definition: grammar.cc:284
@ MATRIX_CMD
Definition: grammar.cc:286
#define idDelete(H)
delete an ideal
Definition: ideals.h:29
BOOLEAN iiCheckTypes(leftv args, const short *type_list, int report)
check a list of arguemys against a given field of types return TRUE if the types match return FALSE (...
Definition: ipshell.cc:6572
STATIC_VAR Poly * h
Definition: janet.cc:971
VAR omBin slists_bin
Definition: lists.cc:23
matrix mpNew(int r, int c)
create a r x c zero-matrix
Definition: matpol.cc:37
matrix mp_Copy(matrix a, const ring r)
copies matrix a (from ring r to r)
Definition: matpol.cc:64
poly mp_DetBareiss(matrix a, const ring r)
returns the determinant of the matrix m; uses Bareiss algorithm
Definition: matpol.cc:1676
#define MATELEM(mat, i, j)
1-based access to matrix
Definition: matpol.h:29
ip_smatrix * matrix
Definition: matpol.h:43
#define MATROWS(i)
Definition: matpol.h:26
#define MATCOLS(i)
Definition: matpol.h:27
#define pNext(p)
Definition: monomials.h:36
static number & pGetCoeff(poly p)
return an alias to the leading coefficient of p assumes that p != NULL NOTE: not copy
Definition: monomials.h:44
slists * lists
Definition: mpr_numeric.h:146
#define nDiv(a, b)
Definition: numbers.h:32
#define nDelete(n)
Definition: numbers.h:16
#define nInpNeg(n)
Definition: numbers.h:21
#define nCopy(n)
Definition: numbers.h:15
#define nGreater(a, b)
Definition: numbers.h:28
#define nGreaterZero(n)
Definition: numbers.h:27
#define omAllocBin(bin)
Definition: omAllocDecl.h:205
#define NULL
Definition: omList.c:12
VAR ring currRing
Widely used global variable which specifies the current polynomial ring for Singular interpreter and ...
Definition: polys.cc:13
Compatibility layer for legacy polynomial operations (over currRing)
#define pDelete(p_ptr)
Definition: polys.h:186
#define pHead(p)
returns newly allocated copy of Lm(p), coef is copied, next=NULL, p might be NULL
Definition: polys.h:67
#define pSetm(p)
Definition: polys.h:271
#define pNSet(n)
Definition: polys.h:313
#define pSub(a, b)
Definition: polys.h:287
#define pGetExp(p, i)
Exponent.
Definition: polys.h:41
#define pNormalize(p)
Definition: polys.h:317
#define pEqualPolys(p1, p2)
Definition: polys.h:399
#define pSetExp(p, i, v)
Definition: polys.h:42
#define pCopy(p)
return a copy of the poly
Definition: polys.h:185
#define pOne()
Definition: polys.h:315
ideal idInit(int idsize, int rank)
initialise an ideal / module
Definition: simpleideals.cc:35
#define IDELEMS(i)
Definition: simpleideals.h:23
#define M
Definition: sirandom.c:25
@ LIST_CMD
Definition: tok.h:118
@ INTVEC_CMD
Definition: tok.h:101
@ INT_CMD
Definition: tok.h:96