source: git/Singular/gms.cc @ 4676d5

spielwiese
Last change on this file since 4676d5 was 4676d5, checked in by Oleksandr Motsak <motsak@…>, 10 years ago
Moved timer.*, feread.cc and febase.* from /kernel/ to /Singular/ (also corrected the option(prot); handling) NOTE: also removes line breakes and timings in the output due to option (prot), which were unaccounted and undocumented
  • Property mode set to 100644
File size: 3.0 KB
Line 
1/*****************************************
2*  Computer Algebra System SINGULAR      *
3*****************************************/
4/*
5* ABSTRACT: Gauss-Manin system normal form
6*/
7
8
9
10
11#include <kernel/mod2.h>
12
13#ifdef HAVE_GMS
14
15#include "gms.h"
16
17#include <coeffs/numbers.h>
18#include <kernel/polys.h>
19#include <Singular/febase.h>
20
21#include "ipid.h"
22
23lists gmsNF(ideal p,ideal g,matrix B,int D,int K)
24{
25  ideal r=idInit(IDELEMS(p),1);
26  ideal q=idInit(IDELEMS(p),1);
27
28  matrix B0=mpNew(MATROWS(B),MATCOLS(B));
29  for(int i=1;i<=MATROWS(B0);i++)
30     for(int j=1;j<=MATCOLS(B0);j++)
31      if(MATELEM(B,i,j)!=NULL)
32        MATELEM(B0,i,j)=pDiff(MATELEM(B,i,j),i+1);
33
34  for(int k=0;k<IDELEMS(p);k++)
35  {
36    while(p->m[k]!=NULL&&pGetExp(p->m[k],1)<=K)
37    {
38      int j=0;
39      while(j<IDELEMS(g)&&!pLmDivisibleBy(g->m[j],p->m[k]))
40        j++;
41
42      if(j<IDELEMS(g))
43      {
44        poly m=pDivideM(pHead(p->m[k]),pHead(g->m[j]));
45        p->m[k]=pSub(p->m[k],ppMult_mm(g->m[j],m));
46        pIncrExp(m,1);
47        pSetm(m);
48        for(int i=0;i<MATROWS(B);i++)
49        {
50          poly m0=pDiff(m,i+2);
51          if(MATELEM(B0,i+1,j+1)!=NULL)
52            p->m[k]=pAdd(p->m[k],ppMult_mm(MATELEM(B0,i+1,j+1),m));
53          if(MATELEM(B,i+1,j+1)!=NULL&&m0!=NULL)
54            p->m[k]=pAdd(p->m[k],ppMult_mm(MATELEM(B,i+1,j+1),m0));
55          pDelete(&m0);
56        }
57        pDelete(&m);
58      }
59      else
60      {
61        poly p0=p->m[k];
62        pIter(p->m[k]);
63        pNext(p0)=NULL;
64        r->m[k]=pAdd(r->m[k],p0);
65      }
66
67      while(p->m[k]!=NULL&&pGetExp(p->m[k],1)<=K&&pWTotaldegree(p->m[k])>D)
68      {
69        int i=pGetExp(p->m[k],1);
70        do
71        {
72          poly p0=p->m[k];
73          pIter(p->m[k]);
74          pNext(p0)=NULL;
75          q->m[k]=pAdd(q->m[k],p0);
76        }while(p->m[k]!=NULL&&pGetExp(p->m[k],1)==i);
77      }
78
79      pNormalize(p->m[k]);
80    }
81
82    q->m[k]=pAdd(q->m[k],p->m[k]);
83    p->m[k]=NULL;
84  }
85  idDelete(&p);
86  idDelete((ideal *)&B0);
87
88  id_Normalize(r, currRing);
89  id_Normalize(q, currRing);
90
91  lists l=(lists)omAllocBin(slists_bin);
92  l->Init(2);
93
94  l->m[0].rtyp=IDEAL_CMD;
95  l->m[0].data=r;
96  l->m[1].rtyp=IDEAL_CMD;
97  l->m[1].data=q;
98
99  return l;
100}
101
102
103BOOLEAN gmsNF(leftv res,leftv h)
104{
105  if(currRingHdl)
106  {
107    if(h&&h->Typ()==IDEAL_CMD)
108    {
109      ideal p=(ideal)h->Data();
110      h=h->next;
111      if(h&&h->Typ()==IDEAL_CMD)
112      {
113        ideal g=(ideal)h->Data();
114        h=h->next;
115        if(h&&h->Typ()==MATRIX_CMD)
116        {
117          matrix B=(matrix)h->Data();
118          h=h->next;
119          if(h&&h->Typ()==INT_CMD)
120          {
121            int D=(int)(long)h->Data();
122            h=h->next;
123            if(h&&h->Typ()==INT_CMD)
124            {
125              int K=(int)(long)h->Data();
126              res->rtyp=LIST_CMD;
127              res->data=(void *)gmsNF(idCopy(p),g,B,D,K);
128              return FALSE;
129            }
130          }
131        }
132      }
133    }
134    WerrorS("<ideal>,<ideal>,<matrix>,<int>,<int> expected");
135    return TRUE;
136  }
137  WerrorS("no ring active");
138  return TRUE;
139}
140
141#endif /* HAVE_GMS */
142
Note: See TracBrowser for help on using the repository browser.