source: git/Singular/eigenval_ip.cc @ 666c90

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