source: git/kernel/kstd2.cc @ b3789a

fieker-DuValspielwiese
Last change on this file since b3789a was 8c36a9, checked in by Hans Schönemann <hannes@…>, 17 years ago
*hannes: option(length) git-svn-id: file:///usr/local/Singular/svn/trunk@9677 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 33.2 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/* $Id: kstd2.cc,v 1.35 2007-01-10 17:01:57 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  //if (h->GetLmTailRing()==NULL) return 0; // HS: SHOULD NOT BE NEEDED!
447  assume(h->FDeg == h->pFDeg());
448
449  poly h_p;
450  int i,j,at,pass, ii;
451  unsigned long not_sev;
452  long reddeg,d;
453
454  pass = j = 0;
455  d = reddeg = h->GetpFDeg();
456  h->SetShortExpVector();
457  int li;
458  h_p = h->GetLmTailRing();
459  not_sev = ~ h->sev;
460  loop
461  {
462    j = kFindDivisibleByInT(strat->T, strat->sevT, strat->tl, h);
463    if (j < 0) return 1;
464
465    li = strat->T[j].pLength;
466    ii = j;
467    /*
468     * the polynomial to reduce with (up to the moment) is;
469     * pi with length li
470     */
471    i = j;
472#if 1
473    if (TEST_OPT_LENGTH)
474    loop
475    {
476      /*- search the shortest possible with respect to length -*/
477      i++;
478      if (i > strat->tl)
479        break;
480      if (li<=1)
481        break;
482      if ((strat->T[i].pLength < li)
483         &&
484          p_LmShortDivisibleBy(strat->T[i].GetLmTailRing(), strat->sevT[i],
485                               h_p, not_sev, strat->tailRing))
486      {
487        /*
488         * the polynomial to reduce with is now;
489         */
490        li = strat->T[i].pLength;
491        ii = i;
492      }
493    }
494#endif
495
496    /*
497     * end of search: have to reduce with pi
498     */
499#ifdef KDEBUG
500    if (TEST_OPT_DEBUG)
501    {
502      PrintS("red:");
503      h->wrp();
504      PrintS(" with ");
505      strat->T[ii].wrp();
506    }
507#endif
508    assume(strat->fromT == FALSE);
509
510    ksReducePoly(h, &(strat->T[ii]), NULL, NULL, strat);
511
512#ifdef KDEBUG
513    if (TEST_OPT_DEBUG)
514    {
515      PrintS("\nto ");
516      h->wrp();
517      PrintLn();
518    }
519#endif
520
521    h_p = h->GetLmTailRing();
522    if (h_p == NULL)
523    {
524      if (h->lcm!=NULL) pLmFree(h->lcm);
525#ifdef KDEBUG
526      h->lcm=NULL;
527#endif
528      return 0;
529    }
530    h->SetShortExpVector();
531    not_sev = ~ h->sev;
532    /*
533     * try to reduce the s-polynomial h
534     *test first whether h should go to the lazyset L
535     *-if the degree jumps
536     *-if the number of pre-defined reductions jumps
537     */
538    pass++;
539    if (!K_TEST_OPT_REDTHROUGH && (strat->Ll >= 0) && (pass > strat->LazyPass))
540    {
541      h->SetLmCurrRing();
542      at = strat->posInL(strat->L,strat->Ll,h,strat);
543      if (at <= strat->Ll)
544      {
545        int dummy=strat->sl;
546        if (kFindDivisibleByInS(strat, &dummy, h) < 0)
547          return 1;
548        enterL(&strat->L,&strat->Ll,&strat->Lmax,*h,at);
549#ifdef KDEBUG
550        if (TEST_OPT_DEBUG)
551          Print(" lazy: -> L%d\n",at);
552#endif
553        h->Clear();
554        return -1;
555      }
556    }
557  }
558}
559
560/*2
561*  reduction procedure for the inhomogeneous case
562*  and not a degree-ordering
563*/
564int redLazy (LObject* h,kStrategy strat)
565{
566  if (strat->tl<0) return 1;
567  int at,d,i,ii,li;
568  int j = 0;
569  int pass = 0;
570  assume(h->pFDeg() == h->FDeg);
571  long reddeg = h->GetpFDeg();
572  unsigned long not_sev;
573
574  h->SetShortExpVector();
575  poly h_p = h->GetLmTailRing();
576  not_sev = ~ h->sev;
577  loop
578  {
579    j = kFindDivisibleByInT(strat->T, strat->sevT, strat->tl, h);
580    if (j < 0) return 1;
581   
582    li = strat->T[j].pLength;
583    #if 0
584    if (li==0)
585    {
586      li=strat->T[j].pLength=pLength(strat->T[j].p);
587    }
588    #endif
589    ii = j;
590    /*
591     * the polynomial to reduce with (up to the moment) is;
592     * pi with length li
593     */
594
595    i = j;
596#if 1
597    if (TEST_OPT_LENGTH)
598    loop
599    {
600      /*- search the shortest possible with respect to length -*/
601      i++;
602      if (i > strat->tl)
603        break;
604      if (li<=1)
605        break;
606    #if 0
607      if (strat->T[i].pLength==0)
608      {
609        PrintS("!");
610        strat->T[i].pLength=pLength(strat->T[i].p);
611      }
612   #endif
613      if ((strat->T[i].pLength < li)
614         &&
615          p_LmShortDivisibleBy(strat->T[i].GetLmTailRing(), strat->sevT[i],
616                               h_p, not_sev, strat->tailRing))
617      {
618        /*
619         * the polynomial to reduce with is now;
620         */
621        PrintS("+");
622        li = strat->T[i].pLength;
623        ii = i;
624      }
625    }
626#endif
627
628    /*
629     * end of search: have to reduce with pi
630     */
631
632
633#ifdef KDEBUG
634    if (TEST_OPT_DEBUG)
635    {
636      PrintS("red:");
637      h->wrp();
638      PrintS(" with ");
639      strat->T[ii].wrp();
640    }
641#endif
642
643    ksReducePoly(h, &(strat->T[ii]), NULL, NULL, strat);
644
645#ifdef KDEBUG
646    if (TEST_OPT_DEBUG)
647    {
648      PrintS("\nto ");
649      h->wrp();
650      PrintLn();
651    }
652#endif
653
654    h_p=h->GetLmTailRing();
655
656    if (h_p == NULL)
657    {
658      if (h->lcm!=NULL) pLmFree(h->lcm);
659#ifdef KDEBUG
660      h->lcm=NULL;
661#endif
662      return 0;
663    }
664    h->SetShortExpVector();
665    not_sev = ~ h->sev;
666    d = h->SetpFDeg();
667    /*- try to reduce the s-polynomial -*/
668    pass++;
669    if (//!K_TEST_OPT_REDTHROUGH &&
670        (strat->Ll >= 0) && ((d > reddeg) || (pass > strat->LazyPass)))
671    {
672      h->SetLmCurrRing();
673      at = strat->posInL(strat->L,strat->Ll,h,strat);
674      if (at <= strat->Ll)
675      {
676#if 1
677        int dummy=strat->sl;
678        if (kFindDivisibleByInS(strat, &dummy, h) < 0)
679          return 1;
680#endif
681#ifdef KDEBUG
682        if (TEST_OPT_DEBUG) Print(" ->L[%d]\n",at);
683#endif
684        enterL(&strat->L,&strat->Ll,&strat->Lmax,*h,at);
685        h->Clear();
686        return -1;
687      }
688    }
689    else if ((TEST_OPT_PROT) && (strat->Ll < 0) && (d != reddeg))
690    {
691      Print(".%d",d);mflush();
692      reddeg = d;
693    }
694  }
695}
696/*2
697*  reduction procedure for the sugar-strategy (honey)
698* reduces h with elements from T choosing first possible
699* element in T with respect to the given ecart
700*/
701int redHoney (LObject* h, kStrategy strat)
702{
703  if (strat->tl<0) return 1;
704  //if (h->GetLmTailRing()==NULL) return 0; // HS: SHOULD NOT BE NEEDED!
705  assume(h->FDeg == h->pFDeg());
706
707  poly h_p;
708  int i,j,at,pass,ei, ii, h_d;
709  unsigned long not_sev;
710  long reddeg,d;
711
712  pass = j = 0;
713  d = reddeg = h->GetpFDeg() + h->ecart;
714  h->SetShortExpVector();
715  int li;
716  h_p = h->GetLmTailRing();
717  not_sev = ~ h->sev;
718  loop
719  {
720    j = kFindDivisibleByInT(strat->T, strat->sevT, strat->tl, h);
721    if (j < 0) return 1;
722
723    ei = strat->T[j].ecart;
724    li = strat->T[j].pLength;
725    #if 0
726    if (li==0)
727    {
728       //PrintS("!");
729       li=strat->T[j].pLength=pLength(strat->T[j].p);
730    }
731    #endif
732    ii = j;
733    /*
734     * the polynomial to reduce with (up to the moment) is;
735     * pi with ecart ei
736     */
737    i = j;
738    if (TEST_OPT_LENGTH)
739    loop
740    {
741      /*- takes the first possible with respect to ecart -*/
742      i++;
743      if (i > strat->tl)
744        break;
745      //if (ei < h->ecart)
746      //  break;
747      if (li<=1)
748        break;
749      if ((((strat->T[i].ecart < ei) && (ei> h->ecart))
750         || ((strat->T[i].ecart <= h->ecart) && (strat->T[i].pLength < li)))
751         &&
752          p_LmShortDivisibleBy(strat->T[i].GetLmTailRing(), strat->sevT[i],
753                               h_p, not_sev, strat->tailRing))
754      {
755        /*
756         * the polynomial to reduce with is now;
757         */
758        ei = strat->T[i].ecart;
759        li = strat->T[i].pLength;
760        ii = i;
761      }
762    }
763
764    /*
765     * end of search: have to reduce with pi
766     */
767    if (!K_TEST_OPT_REDTHROUGH && (pass!=0) && (ei > h->ecart))
768    {
769      h->SetLmCurrRing();
770      /*
771       * It is not possible to reduce h with smaller ecart;
772       * if possible h goes to the lazy-set L,i.e
773       * if its position in L would be not the last one
774       */
775      if (strat->Ll >= 0) /* L is not empty */
776      {
777        at = strat->posInL(strat->L,strat->Ll,h,strat);
778        if(at <= strat->Ll)
779          /*- h will not become the next element to reduce -*/
780        {
781          enterL(&strat->L,&strat->Ll,&strat->Lmax,*h,at);
782#ifdef KDEBUG
783          if (TEST_OPT_DEBUG) Print(" ecart too big: -> L%d\n",at);
784#endif
785          h->Clear();
786          return -1;
787        }
788      }
789    }
790#ifdef KDEBUG
791    if (TEST_OPT_DEBUG)
792    {
793      PrintS("red:");
794      h->wrp();
795      PrintS(" with ");
796      strat->T[ii].wrp();
797    }
798#endif
799    assume(strat->fromT == FALSE);
800
801#if 0 // test poly exchange
802    if (strat->inStdFac==0)
803    {
804      int ll;
805      poly t_p;
806      if (strat->tailRing==currRing)
807        t_p=strat->T[ii].p;
808      else
809        t_p=strat->T[ii].t_p;
810      if ((p_LmCmp(h_p,t_p,strat->tailRing)==0)
811      && ((ll=h->GuessLength()) < strat->T[ii].pLength))
812      {
813        h->GetP();
814        if ((h->pLength=h->GetpLength()) < strat->T[ii].pLength)
815        {
816          if (TEST_OPT_PROT)  PrintS("e");
817          h->GetP();
818          if (h->p!=NULL)
819          {
820            if (strat->T[ii].p!=NULL)
821            {
822              poly swap;
823              omTypeAlloc0Bin(poly,swap,currRing->PolyBin);
824              memcpy(swap,h->p,currRing->PolyBin->sizeW*sizeof(long));
825              memcpy(h->p,strat->T[ii].p,currRing->PolyBin->sizeW*sizeof(long));
826              memcpy(strat->T[ii].p,swap,currRing->PolyBin->sizeW*sizeof(long));
827              omFreeBinAddr(swap);
828            }
829            else
830            {
831              strat->T[ii].p=h->p;
832              h->p=NULL;
833            }
834          }
835          else
836          {
837            if (strat->T[ii].p!=NULL)
838            {
839              h->p=strat->T[ii].p;
840              strat->T[ii].p=NULL;
841            }
842            // else: all NULL
843          }
844          if (h->t_p!=NULL)
845          {
846            if (strat->T[ii].t_p!=NULL)
847            {
848              poly swap;
849              omTypeAlloc0Bin(poly,swap,strat->tailRing->PolyBin);
850              memcpy(swap,h->t_p,strat->tailRing->PolyBin->sizeW*sizeof(long));
851              memcpy(h->t_p,strat->T[ii].t_p,strat->tailRing->PolyBin->sizeW*sizeof(long));
852              memcpy(strat->T[ii].t_p,swap,strat->tailRing->PolyBin->sizeW*sizeof(long));
853              omFreeBinAddr(swap);
854            }
855            else
856            {
857              strat->T[ii].t_p=h->t_p;
858              h->t_p=NULL;
859            }
860          }
861          else
862          {
863            if (strat->T[ii].t_p!=NULL)
864            {
865              h->t_p=strat->T[ii].t_p;
866              strat->T[ii].t_p=NULL;
867            }
868            // else: all NULL
869          }
870          if (strat->tailRing != currRing && (strat->T[ii].p != NULL)
871          && pNext(strat->T[ii].p) != NULL)
872             strat->T[ii].max =p_GetMaxExpP(pNext(strat->T[ii].p), strat->tailRing);
873          else
874             strat->T[ii].max = NULL;
875          h->length=h->pLength=pLength(h->p);
876          strat->T[ii].length=strat->T[ii].pLength=pLength(strat->T[ii].p);
877          if (strat->T[ii].is_normalized)
878          {
879            strat->T[ii].is_normalized=0;
880            strat->T[ii].pNorm();
881          }
882          else
883          {
884            if (TEST_OPT_INTSTRATEGY)
885              strat->T[ii].pCleardenom();
886          }
887          h->PrepareRed(strat->use_buckets);
888        }
889      }
890    }
891#endif // test poly exchange
892    ksReducePoly(h, &(strat->T[ii]), NULL, NULL, strat);
893
894#ifdef KDEBUG
895    if (TEST_OPT_DEBUG)
896    {
897      PrintS("\nto ");
898      h->wrp();
899      PrintLn();
900    }
901#endif
902
903    h_p = h->GetLmTailRing();
904    if (h_p == NULL)
905    {
906      if (h->lcm!=NULL) pLmFree(h->lcm);
907#ifdef KDEBUG
908      h->lcm=NULL;
909#endif
910      return 0;
911    }
912    h->SetShortExpVector();
913    not_sev = ~ h->sev;
914    h_d = h->SetpFDeg();
915    /* compute the ecart */
916    if (ei <= h->ecart)
917      h->ecart = d-h_d;
918    else
919      h->ecart = d-h_d+ei-h->ecart;
920    /*
921     * try to reduce the s-polynomial h
922     *test first whether h should go to the lazyset L
923     *-if the degree jumps
924     *-if the number of pre-defined reductions jumps
925     */
926    pass++;
927    d = h_d + h->ecart;
928    if (!K_TEST_OPT_REDTHROUGH && (strat->Ll >= 0) && ((d > reddeg) || (pass > strat->LazyPass)))
929    {
930      h->SetLmCurrRing();
931      at = strat->posInL(strat->L,strat->Ll,h,strat);
932      if (at <= strat->Ll)
933      {
934        int dummy=strat->sl;
935        if (kFindDivisibleByInS(strat, &dummy, h) < 0)
936          return 1;
937        enterL(&strat->L,&strat->Ll,&strat->Lmax,*h,at);
938#ifdef KDEBUG
939        if (TEST_OPT_DEBUG)
940          Print(" degree jumped: -> L%d\n",at);
941#endif
942        h->Clear();
943        return -1;
944      }
945    }
946    else if (TEST_OPT_PROT && (strat->Ll < 0) && (d > reddeg))
947    {
948      reddeg = d;
949      Print(".%d",d); mflush();
950    }
951  }
952}
953/*2
954*  reduction procedure for the normal form
955*/
956
957poly redNF (poly h,int &max_ind,kStrategy strat)
958{
959  if (h==NULL) return NULL;
960  int j;
961  max_ind=strat->sl;
962
963  if (0 > strat->sl)
964  {
965    return h;
966  }
967  LObject P(h);
968  P.SetShortExpVector();
969  P.bucket = kBucketCreate(currRing);
970  kBucketInit(P.bucket,P.p,pLength(P.p));
971  kbTest(P.bucket);
972  loop
973  {
974/* Obsolete since change in pLmDiv
975#ifdef HAVE_RING2TOM
976    if (currRing->cring == 1)
977    {
978      j=kRingFindDivisibleByInS(strat->S,strat->sevS,strat->sl,&P);
979    }
980    else
981#endif
982*/
983      j=kFindDivisibleByInS(strat,&max_ind,&P);
984    if (j>=0)
985    {
986      nNormalize(pGetCoeff(P.p));
987#ifdef KDEBUG
988      if (TEST_OPT_DEBUG)
989      {
990        PrintS("red:");
991        wrp(h);
992        PrintS(" with ");
993        wrp(strat->S[j]);
994      }
995#endif
996#ifdef HAVE_PLURAL
997      if (rIsPluralRing(currRing))
998      {
999        number coef;
1000        nc_BucketPolyRed(P.bucket,strat->S[j],&coef);
1001        nDelete(&coef);
1002      }
1003      else
1004#endif
1005      {
1006        number coef;
1007        coef=kBucketPolyRed(P.bucket,strat->S[j],pLength(strat->S[j]),strat->kNoether);
1008        nDelete(&coef);
1009      }
1010      h = kBucketGetLm(P.bucket);   // FRAGE OLIVER
1011      if (h==NULL)
1012      {
1013        kBucketDestroy(&P.bucket);
1014        return NULL;
1015      }
1016      kbTest(P.bucket);
1017      P.p=h;
1018      P.t_p=NULL;
1019      P.SetShortExpVector();
1020#ifdef KDEBUG
1021      if (TEST_OPT_DEBUG)
1022      {
1023        PrintS("\nto:");
1024        wrp(h);
1025        PrintLn();
1026      }
1027#endif
1028    }
1029    else
1030    {
1031      P.p=kBucketClear(P.bucket);
1032      kBucketDestroy(&P.bucket);
1033      pNormalize(P.p);
1034      return P.p;
1035    }
1036  }
1037}
1038
1039#ifdef KDEBUG
1040static int bba_count = 0;
1041#endif
1042
1043ideal bba (ideal F, ideal Q,intvec *w,intvec *hilb,kStrategy strat)
1044{
1045#ifdef KDEBUG
1046  bba_count++;
1047  int loop_count = 0;
1048#endif
1049  om_Opts.MinTrack = 5;
1050  int   srmax,lrmax, red_result = 1;
1051  int   olddeg,reduc;
1052  int hilbeledeg=1,hilbcount=0,minimcnt=0;
1053  BOOLEAN withT = FALSE;
1054
1055  initBuchMoraCrit(strat); /*set Gebauer, honey, sugarCrit*/
1056  initBuchMoraPos(strat);
1057  initHilbCrit(F,Q,&hilb,strat);
1058  initBba(F,strat);
1059  /*set enterS, spSpolyShort, reduce, red, initEcart, initEcartPair*/
1060  /*Shdl=*/initBuchMora(F, Q,strat);
1061  if (strat->minim>0) strat->M=idInit(IDELEMS(F),F->rank);
1062  srmax = strat->sl;
1063  reduc = olddeg = lrmax = 0;
1064
1065#ifndef NO_BUCKETS
1066  if (!TEST_OPT_NOT_BUCKETS)
1067    strat->use_buckets = 1;
1068#endif
1069
1070  // redtailBBa against T for inhomogenous input
1071  if (!K_TEST_OPT_OLDSTD)
1072    withT = ! strat->homog;
1073
1074  // strat->posInT = posInT_pLength;
1075  kTest_TS(strat);
1076
1077#ifdef HAVE_TAIL_RING
1078  kStratInitChangeTailRing(strat);
1079#endif
1080
1081  /* compute------------------------------------------------------- */
1082  while (strat->Ll >= 0)
1083  {
1084    if (strat->Ll > lrmax) lrmax =strat->Ll;/*stat.*/
1085#ifdef KDEBUG
1086    loop_count++;
1087#ifdef HAVE_RING2TOM
1088    if (TEST_OPT_DEBUG) PrintS("--- next step ---\n");
1089#endif
1090    if (TEST_OPT_DEBUG) messageSets(strat);
1091#endif
1092    if (strat->Ll== 0) strat->interpt=TRUE;
1093    if (TEST_OPT_DEGBOUND
1094        && ((strat->honey && (strat->L[strat->Ll].ecart+pFDeg(strat->L[strat->Ll].p,currRing)>Kstd1_deg))
1095            || ((!strat->honey) && (pFDeg(strat->L[strat->Ll].p,currRing)>Kstd1_deg))))
1096    {
1097      /*
1098       *stops computation if
1099       * 24 IN test and the degree +ecart of L[strat->Ll] is bigger then
1100       *a predefined number Kstd1_deg
1101       */
1102      while ((strat->Ll >= 0)
1103        && (strat->L[strat->Ll].p1!=NULL) && (strat->L[strat->Ll].p2!=NULL)
1104        && ((strat->honey && (strat->L[strat->Ll].ecart+pFDeg(strat->L[strat->Ll].p,currRing)>Kstd1_deg))
1105            || ((!strat->honey) && (pFDeg(strat->L[strat->Ll].p,currRing)>Kstd1_deg)))
1106        )
1107        deleteInL(strat->L,&strat->Ll,strat->Ll,strat);
1108      if (strat->Ll<0) break;
1109      else strat->noClearS=TRUE;
1110    }
1111    /* picks the last element from the lazyset L */
1112    strat->P = strat->L[strat->Ll];
1113    strat->Ll--;
1114
1115    if (pNext(strat->P.p) == strat->tail)
1116    {
1117      // deletes the short spoly
1118      pLmFree(strat->P.p);
1119      strat->P.p = NULL;
1120      poly m1 = NULL, m2 = NULL;
1121
1122      // check that spoly creation is ok
1123      while (strat->tailRing != currRing &&
1124             !kCheckSpolyCreation(&(strat->P), strat, m1, m2))
1125      {
1126        assume(m1 == NULL && m2 == NULL);
1127        // if not, change to a ring where exponents are at least
1128        // large enough
1129        kStratChangeTailRing(strat);
1130      }
1131      // create the real one
1132      ksCreateSpoly(&(strat->P), NULL, strat->use_buckets,
1133                    strat->tailRing, m1, m2, strat->R);
1134    }
1135    else if (strat->P.p1 == NULL)
1136    {
1137      if (strat->minim > 0)
1138        strat->P.p2=p_Copy(strat->P.p, currRing, strat->tailRing);
1139      // for input polys, prepare reduction
1140      strat->P.PrepareRed(strat->use_buckets);
1141    }
1142
1143    if (strat->P.p == NULL && strat->P.t_p == NULL)
1144    {
1145      red_result = 0;
1146    }
1147    else
1148    {
1149      if (TEST_OPT_PROT)
1150        message((strat->honey ? strat->P.ecart : 0) + strat->P.pFDeg(),
1151                &olddeg,&reduc,strat, red_result);
1152
1153      /* reduction of the element choosen from L */
1154      red_result = strat->red(&strat->P,strat);
1155    }
1156
1157    // reduction to non-zero new poly
1158    if (red_result == 1)
1159    {
1160      /* statistic */
1161      if (TEST_OPT_PROT) PrintS("s");
1162
1163      // get the polynomial (canonicalize bucket, make sure P.p is set)
1164      strat->P.GetP(strat->lmBin);
1165
1166      int pos=posInS(strat,strat->sl,strat->P.p,strat->P.ecart);
1167
1168      // reduce the tail and normalize poly
1169      if (TEST_OPT_INTSTRATEGY)
1170      {
1171        strat->P.pCleardenom();
1172        if ((TEST_OPT_REDSB)||(TEST_OPT_REDTAIL))
1173        {
1174          strat->P.p = redtailBba(&(strat->P),pos-1,strat, withT);
1175          strat->P.pCleardenom();
1176        }
1177      }
1178      else
1179      {
1180        strat->P.pNorm();
1181        if ((TEST_OPT_REDSB)||(TEST_OPT_REDTAIL))
1182          strat->P.p = redtailBba(&(strat->P),pos-1,strat, withT);
1183      }
1184
1185#ifdef KDEBUG
1186      if (TEST_OPT_DEBUG){PrintS("new s:");strat->P.wrp();PrintLn();}
1187#endif
1188
1189      // min_std stuff
1190      if ((strat->P.p1==NULL) && (strat->minim>0))
1191      {
1192        if (strat->minim==1)
1193        {
1194          strat->M->m[minimcnt]=p_Copy(strat->P.p,currRing,strat->tailRing);
1195          p_Delete(&strat->P.p2, currRing, strat->tailRing);
1196        }
1197        else
1198        {
1199          strat->M->m[minimcnt]=strat->P.p2;
1200          strat->P.p2=NULL;
1201        }
1202        if (strat->tailRing!=currRing && pNext(strat->M->m[minimcnt])!=NULL)
1203          pNext(strat->M->m[minimcnt])
1204            = strat->p_shallow_copy_delete(pNext(strat->M->m[minimcnt]),
1205                                           strat->tailRing, currRing,
1206                                           currRing->PolyBin);
1207        minimcnt++;
1208      }
1209
1210      // enter into S, L, and T
1211      //if ((!TEST_OPT_IDLIFT) || (pGetComp(strat->P.p) <= strat->syzComp))
1212        enterT(strat->P, strat);
1213#ifdef HAVE_RING2TOM
1214#ifdef HAVE_VANGB
1215      int at_R = strat->tl;
1216#endif
1217      if (currRing->cring == 1)
1218        superenterpairs(strat->P.p,strat->sl,strat->P.ecart,pos,strat, strat->tl);
1219      else
1220#endif
1221        enterpairs(strat->P.p,strat->sl,strat->P.ecart,pos,strat, strat->tl);
1222      // posInS only depends on the leading term
1223      if ((!TEST_OPT_IDLIFT) || (pGetComp(strat->P.p) <= strat->syzComp))
1224      {
1225#ifdef HAVE_VANGB
1226      strat->enterS(strat->P, pos, strat, at_R);
1227#else
1228      strat->enterS(strat->P, pos, strat, strat->tl);
1229#endif
1230      }
1231      else
1232      {
1233      //  strat->P.Delete(); // syzComp test: it is in T
1234      }
1235#if 0
1236      int pl=pLength(strat->P.p);
1237      if (pl==1)
1238      {
1239        //if (TEST_OPT_PROT)
1240        //PrintS("<1>");
1241      }
1242      else if (pl==2)
1243      {
1244        //if (TEST_OPT_PROT)
1245        //PrintS("<2>");
1246      }
1247#endif
1248      if (hilb!=NULL) khCheck(Q,w,hilb,hilbeledeg,hilbcount,strat);
1249//      Print("[%d]",hilbeledeg);
1250      if (strat->P.lcm!=NULL) pLmFree(strat->P.lcm);
1251      if (strat->sl>srmax) srmax = strat->sl;
1252    }
1253    else if (strat->P.p1 == NULL && strat->minim > 0)
1254    {
1255      p_Delete(&strat->P.p2, currRing, strat->tailRing);
1256    }
1257#ifdef KDEBUG
1258    memset(&(strat->P), 0, sizeof(strat->P));
1259#endif
1260    kTest_TS(strat);
1261  }
1262#ifdef KDEBUG
1263  if (TEST_OPT_DEBUG) messageSets(strat);
1264#endif
1265  /* complete reduction of the standard basis--------- */
1266  if (TEST_OPT_REDSB) completeReduce(strat);
1267  /* release temp data-------------------------------- */
1268  exitBuchMora(strat);
1269  if (TEST_OPT_WEIGHTM)
1270  {
1271    pRestoreDegProcs(pFDegOld, pLDegOld);
1272    if (ecartWeights)
1273    {
1274      omFreeSize((ADDRESS)ecartWeights,(pVariables+1)*sizeof(short));
1275      ecartWeights=NULL;
1276    }
1277  }
1278  if (TEST_OPT_PROT) messageStat(srmax,lrmax,hilbcount,strat);
1279  if (Q!=NULL) updateResult(strat->Shdl,Q,strat);
1280  return (strat->Shdl);
1281}
1282
1283poly kNF2 (ideal F,ideal Q,poly q,kStrategy strat, int lazyReduce)
1284{
1285  poly   p;
1286  int   i;
1287
1288  if ((idIs0(F))&&(Q==NULL))
1289    return pCopy(q); /*F=0*/
1290  strat->ak = idRankFreeModule(F);
1291  /*- creating temp data structures------------------- -*/
1292  BITSET save_test=test;
1293  test|=Sy_bit(OPT_REDTAIL);
1294  initBuchMoraCrit(strat);
1295  strat->initEcart = initEcartBBA;
1296  strat->enterS = enterSBba;
1297  /*- set S -*/
1298  strat->sl = -1;
1299  /*- init local data struct.---------------------------------------- -*/
1300  /*Shdl=*/initS(F,Q,strat);
1301  /*- compute------------------------------------------------------- -*/
1302  //if ((TEST_OPT_INTSTRATEGY)&&(lazyReduce==0))
1303  {
1304    for (i=strat->sl;i>=0;i--)
1305      pNorm(strat->S[i]);
1306  }
1307  kTest(strat);
1308  if (TEST_OPT_PROT) { PrintS("r"); mflush(); }
1309  int max_ind;
1310  p = redNF(pCopy(q),max_ind,strat);
1311  if ((p!=NULL)&&(lazyReduce==0))
1312  {
1313    if (TEST_OPT_PROT) { PrintS("t"); mflush(); }
1314    p = redtailBba(p,max_ind,strat);
1315  }
1316  /*- release temp data------------------------------- -*/
1317  omfree(strat->sevS);
1318  omfree(strat->ecartS);
1319  omfree(strat->T);
1320  omfree(strat->sevT);
1321  omfree(strat->R);
1322  omfree(strat->S_2_R);
1323  omfree(strat->L);
1324  omfree(strat->B);
1325  omfree(strat->fromQ);
1326  idDelete(&strat->Shdl);
1327  test=save_test;
1328  if (TEST_OPT_PROT) PrintLn();
1329  return p;
1330}
1331
1332ideal kNF2 (ideal F,ideal Q,ideal q,kStrategy strat, int lazyReduce)
1333{
1334  poly   p;
1335  int   i;
1336  ideal res;
1337  int max_ind;
1338
1339  if (idIs0(q))
1340    return idInit(IDELEMS(q),q->rank);
1341  if ((idIs0(F))&&(Q==NULL))
1342    return idCopy(q); /*F=0*/
1343  strat->ak = idRankFreeModule(F);
1344  /*- creating temp data structures------------------- -*/
1345  BITSET save_test=test;
1346  test|=Sy_bit(OPT_REDTAIL);
1347  initBuchMoraCrit(strat);
1348  strat->initEcart = initEcartBBA;
1349  strat->enterS = enterSBba;
1350  /*- set S -*/
1351  strat->sl = -1;
1352  /*- init local data struct.---------------------------------------- -*/
1353  /*Shdl=*/initS(F,Q,strat);
1354  /*- compute------------------------------------------------------- -*/
1355  res=idInit(IDELEMS(q),q->rank);
1356  //if ((TEST_OPT_INTSTRATEGY)&&(lazyReduce==0))
1357  {
1358    for (i=strat->sl;i>=0;i--)
1359      pNorm(strat->S[i]);
1360  }
1361  for (i=IDELEMS(q)-1; i>=0; i--)
1362  {
1363    if (q->m[i]!=NULL)
1364    {
1365      if (TEST_OPT_PROT) { PrintS("r");mflush(); }
1366      p = redNF(pCopy(q->m[i]),max_ind,strat);
1367      if ((p!=NULL)&&(lazyReduce==0))
1368      {
1369        if (TEST_OPT_PROT) { PrintS("t"); mflush(); }
1370        p = redtailBba(p,max_ind,strat);
1371      }
1372      res->m[i]=p;
1373    }
1374    //else
1375    //  res->m[i]=NULL;
1376  }
1377  /*- release temp data------------------------------- -*/
1378  omfree(strat->sevS);
1379  omfree(strat->ecartS);
1380  omfree(strat->T);
1381  omfree(strat->sevT);
1382  omfree(strat->R);
1383  omfree(strat->S_2_R);
1384  omfree(strat->L);
1385  omfree(strat->B);
1386  omfree(strat->fromQ);
1387  idDelete(&strat->Shdl);
1388  test=save_test;
1389  if (TEST_OPT_PROT) PrintLn();
1390  return res;
1391}
Note: See TracBrowser for help on using the repository browser.