source: git/kernel/kstd2.cc @ 7f06cca

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