My Project
Loading...
Searching...
No Matches
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"
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
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;
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:57
poly mp_DetBareiss(matrix a, const ring r)
returns the determinant of the matrix m; uses Bareiss algorithm
Definition: matpol.cc:1669
#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