source: git/Singular/gms.cc @ fb1675

spielwiese
Last change on this file since fb1675 was fb1675, checked in by Hans Schoenemann <hannes@…>, 7 years ago
use include ".." for singular related .h, p8
  • 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
20#include "ipid.h"
21
22lists gmsNF(ideal p,ideal g,matrix B,int D,int K)
23{
24  ideal r=idInit(IDELEMS(p),1);
25  ideal q=idInit(IDELEMS(p),1);
26
27  matrix B0=mpNew(MATROWS(B),MATCOLS(B));
28  for(int i=1;i<=MATROWS(B0);i++)
29     for(int j=1;j<=MATCOLS(B0);j++)
30      if(MATELEM(B,i,j)!=NULL)
31        MATELEM(B0,i,j)=pDiff(MATELEM(B,i,j),i+1);
32
33  for(int k=0;k<IDELEMS(p);k++)
34  {
35    while(p->m[k]!=NULL&&pGetExp(p->m[k],1)<=K)
36    {
37      int j=0;
38      while(j<IDELEMS(g)&&!pLmDivisibleBy(g->m[j],p->m[k]))
39        j++;
40
41      if(j<IDELEMS(g))
42      {
43        poly m=pDivideM(pHead(p->m[k]),pHead(g->m[j]));
44        p->m[k]=pSub(p->m[k],ppMult_mm(g->m[j],m));
45        pIncrExp(m,1);
46        pSetm(m);
47        for(int i=0;i<MATROWS(B);i++)
48        {
49          poly m0=pDiff(m,i+2);
50          if(MATELEM(B0,i+1,j+1)!=NULL)
51            p->m[k]=pAdd(p->m[k],ppMult_mm(MATELEM(B0,i+1,j+1),m));
52          if(MATELEM(B,i+1,j+1)!=NULL&&m0!=NULL)
53            p->m[k]=pAdd(p->m[k],ppMult_mm(MATELEM(B,i+1,j+1),m0));
54          pDelete(&m0);
55        }
56        pDelete(&m);
57      }
58      else
59      {
60        poly p0=p->m[k];
61        pIter(p->m[k]);
62        pNext(p0)=NULL;
63        r->m[k]=pAdd(r->m[k],p0);
64      }
65
66      while(p->m[k]!=NULL&&pGetExp(p->m[k],1)<=K&&pWTotaldegree(p->m[k])>D)
67      {
68        int i=pGetExp(p->m[k],1);
69        do
70        {
71          poly p0=p->m[k];
72          pIter(p->m[k]);
73          pNext(p0)=NULL;
74          q->m[k]=pAdd(q->m[k],p0);
75        }while(p->m[k]!=NULL&&pGetExp(p->m[k],1)==i);
76      }
77
78      pNormalize(p->m[k]);
79    }
80
81    q->m[k]=pAdd(q->m[k],p->m[k]);
82    p->m[k]=NULL;
83  }
84  idDelete(&p);
85  idDelete((ideal *)&B0);
86
87  id_Normalize(r, currRing);
88  id_Normalize(q, currRing);
89
90  lists l=(lists)omAllocBin(slists_bin);
91  l->Init(2);
92
93  l->m[0].rtyp=IDEAL_CMD;
94  l->m[0].data=r;
95  l->m[1].rtyp=IDEAL_CMD;
96  l->m[1].data=q;
97
98  return l;
99}
100
101
102BOOLEAN gmsNF(leftv res,leftv h)
103{
104  if(currRingHdl)
105  {
106    if(h&&h->Typ()==IDEAL_CMD)
107    {
108      ideal p=(ideal)h->CopyD();
109      h=h->next;
110      if(h&&h->Typ()==IDEAL_CMD)
111      {
112        ideal g=(ideal)h->Data();
113        h=h->next;
114        if(h&&h->Typ()==MATRIX_CMD)
115        {
116          matrix B=(matrix)h->Data();
117          h=h->next;
118          if(h&&h->Typ()==INT_CMD)
119          {
120            int D=(int)(long)h->Data();
121            h=h->next;
122            if(h&&h->Typ()==INT_CMD)
123            {
124              int K=(int)(long)h->Data();
125              res->rtyp=LIST_CMD;
126              res->data=(void *)gmsNF(p,g,B,D,K);
127              return FALSE;
128            }
129          }
130        }
131      }
132    }
133    WerrorS("<ideal>,<ideal>,<matrix>,<int>,<int> expected");
134    return TRUE;
135  }
136  WerrorS("no ring active");
137  return TRUE;
138}
139
140#endif /* HAVE_GMS */
141
Note: See TracBrowser for help on using the repository browser.