source: git/kernel/kstd2.cc @ 977f94

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