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<idSize(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 = idSize(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  return;
314 }
315 
316 
317 /**
318  * reduces h initially with respect to g,
319  * returns false if h was initially reduced in the first place,
320  * returns true if reductions have taken place.
321  * assumes that h and g are in pReduced form and homogeneous in x of the same degree
322  */
323 bool ppreduceInitially(poly* hStar, const poly g, const ring r)
324 {
325  poly h = *hStar;
326  if (h==NULL || g==NULL)
327  return false;
328  p_Test(h,r);
329  p_Test(g,r);
330  poly hCache;
331  for (hCache=h; hCache; pIter(hCache))
332  if (p_LeadmonomDivisibleBy(g,hCache,r)) break;
333  if (hCache)
334  {
335  number gAlpha = p_GetCoeff(g,r);
336  poly hAlphaT = p_Init(r);
337  p_SetCoeff(hAlphaT,n_Copy(p_GetCoeff(hCache,r),r->cf),r);
338  p_SetExp(hAlphaT,1,p_GetExp(hCache,1,r)-p_GetExp(g,1,r),r);
339  for (int i=2; i<=r->N; i++)
340  p_SetExp(hAlphaT,i,0,r);
341  p_Setm(hAlphaT,r); p_Test(hAlphaT,r);
342  poly q1 = p_Mult_nn(h,gAlpha,r); p_Test(q1,r);
343  poly q2 = p_Mult_q(p_Copy(g,r),hAlphaT,r); p_Test(q2,r);
344  q2 = p_Neg(q2,r); p_Test(q2,r);
345  h = p_Add_q(q1,q2,r);
346  p_Test(h,r);
347  p_Test(g,r);
348  *hStar = h;
349  return true;
350  }
351  p_Test(h,r);
352  p_Test(g,r);
353  return false;
354 }
355 
356 
357 #ifndef NDEBUG
359 {
360  leftv u = args;
361  if ((u != NULL) && (u->Typ() == POLY_CMD))
362  {
363  leftv v = u->next;
364  if ((v != NULL) && (v->Typ() == POLY_CMD))
365  {
366  poly g,h;
367  omUpdateInfo();
368  Print("usedBytesBefore=%ld\n",om_Info.UsedBytes);
369  h = (poly) u->CopyD();
370  g = (poly) v->CopyD();
371  (void)ppreduceInitially(&h,g,currRing);
372  p_Delete(&h,currRing);
373  p_Delete(&g,currRing);
374  omUpdateInfo();
375  Print("usedBytesAfter=%ld\n",om_Info.UsedBytes);
376  h = (poly) u->CopyD();
377  g = (poly) v->CopyD();
378  (void)ppreduceInitially(&h,g,currRing);
379  p_Delete(&g,currRing);
380  res->rtyp = POLY_CMD;
381  res->data = (char*) h;
382  return FALSE;
383  }
384  }
385  return TRUE;
386 }
387 #endif //NDEBUG
388 
389 
390 /***
391  * reduces I initially with respect to itself and with respect to p-t.
392  * also sorts the generators of I with respect to the leading monomials in descending order.
393  * assumes that I is generated by elements which are homogeneous in x of the same degree.
394  **/
395 bool ppreduceInitially(ideal I, const number p, const ring r)
396 {
397  int m=idSize(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)
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  int n=idSize(I);
484  int j;
485  for (j=n-1; j>0; j--)
486  {
487  if (p_LmCmp(I->m[j], I->m[j-1],r)>0)
488  {
489  poly cache = I->m[j];
490  I->m[j] = I->m[j-1];
491  I->m[j-1] = cache;
492  }
493  else
494  break;
495  }
496 
497  /***
498  * the first pass. removing terms with the same monomials in x as lt(g_i) out of g_j for i<j
499  * removing terms with the same monomials in x as lt(g_j) out of g_k for j<k
500  **/
501  for (int i=0; i<j; i++)
502  if (ppreduceInitially(&I->m[j], I->m[i], r))
503  pReduce(I->m[j],p,r);
504  for (int k=j+1; k<n; k++)
505  if (ppreduceInitially(&I->m[k], I->m[j], r))
506  {
507  pReduce(I->m[k],p,r);
508  for (int l=j+1; l<k; l++)
509  if (ppreduceInitially(&I->m[k], I->m[l], r))
510  pReduce(I->m[k],p,r);
511  }
512 
513  /***
514  * the second pass. removing terms divisible by lt(g_j) and lt(g_k) out of g_i for i<j<k
515  * removing terms divisible by lt(g_k) out of g_j for j<k
516  **/
517  for (int i=0; i<j; i++)
518  for (int k=j; k<n; k++)
519  if (ppreduceInitially(&I->m[i], I->m[k], r))
520  pReduce(I->m[i],p,r);
521  for (int k=j; k<n-1; k++)
522  for (int l=k+1; l<n; l++)
523  if (ppreduceInitially(&I->m[k], I->m[l], r))
524  pReduce(I->m[k],p,r);
525 
526  /***
527  * removes the elements of I which have been reduced to 0 in the previous two passes
528  **/
529  idSkipZeroes(I);
530  id_Test(I,r);
531  return j;
532 }
533 
534 
535 #ifndef NDEBUG
537 {
538  leftv u = args;
539  if ((u != NULL) && (u->Typ() == IDEAL_CMD))
540  {
541  leftv v = u->next;
542  if ((v != NULL) && (v->Typ() == NUMBER_CMD))
543  {
544  leftv w = v->next;
545  if ((w != NULL) && (w->Typ() == POLY_CMD))
546  {
547  ideal I; number p; poly g;
548  omUpdateInfo();
549  Print("usedBytesBefore=%ld\n",om_Info.UsedBytes);
550  I = (ideal) u->CopyD();
551  p = (number) v->CopyD();
552  g = (poly) w->CopyD();
553  (void) ppreduceInitially(I,p,g,currRing);
554  id_Delete(&I,currRing);
555  n_Delete(&p,currRing->cf);
556  omUpdateInfo();
557  Print("usedBytesAfter=%ld\n",om_Info.UsedBytes);
558  I = (ideal) u->CopyD();
559  p = (number) v->CopyD();
560  g = (poly) w->CopyD();
561  (void) ppreduceInitially(I,p,g,currRing);
562  n_Delete(&p,currRing->cf);
563  res->rtyp = IDEAL_CMD;
564  res->data = (char*) I;
565  return FALSE;
566  }
567  }
568  }
569  return TRUE;
570 }
571 #endif //NDEBUG
572 
573 
574 static poly ppNext(poly p, int l)
575 {
576  poly q = p;
577  for (int i=0; i<l; i++)
578  {
579  if (q==NULL)
580  break;
581  pIter(q);
582  }
583  return q;
584 }
585 
586 
587 static void sortMarks(const ideal H, const ring r, std::vector<mark> &T)
588 {
589  std::pair<int,int> pointerToTerm;
590  int k=T.size();
591  do
592  {
593  int j=0;
594  for (int i=1; i<k-1; i++)
595  {
596  int generatorA = T[i-1].first;
597  int termA = T[i-1].second;
598  int generatorB = T[i].first;
599  int termB = T[i].second;
600  if (p_LmCmp(ppNext(H->m[generatorA],termA),ppNext(H->m[generatorB],termB),r)<0)
601  {
602  mark cache=T[i-1];
603  T[i-1]=T[i];
604  T[i]=cache;
605  j = i;
606  }
607  }
608  k=j;
609  } while(k);
610  return;
611 }
612 
613 
614 static poly getTerm(const ideal H, const mark ab)
615 {
616  int a = ab.first;
617  int b = ab.second;
618  return ppNext(H->m[a],b);
619 }
620 
621 
622 static void adjustMarks(std::vector<mark> &T, const int newEntry)
623 {
624  for (unsigned i=0; i<T.size(); i++)
625  {
626  if (T[i].first>=newEntry)
627  T[i].first = T[i].first+1;
628  }
629  return;
630 }
631 
632 
633 static void cleanupMarks(const ideal H, std::vector<mark> &T)
634 {
635  for (unsigned i=0; i<T.size();)
636  {
637  if (getTerm(H,T[i])==NULL)
638  T.erase(T.begin()+i);
639  else
640  i++;
641  }
642  return;
643 }
644 
645 
646 /***
647  * reduces H initially with respect to itself, with respect to p-t,
648  * and with respect to G.
649  * assumes that the generators of H are homogeneous in x of the same degree,
650  * assumes that the generators of G are homogeneous in x of lesser degree.
651  **/
652 bool ppreduceInitially(ideal &H, const number p, const ideal G, const ring r)
653 {
654  /***
655  * Step 1: reduce H initially with respect to itself and with respect to p-t
656  **/
657  if (ppreduceInitially(H,p,r)) return true;
658 
659  /***
660  * Step 2: initialize an ideal I in which the reductions will take place-
661  * along the reduction it will be enlarged with elements that will be discarded at the end
662  * initialize a working list T which keeps track.
663  * the working list T is a vector of pairs of integer.
664  * if T contains a pair (i,j) then that means that in the i-th element of H
665  * term j and subsequent terms need to be checked for reduction.
666  * T is sorted by the ordering on the temrs the pairs correspond to.
667  **/
668  int m=idSize(H);
669  ideal I = idInit(m);
670  std::vector<mark> T;
671  for (int i=0; i<m; i++)
672  {
673  I->m[i]=H->m[i];
674  if (pNext(I->m[i])!=NULL)
675  T.push_back(std::make_pair<int,int>(i,1));
676  }
677 
678  /***
679  * Step 3: as long as the working list is not empty, successively reduce terms in it
680  * by adding suitable elements to I and reducing it initially with respect to itself
681  **/
682  int k=idSize(G);
683  while (T.size()>0)
684  {
685  sortMarks(I,r,T);
686  int i=0; for (; i<k; i++)
687  if (p_LeadmonomDivisibleBy(G->m[i],getTerm(I,T[0]),r)) break;
688  if (i<k)
689  {
690  poly g = p_One(r); poly h0 = getTerm(I,T[0]);
691  assume(h0!=NULL);
692  for (int j=2; j<=r->N; j++)
693  p_SetExp(g,j,p_GetExp(h0,j,r)-p_GetExp(G->m[i],j,r),r);
694  p_Setm(g,r);
695  g = p_Mult_q(g,p_Copy(G->m[i],r),r);
696  int newEntry = ppreduceInitially(I,p,g,r);
697  adjustMarks(T,newEntry);
698  }
699  else
700  T[0].second = T[0].second+1;
701  cleanupMarks(I,T);
702  }
703 
704  /***
705  * Step 4: cleanup, delete all polynomials in I which have been added in Step 3
706  **/
707  k=idSize(I);
708  for (int i=0; i<k; i++)
709  {
710  for (int j=0; j<m; j++)
711  {
712  if (p_LeadmonomDivisibleBy(H->m[j],I->m[i],r))
713  {
714  I->m[i]=NULL;
715  break;
716  }
717  }
718  }
719  id_Delete(&I,r);
720  return false;
721 }
722 
723 
724 #ifndef NDEBUG
726 {
727  leftv u = args;
728  if ((u != NULL) && (u->Typ() == IDEAL_CMD))
729  {
730  leftv v = u->next;
731  if ((v != NULL) && (v->Typ() == NUMBER_CMD))
732  {
733  leftv w = v->next;
734  if ((w != NULL) && (w->Typ() == IDEAL_CMD))
735  {
736  ideal H,G; number p;
737  omUpdateInfo();
738  Print("usedBytesBefore=%ld\n",om_Info.UsedBytes);
739  H = (ideal) u->CopyD();
740  p = (number) v->CopyD();
741  G = (ideal) w->CopyD();
742  (void) ppreduceInitially(H,p,G,currRing);
743  n_Delete(&p,currRing->cf);
744  id_Delete(&G,currRing);
745  res->rtyp = IDEAL_CMD;
746  res->data = (char*) H;
747  return FALSE;
748  }
749  }
750  }
751  return TRUE;
752 }
753 #endif //NDEBUG
754 
755 /**
756  * reduces I initially with respect to itself.
757  * assumes that the generators of I are homogeneous in x and that p-t is in I.
758  */
759 bool ppreduceInitially(ideal I, const ring r, const number p)
760 {
761  assume(!n_IsUnit(p,r->cf));
762 
763  /***
764  * Step 1: split up I into components of same degree in x
765  * the lowest component should only contain p-t
766  **/
767  std::map<long,ideal> H; int n = idSize(I);
768  for (int i=0; i<n; i++)
769  {
770  I->m[i] = p_Cleardenom(I->m[i],r);
771  long d = 0;
772  for (int j=2; j<=r->N; j++)
773  d += p_GetExp(I->m[i],j,r);
774  std::map<long,ideal>::iterator it = H.find(d);
775  if (it != H.end())
776  idInsertPoly(it->second,I->m[i]);
777  else
778  {
779  std::pair<long,ideal> Hd(d,idInit(1));
780  Hd.second->m[0] = I->m[i];
781  H.insert(Hd);
782  }
783  }
784 
785  std::map<long,ideal>::iterator it=H.begin();
786  ideal Hi = it->second;
787  idShallowDelete(&Hi);
788  it++;
789  Hi = it->second;
790 
791  /***
792  * Step 2: reduce each component initially with respect to itself
793  * and all lower components
794  **/
795  if (ppreduceInitially(Hi,p,r)) return true;
796  id_Test(Hi,r);
797  id_Test(I,r);
798 
799  ideal G = idInit(n); int m=0;
800  ideal GG = (ideal) omAllocBin(sip_sideal_bin);
801  GG->nrows = 1; GG->rank = 1; GG->m=NULL;
802 
803  for (it++; it!=H.end(); it++)
804  {
805  int l=idSize(Hi); int k=l; poly cache;
806  /**
807  * sorts Hi according to degree in t in descending order
808  * (lowest first, highest last)
809  */
810  do
811  {
812  int j=0;
813  for (int i=1; i<k; i++)
814  {
815  if (p_GetExp(Hi->m[i-1],1,r)<p_GetExp(Hi->m[i],1,r))
816  {
817  cache=Hi->m[i-1];
818  Hi->m[i-1]=Hi->m[i];
819  Hi->m[i]=cache;
820  j = i;
821  }
822  }
823  k=j;
824  } while(k);
825  int kG=n-m, kH=0;
826  for (int i=n-m-l; i<n; i++)
827  {
828  if (kG==n)
829  {
830  memcpy(&(G->m[i]),&(Hi->m[kH]),(n-i)*sizeof(poly));
831  break;
832  }
833  if (kH==l)
834  break;
835  if (p_GetExp(G->m[kG],1,r)>p_GetExp(Hi->m[kH],1,r))
836  G->m[i] = G->m[kG++];
837  else
838  G->m[i] = Hi->m[kH++];
839  }
840  m += l; IDELEMS(GG) = m; GG->m = &G->m[n-m];
841  id_Test(it->second,r);
842  id_Test(GG,r);
843  if (ppreduceInitially(it->second,p,GG,r)) return true;
844  id_Test(it->second,r);
845  id_Test(GG,r);
846  idShallowDelete(&Hi); Hi = it->second;
847  }
848  idShallowDelete(&Hi);
849 
850  ptNormalize(I,p,r);
852  idShallowDelete(&G);
853  return false;
854 }
855 
856 
857 #ifndef NDEBUG
859 {
860  leftv u = args;
861  if ((u != NULL) && (u->Typ() == IDEAL_CMD))
862  {
863  leftv v = u->next;
864  if ((v != NULL) && (v->Typ() == NUMBER_CMD))
865  {
866  omUpdateInfo();
867  Print("usedBytesBefore=%ld\n",om_Info.UsedBytes);
868  ideal I = (ideal) u->CopyD();
869  number p = (number) v->Data();
870  (void) ppreduceInitially(I,currRing,p);
871  res->rtyp = IDEAL_CMD;
872  res->data = (char*) I;
873  return FALSE;
874  }
875  }
876  return TRUE;
877 }
878 #endif
879 
880 
881 // BOOLEAN ppreduceInitially(leftv res, leftv args)
882 // {
883 // leftv u = args;
884 // if ((u != NULL) && (u->Typ() == IDEAL_CMD))
885 // {
886 // ideal I = (ideal) u->CopyD();
887 // (void) ppreduceInitially(I,currRing);
888 // res->rtyp = IDEAL_CMD;
889 // res->data = (char*) I;
890 // return FALSE;
891 // }
892 // return TRUE;
893 // }
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 'a' and 'b' in Z/nZ, Z/2^kZ: computed as in the case Z in Z/pZ...
Definition: coeffs.h:685
const CanonicalForm int s
Definition: facAbsFact.cc:55
Class used for (list of) interpreter objects.
Definition: subexpr.h:83
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:516
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:721
static poly ppNext(poly p, int l)
std::pair< int, int > mark
#define FALSE
Definition: auxiliary.h:140
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 'n' represents the one element.
Definition: coeffs.h:469
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:973
#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:539
const CanonicalForm CFMap CFMap int &both_non_zero int n
Definition: cfEzgcd.cc:52
void id_Delete(ideal *h, ring r)
deletes an ideal/module/matrix
#define TRUE
Definition: auxiliary.h:144
void * ADDRESS
Definition: auxiliary.h:161
g
Definition: cfModGcd.cc:4031
void WerrorS(const char *s)
Definition: feFopen.cc:23
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:955
static number p_SetCoeff(poly p, number n, ring r)
Definition: p_polys.h:401
static poly p_Copy(poly p, const ring r)
returns a copy of p
Definition: p_polys.h:811
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:89
#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 'a' and 'b', i.e., a*b
Definition: coeffs.h:635
ring currRing
Widely used global variable which specifies the current polynomial ring for Singular interpreter and ...
Definition: polys.cc:12
static poly p_Head(poly p, const ring r)
Definition: p_polys.h:819
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:1318
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:465
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:405
static FORCE_INLINE number n_Add(number a, number b, const coeffs r)
return the sum of 'a' and 'b', i.e., a+b
Definition: coeffs.h:655
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 'a' is divisible 'b'; for r encoding a field: TRUE iff 'b' does not represent zero in Z:...
Definition: coeffs.h:771
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:1472
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:902
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 'n' represents the zero element.
Definition: coeffs.h:465
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:87
static void p_Delete(poly *p, const ring r)
Definition: p_polys.h:850
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:692
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:484
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:631
#define NULL
Definition: omList.c:10
static FORCE_INLINE number n_Copy(number n, const coeffs r)
return a copy of 'n'
Definition: coeffs.h:452
BOOLEAN ppreduceInitially0(leftv res, leftv args)
static FORCE_INLINE number n_Div(number a, number b, const coeffs r)
return the quotient of 'a' and 'b', i.e., a/b; raises an error if 'b' is not invertible in r exceptio...
Definition: coeffs.h:615
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:92
#define pNext(p)
Definition: monomials.h:43
void * Data()
Definition: subexpr.cc:1097
static FORCE_INLINE BOOLEAN n_Equal(number a, number b, const coeffs r)
TRUE iff 'a' and 'b' represent the same number; they may have different representations.
Definition: coeffs.h:461
static poly p_LmInit(poly p, const ring r)
Definition: p_polys.h:1263
static void p_Setm(poly p, const ring r)
Definition: p_polys.h:436
static long p_AddExp(poly p, int v, long ee, ring r)
Definition: p_polys.h:602
#define p_GetCoeff(p, r)
Definition: monomials.h:57
#define pPower(p, q)
Definition: polys.h:175
static int idSize(const ideal id)
Count the effective size of an ideal (without the trailing allocated zero-elements) ...
Definition: ideals.h:46
static poly p_Neg(poly p, const ring r)
Definition: p_polys.h:1018
static FORCE_INLINE void n_Delete(number *p, const coeffs r)
delete 'p'
Definition: coeffs.h:456
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:884
#define omFreeBin(addr, bin)
Definition: omAllocDecl.h:259
static Poly * h
Definition: janet.cc:978
int BOOLEAN
Definition: auxiliary.h:131
static poly p_Init(const ring r, omBin bin)
Definition: p_polys.h:1248
const poly b
Definition: syzextra.cc:213
poly p_Cleardenom(poly p, const ring r)
Definition: p_polys.cc:2682
static poly p_Mult_q(poly p, poly q, const ring r)
Definition: p_polys.h:1025
void * CopyD(int t)
Definition: subexpr.cc:662
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