gms.cc
Go to the documentation of this file.
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 
22 lists 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 
89 
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 
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 
#define omAllocBin(bin)
Definition: omAllocDecl.h:205
void id_Normalize(ideal I, const ring r)
normialize all polys in id
sleftv * m
Definition: lists.h:45
#define D(A)
Definition: gentable.cc:119
Class used for (list of) interpreter objects.
Definition: subexpr.h:83
#define pSetm(p)
Definition: polys.h:241
Definition: tok.h:98
#define pAdd(p, q)
Definition: polys.h:174
#define idDelete(H)
delete an ideal
Definition: ideals.h:31
Definition: lists.h:22
#define FALSE
Definition: auxiliary.h:140
Compatiblity layer for legacy polynomial operations (over currRing)
return P p
Definition: myNF.cc:203
#define ppMult_mm(p, m)
Definition: polys.h:172
#define TRUE
Definition: auxiliary.h:144
g
Definition: cfModGcd.cc:4031
void WerrorS(const char *s)
Definition: feFopen.cc:24
int k
Definition: cfEzgcd.cc:93
#define pLmDivisibleBy(a, b)
like pDivisibleBy, except that it is assumed that a!=NULL, b!=NULL
Definition: polys.h:128
int Typ()
Definition: subexpr.cc:976
void * data
Definition: subexpr.h:89
#define pIter(p)
Definition: monomials.h:44
poly res
Definition: myNF.cc:322
ring currRing
Widely used global variable which specifies the current polynomial ring for Singular interpreter and ...
Definition: polys.cc:12
#define pGetExp(p, i)
Exponent.
Definition: polys.h:41
#define pIncrExp(p, i)
Definition: polys.h:43
const ring r
Definition: syzextra.cc:208
#define pSub(a, b)
Definition: polys.h:258
int j
Definition: myNF.cc:70
#define pDivideM(a, b)
Definition: polys.h:265
pNormalize(P.p)
ip_smatrix * matrix
idhdl currRingHdl
Definition: ipid.cc:65
int m
Definition: cfEzgcd.cc:119
int i
Definition: cfEzgcd.cc:123
#define pHead(p)
returns newly allocated copy of Lm(p), coef is copied, next=NULL, p might be NULL ...
Definition: polys.h:67
#define IDELEMS(i)
Definition: simpleideals.h:24
leftv next
Definition: subexpr.h:87
INLINE_THIS void Init(int l=0)
Definition: lists.h:66
matrix mpNew(int r, int c)
create a r x c zero-matrix
Definition: matpol.cc:48
ideal idInit(int idsize, int rank)
initialise an ideal / module
Definition: simpleideals.cc:38
#define MATCOLS(i)
Definition: matpol.h:28
#define NULL
Definition: omList.c:10
slists * lists
Definition: mpr_numeric.h:146
b *CanonicalForm B
Definition: facBivar.cc:51
#define pWTotaldegree(p)
Definition: polys.h:254
#define pDelete(p_ptr)
Definition: polys.h:157
int rtyp
Definition: subexpr.h:92
#define pNext(p)
Definition: monomials.h:43
void * Data()
Definition: subexpr.cc:1118
Definition: tok.h:120
omBin slists_bin
Definition: lists.cc:23
#define pDiff(a, b)
Definition: polys.h:267
#define MATROWS(i)
Definition: matpol.h:27
lists gmsNF(ideal p, ideal g, matrix B, int D, int K)
Definition: gms.cc:22
polyrec * poly
Definition: hilb.h:10
static Poly * h
Definition: janet.cc:978
int BOOLEAN
Definition: auxiliary.h:131
void * CopyD(int t)
Definition: subexpr.cc:676
int l
Definition: cfEzgcd.cc:94
#define MATELEM(mat, i, j)
Definition: matpol.h:29