source: git/kernel/kstd2.cc @ fbc0b6

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