source: git/Singular/gms.cc @ b2aa08

spielwiese
Last change on this file since b2aa08 was 6ce030f, checked in by Oleksandr Motsak <motsak@…>, 12 years ago
removal of the $Id$ svn tag from everywhere NOTE: the git SHA1 may be used instead (only on special places) NOTE: the libraries Singular/LIB/*.lib still contain the marker due to our current use of svn
  • 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#include "config.h"
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
21lists 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
101BOOLEAN 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
Note: See TracBrowser for help on using the repository browser.