source: git/kernel/kstd2.cc @ 019649

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