source: git/kernel/kstd2.cc @ 939847

spielwiese
Last change on this file since 939847 was 939847, checked in by Hans Schönemann <hannes@…>, 18 years ago
*hannes: fix degbound error, deform_l.tst git-svn-id: file:///usr/local/Singular/svn/trunk@9014 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 26.9 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/* $Id: kstd2.cc,v 1.13 2006-03-10 12:53:19 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_ISet((long) 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_ISet((long) 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
890      deleteInL(strat->L,&strat->Ll,strat->Ll,strat);
891      break;
892    }
893    /* picks the last element from the lazyset L */
894    strat->P = strat->L[strat->Ll];
895    strat->Ll--;
896
897    if (pNext(strat->P.p) == strat->tail)
898    {
899      // deletes the short spoly
900      pLmFree(strat->P.p);
901      strat->P.p = NULL;
902      poly m1 = NULL, m2 = NULL;
903
904      // check that spoly creation is ok
905      while (strat->tailRing != currRing &&
906             !kCheckSpolyCreation(&(strat->P), strat, m1, m2))
907      {
908        assume(m1 == NULL && m2 == NULL);
909        // if not, change to a ring where exponents are at least
910        // large enough
911        kStratChangeTailRing(strat);
912      }
913      // create the real one
914      ksCreateSpoly(&(strat->P), NULL, strat->use_buckets,
915                    strat->tailRing, m1, m2, strat->R);
916    }
917    else if (strat->P.p1 == NULL)
918    {
919      if (strat->minim > 0)
920        strat->P.p2=p_Copy(strat->P.p, currRing, strat->tailRing);
921      // for input polys, prepare reduction
922      strat->P.PrepareRed(strat->use_buckets);
923    }
924
925#ifdef HAVE_RING2TOM
926    if (strat->P.p == NULL && strat->P.t_p == NULL) {
927      red_result = 0;
928    }
929    else
930#endif
931    {
932      if (TEST_OPT_PROT)
933        message((strat->honey ? strat->P.ecart : 0) + strat->P.pFDeg(),
934                &olddeg,&reduc,strat, red_result);
935
936      /* reduction of the element choosen from L */
937      red_result = strat->red(&strat->P,strat);
938    }
939
940    // reduction to non-zero new poly
941    if (red_result == 1)
942    {
943      /* statistic */
944      if (TEST_OPT_PROT) PrintS("s");
945
946      // get the polynomial (canonicalize bucket, make sure P.p is set)
947      strat->P.GetP(strat->lmBin);
948
949      int pos=posInS(strat,strat->sl,strat->P.p,strat->P.ecart);
950
951      // reduce the tail and normalize poly
952      if (TEST_OPT_INTSTRATEGY)
953      {
954        strat->P.pCleardenom();
955        if ((TEST_OPT_REDSB)||(TEST_OPT_REDTAIL))
956        {
957          strat->P.p = redtailBba(&(strat->P),pos-1,strat, withT);
958          strat->P.pCleardenom();
959        }
960      }
961      else
962      {
963        strat->P.pNorm();
964        if ((TEST_OPT_REDSB)||(TEST_OPT_REDTAIL))
965          strat->P.p = redtailBba(&(strat->P),pos-1,strat, withT);
966      }
967
968#ifdef KDEBUG
969      if (TEST_OPT_DEBUG){PrintS("new s:");strat->P.wrp();PrintLn();}
970#endif
971
972      // min_std stuff
973      if ((strat->P.p1==NULL) && (strat->minim>0))
974      {
975        if (strat->minim==1)
976        {
977          strat->M->m[minimcnt]=p_Copy(strat->P.p,currRing,strat->tailRing);
978          p_Delete(&strat->P.p2, currRing, strat->tailRing);
979        }
980        else
981        {
982          strat->M->m[minimcnt]=strat->P.p2;
983          strat->P.p2=NULL;
984        }
985        if (strat->tailRing!=currRing && pNext(strat->M->m[minimcnt])!=NULL)
986          pNext(strat->M->m[minimcnt])
987            = strat->p_shallow_copy_delete(pNext(strat->M->m[minimcnt]),
988                                           strat->tailRing, currRing,
989                                           currRing->PolyBin);
990        minimcnt++;
991      }
992
993      // enter into S, L, and T
994      enterT(strat->P, strat);
995      enterpairs(strat->P.p,strat->sl,strat->P.ecart,pos,strat, strat->tl);
996      // posInS only depends on the leading term
997      strat->enterS(strat->P, pos, strat, strat->tl);
998      if (hilb!=NULL) khCheck(Q,w,hilb,hilbeledeg,hilbcount,strat);
999//      Print("[%d]",hilbeledeg);
1000      if (strat->P.lcm!=NULL) pLmFree(strat->P.lcm);
1001      if (strat->sl>srmax) srmax = strat->sl;
1002    }
1003    else if (strat->P.p1 == NULL && strat->minim > 0)
1004    {
1005      p_Delete(&strat->P.p2, currRing, strat->tailRing);
1006    }
1007#ifdef KDEBUG
1008    memset(&(strat->P), 0, sizeof(strat->P));
1009#endif
1010    kTest_TS(strat);
1011  }
1012#ifdef KDEBUG
1013  if (TEST_OPT_DEBUG) messageSets(strat);
1014#endif
1015  /* complete reduction of the standard basis--------- */
1016  if (TEST_OPT_REDSB) completeReduce(strat);
1017  /* release temp data-------------------------------- */
1018  exitBuchMora(strat);
1019  if (TEST_OPT_WEIGHTM)
1020  {
1021    pRestoreDegProcs(pFDegOld, pLDegOld);
1022    if (ecartWeights)
1023    {
1024      omFreeSize((ADDRESS)ecartWeights,(pVariables+1)*sizeof(short));
1025      ecartWeights=NULL;
1026    }
1027  }
1028  if (TEST_OPT_PROT) messageStat(srmax,lrmax,hilbcount,strat);
1029  if (Q!=NULL) updateResult(strat->Shdl,Q,strat);
1030  return (strat->Shdl);
1031}
1032
1033poly kNF2 (ideal F,ideal Q,poly q,kStrategy strat, int lazyReduce)
1034{
1035  poly   p;
1036  int   i;
1037
1038  if ((idIs0(F))&&(Q==NULL))
1039    return pCopy(q); /*F=0*/
1040  strat->ak = idRankFreeModule(F);
1041  /*- creating temp data structures------------------- -*/
1042  BITSET save_test=test;
1043  test|=Sy_bit(OPT_REDTAIL);
1044  initBuchMoraCrit(strat);
1045  strat->initEcart = initEcartBBA;
1046  strat->enterS = enterSBba;
1047  /*- set S -*/
1048  strat->sl = -1;
1049  /*- init local data struct.---------------------------------------- -*/
1050  /*Shdl=*/initS(F,Q,strat);
1051  /*- compute------------------------------------------------------- -*/
1052  //if ((TEST_OPT_INTSTRATEGY)&&(lazyReduce==0))
1053  {
1054    for (i=strat->sl;i>=0;i--)
1055      pNorm(strat->S[i]);
1056  }
1057  kTest(strat);
1058  if (TEST_OPT_PROT) { PrintS("r"); mflush(); }
1059  p = redNF(pCopy(q),strat);
1060  if ((p!=NULL)&&(lazyReduce==0))
1061  {
1062    if (TEST_OPT_PROT) { PrintS("t"); mflush(); }
1063    p = redtailBba(p,strat->sl,strat);
1064  }
1065  /*- release temp data------------------------------- -*/
1066  omfree(strat->sevS);
1067  omfree(strat->ecartS);
1068  omfree(strat->T);
1069  omfree(strat->sevT);
1070  omfree(strat->R);
1071  omfree(strat->S_2_R);
1072  omfree(strat->L);
1073  omfree(strat->B);
1074  omfree(strat->fromQ);
1075  idDelete(&strat->Shdl);
1076  test=save_test;
1077  if (TEST_OPT_PROT) PrintLn();
1078  return p;
1079}
1080
1081ideal kNF2 (ideal F,ideal Q,ideal q,kStrategy strat, int lazyReduce)
1082{
1083  poly   p;
1084  int   i;
1085  ideal res;
1086
1087  if (idIs0(q))
1088    return idInit(IDELEMS(q),q->rank);
1089  if ((idIs0(F))&&(Q==NULL))
1090    return idCopy(q); /*F=0*/
1091  strat->ak = idRankFreeModule(F);
1092  /*- creating temp data structures------------------- -*/
1093  BITSET save_test=test;
1094  test|=Sy_bit(OPT_REDTAIL);
1095  initBuchMoraCrit(strat);
1096  strat->initEcart = initEcartBBA;
1097  strat->enterS = enterSBba;
1098  /*- set S -*/
1099  strat->sl = -1;
1100  /*- init local data struct.---------------------------------------- -*/
1101  /*Shdl=*/initS(F,Q,strat);
1102  /*- compute------------------------------------------------------- -*/
1103  res=idInit(IDELEMS(q),q->rank);
1104  //if ((TEST_OPT_INTSTRATEGY)&&(lazyReduce==0))
1105  {
1106    for (i=strat->sl;i>=0;i--)
1107      pNorm(strat->S[i]);
1108  }
1109  for (i=IDELEMS(q)-1; i>=0; i--)
1110  {
1111    if (q->m[i]!=NULL)
1112    {
1113      if (TEST_OPT_PROT) { PrintS("r");mflush(); }
1114      p = redNF(pCopy(q->m[i]),strat);
1115      if ((p!=NULL)&&(lazyReduce==0))
1116      {
1117        if (TEST_OPT_PROT) { PrintS("t"); mflush(); }
1118        p = redtailBba(p,strat->sl,strat);
1119      }
1120      res->m[i]=p;
1121    }
1122    //else
1123    //  res->m[i]=NULL;
1124  }
1125  /*- release temp data------------------------------- -*/
1126  omfree(strat->sevS);
1127  omfree(strat->ecartS);
1128  omfree(strat->T);
1129  omfree(strat->sevT);
1130  omfree(strat->R);
1131  omfree(strat->S_2_R);
1132  omfree(strat->L);
1133  omfree(strat->B);
1134  omfree(strat->fromQ);
1135  idDelete(&strat->Shdl);
1136  test=save_test;
1137  if (TEST_OPT_PROT) PrintLn();
1138  return res;
1139}
Note: See TracBrowser for help on using the repository browser.