source: git/Singular/eigenval_ip.cc @ 8c1285

spielwiese
Last change on this file since 8c1285 was a9c298, checked in by Hans Schoenemann <hannes@…>, 10 years ago
format stuff
  • Property mode set to 100644
File size: 6.2 KB
Line 
1/*****************************************
2*  Computer Algebra System SINGULAR      *
3*****************************************/
4/*
5* ABSTRACT: eigenvalues of constant square matrices
6*/
7
8#ifdef HAVE_CONFIG_H
9#include "singularconfig.h"
10#endif /* HAVE_CONFIG_H */
11#include <kernel/mod2.h>
12
13#ifdef HAVE_EIGENVAL
14
15#include <kernel/febase.h>
16#include <Singular/tok.h>
17#include <Singular/ipid.h>
18#include <misc/intvec.h>
19#include <coeffs/numbers.h>
20#include <kernel/polys.h>
21#include <kernel/ideals.h>
22#include <Singular/lists.h>
23#include <polys/matpol.h>
24#include <polys/clapsing.h>
25#include <kernel/eigenval.h>
26#include <Singular/eigenval_ip.h>
27
28
29BOOLEAN evSwap(leftv res,leftv h)
30{
31  if(currRing)
32  {
33    if(h&&h->Typ()==MATRIX_CMD)
34    {
35      matrix M=(matrix)h->Data();
36      h=h->next;
37      if(h&&h->Typ()==INT_CMD)
38      {
39        int i=(int)(long)h->Data();
40        h=h->next;
41        if(h&&h->Typ()==INT_CMD)
42        {
43          int j=(int)(long)h->Data();
44          res->rtyp=MATRIX_CMD;
45          res->data=(void *)evSwap(mp_Copy(M, currRing),i,j);
46          return FALSE;
47        }
48      }
49    }
50    WerrorS("<matrix>,<int>,<int> expected");
51    return TRUE;
52  }
53  WerrorS("no ring active");
54  return TRUE;
55}
56
57BOOLEAN evRowElim(leftv res,leftv h)
58{
59  if(currRing)
60  {
61    if(h&&h->Typ()==MATRIX_CMD)
62    {
63      matrix M=(matrix)h->Data();
64      h=h->next;
65      if(h&&h->Typ()==INT_CMD)
66      {
67        int i=(int)(long)h->Data();
68        h=h->next;
69        if(h&&h->Typ()==INT_CMD)
70        {
71          int j=(int)(long)h->Data();
72          h=h->next;
73          if(h&&h->Typ()==INT_CMD)
74          {
75            int k=(int)(long)h->Data();
76            res->rtyp=MATRIX_CMD;
77            res->data=(void *)evRowElim(mp_Copy(M, currRing),i,j,k);
78            return FALSE;
79          }
80        }
81      }
82    }
83    WerrorS("<matrix>,<int>,<int>,<int> expected");
84    return TRUE;
85  }
86  WerrorS("no ring active");
87  return TRUE;
88}
89
90BOOLEAN evColElim(leftv res,leftv h)
91{
92  if(currRing)
93  {
94    if(h&&h->Typ()==MATRIX_CMD)
95    {
96      matrix M=(matrix)h->Data();
97      h=h->next;
98      if(h&&h->Typ()==INT_CMD)
99      {
100        int i=(int)(long)h->Data();
101        h=h->next;
102        if(h&&h->Typ()==INT_CMD)
103        {
104          int j=(int)(long)h->Data();
105          h=h->next;
106          if(h&&h->Typ()==INT_CMD)
107          {
108            int k=(int)(long)h->Data();
109            res->rtyp=MATRIX_CMD;
110            res->data=(void *)evColElim(mp_Copy(M, currRing),i,j,k);
111            return FALSE;
112          }
113        }
114      }
115    }
116    WerrorS("<matrix>,<int>,<int>,<int> expected");
117    return TRUE;
118  }
119  WerrorS("no ring active");
120  return TRUE;
121}
122
123BOOLEAN evHessenberg(leftv res,leftv h)
124{
125  if(currRing)
126  {
127    if(h&&h->Typ()==MATRIX_CMD)
128    {
129      matrix M=(matrix)h->Data();
130      res->rtyp=MATRIX_CMD;
131      res->data=(void *)evHessenberg(mp_Copy(M, currRing));
132      return FALSE;
133    }
134    WerrorS("<matrix> expected");
135    return TRUE;
136  }
137  WerrorS("no ring active");
138  return TRUE;
139}
140
141
142lists evEigenvals(matrix M)
143{
144  lists l=(lists)omAllocBin(slists_bin);
145  if(MATROWS(M)!=MATCOLS(M))
146  {
147    l->Init(0);
148    return(l);
149  }
150
151  M=evHessenberg((matrix)id_Jet((ideal)M,0,currRing));
152
153  int n=MATROWS(M);
154  ideal e=idInit(n,1);
155  intvec *m=new intvec(n);
156
157  poly t=pOne();
158  pSetExp(t,1,1);
159  pSetm(t);
160
161  for(int j0=1,j=2,k=0;j<=n+1;j0=j,j++)
162  {
163    while(j<=n&&MATELEM(M,j,j-1)!=NULL)
164      j++;
165    if(j==j0+1)
166    {
167      e->m[k]=pHead(MATELEM(M,j0,j0));
168      (*m)[k]=1;
169      k++;
170    }
171    else
172    {
173      int n0=j-j0;
174      matrix M0=mpNew(n0,n0);
175
176      j0--;
177      for(int i=1;i<=n0;i++)
178        for(int j=1;j<=n0;j++)
179          MATELEM(M0,i,j)=pCopy(MATELEM(M,j0+i,j0+j));
180      for(int i=1;i<=n0;i++)
181        MATELEM(M0,i,i)=pSub(MATELEM(M0,i,i),pCopy(t));
182
183      intvec *m0;
184      ideal e0=singclap_factorize(mp_DetBareiss(M,currRing),&m0,2, currRing);
185      if (e0==NULL)
186      {
187        l->Init(0);
188        return(l);
189      }
190
191      for(int i=0;i<IDELEMS(e0);i++)
192      {
193        if(pNext(e0->m[i])==NULL)
194        {
195          (*m)[k]=(*m0)[i];
196          k++;
197        }
198        else
199        if(pGetExp(e0->m[i],1)<2&&pGetExp(pNext(e0->m[i]),1)<2&&
200           pNext(pNext(e0->m[i]))==NULL)
201        {
202          number e1=nNeg(nCopy(pGetCoeff(e0->m[i])));
203          if(pGetExp(pNext(e0->m[i]),1)==0)
204            e->m[k]=pNSet(nDiv(pGetCoeff(pNext(e0->m[i])),e1));
205          else
206            e->m[k]=pNSet(nDiv(e1,pGetCoeff(pNext(e0->m[i]))));
207          nDelete(&e1);
208          pNormalize(e->m[k]);
209          (*m)[k]=(*m0)[i];
210          k++;
211        }
212        else
213        {
214          e->m[k]=e0->m[i];
215          pNormalize(e->m[k]);
216          e0->m[i]=NULL;
217          (*m)[k]=(*m0)[i];
218          k++;
219        }
220      }
221
222      delete(m0);
223      idDelete(&e0);
224    }
225  }
226
227  pDelete(&t);
228  idDelete((ideal *)&M);
229
230  for(int i0=0;i0<n-1;i0++)
231  {
232    for(int i1=i0+1;i1<n;i1++)
233    {
234      if(pEqualPolys(e->m[i0],e->m[i1]))
235      {
236        (*m)[i0]+=(*m)[i1];
237        (*m)[i1]=0;
238      }
239      else
240      {
241        if(e->m[i0]==NULL&&!nGreaterZero(pGetCoeff(e->m[i1]))||
242           e->m[i1]==NULL&&
243          (nGreaterZero(pGetCoeff(e->m[i0]))||pNext(e->m[i0])!=NULL)||
244           e->m[i0]!=NULL&&e->m[i1]!=NULL&&
245          (pNext(e->m[i0])!=NULL&&pNext(e->m[i1])==NULL||
246           pNext(e->m[i0])==NULL&&pNext(e->m[i1])==NULL&&
247           nGreater(pGetCoeff(e->m[i0]),pGetCoeff(e->m[i1]))))
248        {
249          poly e1=e->m[i0];
250          e->m[i0]=e->m[i1];
251          e->m[i1]=e1;
252          int m1=(*m)[i0];
253          (*m)[i0]=(*m)[i1];
254          (*m)[i1]=m1;
255        }
256      }
257    }
258  }
259
260  int n0=0;
261  for(int i=0;i<n;i++)
262    if((*m)[i]>0)
263      n0++;
264
265  ideal e0=idInit(n0,1);
266  intvec *m0=new intvec(n0);
267
268  for(int i=0,i0=0;i<n;i++)
269    if((*m)[i]>0)
270    {
271      e0->m[i0]=e->m[i];
272      e->m[i]=NULL;
273      (*m0)[i0]=(*m)[i];
274      i0++;
275    }
276
277  idDelete(&e);
278  delete(m);
279
280  l->Init(2);
281  l->m[0].rtyp=IDEAL_CMD;
282  l->m[0].data=e0;
283  l->m[1].rtyp=INTVEC_CMD;
284  l->m[1].data=m0;
285
286  return(l);
287}
288
289
290BOOLEAN evEigenvals(leftv res,leftv h)
291{
292  if(currRing)
293  {
294    if(h&&h->Typ()==MATRIX_CMD)
295    {
296      matrix M=(matrix)h->Data();
297      res->rtyp=LIST_CMD;
298      res->data=(void *)evEigenvals(mp_Copy(M, currRing));
299      return FALSE;
300    }
301    WerrorS("<matrix> expected");
302    return TRUE;
303  }
304  WerrorS("no ring active");
305  return TRUE;
306}
307
308#endif /* HAVE_EIGENVAL */
Note: See TracBrowser for help on using the repository browser.