source: git/Singular/dyn_modules/gfanlib/ppinitialReduction.cc @ e5618c

spielwiese
Last change on this file since e5618c was e5618c, checked in by Yue Ren <ren@…>, 9 years ago
chg: final version for Singular release
  • Property mode set to 100644
File size: 19.5 KB
Line 
1#include <libpolys/polys/monomials/p_polys.h>
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
13bool 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
26void 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 **/
59void 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
123void ptNormalize(poly* gStar, const number p, const ring r)
124{
125  poly g = *gStar;
126  if (g==NULL || n_DivBy(p_GetCoeff(g,r),p,r->cf))
127    return;
128  p_Test(g,r);
129
130  // create p-t
131  poly pt = p_Init(r);
132  p_SetCoeff(pt,n_Copy(p,r->cf),r);
133
134  pNext(pt) = p_Init(r);
135  p_SetExp(pNext(pt),1,1,r);
136  p_Setm(pNext(pt),r);
137  p_SetCoeff(pNext(pt),n_Init(-1,r->cf),r);
138
139  // make g monic with the help of p-t
140  number a,b;
141  number gcd = n_ExtGcd(p_GetCoeff(g,r),p,&a,&b,r->cf);
142  assume(n_IsUnit(gcd,r->cf));
143  // now a*leadcoef(g)+b*p = gcd with gcd being a unit
144  // so a*g+b*(p-t)*leadmonom(g) should have a unit as leading coefficient
145  // but first check whether b is 0,
146  // since p_Mult_nn doesn't allow 0 as number input
147  if (n_IsZero(b,r->cf))
148  {
149    n_Delete(&a,r->cf);
150    n_Delete(&b,r->cf);
151    n_Delete(&gcd,r->cf);
152    p_Delete(&pt,r);
153    return;
154  }
155  poly m = p_Head(g,r);
156  p_SetCoeff(m,n_Init(1,r->cf),r);
157  g = p_Add_q(p_Mult_nn(g,a,r),p_Mult_nn(p_Mult_mm(pt,m,r),b,r),r);
158  n_Delete(&a,r->cf);
159  n_Delete(&b,r->cf);
160  n_Delete(&gcd,r->cf);
161  p_Delete(&m,r);
162
163  p_Test(g,r);
164  return;
165}
166
167void ptNormalize(ideal I, const number p, const ring r)
168{
169  for (int i=0; i<idSize(I); i++)
170    ptNormalize(&(I->m[i]),p,r);
171  return;
172}
173
174#ifndef NDEBUG
175BOOLEAN ptNormalize(leftv res, leftv args)
176{
177  leftv u = args;
178  if ((u != NULL) && (u->Typ() == IDEAL_CMD))
179  {
180    leftv v = u->next;
181    if ((v!=NULL) && (v->Typ()==NUMBER_CMD))
182    {
183      omUpdateInfo();
184      Print("usedBytesBefore=%ld\n",om_Info.UsedBytes);
185      ideal I = (ideal) u->CopyD();
186      number p = (number) v->CopyD();
187      ptNormalize(I,p,currRing);
188      n_Delete(&p,currRing->cf);
189      res->rtyp = IDEAL_CMD;
190      res->data = (char*) I;
191      return FALSE;
192    }
193  }
194  return TRUE;
195}
196#endif //NDEBUG
197
198#ifndef NDEBUG
199BOOLEAN pReduceDebug(leftv res, leftv args)
200{
201  leftv u = args;
202  if ((u != NULL) && (u->Typ() == POLY_CMD))
203  {
204    poly g; number p = n_Init(3,currRing->cf);
205    omUpdateInfo();
206    Print("usedBytesBefore=%ld\n",om_Info.UsedBytes);
207    g = (poly) u->CopyD();
208    (void) pReduce(g,p,currRing);
209    p_Delete(&g,currRing);
210    omUpdateInfo();
211    Print("usedBytesAfter=%ld\n",om_Info.UsedBytes);
212    g = (poly) u->CopyD();
213    (void) pReduce(g,p,currRing);
214    n_Delete(&p,currRing->cf);
215    res->rtyp = POLY_CMD;
216    res->data = (char*) g;
217    return FALSE;
218  }
219  return TRUE;
220}
221#endif //NDEBUG
222
223void pReduce(ideal &I, const number p, const ring r)
224{
225  int k = idSize(I);
226  for (int i=0; i<k; i++)
227  {
228    if (I->m[i]!=NULL)
229    {
230      number c = p_GetCoeff(I->m[i],r);
231      if (!n_Equal(p,c,r->cf))
232        pReduce(I->m[i],p,r);
233    }
234  }
235  return;
236}
237
238
239/**
240 * reduces h initially with respect to g,
241 * returns false if h was initially reduced in the first place,
242 * returns true if reductions have taken place.
243 * assumes that h and g are in pReduced form and homogeneous in x of the same degree
244 */
245bool ppreduceInitially(poly* hStar, const poly g, const ring r)
246{
247  poly h = *hStar;
248  if (h==NULL || g==NULL)
249    return false;
250  p_Test(h,r);
251  p_Test(g,r);
252  poly hCache;
253  for (hCache=h; hCache; pIter(hCache))
254    if (p_LeadmonomDivisibleBy(g,hCache,r)) break;
255  if (hCache)
256  {
257    number gAlpha = p_GetCoeff(g,r);
258    poly hAlphaT = p_Init(r);
259    p_SetCoeff(hAlphaT,n_Copy(p_GetCoeff(hCache,r),r->cf),r);
260    p_SetExp(hAlphaT,1,p_GetExp(hCache,1,r)-p_GetExp(g,1,r),r);
261    for (int i=2; i<=r->N; i++)
262      p_SetExp(hAlphaT,i,0,r);
263    p_Setm(hAlphaT,r); p_Test(hAlphaT,r);
264    poly q1 = p_Mult_nn(h,gAlpha,r); p_Test(q1,r);
265    poly q2 = p_Mult_q(p_Copy(g,r),hAlphaT,r); p_Test(q2,r);
266    q2 = p_Neg(q2,r); p_Test(q2,r);
267    h = p_Add_q(q1,q2,r);
268    p_Test(h,r);
269    p_Test(g,r);
270    *hStar = h;
271    return true;
272  }
273  p_Test(h,r);
274  p_Test(g,r);
275  return false;
276}
277
278
279#ifndef NDEBUG
280BOOLEAN ppreduceInitially0(leftv res, leftv args)
281{
282  leftv u = args;
283  if ((u != NULL) && (u->Typ() == POLY_CMD))
284  {
285    leftv v = u->next;
286    if ((v != NULL) && (v->Typ() == POLY_CMD))
287    {
288      poly g,h;
289      omUpdateInfo();
290      Print("usedBytesBefore=%ld\n",om_Info.UsedBytes);
291      h = (poly) u->CopyD();
292      g = (poly) v->CopyD();
293      (void)ppreduceInitially(&h,g,currRing);
294      p_Delete(&h,currRing);
295      p_Delete(&g,currRing);
296      omUpdateInfo();
297      Print("usedBytesAfter=%ld\n",om_Info.UsedBytes);
298      h = (poly) u->CopyD();
299      g = (poly) v->CopyD();
300      (void)ppreduceInitially(&h,g,currRing);
301      p_Delete(&g,currRing);
302      res->rtyp = POLY_CMD;
303      res->data = (char*) h;
304      return FALSE;
305    }
306  }
307  return TRUE;
308}
309#endif //NDEBUG
310
311
312/***
313 * reduces I initially with respect to itself and with respect to p-t.
314 * also sorts the generators of I with respect to the leading monomials in descending order.
315 * assumes that I is generated by elements which are homogeneous in x of the same degree.
316 **/
317bool ppreduceInitially(ideal I, const number p, const ring r)
318{
319  int m=idSize(I),n=m; poly cache;
320  do
321  {
322    int j=0;
323    for (int i=1; i<n; i++)
324    {
325      if (p_LmCmp(I->m[i-1],I->m[i],r)<0)
326      {
327        cache=I->m[i-1];
328        I->m[i-1]=I->m[i];
329        I->m[i]=cache;
330        j = i;
331      }
332    }
333    n=j;
334  } while(n);
335  for (int i=0; i<m; i++)
336    pReduce(I->m[i],p,r);
337
338  /***
339   * the first pass. removing terms with the same monomials in x as lt(g_i) out of g_j for i<j
340   **/
341  for (int i=0; i<m-1; i++)
342    for (int j=i+1; j<m; j++)
343      if (ppreduceInitially(&I->m[j], I->m[i], r))
344        pReduce(I->m[j],p,r);
345
346  /***
347   * the second pass. removing terms divisible by lt(g_j) out of g_i for i<j
348   **/
349  for (int i=0; i<m-1; i++)
350    for (int j=i+1; j<m; j++)
351      if (ppreduceInitially(&I->m[i], I->m[j],r))
352        pReduce(I->m[i],p,r);
353
354  /***
355   * removes the elements of I which have been reduced to 0 in the previous two passes
356   **/
357  idSkipZeroes(I);
358  return false;
359}
360
361
362#ifndef NDEBUG
363BOOLEAN ppreduceInitially1(leftv res, leftv args)
364{
365  leftv u = args;
366  if ((u != NULL) && (u->Typ() == IDEAL_CMD))
367  {
368    leftv v = u->next;
369    if ((v != NULL) && (v->Typ() == NUMBER_CMD))
370    {
371      ideal I; number p;
372      omUpdateInfo();
373      Print("usedBytesBefore=%ld\n",om_Info.UsedBytes);
374      I = (ideal) u->CopyD();
375      p = (number) v->CopyD();
376      (void) ppreduceInitially(I,p,currRing);
377      id_Delete(&I,currRing);
378      n_Delete(&p,currRing->cf);
379      omUpdateInfo();
380      Print("usedBytesAfter=%ld\n",om_Info.UsedBytes);
381      I = (ideal) u->CopyD();
382      p = (number) v->CopyD();
383      (void) ppreduceInitially(I,p,currRing);
384      n_Delete(&p,currRing->cf);
385      res->rtyp = IDEAL_CMD;
386      res->data = (char*) I;
387      return FALSE;
388    }
389  }
390  return TRUE;
391}
392#endif //NDEBUG
393
394
395/***
396 * inserts g into I and reduces I with respect to itself and p-t
397 * returns the position in I in which g was inserted
398 * assumes that I was already sorted and initially reduced in the first place
399 **/
400int ppreduceInitially(ideal I, const number p, const poly g, const ring r)
401{
402  id_Test(I,r);
403  p_Test(g,r);
404  idInsertPoly(I,g);
405  int n=idSize(I);
406  int j;
407  for (j=n-1; j>0; j--)
408  {
409    if (p_LmCmp(I->m[j], I->m[j-1],r)>0)
410    {
411      poly cache = I->m[j];
412      I->m[j] = I->m[j-1];
413      I->m[j-1] = cache;
414    }
415    else
416      break;
417  }
418
419  /***
420   * the first pass. removing terms with the same monomials in x as lt(g_i) out of g_j for i<j
421   * removing terms with the same monomials in x as lt(g_j) out of g_k for j<k
422   **/
423  for (int i=0; i<j; i++)
424    if (ppreduceInitially(&I->m[j], I->m[i], r))
425      pReduce(I->m[j],p,r);
426  for (int k=j+1; k<n; k++)
427    if (ppreduceInitially(&I->m[k], I->m[j], r))
428    {
429      pReduce(I->m[k],p,r);
430      for (int l=j+1; l<k; l++)
431        if (ppreduceInitially(&I->m[k], I->m[l], r))
432          pReduce(I->m[k],p,r);
433    }
434
435  /***
436   * the second pass. removing terms divisible by lt(g_j) and lt(g_k) out of g_i for i<j<k
437   * removing terms divisible by lt(g_k) out of g_j for j<k
438   **/
439  for (int i=0; i<j; i++)
440    for (int k=j; k<n; k++)
441      if (ppreduceInitially(&I->m[i], I->m[k], r))
442        pReduce(I->m[i],p,r);
443  for (int k=j; k<n-1; k++)
444    for (int l=k+1; l<n; l++)
445      if (ppreduceInitially(&I->m[k], I->m[l], r))
446        pReduce(I->m[k],p,r);
447
448  /***
449   * removes the elements of I which have been reduced to 0 in the previous two passes
450   **/
451  idSkipZeroes(I);
452  id_Test(I,r);
453  return j;
454}
455
456
457#ifndef NDEBUG
458BOOLEAN ppreduceInitially2(leftv res, leftv args)
459{
460  leftv u = args;
461  if ((u != NULL) && (u->Typ() == IDEAL_CMD))
462  {
463    leftv v = u->next;
464    if ((v != NULL) && (v->Typ() == NUMBER_CMD))
465    {
466      leftv w = v->next;
467      if ((w != NULL) && (w->Typ() == POLY_CMD))
468      {
469        ideal I; number p; poly g;
470        omUpdateInfo();
471        Print("usedBytesBefore=%ld\n",om_Info.UsedBytes);
472        I = (ideal) u->CopyD();
473        p = (number) v->CopyD();
474        g = (poly) w->CopyD();
475        (void) ppreduceInitially(I,p,g,currRing);
476        id_Delete(&I,currRing);
477        n_Delete(&p,currRing->cf);
478        omUpdateInfo();
479        Print("usedBytesAfter=%ld\n",om_Info.UsedBytes);
480        I = (ideal) u->CopyD();
481        p = (number) v->CopyD();
482        g = (poly) w->CopyD();
483        (void) ppreduceInitially(I,p,g,currRing);
484        n_Delete(&p,currRing->cf);
485        res->rtyp = IDEAL_CMD;
486        res->data = (char*) I;
487        return FALSE;
488      }
489    }
490  }
491  return TRUE;
492}
493#endif //NDEBUG
494
495
496static poly ppNext(poly p, int l)
497{
498  poly q = p;
499  for (int i=0; i<l; i++)
500  {
501    if (q==NULL)
502      break;
503    pIter(q);
504  }
505  return q;
506}
507
508
509static void sortMarks(const ideal H, const ring r, std::vector<mark> &T)
510{
511  std::pair<int,int> pointerToTerm;
512  int k=T.size();
513  do
514  {
515    int j=0;
516    for (int i=1; i<k-1; i++)
517    {
518      int generatorA = T[i-1].first;
519      int termA = T[i-1].second;
520      int generatorB = T[i].first;
521      int termB = T[i].second;
522      if (p_LmCmp(ppNext(H->m[generatorA],termA),ppNext(H->m[generatorB],termB),r)<0)
523      {
524        mark cache=T[i-1];
525        T[i-1]=T[i];
526        T[i]=cache;
527        j = i;
528      }
529    }
530    k=j;
531  } while(k);
532  return;
533}
534
535
536static poly getTerm(const ideal H, const mark ab)
537{
538  int a = ab.first;
539  int b = ab.second;
540  return ppNext(H->m[a],b);
541}
542
543
544static void adjustMarks(std::vector<mark> &T, const int newEntry)
545{
546  for (unsigned i=0; i<T.size(); i++)
547  {
548    if (T[i].first>=newEntry)
549      T[i].first = T[i].first+1;
550  }
551  return;
552}
553
554
555static void cleanupMarks(const ideal H, std::vector<mark> &T)
556{
557  for (unsigned i=0; i<T.size();)
558  {
559    if (getTerm(H,T[i])==NULL)
560      T.erase(T.begin()+i);
561    else
562      i++;
563  }
564  return;
565}
566
567
568/***
569 * reduces H initially with respect to itself, with respect to p-t,
570 * and with respect to G.
571 * assumes that the generators of H are homogeneous in x of the same degree,
572 * assumes that the generators of G are homogeneous in x of lesser degree.
573 **/
574bool ppreduceInitially(ideal &H, const number p, const ideal G, const ring r)
575{
576  /***
577   * Step 1: reduce H initially with respect to itself and with respect to p-t
578   **/
579  if (ppreduceInitially(H,p,r)) return true;
580
581  /***
582   * Step 2: initialize an ideal I in which the reductions will take place-
583   *   along the reduction it will be enlarged with elements that will be discarded at the end
584   *   initialize a working list T which keeps track.
585   *   the working list T is a vector of pairs of integer.
586   *   if T contains a pair (i,j) then that means that in the i-th element of H
587   *   term j and subsequent terms need to be checked for reduction.
588   *   T is sorted by the ordering on the temrs the pairs correspond to.
589   **/
590  int m=idSize(H);
591  ideal I = idInit(m);
592  std::vector<mark> T;
593  for (int i=0; i<m; i++)
594  {
595    I->m[i]=H->m[i];
596    if (pNext(I->m[i])!=NULL)
597      T.push_back(std::make_pair<int,int>(i,1));
598  }
599
600  /***
601   * Step 3: as long as the working list is not empty, successively reduce terms in it
602   *   by adding suitable elements to I and reducing it initially with respect to itself
603   **/
604  int k=idSize(G);
605  while (T.size()>0)
606  {
607    sortMarks(I,r,T);
608    int i=0; for (; i<k; i++)
609      if (p_LeadmonomDivisibleBy(G->m[i],getTerm(I,T[0]),r)) break;
610    if (i<k)
611    {
612      poly g = p_One(r); poly h0 = getTerm(I,T[0]);
613      assume(h0!=NULL);
614      for (int j=2; j<=r->N; j++)
615        p_SetExp(g,j,p_GetExp(h0,j,r)-p_GetExp(G->m[i],j,r),r);
616      p_Setm(g,r);
617      g = p_Mult_q(g,p_Copy(G->m[i],r),r);
618      int newEntry = ppreduceInitially(I,p,g,r);
619      adjustMarks(T,newEntry);
620    }
621    else
622      T[0].second = T[0].second+1;
623    cleanupMarks(I,T);
624  }
625
626  /***
627   * Step 4: cleanup, delete all polynomials in I which have been added in Step 3
628   **/
629  k=idSize(I);
630  for (int i=0; i<k; i++)
631  {
632    for (int j=0; j<m; j++)
633    {
634      if (p_LeadmonomDivisibleBy(H->m[j],I->m[i],r))
635      {
636        I->m[i]=NULL;
637        break;
638      }
639    }
640  }
641  id_Delete(&I,r);
642  return false;
643}
644
645
646#ifndef NDEBUG
647BOOLEAN ppreduceInitially3(leftv res, leftv args)
648{
649  leftv u = args;
650  if ((u != NULL) && (u->Typ() == IDEAL_CMD))
651  {
652    leftv v = u->next;
653    if ((v != NULL) && (v->Typ() == NUMBER_CMD))
654    {
655      leftv w = v->next;
656      if ((w != NULL) && (w->Typ() == IDEAL_CMD))
657      {
658        ideal H,G; number p;
659        omUpdateInfo();
660        Print("usedBytesBefore=%ld\n",om_Info.UsedBytes);
661        H = (ideal) u->CopyD();
662        p = (number) v->CopyD();
663        G = (ideal) w->CopyD();
664        (void) ppreduceInitially(H,p,G,currRing);
665        n_Delete(&p,currRing->cf);
666        id_Delete(&G,currRing);
667        res->rtyp = IDEAL_CMD;
668        res->data = (char*) H;
669        return FALSE;
670      }
671    }
672  }
673  return TRUE;
674}
675#endif //NDEBUG
676
677/**
678 * reduces I initially with respect to itself.
679 * assumes that the generators of I are homogeneous in x and that p-t is in I.
680 */
681bool ppreduceInitially(ideal I, const ring r, const number p)
682{
683  assume(!n_IsUnit(p,r->cf));
684
685  /***
686   * Step 1: split up I into components of same degree in x
687   *  the lowest component should only contain p-t
688   **/
689  std::map<long,ideal> H; int n = idSize(I);
690  for (int i=0; i<n; i++)
691  {
692    I->m[i] = p_Cleardenom(I->m[i],r);
693    long d = 0;
694    for (int j=2; j<=r->N; j++)
695      d += p_GetExp(I->m[i],j,r);
696    std::map<long,ideal>::iterator it = H.find(d);
697    if (it != H.end())
698      idInsertPoly(it->second,I->m[i]);
699    else
700    {
701      std::pair<long,ideal> Hd(d,idInit(1));
702      Hd.second->m[0] = I->m[i];
703      H.insert(Hd);
704    }
705  }
706
707  std::map<long,ideal>::iterator it=H.begin();
708  ideal Hi = it->second;
709  idShallowDelete(&Hi);
710  it++;
711  Hi = it->second;
712
713  /***
714   * Step 2: reduce each component initially with respect to itself
715   *  and all lower components
716   **/
717  if (ppreduceInitially(Hi,p,r)) return true;
718  id_Test(Hi,r);
719  id_Test(I,r);
720
721  ideal G = idInit(n); int m=0;
722  ideal GG = (ideal) omAllocBin(sip_sideal_bin);
723  GG->nrows = 1; GG->rank = 1; GG->m=NULL;
724
725  for (it++; it!=H.end(); it++)
726  {
727    int l=idSize(Hi); int k=l; poly cache;
728    /**
729     * sorts Hi according to degree in t in descending order
730     * (lowest first, highest last)
731     */
732    do
733    {
734      int j=0;
735      for (int i=1; i<k; i++)
736      {
737        if (p_GetExp(Hi->m[i-1],1,r)<p_GetExp(Hi->m[i],1,r))
738        {
739          cache=Hi->m[i-1];
740          Hi->m[i-1]=Hi->m[i];
741          Hi->m[i]=cache;
742          j = i;
743        }
744      }
745      k=j;
746    } while(k);
747    int kG=n-m, kH=0;
748    for (int i=n-m-l; i<n; i++)
749    {
750      if (kG==n)
751      {
752        memcpy(&(G->m[i]),&(Hi->m[kH]),(n-i)*sizeof(poly));
753        break;
754      }
755      if (kH==l)
756        break;
757      if (p_GetExp(G->m[kG],1,r)>p_GetExp(Hi->m[kH],1,r))
758        G->m[i] = G->m[kG++];
759      else
760        G->m[i] = Hi->m[kH++];
761    }
762    m += l; IDELEMS(GG) = m; GG->m = &G->m[n-m];
763    id_Test(it->second,r);
764    id_Test(GG,r);
765    if (ppreduceInitially(it->second,p,GG,r)) return true;
766    id_Test(it->second,r);
767    id_Test(GG,r);
768    idShallowDelete(&Hi); Hi = it->second;
769  }
770  idShallowDelete(&Hi);
771
772  ptNormalize(I,p,r);
773  omFreeBin((ADDRESS)GG, sip_sideal_bin);
774  idShallowDelete(&G);
775  return false;
776}
777
778
779#ifndef NDEBUG
780BOOLEAN reduceInitiallyDebug(leftv res, leftv args)
781{
782  leftv u = args;
783  if ((u != NULL) && (u->Typ() == IDEAL_CMD))
784  {
785    leftv v = u->next;
786    if ((v != NULL) && (v->Typ() == NUMBER_CMD))
787    {
788      omUpdateInfo();
789      Print("usedBytesBefore=%ld\n",om_Info.UsedBytes);
790      ideal I = (ideal) u->CopyD();
791      number p = (number) v->Data();
792      (void) ppreduceInitially(I,currRing,p);
793      res->rtyp = IDEAL_CMD;
794      res->data = (char*) I;
795      return FALSE;
796    }
797  }
798  return TRUE;
799}
800#endif
801
802
803// BOOLEAN ppreduceInitially(leftv res, leftv args)
804// {
805//   leftv u = args;
806//   if ((u != NULL) && (u->Typ() == IDEAL_CMD))
807//   {
808//     ideal I = (ideal) u->CopyD();
809//     (void) ppreduceInitially(I,currRing);
810//     res->rtyp = IDEAL_CMD;
811//     res->data = (char*) I;
812//     return FALSE;
813//   }
814//   return TRUE;
815// }
Note: See TracBrowser for help on using the repository browser.