source: git/kernel/kstd2.cc @ 994445

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