source: git/kernel/kstd2.cc @ 6d09f28

spielwiese
Last change on this file since 6d09f28 was 6d09f28, checked in by Oliver Wienand <wienand@…>, 18 years ago
introduced def HAVE_VANGB --> compute GB in the ring of function rather than polynomials kstd2.cc: * moved ind2 and ind_fact_2 to kutil.h/cc * bba: enterpairs changes size of T in VANGB case kutil.h/cc: * ind2/ind_fact_2 see above * kCreateZeroPoly: create a zero polynomial with given exponent * createG0: creates the groebner basis of the vanishing polynomials git-svn-id: file:///usr/local/Singular/svn/trunk@9404 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 26.9 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/* $Id: kstd2.cc,v 1.18 2006-08-14 17:08:36 wienand Exp $ */
5/*
6*  ABSTRACT -  Kernel: alg. of Buchberger
7*/
8
9// #define PDEBUG 2
10// define to enable tailRings
11#define HAVE_TAIL_RING
12// define if no buckets should be used
13// #define NO_BUCKETS
14
15#include "mod2.h"
16#include "kutil.h"
17#include "structs.h"
18#include "omalloc.h"
19#include "polys.h"
20#include "ideals.h"
21#include "febase.h"
22#include "kstd1.h"
23#include "khstd.h"
24#include "kbuckets.h"
25//#include "cntrlc.h"
26#include "weight.h"
27#include "intvec.h"
28#ifdef HAVE_PLURAL
29#include "gring.h"
30#endif
31// #include "timer.h"
32
33// return -1 if no divisor is found
34//        number of first divisor, otherwise
35int kFindDivisibleByInT(const TSet &T, const unsigned long* sevT,
36                        const int tl, const LObject* L, const int start)
37{
38  unsigned long not_sev = ~L->sev;
39  int j = start;
40  poly p;
41  ring r;
42  L->GetLm(p, r);
43
44  pAssume(~not_sev == p_GetShortExpVector(p, r));
45
46  if (r == currRing)
47  {
48    loop
49    {
50      if (j > tl) return -1;
51#if defined(PDEBUG) || defined(PDIV_DEBUG)
52      if (p_LmShortDivisibleBy(T[j].p, sevT[j],
53                               p, not_sev, r))
54        return j;
55#else
56      if (!(sevT[j] & not_sev) &&
57          p_LmDivisibleBy(T[j].p, p, r))
58        return j;
59#endif
60      j++;
61    }
62  }
63  else
64  {
65    loop
66    {
67      if (j > tl) return -1;
68#if defined(PDEBUG) || defined(PDIV_DEBUG)
69      if (p_LmShortDivisibleBy(T[j].t_p, sevT[j],
70                               p, not_sev, r))
71        return j;
72#else
73      if (!(sevT[j] & not_sev) &&
74          p_LmDivisibleBy(T[j].t_p, p, r))
75        return j;
76#endif
77      j++;
78    }
79  }
80}
81
82// same as above, only with set S
83int kFindDivisibleByInS(const polyset &S, const unsigned long* sev, const int sl, LObject* L)
84{
85  unsigned long not_sev = ~L->sev;
86  poly p = L->GetLmCurrRing();
87  int j = 0;
88
89  pAssume(~not_sev == p_GetShortExpVector(p, currRing));
90  loop
91  {
92    if (j > sl) return -1;
93#if defined(PDEBUG) || defined(PDIV_DEBUG)
94    if (p_LmShortDivisibleBy(S[j], sev[j],
95                             p, not_sev, currRing))
96        return j;
97#else
98    if ( !(sev[j] & not_sev) &&
99         p_LmDivisibleBy(S[j], p, currRing))
100      return j;
101#endif
102    j++;
103  }
104}
105
106#ifdef HAVE_RING2TOM
107/*        Obsolute since changes to pLmDiv
108// return -1 if no divisor is found
109//        number of first divisor, otherwise
110int kRingFindDivisibleByInT(const TSet &T, const unsigned long* sevT,
111                        const int tl, const LObject* L, const int start)
112{
113  unsigned long not_sev = ~L->sev;
114  int j = start;
115  poly p;
116  ring r;
117  L->GetLm(p, r);
118
119  pAssume(~not_sev == p_GetShortExpVector(p, r));
120
121  {
122    loop
123    {
124      if (j > tl) return -1;
125#if defined(PDEBUG) || defined(PDIV_DEBUG)
126      if (p_LmRingShortDivisibleBy(T[j].p, sevT[j],
127                               p, not_sev, r))
128        return j;
129#else
130      if ( !(sevT[j] & not_sev) &&
131           p_LmRingDivisibleBy(T[j].p, p, r) )
132        return j;
133#endif
134      j++;
135    }
136  }
137  return -1;
138}
139
140// same as above, only with set S
141int kRingFindDivisibleByInS(const polyset &S, const unsigned long* sev, const int sl, LObject* L)
142{
143  unsigned long not_sev = ~L->sev;
144  poly p = L->GetLmCurrRing();
145  int j = 0;
146  //PrintS("FindDiv: p="); wrp(p); PrintLn();
147  pAssume(~not_sev == p_GetShortExpVector(p, currRing));
148  loop
149  {
150    //PrintS("FindDiv: S[j]="); wrp(S[j]); PrintLn();
151    if (j > sl) return -1;
152#if defined(PDEBUG) || defined(PDIV_DEBUG)
153    if (p_LmRingShortDivisibleBy(S[j], sev[j],
154                             p, not_sev, currRing))
155        return j;
156#else
157    if ( !(sev[j] & not_sev) &&
158         p_LmRingDivisibleBy(S[j], p, currRing) )
159      return j;
160#endif
161    j++;
162  }
163}
164*/
165
166/* now in kutil.cc
167long twoPow(long arg) {
168  long t = arg;
169  long result = 1;
170  while (t > 0) {
171    result = 2 * result;
172    t--;
173  }
174  return result;
175}
176*/
177
178long factorial(long arg)
179{
180   long tmp = 1; arg++;
181   for (int i = 2; i < arg; i++)
182   {
183     tmp *= i;
184   }
185   return tmp;
186}
187
188poly kFindZeroPoly(poly input_p, ring leadRing, ring tailRing)
189{
190  // m = currRing->ch
191
192  if (input_p == NULL) return NULL;
193
194  poly p = input_p;
195  poly zeroPoly = NULL;
196  long a = (long) pGetCoeff(p);
197
198  int k_ind2 = 0;
199  int a_ind2 = ind2(a);
200
201  long k = 1;
202  // of interest is only k_ind2, special routine for improvement ... TODO OLIVER
203  for (int i = 1; i <= leadRing->N; i++)
204  {
205    k_ind2 = k_ind2 + ind_fact_2(p_GetExp(p, i, leadRing));
206  }
207
208  a = (long) pGetCoeff(p);
209
210  number tmp1;
211  poly tmp2, tmp3;
212  poly lead_mult = p_ISet(1, tailRing);
213  if (leadRing->ch <= k_ind2 + a_ind2)
214  {
215    int too_much = k_ind2 + a_ind2 - leadRing->ch;
216    int s_exp;
217    zeroPoly = p_ISet(a, tailRing);
218    for (int i = 1; i <= leadRing->N; i++)
219    {
220      s_exp = p_GetExp(p, i,leadRing);
221      if (s_exp % 2 != 0)
222      {
223        s_exp = s_exp - 1;
224      }
225      while ( (0 < ind2(s_exp)) && (ind2(s_exp) <= too_much) )
226      {
227        too_much = too_much - ind2(s_exp);
228        s_exp = s_exp - 2;
229      }
230      p_SetExp(lead_mult, i, p_GetExp(p, i,leadRing) - s_exp, tailRing);
231      for (long j = 1; j <= s_exp; j++)
232      {
233        tmp1 = nInit(j);
234        tmp2 = p_ISet(1, tailRing);
235        p_SetExp(tmp2, i, 1, tailRing);
236        p_Setm(tmp2, tailRing);
237        if (nIsZero(tmp1))
238        { // should nowbe obsolet, test ! TODO OLIVER
239          zeroPoly = p_Mult_q(zeroPoly, tmp2, tailRing);
240        }
241        else
242        {
243          tmp3 = p_NSet(nCopy(tmp1), tailRing);
244          zeroPoly = p_Mult_q(zeroPoly, p_Add_q(tmp3, tmp2, tailRing), tailRing);
245        }
246      }
247    }
248    p_Setm(lead_mult, tailRing);
249    zeroPoly = p_Mult_mm(zeroPoly, lead_mult, tailRing);
250    tmp2 = p_NSet(nCopy(pGetCoeff(zeroPoly)), leadRing);
251    for (int i = 1; i <= leadRing->N; i++) {
252      pSetExp(tmp2, i, p_GetExp(zeroPoly, i, tailRing));
253    }
254    p_Setm(tmp2, leadRing);
255    zeroPoly = p_LmDeleteAndNext(zeroPoly, tailRing);
256    pNext(tmp2) = zeroPoly;
257    return tmp2;
258  }
259/*  long alpha_k = twoPow(leadRing->ch - k_ind2);
260  if (1 == 0 && alpha_k <= a) {  // Temporarly disabled, reducing coefficients not compatible with std TODO Oliver
261    zeroPoly = p_ISet((a / alpha_k)*alpha_k, tailRing);
262    for (int i = 1; i <= leadRing->N; i++) {
263      for (long j = 1; j <= p_GetExp(p, i, leadRing); j++) {
264        tmp1 = nInit(j);
265        tmp2 = p_ISet(1, tailRing);
266        p_SetExp(tmp2, i, 1, tailRing);
267        p_Setm(tmp2, tailRing);
268        if (nIsZero(tmp1)) {
269          zeroPoly = p_Mult_q(zeroPoly, tmp2, tailRing);
270        }
271        else {
272          tmp3 = p_ISet((long) tmp1, tailRing);
273          zeroPoly = p_Mult_q(zeroPoly, p_Add_q(tmp2, tmp3, tailRing), tailRing);
274        }
275      }
276    }
277    tmp2 = p_ISet((long) pGetCoeff(zeroPoly), leadRing);
278    for (int i = 1; i <= leadRing->N; i++) {
279      pSetExp(tmp2, i, p_GetExp(zeroPoly, i, tailRing));
280    }
281    p_Setm(tmp2, leadRing);
282    zeroPoly = p_LmDeleteAndNext(zeroPoly, tailRing);
283    pNext(tmp2) = zeroPoly;
284    return tmp2;
285  } */
286  return NULL;
287}
288
289poly kFindDivisibleByZeroPoly(LObject* h)
290{
291  return kFindZeroPoly(h->GetLmCurrRing(), currRing, h->tailRing);
292}
293
294#ifdef HAVE_RING2TOM
295/*2
296*  reduction procedure for the ring Z/2^m
297*/
298int redRing2toM (LObject* h,kStrategy strat)
299{
300  if (h->p == NULL && h->t_p == NULL) return 0; // spoly is zero (can only occure with zero divisors)
301
302//  if (strat->tl<0) return 1;
303  int at,d,i;
304  int j = 0;
305  int pass = 0;
306  poly zeroPoly;
307
308// TODO warum SetpFDeg notwendig?
309  h->SetpFDeg();
310  assume(h->pFDeg() == h->FDeg);
311//  if (h->pFDeg() != h->FDeg)
312//  {
313//    Print("h->pFDeg()=%d =!= h->FDeg=%d\n", h->pFDeg(), h->FDeg);
314//  }
315  long reddeg = h->GetpFDeg();
316
317  h->SetShortExpVector();
318  loop
319  {
320#ifdef HAVE_VANGB
321    zeroPoly = kFindDivisibleByZeroPoly(h);
322    if (zeroPoly != NULL)
323    {
324      if (TEST_OPT_PROT)
325      {
326        PrintS("z");
327      }
328#ifdef KDEBUG
329      if (TEST_OPT_DEBUG)
330      {
331        PrintS("zero red: ");
332      }
333#endif
334      LObject tmp_h(zeroPoly, currRing, strat->tailRing);
335      tmp_h.SetShortExpVector();
336      strat->initEcart(&tmp_h);
337      tmp_h.sev = pGetShortExpVector(tmp_h.p);
338      tmp_h.SetpFDeg();
339
340      enterT(tmp_h, strat, strat->tl + 1);
341      j = strat->tl;
342    }
343    else
344#endif
345    {
346      j = kFindDivisibleByInT(strat->T, strat->sevT, strat->tl, h);
347      if (j < 0) return 1;
348#ifdef KDEBUG
349      if (TEST_OPT_DEBUG)
350      {
351        PrintS("T red:");
352      }
353#endif
354    }
355#ifdef KDEBUG
356    if (TEST_OPT_DEBUG)
357    {
358      h->wrp();
359      PrintS(" with ");
360      strat->T[j].wrp();
361    }
362#endif
363
364    ksReducePoly(h, &(strat->T[j]), NULL, NULL, strat);
365    if (zeroPoly != NULL) {
366      // TODO Free memory of zeropoly and last element of L
367      strat->tl--;
368    }
369
370#ifdef KDEBUG
371    if (TEST_OPT_DEBUG)
372    {
373      PrintS("\nto ");
374      h->wrp();
375      PrintLn();
376    }
377#endif
378
379    if (h->GetLmTailRing() == NULL)
380    {
381      if (h->lcm!=NULL) pLmFree(h->lcm);
382#ifdef KDEBUG
383      h->lcm=NULL;
384#endif
385      return 0;
386    }
387    h->SetShortExpVector();
388    d = h->SetpFDeg();
389    /*- try to reduce the s-polynomial -*/
390    pass++;
391    if (!K_TEST_OPT_REDTHROUGH &&
392        (strat->Ll >= 0) && ((d > reddeg) || (pass > strat->LazyPass)))
393    {
394      h->SetLmCurrRing();
395      at = strat->posInL(strat->L,strat->Ll,h,strat);
396      if (at <= strat->Ll)
397      {
398#if 0
399        if (kRingFindDivisibleByInS(strat->S, strat->sevS, strat->sl, h) < 0)
400          return 1;
401#endif
402#ifdef KDEBUG
403        if (TEST_OPT_DEBUG) Print(" ->L[%d]\n",at);
404#endif
405        enterL(&strat->L,&strat->Ll,&strat->Lmax,*h,at);     // NOT RING CHECKED OLIVER
406        h->Clear();
407        return -1;
408      }
409    }
410    else if ((TEST_OPT_PROT) && (strat->Ll < 0) && (d != reddeg))
411    {
412      Print(".%d",d);mflush();
413      reddeg = d;
414    }
415  }
416}
417#endif
418#endif
419
420/*2
421*  reduction procedure for the homogeneous case
422*  and the case of a degree-ordering
423*/
424int redHomog (LObject* h,kStrategy strat)
425{
426//  if (strat->tl<0) return 1;
427#ifdef KDEBUG
428  if (TEST_OPT_DEBUG)
429  {
430    PrintS("red:");
431    h->wrp();
432    PrintS(" ");
433  }
434#endif
435  int j;
436  int pass = 0;
437  loop
438  {
439    // find a poly with which we can reduce
440    h->SetShortExpVector();
441    j = kFindDivisibleByInT(strat->T, strat->sevT, strat->tl, h);
442    if (j < 0)
443    {
444      h->SetpFDeg();
445      return 1;
446    }
447
448    // now we found one which is divisible -- reduce it
449    ksReducePoly(h, &(strat->T[j]), NULL, NULL, strat);
450
451#ifdef KDEBUG
452    if (TEST_OPT_DEBUG)
453    {
454      Print("\nto ", h->t_p);
455      h->wrp();
456      PrintLn();
457    }
458#endif
459    if (h->GetLmTailRing() == NULL)
460    {
461      if (h->lcm!=NULL) pLmFree(h->lcm);
462#ifdef KDEBUG
463      h->lcm=NULL;
464#endif
465      return 0;
466    }
467    if (!K_TEST_OPT_REDTHROUGH &&
468        (strat->Ll >= 0) && (pass > strat->LazyPass))
469    {
470      h->SetLmCurrRing();
471      int at = strat->posInL(strat->L,strat->Ll,h,strat);
472      if (at <= strat->Ll)
473      {
474#ifdef KDEBUG
475        if (TEST_OPT_DEBUG) Print(" ->L[%d]\n",at);
476#endif
477        enterL(&strat->L,&strat->Ll,&strat->Lmax,*h,at);
478        h->Clear();
479        return -1;
480      }
481    }
482  }
483}
484
485/*2
486*  reduction procedure for the inhomogeneous case
487*  and not a degree-ordering
488*/
489int redLazy (LObject* h,kStrategy strat)
490{
491  if (strat->tl<0) return 1;
492  int at,d,i;
493  int j = 0;
494  int pass = 0;
495  assume(h->pFDeg() == h->FDeg);
496  long reddeg = h->GetpFDeg();
497
498  h->SetShortExpVector();
499  loop
500  {
501    j = kFindDivisibleByInT(strat->T, strat->sevT, strat->tl, h);
502    if (j < 0) return 1;
503
504#ifdef KDEBUG
505    if (TEST_OPT_DEBUG)
506    {
507      PrintS("red:");
508      h->wrp();
509      PrintS(" with ");
510      strat->T[j].wrp();
511    }
512#endif
513
514    ksReducePoly(h, &(strat->T[j]), NULL, NULL, strat);
515
516#ifdef KDEBUG
517    if (TEST_OPT_DEBUG)
518    {
519      PrintS("\nto ");
520      h->wrp();
521      PrintLn();
522    }
523#endif
524
525    if (h->GetLmTailRing() == NULL)
526    {
527      if (h->lcm!=NULL) pLmFree(h->lcm);
528#ifdef KDEBUG
529      h->lcm=NULL;
530#endif
531      return 0;
532    }
533    h->SetShortExpVector();
534    d = h->SetpFDeg();
535    /*- try to reduce the s-polynomial -*/
536    pass++;
537    if (!K_TEST_OPT_REDTHROUGH &&
538        (strat->Ll >= 0) && ((d > reddeg) || (pass > strat->LazyPass)))
539    {
540      h->SetLmCurrRing();
541      at = strat->posInL(strat->L,strat->Ll,h,strat);
542      if (at <= strat->Ll)
543      {
544#if 0
545        if (kFindDivisibleByInS(strat->S, strat->sevS, strat->sl, h) < 0)
546          return 1;
547#endif
548#ifdef KDEBUG
549        if (TEST_OPT_DEBUG) Print(" ->L[%d]\n",at);
550#endif
551        enterL(&strat->L,&strat->Ll,&strat->Lmax,*h,at);
552        h->Clear();
553        return -1;
554      }
555    }
556    else if ((TEST_OPT_PROT) && (strat->Ll < 0) && (d != reddeg))
557    {
558      Print(".%d",d);mflush();
559      reddeg = d;
560    }
561  }
562}
563
564/*2
565*  reduction procedure for the sugar-strategy (honey)
566* reduces h with elements from T choosing first possible
567* element in T with respect to the given ecart
568*/
569int redHoney (LObject* h, kStrategy strat)
570{
571  if (strat->tl<0) return 1;
572  //if (h->GetLmTailRing()==NULL) return 0; // HS: SHOULD NOT BE NEEDED!
573  assume(h->FDeg == h->pFDeg());
574
575  poly h_p;
576  int i,j,at,pass,ei, ii, h_d;
577  unsigned long not_sev;
578  long reddeg,d;
579
580  pass = j = 0;
581  d = reddeg = h->GetpFDeg() + h->ecart;
582  h->SetShortExpVector();
583  h_p = h->GetLmTailRing();
584  not_sev = ~ h->sev;
585  loop
586  {
587    j = kFindDivisibleByInT(strat->T, strat->sevT, strat->tl, h);
588    if (j < 0) return 1;
589
590    ei = strat->T[j].ecart;
591    ii = j;
592    /*
593     * the polynomial to reduce with (up to the moment) is;
594     * pi with ecart ei
595     */
596    i = j;
597    loop
598    {
599      /*- takes the first possible with respect to ecart -*/
600      i++;
601      if (i > strat->tl)
602        break;
603      if (ei <= h->ecart)
604        break;
605      if ((strat->T[i].ecart < ei) &&
606          p_LmShortDivisibleBy(strat->T[i].GetLmTailRing(), strat->sevT[i],
607                               h_p, not_sev, strat->tailRing))
608      {
609        /*
610         * the polynomial to reduce with is now;
611         */
612        ei = strat->T[i].ecart;
613        ii = i;
614      }
615    }
616
617    /*
618     * end of search: have to reduce with pi
619     */
620    if (!K_TEST_OPT_REDTHROUGH && (pass!=0) && (ei > h->ecart))
621    {
622      h->SetLmCurrRing();
623      /*
624       * It is not possible to reduce h with smaller ecart;
625       * if possible h goes to the lazy-set L,i.e
626       * if its position in L would be not the last one
627       */
628      if (strat->Ll >= 0) /* L is not empty */
629      {
630        at = strat->posInL(strat->L,strat->Ll,h,strat);
631        if(at <= strat->Ll)
632          /*- h will not become the next element to reduce -*/
633        {
634          enterL(&strat->L,&strat->Ll,&strat->Lmax,*h,at);
635#ifdef KDEBUG
636          if (TEST_OPT_DEBUG) Print(" ecart too big: -> L%d\n",at);
637#endif
638          h->Clear();
639          return -1;
640        }
641      }
642    }
643#ifdef KDEBUG
644    if (TEST_OPT_DEBUG)
645    {
646      PrintS("red:");
647      h->wrp();
648      PrintS(" with ");
649      strat->T[ii].wrp();
650    }
651#endif
652    assume(strat->fromT == FALSE);
653
654    ksReducePoly(h, &(strat->T[ii]), NULL, NULL, strat);
655
656#ifdef KDEBUG
657    if (TEST_OPT_DEBUG)
658    {
659      PrintS("\nto ");
660      h->wrp();
661      PrintLn();
662    }
663#endif
664
665    h_p = h->GetLmTailRing();
666    if (h_p == NULL)
667    {
668      if (h->lcm!=NULL) pLmFree(h->lcm);
669#ifdef KDEBUG
670      h->lcm=NULL;
671#endif
672      return 0;
673    }
674    h->SetShortExpVector();
675    not_sev = ~ h->sev;
676    h_d = h->SetpFDeg();
677    /* compute the ecart */
678    if (ei <= h->ecart)
679      h->ecart = d-h_d;
680    else
681      h->ecart = d-h_d+ei-h->ecart;
682    /*
683     * try to reduce the s-polynomial h
684     *test first whether h should go to the lazyset L
685     *-if the degree jumps
686     *-if the number of pre-defined reductions jumps
687     */
688    pass++;
689    d = h_d + h->ecart;
690    if (!K_TEST_OPT_REDTHROUGH && (strat->Ll >= 0) && ((d > reddeg) || (pass > strat->LazyPass)))
691    {
692      h->SetLmCurrRing();
693      at = strat->posInL(strat->L,strat->Ll,h,strat);
694      if (at <= strat->Ll)
695      {
696        if (kFindDivisibleByInS(strat->S, strat->sevS, strat->sl, h) < 0)
697          return 1;
698        enterL(&strat->L,&strat->Ll,&strat->Lmax,*h,at);
699#ifdef KDEBUG
700        if (TEST_OPT_DEBUG)
701          Print(" degree jumped: -> L%d\n",at);
702#endif
703        h->Clear();
704        return -1;
705      }
706    }
707    else if (TEST_OPT_PROT && (strat->Ll < 0) && (d > reddeg))
708    {
709      reddeg = d;
710      Print(".%d",d); mflush();
711    }
712  }
713}
714/*2
715*  reduction procedure for the normal form
716*/
717
718poly redNF (poly h,kStrategy strat)
719{
720  if (h==NULL) return NULL;
721  int j;
722
723  if (0 > strat->sl)
724  {
725    return h;
726  }
727  LObject P(h);
728  P.SetShortExpVector();
729  P.bucket = kBucketCreate(currRing);
730  kBucketInit(P.bucket,P.p,pLength(P.p));
731  kbTest(P.bucket);
732  loop
733  {
734/* Obsolete since change in pLmDiv
735#ifdef HAVE_RING2TOM
736    if (currRing->cring == 1) {
737      j=kRingFindDivisibleByInS(strat->S,strat->sevS,strat->sl,&P);
738    }
739    else
740#endif
741*/
742      j=kFindDivisibleByInS(strat->S,strat->sevS,strat->sl,&P);
743    if (j>=0)
744    {
745      nNormalize(pGetCoeff(P.p));
746#ifdef KDEBUG
747      if (TEST_OPT_DEBUG)
748      {
749        PrintS("red:");
750        wrp(h);
751        PrintS(" with ");
752        wrp(strat->S[j]);
753      }
754#endif
755#ifdef HAVE_PLURAL
756      if (rIsPluralRing(currRing))
757      {
758        number coef;
759        nc_kBucketPolyRed(P.bucket,strat->S[j],&coef);
760        nDelete(&coef);
761      }
762      else
763#endif
764      {
765        number coef;
766        coef=kBucketPolyRed(P.bucket,strat->S[j],pLength(strat->S[j]),strat->kNoether);
767        nDelete(&coef);
768      }
769      h = kBucketGetLm(P.bucket);   // FRAGE OLIVER
770      if (h==NULL)
771      {
772        kBucketDestroy(&P.bucket);
773        return NULL;
774      }
775      kbTest(P.bucket);
776      P.p=h;
777      P.t_p=NULL;
778      P.SetShortExpVector();
779#ifdef KDEBUG
780      if (TEST_OPT_DEBUG)
781      {
782        PrintS("\nto:");
783        wrp(h);
784        PrintLn();
785      }
786#endif
787    }
788    else
789    {
790      P.p=kBucketClear(P.bucket);
791      kBucketDestroy(&P.bucket);
792      pNormalize(P.p);
793      return P.p;
794    }
795  }
796}
797
798#ifdef KDEBUG
799static int bba_count = 0;
800#endif
801
802ideal bba (ideal F, ideal Q,intvec *w,intvec *hilb,kStrategy strat)
803{
804#ifdef KDEBUG
805  bba_count++;
806  int loop_count = 0;
807#endif
808  om_Opts.MinTrack = 5;
809  int   srmax,lrmax, red_result = 1;
810  int   olddeg,reduc;
811  int hilbeledeg=1,hilbcount=0,minimcnt=0;
812  BOOLEAN withT = FALSE;
813
814  initBuchMoraCrit(strat); /*set Gebauer, honey, sugarCrit*/
815  initBuchMoraPos(strat);
816  initHilbCrit(F,Q,&hilb,strat);
817  initBba(F,strat);
818  /*set enterS, spSpolyShort, reduce, red, initEcart, initEcartPair*/
819  /*Shdl=*/initBuchMora(F, Q,strat);
820  if (strat->minim>0) strat->M=idInit(IDELEMS(F),F->rank);
821  srmax = strat->sl;
822  reduc = olddeg = lrmax = 0;
823
824#ifndef NO_BUCKETS
825  if (!TEST_OPT_NOT_BUCKETS)
826    strat->use_buckets = 1;
827#endif
828
829  // redtailBBa against T for inhomogenous input
830  if (!K_TEST_OPT_OLDSTD)
831    withT = ! strat->homog;
832
833  // strat->posInT = posInT_pLength;
834  kTest_TS(strat);
835
836#ifdef HAVE_TAIL_RING
837  kStratInitChangeTailRing(strat);
838#endif
839
840  /* compute------------------------------------------------------- */
841  while (strat->Ll >= 0)
842  {
843    if (strat->Ll > lrmax) lrmax =strat->Ll;/*stat.*/
844#ifdef KDEBUG
845    loop_count++;
846#ifdef HAVE_RING2TOM
847    if (TEST_OPT_DEBUG) PrintS("--- next step ---\n");
848#endif
849    if (TEST_OPT_DEBUG) messageSets(strat);
850#endif
851    if (strat->Ll== 0) strat->interpt=TRUE;
852    if (TEST_OPT_DEGBOUND
853        && ((strat->honey && (strat->L[strat->Ll].ecart+pFDeg(strat->L[strat->Ll].p,currRing)>Kstd1_deg))
854            || ((!strat->honey) && (pFDeg(strat->L[strat->Ll].p,currRing)>Kstd1_deg))))
855    {
856      /*
857       *stops computation if
858       * 24 IN test and the degree +ecart of L[strat->Ll] is bigger then
859       *a predefined number Kstd1_deg
860       */
861      while ((strat->Ll >= 0)
862        && (strat->L[strat->Ll].p1!=NULL) && (strat->L[strat->Ll].p2!=NULL)
863        && ((strat->honey && (strat->L[strat->Ll].ecart+pFDeg(strat->L[strat->Ll].p,currRing)>Kstd1_deg))
864            || ((!strat->honey) && (pFDeg(strat->L[strat->Ll].p,currRing)>Kstd1_deg)))
865        )
866        deleteInL(strat->L,&strat->Ll,strat->Ll,strat);
867      if (strat->Ll<0) break;
868      else strat->noClearS=TRUE;
869    }
870    /* picks the last element from the lazyset L */
871    strat->P = strat->L[strat->Ll];
872    strat->Ll--;
873
874    if (pNext(strat->P.p) == strat->tail)
875    {
876      // deletes the short spoly
877      pLmFree(strat->P.p);
878      strat->P.p = NULL;
879      poly m1 = NULL, m2 = NULL;
880
881      // check that spoly creation is ok
882      while (strat->tailRing != currRing &&
883             !kCheckSpolyCreation(&(strat->P), strat, m1, m2))
884      {
885        assume(m1 == NULL && m2 == NULL);
886        // if not, change to a ring where exponents are at least
887        // large enough
888        kStratChangeTailRing(strat);
889      }
890      // create the real one
891      ksCreateSpoly(&(strat->P), NULL, strat->use_buckets,
892                    strat->tailRing, m1, m2, strat->R);
893    }
894    else if (strat->P.p1 == NULL)
895    {
896      if (strat->minim > 0)
897        strat->P.p2=p_Copy(strat->P.p, currRing, strat->tailRing);
898      // for input polys, prepare reduction
899      strat->P.PrepareRed(strat->use_buckets);
900    }
901
902#ifdef HAVE_RING2TOM
903    if (strat->P.p == NULL && strat->P.t_p == NULL) {
904      red_result = 0;
905    }
906    else
907#endif
908    {
909      if (TEST_OPT_PROT)
910        message((strat->honey ? strat->P.ecart : 0) + strat->P.pFDeg(),
911                &olddeg,&reduc,strat, red_result);
912
913      /* reduction of the element choosen from L */
914      red_result = strat->red(&strat->P,strat);
915    }
916
917    // reduction to non-zero new poly
918    if (red_result == 1)
919    {
920      /* statistic */
921      if (TEST_OPT_PROT) PrintS("s");
922
923      // get the polynomial (canonicalize bucket, make sure P.p is set)
924      strat->P.GetP(strat->lmBin);
925
926      int pos=posInS(strat,strat->sl,strat->P.p,strat->P.ecart);
927
928      // reduce the tail and normalize poly
929      if (TEST_OPT_INTSTRATEGY)
930      {
931        strat->P.pCleardenom();
932        if ((TEST_OPT_REDSB)||(TEST_OPT_REDTAIL))
933        {
934          strat->P.p = redtailBba(&(strat->P),pos-1,strat, withT);
935          strat->P.pCleardenom();
936        }
937      }
938      else
939      {
940        strat->P.pNorm();
941        if ((TEST_OPT_REDSB)||(TEST_OPT_REDTAIL))
942          strat->P.p = redtailBba(&(strat->P),pos-1,strat, withT);
943      }
944
945#ifdef KDEBUG
946      if (TEST_OPT_DEBUG){PrintS("new s:");strat->P.wrp();PrintLn();}
947#endif
948
949      // min_std stuff
950      if ((strat->P.p1==NULL) && (strat->minim>0))
951      {
952        if (strat->minim==1)
953        {
954          strat->M->m[minimcnt]=p_Copy(strat->P.p,currRing,strat->tailRing);
955          p_Delete(&strat->P.p2, currRing, strat->tailRing);
956        }
957        else
958        {
959          strat->M->m[minimcnt]=strat->P.p2;
960          strat->P.p2=NULL;
961        }
962        if (strat->tailRing!=currRing && pNext(strat->M->m[minimcnt])!=NULL)
963          pNext(strat->M->m[minimcnt])
964            = strat->p_shallow_copy_delete(pNext(strat->M->m[minimcnt]),
965                                           strat->tailRing, currRing,
966                                           currRing->PolyBin);
967        minimcnt++;
968      }
969
970      // enter into S, L, and T
971      enterT(strat->P, strat);
972#ifdef HAVE_RING2TOM
973#ifdef HAVE_VANGB
974      int at_R = strat->tl;
975#endif
976      if (currRing->cring == 1)
977        superenterpairs(strat->P.p,strat->sl,strat->P.ecart,pos,strat, strat->tl);
978      else
979#endif
980        enterpairs(strat->P.p,strat->sl,strat->P.ecart,pos,strat, strat->tl);
981      // posInS only depends on the leading term
982#ifdef HAVE_VANGB
983      strat->enterS(strat->P, pos, strat, at_R);
984#else
985      strat->enterS(strat->P, pos, strat, strat->tl);
986#endif
987      if (hilb!=NULL) khCheck(Q,w,hilb,hilbeledeg,hilbcount,strat);
988//      Print("[%d]",hilbeledeg);
989      if (strat->P.lcm!=NULL) pLmFree(strat->P.lcm);
990      if (strat->sl>srmax) srmax = strat->sl;
991    }
992    else if (strat->P.p1 == NULL && strat->minim > 0)
993    {
994      p_Delete(&strat->P.p2, currRing, strat->tailRing);
995    }
996#ifdef KDEBUG
997    memset(&(strat->P), 0, sizeof(strat->P));
998#endif
999    kTest_TS(strat);
1000  }
1001#ifdef KDEBUG
1002  if (TEST_OPT_DEBUG) messageSets(strat);
1003#endif
1004  /* complete reduction of the standard basis--------- */
1005  if (TEST_OPT_REDSB) completeReduce(strat);
1006  /* release temp data-------------------------------- */
1007  exitBuchMora(strat);
1008  if (TEST_OPT_WEIGHTM)
1009  {
1010    pRestoreDegProcs(pFDegOld, pLDegOld);
1011    if (ecartWeights)
1012    {
1013      omFreeSize((ADDRESS)ecartWeights,(pVariables+1)*sizeof(short));
1014      ecartWeights=NULL;
1015    }
1016  }
1017  if (TEST_OPT_PROT) messageStat(srmax,lrmax,hilbcount,strat);
1018  if (Q!=NULL) updateResult(strat->Shdl,Q,strat);
1019  return (strat->Shdl);
1020}
1021
1022poly kNF2 (ideal F,ideal Q,poly q,kStrategy strat, int lazyReduce)
1023{
1024  poly   p;
1025  int   i;
1026
1027  if ((idIs0(F))&&(Q==NULL))
1028    return pCopy(q); /*F=0*/
1029  strat->ak = idRankFreeModule(F);
1030  /*- creating temp data structures------------------- -*/
1031  BITSET save_test=test;
1032  test|=Sy_bit(OPT_REDTAIL);
1033  initBuchMoraCrit(strat);
1034  strat->initEcart = initEcartBBA;
1035  strat->enterS = enterSBba;
1036  /*- set S -*/
1037  strat->sl = -1;
1038  /*- init local data struct.---------------------------------------- -*/
1039  /*Shdl=*/initS(F,Q,strat);
1040  /*- compute------------------------------------------------------- -*/
1041  //if ((TEST_OPT_INTSTRATEGY)&&(lazyReduce==0))
1042  {
1043    for (i=strat->sl;i>=0;i--)
1044      pNorm(strat->S[i]);
1045  }
1046  kTest(strat);
1047  if (TEST_OPT_PROT) { PrintS("r"); mflush(); }
1048  p = redNF(pCopy(q),strat);
1049  if ((p!=NULL)&&(lazyReduce==0))
1050  {
1051    if (TEST_OPT_PROT) { PrintS("t"); mflush(); }
1052    p = redtailBba(p,strat->sl,strat);
1053  }
1054  /*- release temp data------------------------------- -*/
1055  omfree(strat->sevS);
1056  omfree(strat->ecartS);
1057  omfree(strat->T);
1058  omfree(strat->sevT);
1059  omfree(strat->R);
1060  omfree(strat->S_2_R);
1061  omfree(strat->L);
1062  omfree(strat->B);
1063  omfree(strat->fromQ);
1064  idDelete(&strat->Shdl);
1065  test=save_test;
1066  if (TEST_OPT_PROT) PrintLn();
1067  return p;
1068}
1069
1070ideal kNF2 (ideal F,ideal Q,ideal q,kStrategy strat, int lazyReduce)
1071{
1072  poly   p;
1073  int   i;
1074  ideal res;
1075
1076  if (idIs0(q))
1077    return idInit(IDELEMS(q),q->rank);
1078  if ((idIs0(F))&&(Q==NULL))
1079    return idCopy(q); /*F=0*/
1080  strat->ak = idRankFreeModule(F);
1081  /*- creating temp data structures------------------- -*/
1082  BITSET save_test=test;
1083  test|=Sy_bit(OPT_REDTAIL);
1084  initBuchMoraCrit(strat);
1085  strat->initEcart = initEcartBBA;
1086  strat->enterS = enterSBba;
1087  /*- set S -*/
1088  strat->sl = -1;
1089  /*- init local data struct.---------------------------------------- -*/
1090  /*Shdl=*/initS(F,Q,strat);
1091  /*- compute------------------------------------------------------- -*/
1092  res=idInit(IDELEMS(q),q->rank);
1093  //if ((TEST_OPT_INTSTRATEGY)&&(lazyReduce==0))
1094  {
1095    for (i=strat->sl;i>=0;i--)
1096      pNorm(strat->S[i]);
1097  }
1098  for (i=IDELEMS(q)-1; i>=0; i--)
1099  {
1100    if (q->m[i]!=NULL)
1101    {
1102      if (TEST_OPT_PROT) { PrintS("r");mflush(); }
1103      p = redNF(pCopy(q->m[i]),strat);
1104      if ((p!=NULL)&&(lazyReduce==0))
1105      {
1106        if (TEST_OPT_PROT) { PrintS("t"); mflush(); }
1107        p = redtailBba(p,strat->sl,strat);
1108      }
1109      res->m[i]=p;
1110    }
1111    //else
1112    //  res->m[i]=NULL;
1113  }
1114  /*- release temp data------------------------------- -*/
1115  omfree(strat->sevS);
1116  omfree(strat->ecartS);
1117  omfree(strat->T);
1118  omfree(strat->sevT);
1119  omfree(strat->R);
1120  omfree(strat->S_2_R);
1121  omfree(strat->L);
1122  omfree(strat->B);
1123  omfree(strat->fromQ);
1124  idDelete(&strat->Shdl);
1125  test=save_test;
1126  if (TEST_OPT_PROT) PrintLn();
1127  return res;
1128}
Note: See TracBrowser for help on using the repository browser.