source: git/kernel/kstd2.cc @ b2c236

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