source: git/kernel/kstd2.cc @ 86016d

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