source: git/kernel/kstd2.cc @ 206e158

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