ppinitialReduction.cc
Go to the documentation of this file.
2 #include <Singular/ipid.h>
3 
4 #include "singularWishlist.h"
5 #include "ppinitialReduction.h"
6 
7 #include <map>
8 #include <set>
9 #include <exception>
10 
11 
12 #ifndef NDEBUG
13 bool isOrderingLocalInT(const ring r)
14 {
15  poly one = p_One(r);
16  poly t = p_One(r);
17  p_SetExp(t,1,1,r);
18  p_Setm(t,r);
19  int s = p_LmCmp(one,t,r);
20  p_Delete(&one,r);
21  p_Delete(&t,r);
22  return (s==1);
23 }
24 #endif
25 
26 void divideByCommonGcd(poly &g, const ring r)
27 {
28  number commonGcd = n_Copy(p_GetCoeff(g,r),r->cf);
29  for (poly gCache=pNext(g); gCache; pIter(gCache))
30  {
31  number commonGcdCache = n_Gcd(commonGcd,p_GetCoeff(gCache,r),r->cf);
32  n_Delete(&commonGcd,r->cf);
33  commonGcd = commonGcdCache;
34  if (n_IsOne(commonGcd,r->cf))
35  {
36  n_Delete(&commonGcd,r->cf);
37  return;
38  }
39  }
40  for (poly gCache=g; gCache; pIter(gCache))
41  {
42  number oldCoeff = p_GetCoeff(gCache,r);
43  number newCoeff = n_Div(oldCoeff,commonGcd,r->cf);
44  p_SetCoeff(gCache,newCoeff,r);
45  }
46  p_Test(g,r);
47  n_Delete(&commonGcd,r->cf);
48  return;
49 }
50 
51 /***
52  * changes a polynomial g with the help p-t such that
53  * 1) each term of g has a distinct monomial in x
54  * 2) no term of g has a coefficient divisible by p
55  * in particular, this means that all g_\alpha can be obtained
56  * by reading the coefficients and that g is initially reduced
57  * with respect to p-t
58  **/
59 void pReduce(poly &g, const number p, const ring r)
60 {
61  if (g==NULL)
62  return;
63  p_Test(g,r);
64 
65  poly toBeChecked = pNext(g);
66  pNext(g) = NULL; poly gEnd = g;
67  poly gCache;
68 
69  number coeff, pPower; int power; poly subst;
70  while(toBeChecked)
71  {
72  for (gCache = g; gCache; pIter(gCache))
73  if (p_LeadmonomDivisibleBy(gCache,toBeChecked,r)) break;
74  if (gCache)
75  {
76  n_Power(p,p_GetExp(toBeChecked,1,r)-p_GetExp(gCache,1,r),&pPower,r->cf);
77  coeff = n_Mult(p_GetCoeff(toBeChecked,r),pPower,r->cf);
78  p_SetCoeff(gCache,n_Add(p_GetCoeff(gCache,r),coeff,r->cf),r);
79  n_Delete(&pPower,r->cf); n_Delete(&coeff,r->cf);
80  toBeChecked=p_LmDeleteAndNext(toBeChecked,r);
81  }
82  else
83  {
84  if (n_DivBy(p_GetCoeff(toBeChecked,r),p,r->cf))
85  {
86  power=1;
87  coeff=n_Div(p_GetCoeff(toBeChecked,r),p,r->cf);
88  while (n_DivBy(coeff,p,r->cf))
89  {
90  power++;
91  number coeff0 = n_Div(coeff,p,r->cf);
92  n_Delete(&coeff,r->cf);
93  coeff = coeff0;
94  coeff0 = NULL;
95  if (power<1)
96  {
97  WerrorS("pReduce: overflow in exponent");
98  throw 0;
99  }
100  }
101  subst=p_LmInit(toBeChecked,r);
102  p_AddExp(subst,1,power,r);
103  p_SetCoeff(subst,coeff,r);
104  p_Setm(subst,r); p_Test(subst,r);
105  toBeChecked=p_LmDeleteAndNext(toBeChecked,r);
106  toBeChecked=p_Add_q(toBeChecked,subst,r);
107  p_Test(toBeChecked,r);
108  }
109  else
110  {
111  pNext(gEnd)=toBeChecked;
112  pIter(gEnd); pIter(toBeChecked);
113  pNext(gEnd)=NULL;
114  p_Test(g,r);
115  }
116  }
117  }
118  p_Test(g,r);
119  divideByCommonGcd(g,r);
120  return;
121 }
122 
123 bool p_xLeadmonomDivisibleBy(const poly g, const poly f, const ring r)
124 {
125  poly gx = p_Head(g,r);
126  poly fx = p_Head(f,r);
127  p_SetExp(gx,1,0,r);
128  p_SetExp(fx,1,0,r);
129  p_Setm(gx,r);
130  p_Setm(fx,r);
131  bool b = p_LeadmonomDivisibleBy(gx,fx,r);
132  p_Delete(&gx,r);
133  p_Delete(&fx,r);
134  return b;
135 }
136 
137 void pReduceInhomogeneous(poly &g, const number p, const ring r)
138 {
139  if (g==NULL)
140  return;
141  p_Test(g,r);
142 
143  poly toBeChecked = pNext(g);
144  pNext(g) = NULL; poly gEnd = g;
145  poly gCache;
146 
147  number coeff, pPower; int power; poly subst;
148  while(toBeChecked)
149  {
150  for (gCache = g; gCache; pIter(gCache))
151  if (p_xLeadmonomDivisibleBy(gCache,toBeChecked,r)) break;
152  if (gCache)
153  {
154  n_Power(p,p_GetExp(toBeChecked,1,r)-p_GetExp(gCache,1,r),&pPower,r->cf);
155  coeff = n_Mult(p_GetCoeff(toBeChecked,r),pPower,r->cf);
156  p_SetCoeff(gCache,n_Add(p_GetCoeff(gCache,r),coeff,r->cf),r);
157  n_Delete(&pPower,r->cf); n_Delete(&coeff,r->cf);
158  toBeChecked=p_LmDeleteAndNext(toBeChecked,r);
159  }
160  else
161  {
162  if (n_DivBy(p_GetCoeff(toBeChecked,r),p,r->cf))
163  {
164  power=1;
165  coeff=n_Div(p_GetCoeff(toBeChecked,r),p,r->cf);
166  while (n_DivBy(coeff,p,r->cf))
167  {
168  power++;
169  number coeff0 = n_Div(coeff,p,r->cf);
170  n_Delete(&coeff,r->cf);
171  coeff = coeff0;
172  coeff0 = NULL;
173  if (power<1)
174  {
175  WerrorS("pReduce: overflow in exponent");
176  throw 0;
177  }
178  }
179  subst=p_LmInit(toBeChecked,r);
180  p_AddExp(subst,1,power,r);
181  p_SetCoeff(subst,coeff,r);
182  p_Setm(subst,r); p_Test(subst,r);
183  toBeChecked=p_LmDeleteAndNext(toBeChecked,r);
184  toBeChecked=p_Add_q(toBeChecked,subst,r);
185  p_Test(toBeChecked,r);
186  }
187  else
188  {
189  pNext(gEnd)=toBeChecked;
190  pIter(gEnd); pIter(toBeChecked);
191  pNext(gEnd)=NULL;
192  p_Test(g,r);
193  }
194  }
195  }
196  p_Test(g,r);
197  divideByCommonGcd(g,r);
198  return;
199 }
200 
201 void ptNormalize(poly* gStar, const number p, const ring r)
202 {
203  poly g = *gStar;
204  if (g==NULL || n_DivBy(p_GetCoeff(g,r),p,r->cf))
205  return;
206  p_Test(g,r);
207 
208  // create p-t
209  poly pt = p_Init(r);
210  p_SetCoeff(pt,n_Copy(p,r->cf),r);
211 
212  pNext(pt) = p_Init(r);
213  p_SetExp(pNext(pt),1,1,r);
214  p_Setm(pNext(pt),r);
215  p_SetCoeff(pNext(pt),n_Init(-1,r->cf),r);
216 
217  // make g monic with the help of p-t
218  number a,b;
219  number gcd = n_ExtGcd(p_GetCoeff(g,r),p,&a,&b,r->cf);
220  assume(n_IsUnit(gcd,r->cf));
221  // now a*leadcoef(g)+b*p = gcd with gcd being a unit
222  // so a*g+b*(p-t)*leadmonom(g) should have a unit as leading coefficient
223  // but first check whether b is 0,
224  // since p_Mult_nn doesn't allow 0 as number input
225  if (n_IsZero(b,r->cf))
226  {
227  n_Delete(&a,r->cf);
228  n_Delete(&b,r->cf);
229  n_Delete(&gcd,r->cf);
230  p_Delete(&pt,r);
231  return;
232  }
233  poly m = p_Head(g,r);
234  p_SetCoeff(m,n_Init(1,r->cf),r);
235  g = p_Add_q(p_Mult_nn(g,a,r),p_Mult_nn(p_Mult_mm(pt,m,r),b,r),r);
236  n_Delete(&a,r->cf);
237  n_Delete(&b,r->cf);
238  n_Delete(&gcd,r->cf);
239  p_Delete(&m,r);
240 
241  p_Test(g,r);
242  return;
243 }
244 
245 void ptNormalize(ideal I, const number p, const ring r)
246 {
247  for (int i=0; i<IDELEMS(I); i++)
248  ptNormalize(&(I->m[i]),p,r);
249  return;
250 }
251 
252 #ifndef NDEBUG
254 {
255  leftv u = args;
256  if ((u != NULL) && (u->Typ() == IDEAL_CMD))
257  {
258  leftv v = u->next;
259  if ((v!=NULL) && (v->Typ()==NUMBER_CMD))
260  {
261  omUpdateInfo();
262  Print("usedBytesBefore=%ld\n",om_Info.UsedBytes);
263  ideal I = (ideal) u->CopyD();
264  number p = (number) v->CopyD();
265  ptNormalize(I,p,currRing);
266  n_Delete(&p,currRing->cf);
267  res->rtyp = IDEAL_CMD;
268  res->data = (char*) I;
269  return FALSE;
270  }
271  }
272  return TRUE;
273 }
274 #endif //NDEBUG
275 
276 #ifndef NDEBUG
278 {
279  leftv u = args;
280  if ((u != NULL) && (u->Typ() == POLY_CMD))
281  {
282  poly g; number p = n_Init(3,currRing->cf);
283  omUpdateInfo();
284  Print("usedBytesBefore=%ld\n",om_Info.UsedBytes);
285  g = (poly) u->CopyD();
286  (void) pReduce(g,p,currRing);
287  p_Delete(&g,currRing);
288  omUpdateInfo();
289  Print("usedBytesAfter=%ld\n",om_Info.UsedBytes);
290  g = (poly) u->CopyD();
291  (void) pReduce(g,p,currRing);
292  n_Delete(&p,currRing->cf);
293  res->rtyp = POLY_CMD;
294  res->data = (char*) g;
295  return FALSE;
296  }
297  return TRUE;
298 }
299 #endif //NDEBUG
300 
301 void pReduce(ideal &I, const number p, const ring r)
302 {
303  int k = IDELEMS(I);
304  for (int i=0; i<k; i++)
305  {
306  if (I->m[i]!=NULL)
307  {
308  number c = p_GetCoeff(I->m[i],r);
309  if (!n_Equal(p,c,r->cf))
310  pReduce(I->m[i],p,r);
311  }
312  }
313 }
314 
315 
316 /**
317  * reduces h initially with respect to g,
318  * returns false if h was initially reduced in the first place,
319  * returns true if reductions have taken place.
320  * assumes that h and g are in pReduced form and homogeneous in x of the same degree
321  */
322 bool ppreduceInitially(poly* hStar, const poly g, const ring r)
323 {
324  poly h = *hStar;
325  if (h==NULL || g==NULL)
326  return false;
327  p_Test(h,r);
328  p_Test(g,r);
329  poly hCache;
330  for (hCache=h; hCache; pIter(hCache))
331  if (p_LeadmonomDivisibleBy(g,hCache,r)) break;
332  if (hCache)
333  {
334  number gAlpha = p_GetCoeff(g,r);
335  poly hAlphaT = p_Init(r);
336  p_SetCoeff(hAlphaT,n_Copy(p_GetCoeff(hCache,r),r->cf),r);
337  p_SetExp(hAlphaT,1,p_GetExp(hCache,1,r)-p_GetExp(g,1,r),r);
338  for (int i=2; i<=r->N; i++)
339  p_SetExp(hAlphaT,i,0,r);
340  p_Setm(hAlphaT,r); p_Test(hAlphaT,r);
341  poly q1 = p_Mult_nn(h,gAlpha,r); p_Test(q1,r);
342  poly q2 = p_Mult_q(p_Copy(g,r),hAlphaT,r); p_Test(q2,r);
343  q2 = p_Neg(q2,r); p_Test(q2,r);
344  h = p_Add_q(q1,q2,r);
345  p_Test(h,r);
346  p_Test(g,r);
347  *hStar = h;
348  return true;
349  }
350  p_Test(h,r);
351  p_Test(g,r);
352  return false;
353 }
354 
355 
356 #ifndef NDEBUG
358 {
359  leftv u = args;
360  if ((u != NULL) && (u->Typ() == POLY_CMD))
361  {
362  leftv v = u->next;
363  if ((v != NULL) && (v->Typ() == POLY_CMD))
364  {
365  poly g,h;
366  omUpdateInfo();
367  Print("usedBytesBefore=%ld\n",om_Info.UsedBytes);
368  h = (poly) u->CopyD();
369  g = (poly) v->CopyD();
370  (void)ppreduceInitially(&h,g,currRing);
371  p_Delete(&h,currRing);
372  p_Delete(&g,currRing);
373  omUpdateInfo();
374  Print("usedBytesAfter=%ld\n",om_Info.UsedBytes);
375  h = (poly) u->CopyD();
376  g = (poly) v->CopyD();
377  (void)ppreduceInitially(&h,g,currRing);
378  p_Delete(&g,currRing);
379  res->rtyp = POLY_CMD;
380  res->data = (char*) h;
381  return FALSE;
382  }
383  }
384  return TRUE;
385 }
386 #endif //NDEBUG
387 
388 
389 /***
390  * reduces I initially with respect to itself and with respect to p-t.
391  * also sorts the generators of I with respect to the leading monomials in descending order.
392  * assumes that I is generated by elements which are homogeneous in x of the same degree.
393  **/
394 bool ppreduceInitially(ideal I, const number p, const ring r)
395 {
396  idSkipZeroes(I);
397  int m=IDELEMS(I),n=m; poly cache;
398  do
399  {
400  int j=0;
401  for (int i=1; i<n; i++)
402  {
403  if (p_LmCmp(I->m[i-1],I->m[i],r)<0) /*p_LmCmp(p,q): requires: p,q!=NULL*/
404  {
405  cache=I->m[i-1];
406  I->m[i-1]=I->m[i];
407  I->m[i]=cache;
408  j = i;
409  }
410  }
411  n=j;
412  } while(n);
413  for (int i=0; i<m; i++)
414  pReduce(I->m[i],p,r);
415 
416  /***
417  * the first pass. removing terms with the same monomials in x as lt(g_i) out of g_j for i<j
418  **/
419  for (int i=0; i<m-1; i++)
420  for (int j=i+1; j<m; j++)
421  if (ppreduceInitially(&I->m[j], I->m[i], r))
422  pReduce(I->m[j],p,r);
423 
424  /***
425  * the second pass. removing terms divisible by lt(g_j) out of g_i for i<j
426  **/
427  for (int i=0; i<m-1; i++)
428  for (int j=i+1; j<m; j++)
429  if (ppreduceInitially(&I->m[i], I->m[j],r))
430  pReduce(I->m[i],p,r);
431 
432  /***
433  * removes the elements of I which have been reduced to 0 in the previous two passes
434  **/
435  idSkipZeroes(I);
436  return false;
437 }
438 
439 
440 #ifndef NDEBUG
442 {
443  leftv u = args;
444  if ((u != NULL) && (u->Typ() == IDEAL_CMD))
445  {
446  leftv v = u->next;
447  if ((v != NULL) && (v->Typ() == NUMBER_CMD))
448  {
449  ideal I; number p;
450  omUpdateInfo();
451  Print("usedBytesBefore=%ld\n",om_Info.UsedBytes);
452  I = (ideal) u->CopyD();
453  p = (number) v->CopyD();
454  (void) ppreduceInitially(I,p,currRing);
455  id_Delete(&I,currRing);
456  n_Delete(&p,currRing->cf);
457  omUpdateInfo();
458  Print("usedBytesAfter=%ld\n",om_Info.UsedBytes);
459  I = (ideal) u->CopyD();
460  p = (number) v->CopyD();
461  (void) ppreduceInitially(I,p,currRing);
462  n_Delete(&p,currRing->cf);
463  res->rtyp = IDEAL_CMD;
464  res->data = (char*) I;
465  return FALSE;
466  }
467  }
468  return TRUE;
469 }
470 #endif //NDEBUG
471 
472 
473 /***
474  * inserts g into I and reduces I with respect to itself and p-t
475  * returns the position in I in which g was inserted
476  * assumes that I was already sorted and initially reduced in the first place
477  **/
478 int ppreduceInitially(ideal I, const number p, const poly g, const ring r)
479 {
480  id_Test(I,r);
481  p_Test(g,r);
482  idInsertPoly(I,g);
483  idSkipZeroes(I);
484  int n=IDELEMS(I);
485  int j;
486  for (j=n-1; j>0; j--)
487  {
488  if (p_LmCmp(I->m[j], I->m[j-1],r)>0) /*p_LmCmp(p,q) requires: p,q!=NULL */
489  {
490  poly cache = I->m[j];
491  I->m[j] = I->m[j-1];
492  I->m[j-1] = cache;
493  }
494  else
495  break;
496  }
497 
498  /***
499  * the first pass. removing terms with the same monomials in x as lt(g_i) out of g_j for i<j
500  * removing terms with the same monomials in x as lt(g_j) out of g_k for j<k
501  **/
502  for (int i=0; i<j; i++)
503  if (ppreduceInitially(&I->m[j], I->m[i], r))
504  pReduce(I->m[j],p,r);
505  for (int k=j+1; k<n; k++)
506  if (ppreduceInitially(&I->m[k], I->m[j], r))
507  {
508  pReduce(I->m[k],p,r);
509  for (int l=j+1; l<k; l++)
510  if (ppreduceInitially(&I->m[k], I->m[l], r))
511  pReduce(I->m[k],p,r);
512  }
513 
514  /***
515  * the second pass. removing terms divisible by lt(g_j) and lt(g_k) out of g_i for i<j<k
516  * removing terms divisible by lt(g_k) out of g_j for j<k
517  **/
518  for (int i=0; i<j; i++)
519  for (int k=j; k<n; k++)
520  if (ppreduceInitially(&I->m[i], I->m[k], r))
521  pReduce(I->m[i],p,r);
522  for (int k=j; k<n-1; k++)
523  for (int l=k+1; l<n; l++)
524  if (ppreduceInitially(&I->m[k], I->m[l], r))
525  pReduce(I->m[k],p,r);
526 
527  /***
528  * removes the elements of I which have been reduced to 0 in the previous two passes
529  **/
530  idSkipZeroes(I);
531  id_Test(I,r);
532  return j;
533 }
534 
535 
536 #ifndef NDEBUG
538 {
539  leftv u = args;
540  if ((u != NULL) && (u->Typ() == IDEAL_CMD))
541  {
542  leftv v = u->next;
543  if ((v != NULL) && (v->Typ() == NUMBER_CMD))
544  {
545  leftv w = v->next;
546  if ((w != NULL) && (w->Typ() == POLY_CMD))
547  {
548  ideal I; number p; poly g;
549  omUpdateInfo();
550  Print("usedBytesBefore=%ld\n",om_Info.UsedBytes);
551  I = (ideal) u->CopyD();
552  p = (number) v->CopyD();
553  g = (poly) w->CopyD();
554  (void) ppreduceInitially(I,p,g,currRing);
555  id_Delete(&I,currRing);
556  n_Delete(&p,currRing->cf);
557  omUpdateInfo();
558  Print("usedBytesAfter=%ld\n",om_Info.UsedBytes);
559  I = (ideal) u->CopyD();
560  p = (number) v->CopyD();
561  g = (poly) w->CopyD();
562  (void) ppreduceInitially(I,p,g,currRing);
563  n_Delete(&p,currRing->cf);
564  res->rtyp = IDEAL_CMD;
565  res->data = (char*) I;
566  return FALSE;
567  }
568  }
569  }
570  return TRUE;
571 }
572 #endif //NDEBUG
573 
574 
575 static poly ppNext(poly p, int l)
576 {
577  poly q = p;
578  for (int i=0; i<l; i++)
579  {
580  if (q==NULL)
581  break;
582  pIter(q);
583  }
584  return q;
585 }
586 
587 
588 static void sortMarks(const ideal H, const ring r, std::vector<mark> &T)
589 {
590  std::pair<int,int> pointerToTerm;
591  int k=T.size();
592  do
593  {
594  int j=0;
595  for (int i=1; i<k-1; i++)
596  {
597  int generatorA = T[i-1].first;
598  int termA = T[i-1].second;
599  int generatorB = T[i].first;
600  int termB = T[i].second;
601  if (p_LmCmp(ppNext(H->m[generatorA],termA),ppNext(H->m[generatorB],termB),r)<0)
602  {
603  mark cache=T[i-1];
604  T[i-1]=T[i];
605  T[i]=cache;
606  j = i;
607  }
608  }
609  k=j;
610  } while(k);
611  return;
612 }
613 
614 
615 static poly getTerm(const ideal H, const mark ab)
616 {
617  int a = ab.first;
618  int b = ab.second;
619  return ppNext(H->m[a],b);
620 }
621 
622 
623 static void adjustMarks(std::vector<mark> &T, const int newEntry)
624 {
625  for (unsigned i=0; i<T.size(); i++)
626  {
627  if (T[i].first>=newEntry)
628  T[i].first = T[i].first+1;
629  }
630  return;
631 }
632 
633 
634 static void cleanupMarks(const ideal H, std::vector<mark> &T)
635 {
636  for (unsigned i=0; i<T.size();)
637  {
638  if (getTerm(H,T[i])==NULL)
639  T.erase(T.begin()+i);
640  else
641  i++;
642  }
643  return;
644 }
645 
646 
647 /***
648  * reduces H initially with respect to itself, with respect to p-t,
649  * and with respect to G.
650  * assumes that the generators of H are homogeneous in x of the same degree,
651  * assumes that the generators of G are homogeneous in x of lesser degree.
652  **/
653 bool ppreduceInitially(ideal &H, const number p, const ideal G, const ring r)
654 {
655  /***
656  * Step 1: reduce H initially with respect to itself and with respect to p-t
657  **/
658  if (ppreduceInitially(H,p,r)) return true;
659 
660  /***
661  * Step 2: initialize an ideal I in which the reductions will take place-
662  * along the reduction it will be enlarged with elements that will be discarded at the end
663  * initialize a working list T which keeps track.
664  * the working list T is a vector of pairs of integer.
665  * if T contains a pair (i,j) then that means that in the i-th element of H
666  * term j and subsequent terms need to be checked for reduction.
667  * T is sorted by the ordering on the temrs the pairs correspond to.
668  **/
669  int m=IDELEMS(H);
670  ideal I = idInit(m);
671  std::vector<mark> T;
672  for (int i=0; i<m; i++)
673  {
674  if(H->m[i]!=NULL)
675  {
676  I->m[i]=H->m[i];
677  if (pNext(I->m[i])!=NULL)
678  T.push_back(std::pair<int,int>(i,1));
679  }
680  }
681 
682  /***
683  * Step 3: as long as the working list is not empty, successively reduce terms in it
684  * by adding suitable elements to I and reducing it initially with respect to itself
685  **/
686  int k=IDELEMS(G);
687  while (T.size()>0)
688  {
689  sortMarks(I,r,T);
690  int i=0;
691  for (; i<k; i++)
692  {
693  if(G->m[i]!=NULL)
694  {
695  if (p_LeadmonomDivisibleBy(G->m[i],getTerm(I,T[0]),r)) break;
696  }
697  }
698  if (i<k)
699  {
700  poly g = p_One(r); poly h0 = getTerm(I,T[0]);
701  assume(h0!=NULL);
702  for (int j=2; j<=r->N; j++)
703  p_SetExp(g,j,p_GetExp(h0,j,r)-p_GetExp(G->m[i],j,r),r);
704  p_Setm(g,r);
705  g = p_Mult_q(g,p_Copy(G->m[i],r),r);
706  int newEntry = ppreduceInitially(I,p,g,r);
707  adjustMarks(T,newEntry);
708  }
709  else
710  T[0].second = T[0].second+1;
711  cleanupMarks(I,T);
712  }
713 
714  /***
715  * Step 4: cleanup, delete all polynomials in I which have been added in Step 3
716  **/
717  k=IDELEMS(I);
718  for (int i=0; i<k; i++)
719  {
720  if(I->m[i]!=NULL)
721  {
722  for (int j=0; j<m; j++)
723  {
724  if ((H->m[j]!=NULL)
725  && (p_LeadmonomDivisibleBy(H->m[j],I->m[i],r)))
726  {
727  I->m[i]=NULL;
728  break;
729  }
730  }
731  }
732  }
733  id_Delete(&I,r);
734  return false;
735 }
736 
737 
738 #ifndef NDEBUG
740 {
741  leftv u = args;
742  if ((u != NULL) && (u->Typ() == IDEAL_CMD))
743  {
744  leftv v = u->next;
745  if ((v != NULL) && (v->Typ() == NUMBER_CMD))
746  {
747  leftv w = v->next;
748  if ((w != NULL) && (w->Typ() == IDEAL_CMD))
749  {
750  ideal H,G; number p;
751  omUpdateInfo();
752  Print("usedBytesBefore=%ld\n",om_Info.UsedBytes);
753  H = (ideal) u->CopyD();
754  p = (number) v->CopyD();
755  G = (ideal) w->CopyD();
756  (void) ppreduceInitially(H,p,G,currRing);
757  n_Delete(&p,currRing->cf);
758  id_Delete(&G,currRing);
759  res->rtyp = IDEAL_CMD;
760  res->data = (char*) H;
761  return FALSE;
762  }
763  }
764  }
765  return TRUE;
766 }
767 #endif //NDEBUG
768 
769 /**
770  * reduces I initially with respect to itself.
771  * assumes that the generators of I are homogeneous in x and that p-t is in I.
772  */
773 bool ppreduceInitially(ideal I, const ring r, const number p)
774 {
775  assume(!n_IsUnit(p,r->cf));
776 
777  /***
778  * Step 1: split up I into components of same degree in x
779  * the lowest component should only contain p-t
780  **/
781  std::map<long,ideal> H; int n = IDELEMS(I);
782  for (int i=0; i<n; i++)
783  {
784  if(I->m[i]!=NULL)
785  {
786  I->m[i] = p_Cleardenom(I->m[i],r);
787  long d = 0;
788  for (int j=2; j<=r->N; j++)
789  d += p_GetExp(I->m[i],j,r);
790  std::map<long,ideal>::iterator it = H.find(d);
791  if (it != H.end())
792  idInsertPoly(it->second,I->m[i]);
793  else
794  {
795  std::pair<long,ideal> Hd(d,idInit(1));
796  Hd.second->m[0] = I->m[i];
797  H.insert(Hd);
798  }
799  }
800  }
801 
802  std::map<long,ideal>::iterator it=H.begin();
803  ideal Hi = it->second;
804  idShallowDelete(&Hi);
805  it++;
806  Hi = it->second;
807 
808  /***
809  * Step 2: reduce each component initially with respect to itself
810  * and all lower components
811  **/
812  if (ppreduceInitially(Hi,p,r)) return true;
813  idSkipZeroes(Hi);
814  id_Test(Hi,r);
815  id_Test(I,r);
816 
817  ideal G = idInit(n); int m=0;
818  ideal GG = (ideal) omAllocBin(sip_sideal_bin);
819  GG->nrows = 1; GG->rank = 1; GG->m=NULL;
820 
821  for (it++; it!=H.end(); it++)
822  {
823  int l=IDELEMS(Hi); int k=l; poly cache;
824  /**
825  * sorts Hi according to degree in t in descending order
826  * (lowest first, highest last)
827  */
828  do
829  {
830  int j=0;
831  for (int i=1; i<k; i++)
832  {
833  if (p_GetExp(Hi->m[i-1],1,r)<p_GetExp(Hi->m[i],1,r))
834  {
835  cache=Hi->m[i-1];
836  Hi->m[i-1]=Hi->m[i];
837  Hi->m[i]=cache;
838  j = i;
839  }
840  }
841  k=j;
842  } while(k);
843  int kG=n-m, kH=0;
844  for (int i=n-m-l; i<n; i++)
845  {
846  if (kG==n)
847  {
848  memcpy(&(G->m[i]),&(Hi->m[kH]),(n-i)*sizeof(poly));
849  break;
850  }
851  if (kH==l)
852  break;
853  if (p_GetExp(G->m[kG],1,r)>p_GetExp(Hi->m[kH],1,r))
854  G->m[i] = G->m[kG++];
855  else
856  G->m[i] = Hi->m[kH++];
857  }
858  m += l; IDELEMS(GG) = m; GG->m = &G->m[n-m];
859  id_Test(it->second,r);
860  id_Test(GG,r);
861  if (ppreduceInitially(it->second,p,GG,r)) return true;
862  id_Test(it->second,r);
863  id_Test(GG,r);
864  idShallowDelete(&Hi); Hi = it->second;
865  }
866  idShallowDelete(&Hi);
867 
868  ptNormalize(I,p,r);
870  idShallowDelete(&G);
871  return false;
872 }
873 
874 
875 #ifndef NDEBUG
877 {
878  leftv u = args;
879  if ((u != NULL) && (u->Typ() == IDEAL_CMD))
880  {
881  leftv v = u->next;
882  if ((v != NULL) && (v->Typ() == NUMBER_CMD))
883  {
884  omUpdateInfo();
885  Print("usedBytesBefore=%ld\n",om_Info.UsedBytes);
886  ideal I = (ideal) u->CopyD();
887  number p = (number) v->Data();
888  (void) ppreduceInitially(I,currRing,p);
889  res->rtyp = IDEAL_CMD;
890  res->data = (char*) I;
891  return FALSE;
892  }
893  }
894  return TRUE;
895 }
896 #endif
897 
898 
899 // BOOLEAN ppreduceInitially(leftv res, leftv args)
900 // {
901 // leftv u = args;
902 // if ((u != NULL) && (u->Typ() == IDEAL_CMD))
903 // {
904 // ideal I = (ideal) u->CopyD();
905 // (void) ppreduceInitially(I,currRing);
906 // res->rtyp = IDEAL_CMD;
907 // res->data = (char*) I;
908 // return FALSE;
909 // }
910 // return TRUE;
911 // }
BOOLEAN ppreduceInitially2(leftv res, leftv args)
CanonicalForm power(const CanonicalForm &f, int n)
exponentiation
#define omAllocBin(bin)
Definition: omAllocDecl.h:205
static FORCE_INLINE number n_Gcd(number a, number b, const coeffs r)
in Z: return the gcd of &#39;a&#39; and &#39;b&#39; in Z/nZ, Z/2^kZ: computed as in the case Z in Z/pZ...
Definition: coeffs.h:690
const CanonicalForm int s
Definition: facAbsFact.cc:55
Class used for (list of) interpreter objects.
Definition: subexpr.h:84
static FORCE_INLINE BOOLEAN n_IsUnit(number n, const coeffs r)
TRUE iff n has a multiplicative inverse in the given coeff field/ring r.
Definition: coeffs.h:519
const poly a
Definition: syzextra.cc:212
BOOLEAN reduceInitiallyDebug(leftv res, leftv args)
#define Print
Definition: emacs.cc:83
static poly p_LmDeleteAndNext(poly p, const ring r)
Definition: p_polys.h:720
static poly ppNext(poly p, int l)
std::pair< int, int > mark
#define FALSE
Definition: auxiliary.h:97
bool isOrderingLocalInT(const ring r)
return P p
Definition: myNF.cc:203
static FORCE_INLINE BOOLEAN n_IsOne(number n, const coeffs r)
TRUE iff &#39;n&#39; represents the one element.
Definition: coeffs.h:472
omBin sip_sideal_bin
Definition: simpleideals.cc:30
static void sortMarks(const ideal H, const ring r, std::vector< mark > &T)
bool ppreduceInitially(poly *hStar, const poly g, const ring r)
reduces h initially with respect to g, returns false if h was initially reduced in the first place...
static poly p_Mult_mm(poly p, poly m, const ring r)
Definition: p_polys.h:968
#define id_Test(A, lR)
Definition: simpleideals.h:80
static poly getTerm(const ideal H, const mark ab)
static FORCE_INLINE number n_Init(long i, const coeffs r)
a number representing i in the given coeff field/ring r
Definition: coeffs.h:542
void id_Delete(ideal *h, ring r)
deletes an ideal/module/matrix
#define TRUE
Definition: auxiliary.h:101
void * ADDRESS
Definition: auxiliary.h:118
g
Definition: cfModGcd.cc:4031
void WerrorS(const char *s)
Definition: feFopen.cc:24
int k
Definition: cfEzgcd.cc:93
static void cleanupMarks(const ideal H, std::vector< mark > &T)
static TreeM * G
Definition: janet.cc:38
int Typ()
Definition: subexpr.cc:979
static number p_SetCoeff(poly p, number n, ring r)
Definition: p_polys.h:407
static poly p_Copy(poly p, const ring r)
returns a copy of p
Definition: p_polys.h:804
CanonicalForm subst(const CanonicalForm &f, const CFList &a, const CFList &b, const CanonicalForm &Rstar, bool isFunctionField)
const CanonicalForm & GG
Definition: cfModGcd.cc:4017
void * data
Definition: subexpr.h:90
#define pIter(p)
Definition: monomials.h:44
poly res
Definition: myNF.cc:322
static FORCE_INLINE number n_Mult(number a, number b, const coeffs r)
return the product of &#39;a&#39; and &#39;b&#39;, i.e., a*b
Definition: coeffs.h:640
ring currRing
Widely used global variable which specifies the current polynomial ring for Singular interpreter and ...
Definition: polys.cc:10
static poly p_Head(poly p, const ring r)
Definition: p_polys.h:812
const ring r
Definition: syzextra.cc:208
bool p_xLeadmonomDivisibleBy(const poly g, const poly f, const ring r)
poly p_One(const ring r)
Definition: p_polys.cc:1313
static long p_GetExp(const poly p, const unsigned long iBitmask, const int VarOffset)
get a single variable exponent : the integer VarOffset encodes:
Definition: p_polys.h:464
int j
Definition: myNF.cc:70
BOOLEAN pReduceDebug(leftv res, leftv args)
omInfo_t om_Info
Definition: omStats.c:13
#define assume(x)
Definition: mod2.h:403
static FORCE_INLINE number n_Add(number a, number b, const coeffs r)
return the sum of &#39;a&#39; and &#39;b&#39;, i.e., a+b
Definition: coeffs.h:660
BOOLEAN ppreduceInitially3(leftv res, leftv args)
void pReduce(poly &g, const number p, const ring r)
static FORCE_INLINE BOOLEAN n_DivBy(number a, number b, const coeffs r)
test whether &#39;a&#39; is divisible &#39;b&#39;; for r encoding a field: TRUE iff &#39;b&#39; does not represent zero in Z:...
Definition: coeffs.h:787
BOOLEAN idInsertPoly(ideal h1, poly h2)
insert h2 into h1 (if h2 is not the zero polynomial) return TRUE iff h2 was indeed inserted ...
static int p_LmCmp(poly p, poly q, const ring r)
Definition: p_polys.h:1467
int m
Definition: cfEzgcd.cc:119
FILE * f
Definition: checklibs.c:7
int i
Definition: cfEzgcd.cc:123
static poly p_Mult_nn(poly p, number n, const ring r)
Definition: p_polys.h:895
CanonicalForm H
Definition: facAbsFact.cc:64
#define IDELEMS(i)
Definition: simpleideals.h:24
static FORCE_INLINE BOOLEAN n_IsZero(number n, const coeffs r)
TRUE iff &#39;n&#39; represents the zero element.
Definition: coeffs.h:468
void idSkipZeroes(ideal ide)
gives an ideal/module the minimal possible size
#define p_Test(p, r)
Definition: p_polys.h:160
void ptNormalize(poly *gStar, const number p, const ring r)
leftv next
Definition: subexpr.h:88
static void p_Delete(poly *p, const ring r)
Definition: p_polys.h:843
static FORCE_INLINE number n_ExtGcd(number a, number b, number *s, number *t, const coeffs r)
beware that ExtGCD is only relevant for a few chosen coeff. domains and may perform something unexpec...
Definition: coeffs.h:697
ideal idInit(int idsize, int rank)
initialise an ideal / module
Definition: simpleideals.cc:38
const Variable & v
< [in] a sqrfree bivariate poly
Definition: facBivar.h:37
static unsigned long p_SetExp(poly p, const unsigned long e, const unsigned long iBitmask, const int VarOffset)
set a single variable exponent : VarOffset encodes the position in p->exp
Definition: p_polys.h:483
static FORCE_INLINE void n_Power(number a, int b, number *res, const coeffs r)
fill res with the power a^b
Definition: coeffs.h:636
#define NULL
Definition: omList.c:10
static FORCE_INLINE number n_Copy(number n, const coeffs r)
return a copy of &#39;n&#39;
Definition: coeffs.h:455
BOOLEAN ppreduceInitially0(leftv res, leftv args)
static FORCE_INLINE number n_Div(number a, number b, const coeffs r)
return the quotient of &#39;a&#39; and &#39;b&#39;, i.e., a/b; raises an error if &#39;b&#39; is not invertible in r exceptio...
Definition: coeffs.h:619
int gcd(int a, int b)
Definition: walkSupport.cc:839
void omUpdateInfo()
Definition: omStats.c:24
void divideByCommonGcd(poly &g, const ring r)
const CanonicalForm & w
Definition: facAbsFact.cc:55
BOOLEAN ppreduceInitially1(leftv res, leftv args)
static void adjustMarks(std::vector< mark > &T, const int newEntry)
int rtyp
Definition: subexpr.h:93
#define pNext(p)
Definition: monomials.h:43
void * Data()
Definition: subexpr.cc:1121
static FORCE_INLINE BOOLEAN n_Equal(number a, number b, const coeffs r)
TRUE iff &#39;a&#39; and &#39;b&#39; represent the same number; they may have different representations.
Definition: coeffs.h:464
static poly p_LmInit(poly p, const ring r)
Definition: p_polys.h:1258
static void p_Setm(poly p, const ring r)
Definition: p_polys.h:228
static long p_AddExp(poly p, int v, long ee, ring r)
Definition: p_polys.h:601
#define p_GetCoeff(p, r)
Definition: monomials.h:57
#define pPower(p, q)
Definition: polys.h:187
static poly p_Neg(poly p, const ring r)
Definition: p_polys.h:1013
static FORCE_INLINE void n_Delete(number *p, const coeffs r)
delete &#39;p&#39;
Definition: coeffs.h:459
static jList * T
Definition: janet.cc:37
polyrec * poly
Definition: hilb.h:10
static poly p_Add_q(poly p, poly q, const ring r)
Definition: p_polys.h:877
#define omFreeBin(addr, bin)
Definition: omAllocDecl.h:259
static Poly * h
Definition: janet.cc:978
int BOOLEAN
Definition: auxiliary.h:88
static poly p_Init(const ring r, omBin bin)
Definition: p_polys.h:1243
const poly b
Definition: syzextra.cc:213
poly p_Cleardenom(poly p, const ring r)
Definition: p_polys.cc:2716
static poly p_Mult_q(poly p, poly q, const ring r)
Definition: p_polys.h:1020
void * CopyD(int t)
Definition: subexpr.cc:679
static BOOLEAN p_LeadmonomDivisibleBy(poly a, poly b, const ring r)
p_LmDivisibleBy checks also the divisibility of coefficients
idhdl h0
Definition: libparse.cc:1141
int l
Definition: cfEzgcd.cc:94
void pReduceInhomogeneous(poly &g, const number p, const ring r)
void idShallowDelete(ideal *h)
id_ShallowDelete deletes the monomials of the polynomials stored inside of it