1 | /***************************************** |
---|
2 | * Computer Algebra System SINGULAR * |
---|
3 | *****************************************/ |
---|
4 | /* $Id$ */ |
---|
5 | /* |
---|
6 | * ABSTRACT: Gauss-Manin system normal form |
---|
7 | */ |
---|
8 | |
---|
9 | #include <kernel/mod2.h> |
---|
10 | |
---|
11 | #ifdef HAVE_GMS |
---|
12 | |
---|
13 | #include "gms.h" |
---|
14 | |
---|
15 | #include <coeffs/numbers.h> |
---|
16 | #include <kernel/polys.h> |
---|
17 | #include <kernel/febase.h> |
---|
18 | |
---|
19 | #include "ipid.h" |
---|
20 | |
---|
21 | lists gmsNF(ideal p,ideal g,matrix B,int D,int K) |
---|
22 | { |
---|
23 | ideal r=idInit(IDELEMS(p),1); |
---|
24 | ideal q=idInit(IDELEMS(p),1); |
---|
25 | |
---|
26 | matrix B0=mpNew(MATROWS(B),MATCOLS(B)); |
---|
27 | for(int i=1;i<=MATROWS(B0);i++) |
---|
28 | for(int j=1;j<=MATCOLS(B0);j++) |
---|
29 | if(MATELEM(B,i,j)!=NULL) |
---|
30 | MATELEM(B0,i,j)=pDiff(MATELEM(B,i,j),i+1); |
---|
31 | |
---|
32 | for(int k=0;k<IDELEMS(p);k++) |
---|
33 | { |
---|
34 | while(p->m[k]!=NULL&&pGetExp(p->m[k],1)<=K) |
---|
35 | { |
---|
36 | int j=0; |
---|
37 | while(j<IDELEMS(g)&&!pLmDivisibleBy(g->m[j],p->m[k])) |
---|
38 | j++; |
---|
39 | |
---|
40 | if(j<IDELEMS(g)) |
---|
41 | { |
---|
42 | poly m=pDivideM(pHead(p->m[k]),pHead(g->m[j])); |
---|
43 | p->m[k]=pSub(p->m[k],ppMult_mm(g->m[j],m)); |
---|
44 | pIncrExp(m,1); |
---|
45 | pSetm(m); |
---|
46 | for(int i=0;i<MATROWS(B);i++) |
---|
47 | { |
---|
48 | poly m0=pDiff(m,i+2); |
---|
49 | if(MATELEM(B0,i+1,j+1)!=NULL) |
---|
50 | p->m[k]=pAdd(p->m[k],ppMult_mm(MATELEM(B0,i+1,j+1),m)); |
---|
51 | if(MATELEM(B,i+1,j+1)!=NULL&&m0!=NULL) |
---|
52 | p->m[k]=pAdd(p->m[k],ppMult_mm(MATELEM(B,i+1,j+1),m0)); |
---|
53 | pDelete(&m0); |
---|
54 | } |
---|
55 | pDelete(&m); |
---|
56 | } |
---|
57 | else |
---|
58 | { |
---|
59 | poly p0=p->m[k]; |
---|
60 | pIter(p->m[k]); |
---|
61 | pNext(p0)=NULL; |
---|
62 | r->m[k]=pAdd(r->m[k],p0); |
---|
63 | } |
---|
64 | |
---|
65 | while(p->m[k]!=NULL&&pGetExp(p->m[k],1)<=K&&pWTotaldegree(p->m[k])>D) |
---|
66 | { |
---|
67 | int i=pGetExp(p->m[k],1); |
---|
68 | do |
---|
69 | { |
---|
70 | poly p0=p->m[k]; |
---|
71 | pIter(p->m[k]); |
---|
72 | pNext(p0)=NULL; |
---|
73 | q->m[k]=pAdd(q->m[k],p0); |
---|
74 | }while(p->m[k]!=NULL&&pGetExp(p->m[k],1)==i); |
---|
75 | } |
---|
76 | |
---|
77 | pNormalize(p->m[k]); |
---|
78 | } |
---|
79 | |
---|
80 | q->m[k]=pAdd(q->m[k],p->m[k]); |
---|
81 | p->m[k]=NULL; |
---|
82 | } |
---|
83 | idDelete(&p); |
---|
84 | idDelete((ideal *)&B0); |
---|
85 | |
---|
86 | id_Normalize(r, currRing); |
---|
87 | id_Normalize(q, currRing); |
---|
88 | |
---|
89 | lists l=(lists)omAllocBin(slists_bin); |
---|
90 | l->Init(2); |
---|
91 | |
---|
92 | l->m[0].rtyp=IDEAL_CMD; |
---|
93 | l->m[0].data=r; |
---|
94 | l->m[1].rtyp=IDEAL_CMD; |
---|
95 | l->m[1].data=q; |
---|
96 | |
---|
97 | return l; |
---|
98 | } |
---|
99 | |
---|
100 | |
---|
101 | BOOLEAN gmsNF(leftv res,leftv h) |
---|
102 | { |
---|
103 | if(currRingHdl) |
---|
104 | { |
---|
105 | if(h&&h->Typ()==IDEAL_CMD) |
---|
106 | { |
---|
107 | ideal p=(ideal)h->Data(); |
---|
108 | h=h->next; |
---|
109 | if(h&&h->Typ()==IDEAL_CMD) |
---|
110 | { |
---|
111 | ideal g=(ideal)h->Data(); |
---|
112 | h=h->next; |
---|
113 | if(h&&h->Typ()==MATRIX_CMD) |
---|
114 | { |
---|
115 | matrix B=(matrix)h->Data(); |
---|
116 | h=h->next; |
---|
117 | if(h&&h->Typ()==INT_CMD) |
---|
118 | { |
---|
119 | int D=(int)(long)h->Data(); |
---|
120 | h=h->next; |
---|
121 | if(h&&h->Typ()==INT_CMD) |
---|
122 | { |
---|
123 | int K=(int)(long)h->Data(); |
---|
124 | res->rtyp=LIST_CMD; |
---|
125 | res->data=(void *)gmsNF(idCopy(p),g,B,D,K); |
---|
126 | return FALSE; |
---|
127 | } |
---|
128 | } |
---|
129 | } |
---|
130 | } |
---|
131 | } |
---|
132 | WerrorS("<ideal>,<ideal>,<matrix>,<int>,<int> expected"); |
---|
133 | return TRUE; |
---|
134 | } |
---|
135 | WerrorS("no ring active"); |
---|
136 | return TRUE; |
---|
137 | } |
---|
138 | |
---|
139 | #endif /* HAVE_GMS */ |
---|
140 | |
---|