source: git/kernel/kstd2.cc @ a725dae

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