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