[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 | |
---|
| 27 | BOOLEAN 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 | |
---|
| 55 | BOOLEAN 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] | 88 | BOOLEAN 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 | |
---|
| 121 | BOOLEAN 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] | 140 | lists 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] | 288 | BOOLEAN 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 */ |
---|