source: git/kernel/kutil.cc @ 6d2ef2

spielwiese
Last change on this file since 6d2ef2 was 6d2ef2, checked in by Hans Schoenemann <hannes@…>, 12 years ago
fix asujmes: strat->last git-svn-id: file:///usr/local/Singular/svn/trunk@14050 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 192.0 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/* $Id$ */
5/*
6* ABSTRACT: kernel: utils for kStd
7*/
8
9// #define PDEBUG 2
10// #define PDIV_DEBUG
11#define KUTIL_CC
12#include <stdlib.h>
13#include <string.h>
14#include <kernel/mod2.h>
15
16#ifndef NDEBUG
17# define MYTEST 0
18#else /* ifndef NDEBUG */
19# define MYTEST 0
20#endif /* ifndef NDEBUG */
21
22
23#include <omalloc/mylimits.h>
24#include <kernel/options.h>
25#include <kernel/gring.h>
26#include <kernel/sca.h>
27#ifdef KDEBUG
28#undef KDEBUG
29#define KDEBUG 2
30#endif
31
32#ifdef HAVE_RINGS
33#include <kernel/ideals.h>
34#endif
35
36// define if enterL, enterT should use memmove instead of doing it manually
37// on topgun, this is slightly faster (see monodromy_l.tst, homog_gonnet.sing)
38#ifndef SunOS_4
39#define ENTER_USE_MEMMOVE
40#endif
41
42// define, if the my_memmove inlines should be used instead of
43// system memmove -- it does not seem to pay off, though
44// #define ENTER_USE_MYMEMMOVE
45
46#include <kernel/kutil.h>
47#include <kernel/kbuckets.h>
48#include <kernel/febase.h>
49#include <omalloc/omalloc.h>
50#include <kernel/numbers.h>
51#include <kernel/polys.h>
52#include <kernel/ring.h>
53#include <kernel/ideals.h>
54#include <kernel/timer.h>
55//#include "cntrlc.h"
56#include <kernel/stairc.h>
57#include <kernel/kstd1.h>
58#include <kernel/pShallowCopyDelete.h>
59
60/* shiftgb stuff */
61#include <kernel/shiftgb.h>
62#include <kernel/prCopy.h>
63
64#ifdef HAVE_RATGRING
65#include <kernel/ratgring.h>
66#endif
67
68#ifdef KDEBUG
69#undef KDEBUG
70#define KDEBUG 2
71#endif
72
73
74#ifdef ENTER_USE_MYMEMMOVE
75inline void _my_memmove_d_gt_s(unsigned long* d, unsigned long* s, long l)
76{
77  register unsigned long* _dl = (unsigned long*) d;
78  register unsigned long* _sl = (unsigned long*) s;
79  register long _i = l - 1;
80
81  do
82  {
83    _dl[_i] = _sl[_i];
84    _i--;
85  }
86  while (_i >= 0);
87}
88
89inline void _my_memmove_d_lt_s(unsigned long* d, unsigned long* s, long l)
90{
91  register long _ll = l;
92  register unsigned long* _dl = (unsigned long*) d;
93  register unsigned long* _sl = (unsigned long*) s;
94  register long _i = 0;
95
96  do
97  {
98    _dl[_i] = _sl[_i];
99    _i++;
100  }
101  while (_i < _ll);
102}
103
104inline void _my_memmove(void* d, void* s, long l)
105{
106  unsigned long _d = (unsigned long) d;
107  unsigned long _s = (unsigned long) s;
108  unsigned long _l = ((l) + SIZEOF_LONG - 1) >> LOG_SIZEOF_LONG;
109
110  if (_d > _s) _my_memmove_d_gt_s(_d, _s, _l);
111  else _my_memmove_d_lt_s(_d, _s, _l);
112}
113
114#undef memmove
115#define memmove(d,s,l) _my_memmove(d, s, l)
116#endif
117
118static poly redMora (poly h,int maxIndex,kStrategy strat);
119static poly redBba (poly h,int maxIndex,kStrategy strat);
120
121#ifdef HAVE_RINGS
122#define pDivComp_EQUAL 2
123#define pDivComp_LESS 1
124#define pDivComp_GREATER -1
125#define pDivComp_INCOMP 0
126/* Checks the relation of LM(p) and LM(q)
127     LM(p) = LM(q) => return pDivComp_EQUAL
128     LM(p) | LM(q) => return pDivComp_LESS
129     LM(q) | LM(p) => return pDivComp_GREATER
130     else return pDivComp_INCOMP */
131static inline int pDivCompRing(poly p, poly q)
132{
133  if (pGetComp(p) == pGetComp(q))
134  {
135    BOOLEAN a=FALSE, b=FALSE;
136    int i;
137    unsigned long la, lb;
138    unsigned long divmask = currRing->divmask;
139    for (i=0; i<currRing->VarL_Size; i++)
140    {
141      la = p->exp[currRing->VarL_Offset[i]];
142      lb = q->exp[currRing->VarL_Offset[i]];
143      if (la != lb)
144      {
145        if (la < lb)
146        {
147          if (b) return pDivComp_INCOMP;
148          if (((la & divmask) ^ (lb & divmask)) != ((lb - la) & divmask))
149            return pDivComp_INCOMP;
150          a = TRUE;
151        }
152        else
153        {
154          if (a) return pDivComp_INCOMP;
155          if (((la & divmask) ^ (lb & divmask)) != ((la - lb) & divmask))
156            return pDivComp_INCOMP;
157          b = TRUE;
158        }
159      }
160    }
161    if (a) return pDivComp_LESS;
162    if (b) return pDivComp_GREATER;
163    if (!a & !b) return pDivComp_EQUAL;
164  }
165  return pDivComp_INCOMP;
166}
167#endif
168
169static inline int pDivComp(poly p, poly q)
170{
171  if (pGetComp(p) == pGetComp(q))
172  {
173#ifdef HAVE_RATGRING
174    if (rIsRatGRing(currRing))
175    {
176      if (_p_LmDivisibleByPart(p,currRing,
177                           q,currRing,
178                           currRing->real_var_start, currRing->real_var_end))
179        return 0;
180      return pLmCmp(q,p); // ONLY FOR GLOBAL ORDER!
181    }
182#endif
183    BOOLEAN a=FALSE, b=FALSE;
184    int i;
185    unsigned long la, lb;
186    unsigned long divmask = currRing->divmask;
187    for (i=0; i<currRing->VarL_Size; i++)
188    {
189      la = p->exp[currRing->VarL_Offset[i]];
190      lb = q->exp[currRing->VarL_Offset[i]];
191      if (la != lb)
192      {
193        if (la < lb)
194        {
195          if (b) return 0;
196          if (((la & divmask) ^ (lb & divmask)) != ((lb - la) & divmask))
197            return 0;
198          a = TRUE;
199        }
200        else
201        {
202          if (a) return 0;
203          if (((la & divmask) ^ (lb & divmask)) != ((la - lb) & divmask))
204            return 0;
205          b = TRUE;
206        }
207      }
208    }
209    if (a) { /*assume(pLmCmp(q,p)==1);*/ return 1; }
210    if (b) { /*assume(pLmCmp(q,p)==-1);*/return -1; }
211    /*assume(pLmCmp(q,p)==0);*/
212  }
213  return 0;
214}
215
216
217int     HCord;
218int     Kstd1_deg;
219int     Kstd1_mu=32000;
220
221/*2
222*deletes higher monomial of p, re-compute ecart and length
223*works only for orderings with ecart =pFDeg(end)-pFDeg(start)
224*/
225void deleteHC(LObject *L, kStrategy strat, BOOLEAN fromNext)
226{
227  if (strat->kHEdgeFound)
228  {
229    kTest_L(L);
230    poly p1;
231    poly p = L->GetLmTailRing();
232    int l = 1;
233    kBucket_pt bucket = NULL;
234    if (L->bucket != NULL)
235    {
236      kBucketClear(L->bucket, &pNext(p), &L->pLength);
237      L->pLength++;
238      bucket = L->bucket;
239      L->bucket = NULL;
240      L->last = NULL;
241    }
242
243    if (!fromNext && p_Cmp(p,strat->kNoetherTail(), L->tailRing) == -1)
244    {
245      L->Delete();
246      L->Clear();
247      L->ecart = -1;
248      if (bucket != NULL) kBucketDestroy(&bucket);
249      return;
250    }
251    p1 = p;
252    while (pNext(p1)!=NULL)
253    {
254      if (p_LmCmp(pNext(p1), strat->kNoetherTail(), L->tailRing) == -1)
255      {
256        L->last = p1;
257        p_Delete(&pNext(p1), L->tailRing);
258        if (p1 == p)
259        {
260          if (L->t_p != NULL)
261          {
262            assume(L->p != NULL && p == L->t_p);
263            pNext(L->p) = NULL;
264          }
265          L->max  = NULL;
266        }
267        else if (fromNext)
268          L->max  = p_GetMaxExpP(pNext(L->p), L->tailRing ); // p1;
269        //if (L->pLength != 0)
270        L->pLength = l;
271        // Hmmm when called from updateT, then only
272        // reset ecart when cut
273        if (fromNext)
274          L->ecart = L->pLDeg() - L->GetpFDeg();
275        break;
276      }
277      l++;
278      pIter(p1);
279    }
280    if (! fromNext)
281    {
282      L->SetpFDeg();
283      L->ecart = L->pLDeg(strat->LDegLast) - L->GetpFDeg();
284    }
285    if (bucket != NULL)
286    {
287      if (L->pLength > 1)
288      {
289        kBucketInit(bucket, pNext(p), L->pLength - 1);
290        pNext(p) = NULL;
291        if (L->t_p != NULL) pNext(L->t_p) = NULL;
292        L->pLength = 0;
293        L->bucket = bucket;
294        L->last = NULL;
295      }
296      else
297        kBucketDestroy(&bucket);
298    }
299    kTest_L(L);
300  }
301}
302
303void deleteHC(poly* p, int* e, int* l,kStrategy strat)
304{
305  LObject L(*p, currRing, strat->tailRing);
306
307  deleteHC(&L, strat);
308  *p = L.p;
309  *e = L.ecart;
310  *l = L.length;
311  if (L.t_p != NULL) p_LmFree(L.t_p, strat->tailRing);
312}
313
314/*2
315*tests if p.p=monomial*unit and cancels the unit
316*/
317void cancelunit (LObject* L,BOOLEAN inNF)
318{
319  int  i;
320  poly h;
321
322  if(rHasGlobalOrdering_currRing()) return;
323  if(TEST_OPT_CANCELUNIT) return;
324
325  ring r = L->tailRing;
326  poly p = L->GetLmTailRing();
327
328#ifdef HAVE_RINGS_LOC
329  // Leading coef have to be a unit
330  if ( !(nIsUnit(p_GetCoeff(p, r))) ) return;
331#endif
332
333  if(p_GetComp(p, r) != 0 && !p_OneComp(p, r)) return;
334
335//    for(i=r->N;i>0;i--)
336//    {
337//      if ((p_GetExp(p,i,r)>0) && (rIsPolyVar(i, r)==TRUE)) return;
338//    }
339  h = pNext(p);
340  loop
341  {
342    if (h==NULL)
343    {
344      p_Delete(&pNext(p), r);
345      if (!inNF)
346      {
347        number eins=nInit(1);
348        if (L->p != NULL)  pSetCoeff(L->p,eins);
349        else if (L->t_p != NULL) nDelete(&pGetCoeff(L->t_p));
350        if (L->t_p != NULL) pSetCoeff0(L->t_p,eins);
351      }
352      L->ecart = 0;
353      L->length = 1;
354      //if (L->pLength > 0)
355      L->pLength = 1;
356      if (L->last != NULL) L->last = p;
357      L->max = NULL;
358
359      if (L->t_p != NULL && pNext(L->t_p) != NULL)
360        pNext(L->t_p) = NULL;
361      if (L->p != NULL && pNext(L->p) != NULL)
362        pNext(L->p) = NULL;
363      return;
364    }
365    i = 0;
366    loop
367    {
368      i++;
369      if (p_GetExp(p,i,r) > p_GetExp(h,i,r)) return ; // does not divide
370      if (i == r->N) break; // does divide, try next monom
371    }
372    pIter(h);
373  }
374}
375
376/*2
377*pp is the new element in s
378*returns TRUE (in strat->kHEdgeFound) if
379*-HEcke is allowed
380*-we are in the last componente of the vector
381*-on all axis are monomials (all elements in NotUsedAxis are FALSE)
382*returns FALSE for pLexOrderings,
383*assumes in module case an ordering of type c* !!
384* HEckeTest is only called with strat->kHEdgeFound==FALSE !
385*/
386void HEckeTest (poly pp,kStrategy strat)
387{
388  int   j,k,p;
389
390  strat->kHEdgeFound=FALSE;
391  if (pLexOrder || currRing->MixedOrder)
392  {
393    return;
394  }
395  if (strat->ak > 1)           /*we are in the module case*/
396  {
397    return; // until ....
398    //if (!pVectorOut)     /*pVectorOut <=> order = c,* */
399    //  return FALSE;
400    //if (pGetComp(pp) < strat->ak) /* ak is the number of the last component */
401    //  return FALSE;
402  }
403  k = 0;
404  p=pIsPurePower(pp);
405  if (p!=0) strat->NotUsedAxis[p] = FALSE;
406  /*- the leading term of pp is a power of the p-th variable -*/
407  for (j=pVariables;j>0; j--)
408  {
409    if (strat->NotUsedAxis[j])
410    {
411      return;
412    }
413  }
414  strat->kHEdgeFound=TRUE;
415}
416
417/*2
418*utilities for TSet, LSet
419*/
420inline static intset initec (const int maxnr)
421{
422  return (intset)omAlloc(maxnr*sizeof(int));
423}
424
425inline static unsigned long* initsevS (const int maxnr)
426{
427  return (unsigned long*)omAlloc0(maxnr*sizeof(unsigned long));
428}
429inline static int* initS_2_R (const int maxnr)
430{
431  return (int*)omAlloc0(maxnr*sizeof(int));
432}
433
434static inline void enlargeT (TSet &T, TObject** &R, unsigned long* &sevT,
435                             int &length, const int incr)
436{
437  assume(T!=NULL);
438  assume(sevT!=NULL);
439  assume(R!=NULL);
440  assume((length+incr) > 0);
441
442  int i;
443  T = (TSet)omRealloc0Size(T, length*sizeof(TObject),
444                           (length+incr)*sizeof(TObject));
445
446  sevT = (unsigned long*) omReallocSize(sevT, length*sizeof(long*),
447                           (length+incr)*sizeof(long*));
448
449  R = (TObject**)omRealloc0Size(R,length*sizeof(TObject*),
450                                (length+incr)*sizeof(TObject*));
451  for (i=length-1;i>=0;i--) R[T[i].i_r] = &(T[i]);
452  length += incr;
453}
454
455void cleanT (kStrategy strat)
456{
457  int i,j;
458  poly  p;
459  assume(currRing == strat->tailRing || strat->tailRing != NULL);
460
461  pShallowCopyDeleteProc p_shallow_copy_delete =
462    (strat->tailRing != currRing ?
463     pGetShallowCopyDeleteProc(strat->tailRing, currRing) :
464     NULL);
465
466  for (j=0; j<=strat->tl; j++)
467  {
468    p = strat->T[j].p;
469    strat->T[j].p=NULL;
470    if (strat->T[j].max != NULL)
471    {
472      p_LmFree(strat->T[j].max, strat->tailRing);
473    }
474    i = -1;
475    loop
476    {
477      i++;
478      if (i>strat->sl)
479      {
480        if (strat->T[j].t_p != NULL)
481        {
482          p_Delete(&(strat->T[j].t_p), strat->tailRing);
483          p_LmFree(p, currRing);
484        }
485        else
486          pDelete(&p);
487        break;
488      }
489      if (p == strat->S[i])
490      {
491        if (strat->T[j].t_p != NULL)
492        {
493          assume(p_shallow_copy_delete != NULL);
494          pNext(p) = p_shallow_copy_delete(pNext(p),strat->tailRing,currRing,
495                                           currRing->PolyBin);
496          p_LmFree(strat->T[j].t_p, strat->tailRing);
497        }
498        break;
499      }
500    }
501  }
502  strat->tl=-1;
503}
504
505//LSet initL ()
506//{
507//  int i;
508//  LSet l = (LSet)omAlloc(setmaxL*sizeof(LObject));
509//  return l;
510//}
511
512static inline void enlargeL (LSet* L,int* length,const int incr)
513{
514  assume((*L)!=NULL);
515  assume((length+incr)>0);
516
517  *L = (LSet)omReallocSize((*L),(*length)*sizeof(LObject),
518                                   ((*length)+incr)*sizeof(LObject));
519  (*length) += incr;
520}
521
522void initPairtest(kStrategy strat)
523{
524  strat->pairtest = (BOOLEAN *)omAlloc0((strat->sl+2)*sizeof(BOOLEAN));
525}
526
527/*2
528*test whether (p1,p2) or (p2,p1) is in L up position length
529*it returns TRUE if yes and the position k
530*/
531BOOLEAN isInPairsetL(int length,poly p1,poly p2,int*  k,kStrategy strat)
532{
533  LObject *p=&(strat->L[length]);
534
535  *k = length;
536  loop
537  {
538    if ((*k) < 0) return FALSE;
539    if (((p1 == (*p).p1) && (p2 == (*p).p2))
540    ||  ((p1 == (*p).p2) && (p2 == (*p).p1)))
541      return TRUE;
542    (*k)--;
543    p--;
544  }
545}
546
547/*2
548*in B all pairs have the same element p on the right
549*it tests whether (q,p) is in B and returns TRUE if yes
550*and the position k
551*/
552BOOLEAN isInPairsetB(poly q,int*  k,kStrategy strat)
553{
554  LObject *p=&(strat->B[strat->Bl]);
555
556  *k = strat->Bl;
557  loop
558  {
559    if ((*k) < 0) return FALSE;
560    if (q == (*p).p1)
561      return TRUE;
562    (*k)--;
563    p--;
564  }
565}
566
567int kFindInT(poly p, TSet T, int tlength)
568{
569  int i;
570
571  for (i=0; i<=tlength; i++)
572  {
573    if (T[i].p == p) return i;
574  }
575  return -1;
576}
577
578int kFindInT(poly p, kStrategy strat)
579{
580  int i;
581  do
582  {
583    i = kFindInT(p, strat->T, strat->tl);
584    if (i >= 0) return i;
585    strat = strat->next;
586  }
587  while (strat != NULL);
588  return -1;
589}
590
591#ifdef KDEBUG
592
593void sTObject::wrp()
594{
595  if (t_p != NULL) p_wrp(t_p, tailRing);
596  else if (p != NULL) p_wrp(p, currRing, tailRing);
597  else ::wrp(NULL);
598}
599
600#define kFalseReturn(x) do { if (!x) return FALSE;} while (0)
601
602// check that Lm's of a poly from T are "equal"
603static const char* kTest_LmEqual(poly p, poly t_p, ring tailRing)
604{
605  int i;
606  for (i=1; i<=tailRing->N; i++)
607  {
608    if (p_GetExp(p, i, currRing) != p_GetExp(t_p, i, tailRing))
609      return "Lm[i] different";
610  }
611  if (p_GetComp(p, currRing) != p_GetComp(t_p, tailRing))
612    return "Lm[0] different";
613  if (pNext(p) != pNext(t_p))
614    return "Lm.next different";
615  if (pGetCoeff(p) != pGetCoeff(t_p))
616    return "Lm.coeff different";
617  return NULL;
618}
619
620static BOOLEAN sloppy_max = FALSE;
621BOOLEAN kTest_T(TObject * T, ring strat_tailRing, int i, char TN)
622{
623  ring tailRing = T->tailRing;
624  if (strat_tailRing == NULL) strat_tailRing = tailRing;
625  r_assume(strat_tailRing == tailRing);
626
627  poly p = T->p;
628  ring r = currRing;
629
630  if (T->p == NULL && T->t_p == NULL && i >= 0)
631    return dReportError("%c[%d].poly is NULL", TN, i);
632
633  if (T->tailRing != currRing)
634  {
635    if (T->t_p == NULL && i > 0)
636      return dReportError("%c[%d].t_p is NULL", TN, i);
637    pFalseReturn(p_Test(T->t_p, T->tailRing));
638    if (T->p != NULL) pFalseReturn(p_LmTest(T->p, currRing));
639    if (T->p != NULL && T->t_p != NULL)
640    {
641      const char* msg = kTest_LmEqual(T->p, T->t_p, T->tailRing);
642      if (msg != NULL)
643        return dReportError("%c[%d] %s", TN, i, msg);
644      r = T->tailRing;
645      p = T->t_p;
646    }
647    if (T->p == NULL)
648    {
649      p = T->t_p;
650      r = T->tailRing;
651    }
652    if (T->t_p != NULL && i >= 0 && TN == 'T')
653    {
654      if (pNext(T->t_p) == NULL)
655      {
656        if (T->max != NULL)
657          return dReportError("%c[%d].max is not NULL as it should be", TN, i);
658      }
659      else
660      {
661        if (T->max == NULL)
662          return dReportError("%c[%d].max is NULL", TN, i);
663        if (pNext(T->max) != NULL)
664          return dReportError("pNext(%c[%d].max) != NULL", TN, i);
665
666        pFalseReturn(p_CheckPolyRing(T->max, tailRing));
667        omCheckBinAddrSize(T->max, (tailRing->PolyBin->sizeW)*SIZEOF_LONG);
668#if KDEBUG > 0
669        if (! sloppy_max)
670        {
671          poly test_max = p_GetMaxExpP(pNext(T->t_p), tailRing);
672          p_Setm(T->max, tailRing);
673          p_Setm(test_max, tailRing);
674          BOOLEAN equal = p_ExpVectorEqual(T->max, test_max, tailRing);
675          if (! equal)
676            return dReportError("%c[%d].max out of sync", TN, i);
677          p_LmFree(test_max, tailRing);
678        }
679#endif
680      }
681    }
682  }
683  else
684  {
685    if (T->max != NULL)
686      return dReportError("%c[%d].max != NULL but tailRing == currRing",TN,i);
687    if (T->t_p != NULL)
688      return dReportError("%c[%d].t_p != NULL but tailRing == currRing",TN,i);
689    if (T->p == NULL && i > 0)
690      return dReportError("%c[%d].p is NULL", TN, i);
691    pFalseReturn(p_Test(T->p, currRing));
692  }
693
694  if (i >= 0 && T->pLength != 0
695  && ! rIsSyzIndexRing(currRing) && T->pLength != pLength(p))
696  {
697    int l=T->pLength;
698    T->pLength=pLength(p);
699    return dReportError("%c[%d] pLength error: has %d, specified to have %d",
700                        TN, i , pLength(p), l);
701  }
702
703  // check FDeg,  for elements in L and T
704  if (i >= 0 && (TN == 'T' || TN == 'L'))
705  {
706    // FDeg has ir element from T of L set
707    if (T->FDeg  != T->pFDeg())
708    {
709      int d=T->FDeg;
710      T->FDeg=T->pFDeg();
711      return dReportError("%c[%d] FDeg error: has %d, specified to have %d",
712                          TN, i , T->pFDeg(), d);
713    }
714  }
715
716  // check is_normalized for elements in T
717  if (i >= 0 && TN == 'T')
718  {
719    if (T->is_normalized && ! nIsOne(pGetCoeff(p)))
720      return dReportError("T[%d] is_normalized error", i);
721
722  }
723  return TRUE;
724}
725
726BOOLEAN kTest_L(LObject *L, ring strat_tailRing,
727                BOOLEAN testp, int lpos, TSet T, int tlength)
728{
729  if (testp)
730  {
731    poly pn = NULL;
732    if (L->bucket != NULL)
733    {
734      kFalseReturn(kbTest(L->bucket));
735      r_assume(L->bucket->bucket_ring == L->tailRing);
736      if (L->p != NULL && pNext(L->p) != NULL)
737      {
738        pn = pNext(L->p);
739        pNext(L->p) = NULL;
740      }
741    }
742    kFalseReturn(kTest_T(L, strat_tailRing, lpos, 'L'));
743    if (pn != NULL)
744      pNext(L->p) = pn;
745
746    ring r;
747    poly p;
748    L->GetLm(p, r);
749    if (L->sev != 0 && p_GetShortExpVector(p, r) != L->sev)
750    {
751      return dReportError("L[%d] wrong sev: has %o, specified to have %o",
752                          lpos, p_GetShortExpVector(p, r), L->sev);
753    }
754    if (lpos > 0 && L->last != NULL && pLast(p) != L->last)
755    {
756      return dReportError("L[%d] last wrong: has %p specified to have %p",
757                          lpos, pLast(p), L->last);
758    }
759  }
760  if (L->p1 == NULL)
761  {
762    // L->p2 either NULL or "normal" poly
763    pFalseReturn(pp_Test(L->p2, currRing, L->tailRing));
764  }
765  else if (tlength > 0 && T != NULL && (lpos >=0))
766  {
767    // now p1 and p2 must be != NULL and must be contained in T
768    int i;
769    i = kFindInT(L->p1, T, tlength);
770    if (i < 0)
771      return dReportError("L[%d].p1 not in T",lpos);
772    i = kFindInT(L->p2, T, tlength);
773    if (i < 0)
774      return dReportError("L[%d].p2 not in T",lpos);
775  }
776  return TRUE;
777}
778
779BOOLEAN kTest (kStrategy strat)
780{
781  int i;
782
783  // test P
784  kFalseReturn(kTest_L(&(strat->P), strat->tailRing,
785                       (strat->P.p != NULL && pNext(strat->P.p)!=strat->tail),
786                       -1, strat->T, strat->tl));
787
788  // test T
789  if (strat->T != NULL)
790  {
791    for (i=0; i<=strat->tl; i++)
792    {
793      kFalseReturn(kTest_T(&(strat->T[i]), strat->tailRing, i, 'T'));
794      if (strat->sevT[i] != pGetShortExpVector(strat->T[i].p))
795        return dReportError("strat->sevT[%d] out of sync", i);
796    }
797  }
798
799  // test L
800  if (strat->L != NULL)
801  {
802    for (i=0; i<=strat->Ll; i++)
803    {
804      kFalseReturn(kTest_L(&(strat->L[i]), strat->tailRing,
805                           strat->L[i].Next() != strat->tail, i,
806                           strat->T, strat->tl));
807      // may be unused
808      //if (strat->use_buckets && strat->L[i].Next() != strat->tail &&
809      //    strat->L[i].Next() != NULL && strat->L[i].p1 != NULL)
810      //{
811      //  assume(strat->L[i].bucket != NULL);
812      //}
813    }
814  }
815
816  // test S
817  if (strat->S != NULL)
818    kFalseReturn(kTest_S(strat));
819
820  return TRUE;
821}
822
823BOOLEAN kTest_S(kStrategy strat)
824{
825  int i;
826  BOOLEAN ret = TRUE;
827  for (i=0; i<=strat->sl; i++)
828  {
829    if (strat->S[i] != NULL &&
830        strat->sevS[i] != pGetShortExpVector(strat->S[i]))
831    {
832      return dReportError("S[%d] wrong sev: has %o, specified to have %o",
833                          i , pGetShortExpVector(strat->S[i]), strat->sevS[i]);
834    }
835  }
836  return ret;
837}
838
839
840
841BOOLEAN kTest_TS(kStrategy strat)
842{
843  int i, j;
844  BOOLEAN ret = TRUE;
845  kFalseReturn(kTest(strat));
846
847  // test strat->R, strat->T[i].i_r
848  for (i=0; i<=strat->tl; i++)
849  {
850    if (strat->T[i].i_r < 0 || strat->T[i].i_r > strat->tl)
851      return dReportError("strat->T[%d].i_r == %d out of bounds", i,
852                          strat->T[i].i_r);
853    if (strat->R[strat->T[i].i_r] != &(strat->T[i]))
854      return dReportError("T[%d].i_r with R out of sync", i);
855  }
856  // test containment of S inT
857  if (strat->S != NULL)
858  {
859    for (i=0; i<=strat->sl; i++)
860    {
861      j = kFindInT(strat->S[i], strat->T, strat->tl);
862      if (j < 0)
863        return dReportError("S[%d] not in T", i);
864      if (strat->S_2_R[i] != strat->T[j].i_r)
865        return dReportError("S_2_R[%d]=%d != T[%d].i_r=%d\n",
866                            i, strat->S_2_R[i], j, strat->T[j].i_r);
867    }
868  }
869  // test strat->L[i].i_r1
870  for (i=0; i<=strat->Ll; i++)
871  {
872    if (strat->L[i].p1 != NULL && strat->L[i].p2)
873    {
874      if (strat->L[i].i_r1 < 0 ||
875          strat->L[i].i_r1 > strat->tl ||
876          strat->L[i].T_1(strat)->p != strat->L[i].p1)
877        return dReportError("L[%d].i_r1 out of sync", i);
878      if (strat->L[i].i_r2 < 0 ||
879          strat->L[i].i_r2 > strat->tl ||
880          strat->L[i].T_2(strat)->p != strat->L[i].p2);
881    }
882    else
883    {
884      if (strat->L[i].i_r1 != -1)
885        return dReportError("L[%d].i_r1 out of sync", i);
886      if (strat->L[i].i_r2 != -1)
887        return dReportError("L[%d].i_r2 out of sync", i);
888    }
889    if (strat->L[i].i_r != -1)
890      return dReportError("L[%d].i_r out of sync", i);
891  }
892  return TRUE;
893}
894
895#endif // KDEBUG
896
897/*2
898*cancels the i-th polynomial in the standardbase s
899*/
900void deleteInS (int i,kStrategy strat)
901{
902#ifdef ENTER_USE_MEMMOVE
903  memmove(&(strat->S[i]), &(strat->S[i+1]), (strat->sl - i)*sizeof(poly));
904  memmove(&(strat->ecartS[i]),&(strat->ecartS[i+1]),(strat->sl - i)*sizeof(int));
905  memmove(&(strat->sevS[i]),&(strat->sevS[i+1]),(strat->sl - i)*sizeof(long));
906  memmove(&(strat->S_2_R[i]),&(strat->S_2_R[i+1]),(strat->sl - i)*sizeof(int));
907#else
908  int j;
909  for (j=i; j<strat->sl; j++)
910  {
911    strat->S[j] = strat->S[j+1];
912    strat->ecartS[j] = strat->ecartS[j+1];
913    strat->sevS[j] = strat->sevS[j+1];
914    strat->S_2_R[j] = strat->S_2_R[j+1];
915  }
916#endif
917  if (strat->lenS!=NULL)
918  {
919#ifdef ENTER_USE_MEMMOVE
920    memmove(&(strat->lenS[i]),&(strat->lenS[i+1]),(strat->sl - i)*sizeof(int));
921#else
922    for (j=i; j<strat->sl; j++) strat->lenS[j] = strat->lenS[j+1];
923#endif
924  }
925  if (strat->lenSw!=NULL)
926  {
927#ifdef ENTER_USE_MEMMOVE
928    memmove(&(strat->lenSw[i]),&(strat->lenSw[i+1]),(strat->sl - i)*sizeof(wlen_type));
929#else
930    for (j=i; j<strat->sl; j++) strat->lenSw[j] = strat->lenSw[j+1];
931#endif
932  }
933  if (strat->fromQ!=NULL)
934  {
935#ifdef ENTER_USE_MEMMOVE
936    memmove(&(strat->fromQ[i]),&(strat->fromQ[i+1]),(strat->sl - i)*sizeof(int));
937#else
938    for (j=i; j<strat->sl; j++)
939    {
940      strat->fromQ[j] = strat->fromQ[j+1];
941    }
942#endif
943  }
944  strat->S[strat->sl] = NULL;
945  strat->sl--;
946}
947
948/*2
949*cancels the j-th polynomial in the set
950*/
951void deleteInL (LSet set, int *length, int j,kStrategy strat)
952{
953  if (set[j].lcm!=NULL)
954  {
955#ifdef HAVE_RINGS
956    if (pGetCoeff(set[j].lcm) != NULL)
957      pLmDelete(set[j].lcm);
958    else
959#endif
960      pLmFree(set[j].lcm);
961  }
962  if (set[j].p!=NULL)
963  {
964    if (pNext(set[j].p) == strat->tail)
965    {
966#ifdef HAVE_RINGS
967      if (pGetCoeff(set[j].p) != NULL)
968        pLmDelete(set[j].p);
969      else
970#endif
971        pLmFree(set[j].p);
972      /*- tail belongs to several int spolys -*/
973    }
974    else
975    {
976      // search p in T, if it is there, do not delete it
977      if (pOrdSgn != -1 || kFindInT(set[j].p, strat) < 0)
978      {
979        // assure that for global orderings kFindInT fails
980        assume(pOrdSgn == -1 || kFindInT(set[j].p, strat) < 0);
981        set[j].Delete();
982      }
983    }
984  }
985  if (*length > 0 && j < *length)
986  {
987#ifdef ENTER_USE_MEMMOVE
988    memmove(&(set[j]), &(set[j+1]), (*length - j)*sizeof(LObject));
989#else
990    int i;
991    for (i=j; i < (*length); i++)
992      set[i] = set[i+1];
993#endif
994  }
995#ifdef KDEBUG
996  memset(&(set[*length]),0,sizeof(LObject));
997#endif
998  (*length)--;
999}
1000
1001/*2
1002*enters p at position at in L
1003*/
1004void enterL (LSet *set,int *length, int *LSetmax, LObject p,int at)
1005{
1006#ifdef PDEBUG
1007  /*  zaehler++; */
1008#endif /*PDEBUG*/
1009  int i;
1010  // this should be corrected
1011  assume(p.FDeg == p.pFDeg());
1012
1013  if ((*length)>=0)
1014  {
1015    if ((*length) == (*LSetmax)-1) enlargeL(set,LSetmax,setmaxLinc);
1016    if (at <= (*length))
1017#ifdef ENTER_USE_MEMMOVE
1018      memmove(&((*set)[at+1]), &((*set)[at]), ((*length)-at+1)*sizeof(LObject));
1019#else
1020    for (i=(*length)+1; i>=at+1; i--) (*set)[i] = (*set)[i-1];
1021#endif
1022  }
1023  else at = 0;
1024  (*set)[at] = p;
1025  (*length)++;
1026}
1027
1028/*2
1029* computes the normal ecart;
1030* used in mora case and if pLexOrder & sugar in bba case
1031*/
1032void initEcartNormal (LObject* h)
1033{
1034  h->FDeg = h->pFDeg();
1035  h->ecart = h->pLDeg() - h->FDeg;
1036  // h->length is set by h->pLDeg
1037  h->length=h->pLength=pLength(h->p);
1038}
1039
1040void initEcartBBA (LObject* h)
1041{
1042  h->FDeg = h->pFDeg();
1043  (*h).ecart = 0;
1044  h->length=h->pLength=pLength(h->p);
1045}
1046
1047void initEcartPairBba (LObject* Lp,poly f,poly g,int ecartF,int ecartG)
1048{
1049  Lp->FDeg = Lp->pFDeg();
1050  (*Lp).ecart = 0;
1051  (*Lp).length = 0;
1052}
1053
1054void initEcartPairMora (LObject* Lp,poly f,poly g,int ecartF,int ecartG)
1055{
1056  Lp->FDeg = Lp->pFDeg();
1057  (*Lp).ecart = si_max(ecartF,ecartG);
1058  (*Lp).ecart = (*Lp).ecart- (Lp->FDeg -pFDeg((*Lp).lcm,currRing));
1059  (*Lp).length = 0;
1060}
1061
1062/*2
1063*if ecart1<=ecart2 it returns TRUE
1064*/
1065static inline BOOLEAN sugarDivisibleBy(int ecart1, int ecart2)
1066{
1067  return (ecart1 <= ecart2);
1068}
1069
1070#ifdef HAVE_RINGS
1071/*2
1072* put the pair (s[i],p)  into the set B, ecart=ecart(p) (ring case)
1073*/
1074void enterOnePairRing (int i,poly p,int ecart, int isFromQ,kStrategy strat, int atR = -1)
1075{
1076  assume(i<=strat->sl);
1077  int      l,j,compare,compareCoeff;
1078  LObject  Lp;
1079
1080  if (strat->interred_flag) return;
1081#ifdef KDEBUG
1082  Lp.ecart=0; Lp.length=0;
1083#endif
1084  /*- computes the lcm(s[i],p) -*/
1085  Lp.lcm = pInit();
1086  pSetCoeff0(Lp.lcm, nLcm(pGetCoeff(p), pGetCoeff(strat->S[i]), currRing));
1087  // Lp.lcm == 0
1088  if (nIsZero(pGetCoeff(Lp.lcm)))
1089  {
1090#ifdef KDEBUG
1091      if (TEST_OPT_DEBUG)
1092      {
1093        PrintS("--- Lp.lcm == 0\n");
1094        PrintS("p:");
1095        wrp(p);
1096        Print("  strat->S[%d]:", i);
1097        wrp(strat->S[i]);
1098        PrintLn();
1099      }
1100#endif
1101      strat->cp++;
1102      pLmDelete(Lp.lcm);
1103      return;
1104  }
1105  // basic product criterion
1106  pLcm(p,strat->S[i],Lp.lcm);
1107  pSetm(Lp.lcm);
1108  assume(!strat->sugarCrit);
1109  if (pHasNotCF(p,strat->S[i]) && nIsUnit(pGetCoeff(p)) && nIsUnit(pGetCoeff(strat->S[i])))
1110  {
1111#ifdef KDEBUG
1112      if (TEST_OPT_DEBUG)
1113      {
1114        PrintS("--- product criterion func enterOnePairRing type 1\n");
1115        PrintS("p:");
1116        wrp(p);
1117        Print("  strat->S[%d]:", i);
1118        wrp(strat->S[i]);
1119        PrintLn();
1120      }
1121#endif
1122      strat->cp++;
1123      pLmDelete(Lp.lcm);
1124      return;
1125  }
1126  assume(!strat->fromT);
1127  /*
1128  *the set B collects the pairs of type (S[j],p)
1129  *suppose (r,p) is in B and (s,p) is the new pair and lcm(s,p) != lcm(r,p)
1130  *if the leading term of s devides lcm(r,p) then (r,p) will be canceled
1131  *if the leading term of r devides lcm(s,p) then (s,p) will not enter B
1132  */
1133  for(j = strat->Bl;j>=0;j--)
1134  {
1135    compare=pDivCompRing(strat->B[j].lcm,Lp.lcm);
1136    compareCoeff = nDivComp(pGetCoeff(strat->B[j].lcm), pGetCoeff(Lp.lcm));
1137    if (compareCoeff == pDivComp_EQUAL || compare == compareCoeff)
1138    {
1139      if (compare == 1)
1140      {
1141        strat->c3++;
1142#ifdef KDEBUG
1143        if (TEST_OPT_DEBUG)
1144        {
1145          PrintS("--- chain criterion type 1\n");
1146          PrintS("strat->B[j]:");
1147          wrp(strat->B[j].lcm);
1148          PrintS("  Lp.lcm:");
1149          wrp(Lp.lcm);
1150          PrintLn();
1151        }
1152#endif
1153        if ((strat->fromQ==NULL) || (isFromQ==0) || (strat->fromQ[i]==0))
1154        {
1155          pLmDelete(Lp.lcm);
1156          return;
1157        }
1158        break;
1159      }
1160      else
1161      if (compare == -1)
1162      {
1163#ifdef KDEBUG
1164        if (TEST_OPT_DEBUG)
1165        {
1166          PrintS("--- chain criterion type 2\n");
1167          Print("strat->B[%d].lcm:",j);
1168          wrp(strat->B[j].lcm);
1169          PrintS("  Lp.lcm:");
1170          wrp(Lp.lcm);
1171          PrintLn();
1172        }
1173#endif
1174        deleteInL(strat->B,&strat->Bl,j,strat);
1175        strat->c3++;
1176      }
1177    }
1178    if ((compare == pDivComp_EQUAL) && (compareCoeff != 2))
1179    {
1180      if (compareCoeff == pDivComp_LESS)
1181      {
1182#ifdef KDEBUG
1183        if (TEST_OPT_DEBUG)
1184        {
1185          PrintS("--- chain criterion type 3\n");
1186          Print("strat->B[%d].lcm:", j);
1187          wrp(strat->B[j].lcm);
1188          PrintS("  Lp.lcm:");
1189          wrp(Lp.lcm);
1190          PrintLn();
1191        }
1192#endif
1193        strat->c3++;
1194        if ((strat->fromQ==NULL) || (isFromQ==0) || (strat->fromQ[i]==0))
1195        {
1196          pLmDelete(Lp.lcm);
1197          return;
1198        }
1199        break;
1200      }
1201      else
1202      // Add hint for same LM and LC (later) (TODO Oliver)
1203      // if (compareCoeff == pDivComp_GREATER)
1204      {
1205#ifdef KDEBUG
1206        if (TEST_OPT_DEBUG)
1207        {
1208          PrintS("--- chain criterion type 4\n");
1209          Print("strat->B[%d].lcm:", j);
1210          wrp(strat->B[j].lcm);
1211          PrintS("  Lp.lcm:");
1212          wrp(Lp.lcm);
1213          PrintLn();
1214        }
1215#endif
1216        deleteInL(strat->B,&strat->Bl,j,strat);
1217        strat->c3++;
1218      }
1219    }
1220  }
1221  /*
1222  *the pair (S[i],p) enters B if the spoly != 0
1223  */
1224  /*-  compute the short s-polynomial -*/
1225  if ((strat->S[i]==NULL) || (p==NULL)) {
1226#ifdef KDEBUG
1227    if (TEST_OPT_DEBUG)
1228    {
1229      PrintS("--- spoly = NULL\n");
1230    }
1231#endif
1232    pLmDelete(Lp.lcm);
1233    return;
1234  }
1235  if ((strat->fromQ!=NULL) && (isFromQ!=0) && (strat->fromQ[i]!=0))
1236  {
1237    // Is from a previous computed GB, therefore we know that spoly will
1238    // reduce to zero. Oliver.
1239    WarnS("Could we come here? 8738947389");
1240    Lp.p=NULL;
1241  }
1242  else
1243  {
1244    Lp.p = ksCreateShortSpoly(strat->S[i], p, strat->tailRing);
1245  }
1246  if (Lp.p == NULL)
1247  {
1248#ifdef KDEBUG
1249    if (TEST_OPT_DEBUG)
1250    {
1251      PrintS("--- spoly = NULL\n");
1252    }
1253#endif
1254    /*- the case that the s-poly is 0 -*/
1255    if (strat->pairtest==NULL) initPairtest(strat);
1256    strat->pairtest[i] = TRUE;/*- hint for spoly(S^[i],p)=0 -*/
1257    strat->pairtest[strat->sl+1] = TRUE;
1258    /*hint for spoly(S[i],p) == 0 for some i,0 <= i <= sl*/
1259    /*
1260    *suppose we have (s,r),(r,p),(s,p) and spoly(s,p) == 0 and (r,p) is
1261    *still in B (i.e. lcm(r,p) == lcm(s,p) or the leading term of s does not
1262    *devide lcm(r,p)). In the last case (s,r) can be canceled if the leading
1263    *term of p devides the lcm(s,r)
1264    *(this canceling should be done here because
1265    *the case lcm(s,p) == lcm(s,r) is not covered in chainCrit)
1266    *the first case is handeled in chainCrit
1267    */
1268    pLmDelete(Lp.lcm);
1269  }
1270  else
1271  {
1272    /*- the pair (S[i],p) enters B -*/
1273    Lp.p1 = strat->S[i];
1274    Lp.p2 = p;
1275
1276    pNext(Lp.p) = strat->tail;
1277
1278    if (atR >= 0)
1279    {
1280      Lp.i_r2 = atR;
1281      Lp.i_r1 = strat->S_2_R[i];
1282    }
1283    strat->initEcartPair(&Lp,strat->S[i],p,strat->ecartS[i],ecart);
1284    l = strat->posInL(strat->B,strat->Bl,&Lp,strat);
1285    enterL(&strat->B,&strat->Bl,&strat->Bmax,Lp,l);
1286  }
1287}
1288
1289
1290/*2
1291* put the  lcm(s[i],p)  into the set B
1292*/
1293
1294BOOLEAN enterOneStrongPoly (int i,poly p,int ecart, int isFromQ,kStrategy strat, int atR = -1)
1295{
1296  number d, s, t;
1297  assume(i<=strat->sl);
1298  assume(atR >= 0);
1299  poly m1, m2, gcd;
1300
1301  d = nExtGcd(pGetCoeff(p), pGetCoeff(strat->S[i]), &s, &t);
1302
1303  if (nIsZero(s) || nIsZero(t))  // evtl. durch divBy tests ersetzen
1304  {
1305    nDelete(&d);
1306    nDelete(&s);
1307    nDelete(&t);
1308    return FALSE;
1309  }
1310
1311  k_GetStrongLeadTerms(p, strat->S[i], currRing, m1, m2, gcd, strat->tailRing);
1312  //p_Test(m1,strat->tailRing);
1313  //p_Test(m2,strat->tailRing);
1314  while (! kCheckStrongCreation(atR, m1, i, m2, strat) )
1315  {
1316    memset(&(strat->P), 0, sizeof(strat->P));
1317    kStratChangeTailRing(strat);
1318    strat->P = *(strat->R[atR]);
1319    p_LmFree(m1, strat->tailRing);
1320    p_LmFree(m2, strat->tailRing);
1321    p_LmFree(gcd, currRing);
1322    k_GetStrongLeadTerms(p, strat->S[i], currRing, m1, m2, gcd, strat->tailRing);
1323  }
1324  pSetCoeff0(m1, s);
1325  pSetCoeff0(m2, t);
1326  pSetCoeff0(gcd, d);
1327  p_Test(m1,strat->tailRing);
1328  p_Test(m2,strat->tailRing);
1329
1330#ifdef KDEBUG
1331  if (TEST_OPT_DEBUG)
1332  {
1333    // Print("t = %d; s = %d; d = %d\n", nInt(t), nInt(s), nInt(d));
1334    PrintS("m1 = ");
1335    p_wrp(m1, strat->tailRing);
1336    PrintS(" ; m2 = ");
1337    p_wrp(m2, strat->tailRing);
1338    PrintS(" ; gcd = ");
1339    wrp(gcd);
1340    PrintS("\n--- create strong gcd poly: ");
1341    Print("\n p: ", i);
1342    wrp(p);
1343    Print("\n strat->S[%d]: ", i);
1344    wrp(strat->S[i]);
1345    PrintS(" ---> ");
1346  }
1347#endif
1348
1349  pNext(gcd) = p_Add_q(pp_Mult_mm(pNext(p), m1, strat->tailRing), pp_Mult_mm(pNext(strat->S[i]), m2, strat->tailRing), strat->tailRing);
1350  p_LmDelete(m1, strat->tailRing);
1351  p_LmDelete(m2, strat->tailRing);
1352
1353#ifdef KDEBUG
1354  if (TEST_OPT_DEBUG)
1355  {
1356    wrp(gcd);
1357    PrintLn();
1358  }
1359#endif
1360
1361  LObject h;
1362  h.p = gcd;
1363  h.tailRing = strat->tailRing;
1364  int posx;
1365  h.pCleardenom();
1366  strat->initEcart(&h);
1367  if (strat->Ll==-1)
1368    posx =0;
1369  else
1370    posx = strat->posInL(strat->L,strat->Ll,&h,strat);
1371  h.sev = pGetShortExpVector(h.p);
1372  if (currRing!=strat->tailRing)
1373    h.t_p = k_LmInit_currRing_2_tailRing(h.p, strat->tailRing);
1374  enterL(&strat->L,&strat->Ll,&strat->Lmax,h,posx);
1375  return TRUE;
1376}
1377#endif
1378
1379/*2
1380* put the pair (s[i],p)  into the set B, ecart=ecart(p)
1381*/
1382
1383void enterOnePairNormal (int i,poly p,int ecart, int isFromQ,kStrategy strat, int atR = -1)
1384{
1385  assume(i<=strat->sl);
1386  if (strat->interred_flag) return;
1387
1388  int      l,j,compare;
1389  LObject  Lp;
1390  Lp.i_r = -1;
1391
1392#ifdef KDEBUG
1393  Lp.ecart=0; Lp.length=0;
1394#endif
1395  /*- computes the lcm(s[i],p) -*/
1396  Lp.lcm = pInit();
1397
1398#ifndef HAVE_RATGRING
1399  pLcm(p,strat->S[i],Lp.lcm);
1400#elif defined(HAVE_RATGRING)
1401  //  if (rIsRatGRing(currRing))
1402  pLcmRat(p,strat->S[i],Lp.lcm, currRing->real_var_start); // int rat_shift
1403#endif
1404  pSetm(Lp.lcm);
1405
1406
1407  if (strat->sugarCrit && ALLOW_PROD_CRIT(strat))
1408  {
1409    if((!((strat->ecartS[i]>0)&&(ecart>0)))
1410    && pHasNotCF(p,strat->S[i]))
1411    {
1412    /*
1413    *the product criterion has applied for (s,p),
1414    *i.e. lcm(s,p)=product of the leading terms of s and p.
1415    *Suppose (s,r) is in L and the leading term
1416    *of p divides lcm(s,r)
1417    *(==> the leading term of p divides the leading term of r)
1418    *but the leading term of s does not divide the leading term of r
1419    *(notice that tis condition is automatically satisfied if r is still
1420    *in S), then (s,r) can be cancelled.
1421    *This should be done here because the
1422    *case lcm(s,r)=lcm(s,p) is not covered by chainCrit.
1423    *
1424    *Moreover, skipping (s,r) holds also for the noncommutative case.
1425    */
1426      strat->cp++;
1427      pLmFree(Lp.lcm);
1428      Lp.lcm=NULL;
1429      return;
1430    }
1431    else
1432      Lp.ecart = si_max(ecart,strat->ecartS[i]);
1433    if (strat->fromT && (strat->ecartS[i]>ecart))
1434    {
1435      pLmFree(Lp.lcm);
1436      Lp.lcm=NULL;
1437      return;
1438      /*the pair is (s[i],t[.]), discard it if the ecart is too big*/
1439    }
1440    /*
1441    *the set B collects the pairs of type (S[j],p)
1442    *suppose (r,p) is in B and (s,p) is the new pair and lcm(s,p)#lcm(r,p)
1443    *if the leading term of s devides lcm(r,p) then (r,p) will be canceled
1444    *if the leading term of r devides lcm(s,p) then (s,p) will not enter B
1445    */
1446    {
1447      j = strat->Bl;
1448      loop
1449      {
1450        if (j < 0)  break;
1451        compare=pDivComp(strat->B[j].lcm,Lp.lcm);
1452        if ((compare==1)
1453        &&(sugarDivisibleBy(strat->B[j].ecart,Lp.ecart)))
1454        {
1455          strat->c3++;
1456          if ((strat->fromQ==NULL) || (isFromQ==0) || (strat->fromQ[i]==0))
1457          {
1458            pLmFree(Lp.lcm);
1459            return;
1460          }
1461          break;
1462        }
1463        else
1464        if ((compare ==-1)
1465        && sugarDivisibleBy(Lp.ecart,strat->B[j].ecart))
1466        {
1467          deleteInL(strat->B,&strat->Bl,j,strat);
1468          strat->c3++;
1469        }
1470        j--;
1471      }
1472    }
1473  }
1474  else /*sugarcrit*/
1475  {
1476    if (ALLOW_PROD_CRIT(strat))
1477    {
1478      // if currRing->nc_type!=quasi (or skew)
1479      // TODO: enable productCrit for super commutative algebras...
1480      if(/*(strat->ak==0) && productCrit(p,strat->S[i])*/
1481      pHasNotCF(p,strat->S[i]))
1482      {
1483      /*
1484      *the product criterion has applied for (s,p),
1485      *i.e. lcm(s,p)=product of the leading terms of s and p.
1486      *Suppose (s,r) is in L and the leading term
1487      *of p devides lcm(s,r)
1488      *(==> the leading term of p devides the leading term of r)
1489      *but the leading term of s does not devide the leading term of r
1490      *(notice that tis condition is automatically satisfied if r is still
1491      *in S), then (s,r) can be canceled.
1492      *This should be done here because the
1493      *case lcm(s,r)=lcm(s,p) is not covered by chainCrit.
1494      */
1495          strat->cp++;
1496          pLmFree(Lp.lcm);
1497          Lp.lcm=NULL;
1498          return;
1499      }
1500      if (strat->fromT && (strat->ecartS[i]>ecart))
1501      {
1502        pLmFree(Lp.lcm);
1503        Lp.lcm=NULL;
1504        return;
1505        /*the pair is (s[i],t[.]), discard it if the ecart is too big*/
1506      }
1507      /*
1508      *the set B collects the pairs of type (S[j],p)
1509      *suppose (r,p) is in B and (s,p) is the new pair and lcm(s,p)#lcm(r,p)
1510      *if the leading term of s devides lcm(r,p) then (r,p) will be canceled
1511      *if the leading term of r devides lcm(s,p) then (s,p) will not enter B
1512      */
1513      for(j = strat->Bl;j>=0;j--)
1514      {
1515        compare=pDivComp(strat->B[j].lcm,Lp.lcm);
1516        if (compare==1)
1517        {
1518          strat->c3++;
1519          if ((strat->fromQ==NULL) || (isFromQ==0) || (strat->fromQ[i]==0))
1520          {
1521            pLmFree(Lp.lcm);
1522            return;
1523          }
1524          break;
1525        }
1526        else
1527        if (compare ==-1)
1528        {
1529          deleteInL(strat->B,&strat->Bl,j,strat);
1530          strat->c3++;
1531        }
1532      }
1533    }
1534  }
1535  /*
1536  *the pair (S[i],p) enters B if the spoly != 0
1537  */
1538  /*-  compute the short s-polynomial -*/
1539  if (strat->fromT && !TEST_OPT_INTSTRATEGY)
1540    pNorm(p);
1541
1542  if ((strat->S[i]==NULL) || (p==NULL))
1543    return;
1544
1545  if ((strat->fromQ!=NULL) && (isFromQ!=0) && (strat->fromQ[i]!=0))
1546    Lp.p=NULL;
1547  else
1548  {
1549    #ifdef HAVE_PLURAL
1550    if ( rIsPluralRing(currRing) )
1551    {
1552      if(pHasNotCF(p, strat->S[i]))
1553      {
1554         if(ncRingType(currRing) == nc_lie)
1555         {
1556             // generalized prod-crit for lie-type
1557             strat->cp++;
1558             Lp.p = nc_p_Bracket_qq(pCopy(p),strat->S[i]);
1559         }
1560         else
1561        if( ALLOW_PROD_CRIT(strat) )
1562        {
1563            // product criterion for homogeneous case in SCA
1564            strat->cp++;
1565            Lp.p = NULL;
1566        }
1567        else
1568        {
1569          Lp.p = // nc_CreateSpoly(strat->S[i],p,currRing);
1570                nc_CreateShortSpoly(strat->S[i], p, currRing);
1571
1572          assume(pNext(Lp.p)==NULL); // TODO: this may be violated whenever ext.prod.crit. for Lie alg. is used
1573          pNext(Lp.p) = strat->tail; // !!!
1574        }
1575      }
1576      else
1577      {
1578        Lp.p = // nc_CreateSpoly(strat->S[i],p,currRing);
1579              nc_CreateShortSpoly(strat->S[i], p, currRing);
1580
1581        assume(pNext(Lp.p)==NULL); // TODO: this may be violated whenever ext.prod.crit. for Lie alg. is used
1582        pNext(Lp.p) = strat->tail; // !!!
1583
1584      }
1585
1586
1587#if MYTEST
1588      if (TEST_OPT_DEBUG)
1589      {
1590        PrintS("enterOnePairNormal::\n strat->S[i]: "); pWrite(strat->S[i]);
1591        PrintS("p: "); pWrite(p);
1592        PrintS("SPoly: "); pWrite(Lp.p);
1593      }
1594#endif
1595
1596    }
1597    else
1598    #endif
1599    {
1600      assume(!rIsPluralRing(currRing));
1601      Lp.p = ksCreateShortSpoly(strat->S[i], p, strat->tailRing);
1602#if MYTEST
1603      if (TEST_OPT_DEBUG)
1604      {
1605        PrintS("enterOnePairNormal::\n strat->S[i]: "); pWrite(strat->S[i]);
1606        PrintS("p: "); pWrite(p);
1607        PrintS("commutative SPoly: "); pWrite(Lp.p);
1608      }
1609#endif
1610
1611      }
1612  }
1613  if (Lp.p == NULL)
1614  {
1615    /*- the case that the s-poly is 0 -*/
1616    if (strat->pairtest==NULL) initPairtest(strat);
1617    strat->pairtest[i] = TRUE;/*- hint for spoly(S^[i],p)=0 -*/
1618    strat->pairtest[strat->sl+1] = TRUE;
1619    /*hint for spoly(S[i],p) == 0 for some i,0 <= i <= sl*/
1620    /*
1621    *suppose we have (s,r),(r,p),(s,p) and spoly(s,p) == 0 and (r,p) is
1622    *still in B (i.e. lcm(r,p) == lcm(s,p) or the leading term of s does not
1623    *devide lcm(r,p)). In the last case (s,r) can be canceled if the leading
1624    *term of p devides the lcm(s,r)
1625    *(this canceling should be done here because
1626    *the case lcm(s,p) == lcm(s,r) is not covered in chainCrit)
1627    *the first case is handeled in chainCrit
1628    */
1629    if (Lp.lcm!=NULL) pLmFree(Lp.lcm);
1630  }
1631  else
1632  {
1633    /*- the pair (S[i],p) enters B -*/
1634    Lp.p1 = strat->S[i];
1635    Lp.p2 = p;
1636
1637    if (
1638        (!rIsPluralRing(currRing))
1639//      ||  (rIsPluralRing(currRing) && (ncRingType(currRing) != nc_lie))
1640       )
1641    {
1642      assume(pNext(Lp.p)==NULL); // TODO: this may be violated whenever ext.prod.crit. for Lie alg. is used
1643      pNext(Lp.p) = strat->tail; // !!!
1644    }
1645
1646    if (atR >= 0)
1647    {
1648      Lp.i_r1 = strat->S_2_R[i];
1649      Lp.i_r2 = atR;
1650    }
1651    else
1652    {
1653      Lp.i_r1 = -1;
1654      Lp.i_r2 = -1;
1655    }
1656    strat->initEcartPair(&Lp,strat->S[i],p,strat->ecartS[i],ecart);
1657
1658    if (TEST_OPT_INTSTRATEGY)
1659    {
1660      if (!rIsPluralRing(currRing))
1661        nDelete(&(Lp.p->coef));
1662    }
1663
1664    l = strat->posInL(strat->B,strat->Bl,&Lp,strat);
1665    enterL(&strat->B,&strat->Bl,&strat->Bmax,Lp,l);
1666  }
1667}
1668
1669/*2
1670* put the pair (s[i],p) into the set L, ecart=ecart(p)
1671* in the case that s forms a SB of (s)
1672*/
1673void enterOnePairSpecial (int i,poly p,int ecart,kStrategy strat, int atR = -1)
1674{
1675  //PrintS("try ");wrp(strat->S[i]);PrintS(" and ");wrp(p);PrintLn();
1676  if(pHasNotCF(p,strat->S[i]))
1677  {
1678    //PrintS("prod-crit\n");
1679    if(ALLOW_PROD_CRIT(strat))
1680    {
1681      //PrintS("prod-crit\n");
1682      strat->cp++;
1683      return;
1684    }
1685  }
1686
1687  int      l,j,compare;
1688  LObject  Lp;
1689  Lp.i_r = -1;
1690
1691  Lp.lcm = pInit();
1692  pLcm(p,strat->S[i],Lp.lcm);
1693  pSetm(Lp.lcm);
1694  for(j = strat->Ll;j>=0;j--)
1695  {
1696    compare=pDivComp(strat->L[j].lcm,Lp.lcm);
1697    if ((compare==1) || (pLmEqual(strat->L[j].lcm,Lp.lcm)))
1698    {
1699      //PrintS("c3-crit\n");
1700      strat->c3++;
1701      pLmFree(Lp.lcm);
1702      return;
1703    }
1704    else if (compare ==-1)
1705    {
1706      //Print("c3-crit with L[%d]\n",j);
1707      deleteInL(strat->L,&strat->Ll,j,strat);
1708      strat->c3++;
1709    }
1710  }
1711  /*-  compute the short s-polynomial -*/
1712
1713  #ifdef HAVE_PLURAL
1714  if (rIsPluralRing(currRing))
1715  {
1716    Lp.p = nc_CreateShortSpoly(strat->S[i],p); // ??? strat->tailRing?
1717  }
1718  else
1719  #endif
1720    Lp.p = ksCreateShortSpoly(strat->S[i],p,strat->tailRing);
1721
1722  if (Lp.p == NULL)
1723  {
1724     //PrintS("short spoly==NULL\n");
1725     pLmFree(Lp.lcm);
1726  }
1727  else
1728  {
1729    /*- the pair (S[i],p) enters L -*/
1730    Lp.p1 = strat->S[i];
1731    Lp.p2 = p;
1732    if (atR >= 0)
1733    {
1734      Lp.i_r1 = strat->S_2_R[i];
1735      Lp.i_r2 = atR;
1736    }
1737    else
1738    {
1739      Lp.i_r1 = -1;
1740      Lp.i_r2 = -1;
1741    }
1742    assume(pNext(Lp.p) == NULL);
1743    pNext(Lp.p) = strat->tail;
1744    strat->initEcartPair(&Lp,strat->S[i],p,strat->ecartS[i],ecart);
1745    if (TEST_OPT_INTSTRATEGY)
1746    {
1747      nDelete(&(Lp.p->coef));
1748    }
1749    l = strat->posInL(strat->L,strat->Ll,&Lp,strat);
1750    //Print("-> L[%d]\n",l);
1751    enterL(&strat->L,&strat->Ll,&strat->Lmax,Lp,l);
1752  }
1753}
1754
1755/*2
1756* merge set B into L
1757*/
1758void kMergeBintoL(kStrategy strat)
1759{
1760  int j=strat->Ll+strat->Bl+1;
1761  if (j>strat->Lmax)
1762  {
1763    j=((j+setmaxLinc-1)/setmaxLinc)*setmaxLinc;
1764    strat->L = (LSet)omReallocSize(strat->L,strat->Lmax*sizeof(LObject),
1765                                 j*sizeof(LObject));
1766    strat->Lmax=j;
1767  }
1768  j = strat->Ll;
1769  int i;
1770  for (i=strat->Bl; i>=0; i--)
1771  {
1772    j = strat->posInL(strat->L,j,&(strat->B[i]),strat);
1773    enterL(&strat->L,&strat->Ll,&strat->Lmax,strat->B[i],j);
1774  }
1775  strat->Bl = -1;
1776}
1777/*2
1778*the pairset B of pairs of type (s[i],p) is complete now. It will be updated
1779*using the chain-criterion in B and L and enters B to L
1780*/
1781void chainCritNormal (poly p,int ecart,kStrategy strat)
1782{
1783  int i,j,l;
1784
1785  /*
1786  *pairtest[i] is TRUE if spoly(S[i],p) == 0.
1787  *In this case all elements in B such
1788  *that their lcm is divisible by the leading term of S[i] can be canceled
1789  */
1790  if (strat->pairtest!=NULL)
1791  {
1792    {
1793      /*- i.e. there is an i with pairtest[i]==TRUE -*/
1794      for (j=0; j<=strat->sl; j++)
1795      {
1796        if (strat->pairtest[j])
1797        {
1798          for (i=strat->Bl; i>=0; i--)
1799          {
1800            if (pDivisibleBy(strat->S[j],strat->B[i].lcm))
1801            {
1802              deleteInL(strat->B,&strat->Bl,i,strat);
1803              strat->c3++;
1804            }
1805          }
1806        }
1807      }
1808    }
1809    omFreeSize(strat->pairtest,(strat->sl+2)*sizeof(BOOLEAN));
1810    strat->pairtest=NULL;
1811  }
1812  if (strat->Gebauer || strat->fromT)
1813  {
1814    if (strat->sugarCrit)
1815    {
1816    /*
1817    *suppose L[j] == (s,r) and p/lcm(s,r)
1818    *and lcm(s,r)#lcm(s,p) and lcm(s,r)#lcm(r,p)
1819    *and in case the sugar is o.k. then L[j] can be canceled
1820    */
1821      for (j=strat->Ll; j>=0; j--)
1822      {
1823        if (sugarDivisibleBy(ecart,strat->L[j].ecart)
1824        && ((pNext(strat->L[j].p) == strat->tail) || (pOrdSgn==1))
1825        && pCompareChain(p,strat->L[j].p1,strat->L[j].p2,strat->L[j].lcm))
1826        {
1827          if (strat->L[j].p == strat->tail)
1828          {
1829              deleteInL(strat->L,&strat->Ll,j,strat);
1830              strat->c3++;
1831          }
1832        }
1833      }
1834      /*
1835      *this is GEBAUER-MOELLER:
1836      *in B all elements with the same lcm except the "best"
1837      *(i.e. the last one in B with this property) will be canceled
1838      */
1839      j = strat->Bl;
1840      loop /*cannot be changed into a for !!! */
1841      {
1842        if (j <= 0) break;
1843        i = j-1;
1844        loop
1845        {
1846          if (i <  0) break;
1847          if (pLmEqual(strat->B[j].lcm,strat->B[i].lcm))
1848          {
1849            strat->c3++;
1850            if (sugarDivisibleBy(strat->B[j].ecart,strat->B[i].ecart))
1851            {
1852              deleteInL(strat->B,&strat->Bl,i,strat);
1853              j--;
1854            }
1855            else
1856            {
1857              deleteInL(strat->B,&strat->Bl,j,strat);
1858              break;
1859            }
1860          }
1861          i--;
1862        }
1863        j--;
1864      }
1865    }
1866    else /*sugarCrit*/
1867    {
1868      /*
1869      *suppose L[j] == (s,r) and p/lcm(s,r)
1870      *and lcm(s,r)#lcm(s,p) and lcm(s,r)#lcm(r,p)
1871      *and in case the sugar is o.k. then L[j] can be canceled
1872      */
1873      for (j=strat->Ll; j>=0; j--)
1874      {
1875        if (pCompareChain(p,strat->L[j].p1,strat->L[j].p2,strat->L[j].lcm))
1876        {
1877          if ((pNext(strat->L[j].p) == strat->tail)||(pOrdSgn==1))
1878          {
1879            deleteInL(strat->L,&strat->Ll,j,strat);
1880            strat->c3++;
1881          }
1882        }
1883      }
1884      /*
1885      *this is GEBAUER-MOELLER:
1886      *in B all elements with the same lcm except the "best"
1887      *(i.e. the last one in B with this property) will be canceled
1888      */
1889      j = strat->Bl;
1890      loop   /*cannot be changed into a for !!! */
1891      {
1892        if (j <= 0) break;
1893        for(i=j-1; i>=0; i--)
1894        {
1895          if (pLmEqual(strat->B[j].lcm,strat->B[i].lcm))
1896          {
1897            strat->c3++;
1898            deleteInL(strat->B,&strat->Bl,i,strat);
1899            j--;
1900          }
1901        }
1902        j--;
1903      }
1904    }
1905    /*
1906    *the elements of B enter L
1907    */
1908    kMergeBintoL(strat);
1909  }
1910  else
1911  {
1912    for (j=strat->Ll; j>=0; j--)
1913    {
1914      if (pCompareChain(p,strat->L[j].p1,strat->L[j].p2,strat->L[j].lcm))
1915      {
1916        if ((pNext(strat->L[j].p) == strat->tail)||(pOrdSgn==1))
1917        {
1918          deleteInL(strat->L,&strat->Ll,j,strat);
1919          strat->c3++;
1920        }
1921      }
1922    }
1923    /*
1924    *this is our MODIFICATION of GEBAUER-MOELLER:
1925    *First the elements of B enter L,
1926    *then we fix a lcm and the "best" element in L
1927    *(i.e the last in L with this lcm and of type (s,p))
1928    *and cancel all the other elements of type (r,p) with this lcm
1929    *except the case the element (s,r) has also the same lcm
1930    *and is on the worst position with respect to (s,p) and (r,p)
1931    */
1932    /*
1933    *B enters to L/their order with respect to B is permutated for elements
1934    *B[i].p with the same leading term
1935    */
1936    kMergeBintoL(strat);
1937    j = strat->Ll;
1938    loop  /*cannot be changed into a for !!! */
1939    {
1940      if (j <= 0)
1941      {
1942        /*now L[0] cannot be canceled any more and the tail can be removed*/
1943        if (strat->L[0].p2 == strat->tail) strat->L[0].p2 = p;
1944        break;
1945      }
1946      if (strat->L[j].p2 == p)
1947      {
1948        i = j-1;
1949        loop
1950        {
1951          if (i < 0)  break;
1952          if ((strat->L[i].p2 == p) && pLmEqual(strat->L[j].lcm,strat->L[i].lcm))
1953          {
1954            /*L[i] could be canceled but we search for a better one to cancel*/
1955            strat->c3++;
1956            if (isInPairsetL(i-1,strat->L[j].p1,strat->L[i].p1,&l,strat)
1957            && (pNext(strat->L[l].p) == strat->tail)
1958            && (!pLmEqual(strat->L[i].p,strat->L[l].p))
1959            && pDivisibleBy(p,strat->L[l].lcm))
1960            {
1961              /*
1962              *"NOT equal(...)" because in case of "equal" the element L[l]
1963              *is "older" and has to be from theoretical point of view behind
1964              *L[i], but we do not want to reorder L
1965              */
1966              strat->L[i].p2 = strat->tail;
1967              /*
1968              *L[l] will be canceled, we cannot cancel L[i] later on,
1969              *so we mark it with "tail"
1970              */
1971              deleteInL(strat->L,&strat->Ll,l,strat);
1972              i--;
1973            }
1974            else
1975            {
1976              deleteInL(strat->L,&strat->Ll,i,strat);
1977            }
1978            j--;
1979          }
1980          i--;
1981        }
1982      }
1983      else if (strat->L[j].p2 == strat->tail)
1984      {
1985        /*now L[j] cannot be canceled any more and the tail can be removed*/
1986        strat->L[j].p2 = p;
1987      }
1988      j--;
1989    }
1990  }
1991}
1992#ifdef HAVE_RATGRING
1993void chainCritPart (poly p,int ecart,kStrategy strat)
1994{
1995  int i,j,l;
1996
1997  /*
1998  *pairtest[i] is TRUE if spoly(S[i],p) == 0.
1999  *In this case all elements in B such
2000  *that their lcm is divisible by the leading term of S[i] can be canceled
2001  */
2002  if (strat->pairtest!=NULL)
2003  {
2004    {
2005      /*- i.e. there is an i with pairtest[i]==TRUE -*/
2006      for (j=0; j<=strat->sl; j++)
2007      {
2008        if (strat->pairtest[j])
2009        {
2010          for (i=strat->Bl; i>=0; i--)
2011          {
2012            if (_p_LmDivisibleByPart(strat->S[j],currRing,
2013               strat->B[i].lcm,currRing,
2014               currRing->real_var_start,currRing->real_var_end))
2015            {
2016              if(TEST_OPT_DEBUG)
2017              {
2018                 Print("chain-crit-part: S[%d]=",j);
2019                 p_wrp(strat->S[j],currRing);
2020                 Print(" divide B[%d].lcm=",i);
2021                 p_wrp(strat->B[i].lcm,currRing);
2022                 PrintLn();
2023              }
2024              deleteInL(strat->B,&strat->Bl,i,strat);
2025              strat->c3++;
2026            }
2027          }
2028        }
2029      }
2030    }
2031    omFreeSize(strat->pairtest,(strat->sl+2)*sizeof(BOOLEAN));
2032    strat->pairtest=NULL;
2033  }
2034  if (strat->Gebauer || strat->fromT)
2035  {
2036    if (strat->sugarCrit)
2037    {
2038    /*
2039    *suppose L[j] == (s,r) and p/lcm(s,r)
2040    *and lcm(s,r)#lcm(s,p) and lcm(s,r)#lcm(r,p)
2041    *and in case the sugar is o.k. then L[j] can be canceled
2042    */
2043      for (j=strat->Ll; j>=0; j--)
2044      {
2045        if (sugarDivisibleBy(ecart,strat->L[j].ecart)
2046        && ((pNext(strat->L[j].p) == strat->tail) || (pOrdSgn==1))
2047        && pCompareChainPart(p,strat->L[j].p1,strat->L[j].p2,strat->L[j].lcm))
2048        {
2049          if (strat->L[j].p == strat->tail)
2050          {
2051              if(TEST_OPT_DEBUG)
2052              {
2053                 PrintS("chain-crit-part: pCompareChainPart p=");
2054                 p_wrp(p,currRing);
2055                 Print(" delete L[%d]",j);
2056                 p_wrp(strat->L[j].lcm,currRing);
2057                 PrintLn();
2058              }
2059              deleteInL(strat->L,&strat->Ll,j,strat);
2060              strat->c3++;
2061          }
2062        }
2063      }
2064      /*
2065      *this is GEBAUER-MOELLER:
2066      *in B all elements with the same lcm except the "best"
2067      *(i.e. the last one in B with this property) will be canceled
2068      */
2069      j = strat->Bl;
2070      loop /*cannot be changed into a for !!! */
2071      {
2072        if (j <= 0) break;
2073        i = j-1;
2074        loop
2075        {
2076          if (i <  0) break;
2077          if (pLmEqual(strat->B[j].lcm,strat->B[i].lcm))
2078          {
2079            strat->c3++;
2080            if (sugarDivisibleBy(strat->B[j].ecart,strat->B[i].ecart))
2081            {
2082              if(TEST_OPT_DEBUG)
2083              {
2084                 Print("chain-crit-part: sugar B[%d].lcm=",j);
2085                 p_wrp(strat->B[j].lcm,currRing);
2086                 Print(" delete B[%d]",i);
2087                 p_wrp(strat->B[i].lcm,currRing);
2088                 PrintLn();
2089              }
2090              deleteInL(strat->B,&strat->Bl,i,strat);
2091              j--;
2092            }
2093            else
2094            {
2095              if(TEST_OPT_DEBUG)
2096              {
2097                 Print("chain-crit-part: sugar B[%d].lcm=",i);
2098                 p_wrp(strat->B[i].lcm,currRing);
2099                 Print(" delete B[%d]",j);
2100                 p_wrp(strat->B[j].lcm,currRing);
2101                 PrintLn();
2102              }
2103              deleteInL(strat->B,&strat->Bl,j,strat);
2104              break;
2105            }
2106          }
2107          i--;
2108        }
2109        j--;
2110      }
2111    }
2112    else /*sugarCrit*/
2113    {
2114      /*
2115      *suppose L[j] == (s,r) and p/lcm(s,r)
2116      *and lcm(s,r)#lcm(s,p) and lcm(s,r)#lcm(r,p)
2117      *and in case the sugar is o.k. then L[j] can be canceled
2118      */
2119      for (j=strat->Ll; j>=0; j--)
2120      {
2121        if (pCompareChainPart(p,strat->L[j].p1,strat->L[j].p2,strat->L[j].lcm))
2122        {
2123          if ((pNext(strat->L[j].p) == strat->tail)||(pOrdSgn==1))
2124          {
2125              if(TEST_OPT_DEBUG)
2126              {
2127                 PrintS("chain-crit-part: sugar:pCompareChainPart p=");
2128                 p_wrp(p,currRing);
2129                 Print(" delete L[%d]",j);
2130                 p_wrp(strat->L[j].lcm,currRing);
2131                 PrintLn();
2132              }
2133            deleteInL(strat->L,&strat->Ll,j,strat);
2134            strat->c3++;
2135          }
2136        }
2137      }
2138      /*
2139      *this is GEBAUER-MOELLER:
2140      *in B all elements with the same lcm except the "best"
2141      *(i.e. the last one in B with this property) will be canceled
2142      */
2143      j = strat->Bl;
2144      loop   /*cannot be changed into a for !!! */
2145      {
2146        if (j <= 0) break;
2147        for(i=j-1; i>=0; i--)
2148        {
2149          if (pLmEqual(strat->B[j].lcm,strat->B[i].lcm))
2150          {
2151              if(TEST_OPT_DEBUG)
2152              {
2153                 Print("chain-crit-part: equal lcm B[%d].lcm=",j);
2154                 p_wrp(strat->B[j].lcm,currRing);
2155                 Print(" delete B[%d]\n",i);
2156              }
2157            strat->c3++;
2158            deleteInL(strat->B,&strat->Bl,i,strat);
2159            j--;
2160          }
2161        }
2162        j--;
2163      }
2164    }
2165    /*
2166    *the elements of B enter L
2167    */
2168    kMergeBintoL(strat);
2169  }
2170  else
2171  {
2172    for (j=strat->Ll; j>=0; j--)
2173    {
2174      if (pCompareChainPart(p,strat->L[j].p1,strat->L[j].p2,strat->L[j].lcm))
2175      {
2176        if ((pNext(strat->L[j].p) == strat->tail)||(pOrdSgn==1))
2177        {
2178              if(TEST_OPT_DEBUG)
2179              {
2180                 PrintS("chain-crit-part: pCompareChainPart p=");
2181                 p_wrp(p,currRing);
2182                 Print(" delete L[%d]",j);
2183                 p_wrp(strat->L[j].lcm,currRing);
2184                 PrintLn();
2185              }
2186          deleteInL(strat->L,&strat->Ll,j,strat);
2187          strat->c3++;
2188        }
2189      }
2190    }
2191    /*
2192    *this is our MODIFICATION of GEBAUER-MOELLER:
2193    *First the elements of B enter L,
2194    *then we fix a lcm and the "best" element in L
2195    *(i.e the last in L with this lcm and of type (s,p))
2196    *and cancel all the other elements of type (r,p) with this lcm
2197    *except the case the element (s,r) has also the same lcm
2198    *and is on the worst position with respect to (s,p) and (r,p)
2199    */
2200    /*
2201    *B enters to L/their order with respect to B is permutated for elements
2202    *B[i].p with the same leading term
2203    */
2204    kMergeBintoL(strat);
2205    j = strat->Ll;
2206    loop  /*cannot be changed into a for !!! */
2207    {
2208      if (j <= 0)
2209      {
2210        /*now L[0] cannot be canceled any more and the tail can be removed*/
2211        if (strat->L[0].p2 == strat->tail) strat->L[0].p2 = p;
2212        break;
2213      }
2214      if (strat->L[j].p2 == p)
2215      {
2216        i = j-1;
2217        loop
2218        {
2219          if (i < 0)  break;
2220          if ((strat->L[i].p2 == p) && pLmEqual(strat->L[j].lcm,strat->L[i].lcm))
2221          {
2222            /*L[i] could be canceled but we search for a better one to cancel*/
2223            strat->c3++;
2224            if (isInPairsetL(i-1,strat->L[j].p1,strat->L[i].p1,&l,strat)
2225            && (pNext(strat->L[l].p) == strat->tail)
2226            && (!pLmEqual(strat->L[i].p,strat->L[l].p))
2227            && _p_LmDivisibleByPart(p,currRing,
2228                           strat->L[l].lcm,currRing,
2229                           currRing->real_var_start, currRing->real_var_end))
2230
2231            {
2232              /*
2233              *"NOT equal(...)" because in case of "equal" the element L[l]
2234              *is "older" and has to be from theoretical point of view behind
2235              *L[i], but we do not want to reorder L
2236              */
2237              strat->L[i].p2 = strat->tail;
2238              /*
2239              *L[l] will be canceled, we cannot cancel L[i] later on,
2240              *so we mark it with "tail"
2241              */
2242              if(TEST_OPT_DEBUG)
2243              {
2244                 PrintS("chain-crit-part: divisible_by p=");
2245                 p_wrp(p,currRing);
2246                 Print(" delete L[%d]",l);
2247                 p_wrp(strat->L[l].lcm,currRing);
2248                 PrintLn();
2249              }
2250              deleteInL(strat->L,&strat->Ll,l,strat);
2251              i--;
2252            }
2253            else
2254            {
2255              if(TEST_OPT_DEBUG)
2256              {
2257                 PrintS("chain-crit-part: divisible_by(2) p=");
2258                 p_wrp(p,currRing);
2259                 Print(" delete L[%d]",i);
2260                 p_wrp(strat->L[i].lcm,currRing);
2261                 PrintLn();
2262              }
2263              deleteInL(strat->L,&strat->Ll,i,strat);
2264            }
2265            j--;
2266          }
2267          i--;
2268        }
2269      }
2270      else if (strat->L[j].p2 == strat->tail)
2271      {
2272        /*now L[j] cannot be canceled any more and the tail can be removed*/
2273        strat->L[j].p2 = p;
2274      }
2275      j--;
2276    }
2277  }
2278}
2279#endif
2280
2281/*2
2282*(s[0],h),...,(s[k],h) will be put to the pairset L
2283*/
2284void initenterpairs (poly h,int k,int ecart,int isFromQ,kStrategy strat, int atR = -1)
2285{
2286
2287  if ((strat->syzComp==0)
2288  || (pGetComp(h)<=strat->syzComp))
2289  {
2290    int j;
2291    BOOLEAN new_pair=FALSE;
2292
2293    if (pGetComp(h)==0)
2294    {
2295      /* for Q!=NULL: build pairs (f,q),(f1,f2), but not (q1,q2)*/
2296      if ((isFromQ)&&(strat->fromQ!=NULL))
2297      {
2298        for (j=0; j<=k; j++)
2299        {
2300          if (!strat->fromQ[j])
2301          {
2302            new_pair=TRUE;
2303            strat->enterOnePair(j,h,ecart,isFromQ,strat, atR);
2304          //Print("j:%d, Ll:%d\n",j,strat->Ll);
2305          }
2306        }
2307      }
2308      else
2309      {
2310        new_pair=TRUE;
2311        for (j=0; j<=k; j++)
2312        {
2313          strat->enterOnePair(j,h,ecart,isFromQ,strat, atR);
2314          //Print("j:%d, Ll:%d\n",j,strat->Ll);
2315        }
2316      }
2317    }
2318    else
2319    {
2320      for (j=0; j<=k; j++)
2321      {
2322        if ((pGetComp(h)==pGetComp(strat->S[j]))
2323        || (pGetComp(strat->S[j])==0))
2324        {
2325          new_pair=TRUE;
2326          strat->enterOnePair(j,h,ecart,isFromQ,strat, atR);
2327        //Print("j:%d, Ll:%d\n",j,strat->Ll);
2328        }
2329      }
2330    }
2331
2332    if (new_pair)
2333    {
2334#ifdef HAVE_RATGRING
2335      if (currRing->real_var_start>0)
2336        chainCritPart(h,ecart,strat);
2337      else
2338#endif
2339      strat->chainCrit(h,ecart,strat);
2340    }
2341  }
2342}
2343
2344#ifdef HAVE_RINGS
2345/*2
2346*the pairset B of pairs of type (s[i],p) is complete now. It will be updated
2347*using the chain-criterion in B and L and enters B to L
2348*/
2349void chainCritRing (poly p,int ecart,kStrategy strat)
2350{
2351  int i,j,l;
2352  /*
2353  *pairtest[i] is TRUE if spoly(S[i],p) == 0.
2354  *In this case all elements in B such
2355  *that their lcm is divisible by the leading term of S[i] can be canceled
2356  */
2357  if (strat->pairtest!=NULL)
2358  {
2359    {
2360      /*- i.e. there is an i with pairtest[i]==TRUE -*/
2361      for (j=0; j<=strat->sl; j++)
2362      {
2363        if (strat->pairtest[j])
2364        {
2365          for (i=strat->Bl; i>=0; i--)
2366          {
2367            if (pDivisibleBy(strat->S[j],strat->B[i].lcm))
2368            {
2369#ifdef KDEBUG
2370              if (TEST_OPT_DEBUG)
2371              {
2372                PrintS("--- chain criterion func chainCritRing type 1\n");
2373                PrintS("strat->S[j]:");
2374                wrp(strat->S[j]);
2375                PrintS("  strat->B[i].lcm:");
2376                wrp(strat->B[i].lcm);
2377                PrintLn();
2378              }
2379#endif
2380              deleteInL(strat->B,&strat->Bl,i,strat);
2381              strat->c3++;
2382            }
2383          }
2384        }
2385      }
2386    }
2387    omFreeSize(strat->pairtest,(strat->sl+2)*sizeof(BOOLEAN));
2388    strat->pairtest=NULL;
2389  }
2390  assume(!(strat->Gebauer || strat->fromT));
2391  for (j=strat->Ll; j>=0; j--)
2392  {
2393    if (strat->L[j].lcm != NULL && nDivBy(pGetCoeff(strat->L[j].lcm), pGetCoeff(p)))
2394    {
2395      if (pCompareChain(p,strat->L[j].p1,strat->L[j].p2,strat->L[j].lcm))
2396      {
2397        if ((pNext(strat->L[j].p) == strat->tail) || (pOrdSgn==1))
2398        {
2399          deleteInL(strat->L,&strat->Ll,j,strat);
2400          strat->c3++;
2401#ifdef KDEBUG
2402              if (TEST_OPT_DEBUG)
2403              {
2404                PrintS("--- chain criterion func chainCritRing type 2\n");
2405                PrintS("strat->L[j].p:");
2406                wrp(strat->L[j].p);
2407                PrintS("  p:");
2408                wrp(p);
2409                PrintLn();
2410              }
2411#endif
2412        }
2413      }
2414    }
2415  }
2416  /*
2417  *this is our MODIFICATION of GEBAUER-MOELLER:
2418  *First the elements of B enter L,
2419  *then we fix a lcm and the "best" element in L
2420  *(i.e the last in L with this lcm and of type (s,p))
2421  *and cancel all the other elements of type (r,p) with this lcm
2422  *except the case the element (s,r) has also the same lcm
2423  *and is on the worst position with respect to (s,p) and (r,p)
2424  */
2425  /*
2426  *B enters to L/their order with respect to B is permutated for elements
2427  *B[i].p with the same leading term
2428  */
2429  kMergeBintoL(strat);
2430  j = strat->Ll;
2431  loop  /*cannot be changed into a for !!! */
2432  {
2433    if (j <= 0)
2434    {
2435      /*now L[0] cannot be canceled any more and the tail can be removed*/
2436      if (strat->L[0].p2 == strat->tail) strat->L[0].p2 = p;
2437      break;
2438    }
2439    if (strat->L[j].p2 == p) // Was the element added from B?
2440    {
2441      i = j-1;
2442      loop
2443      {
2444        if (i < 0)  break;
2445        // Element is from B and has the same lcm as L[j]
2446        if ((strat->L[i].p2 == p) && nDivBy(pGetCoeff(strat->L[j].lcm), pGetCoeff(strat->L[i].lcm))
2447             && pLmEqual(strat->L[j].lcm,strat->L[i].lcm))
2448        {
2449          /*L[i] could be canceled but we search for a better one to cancel*/
2450          strat->c3++;
2451#ifdef KDEBUG
2452          if (TEST_OPT_DEBUG)
2453          {
2454            PrintS("--- chain criterion func chainCritRing type 3\n");
2455            PrintS("strat->L[j].lcm:");
2456            wrp(strat->L[j].lcm);
2457            PrintS("  strat->L[i].lcm:");
2458            wrp(strat->L[i].lcm);
2459            PrintLn();
2460          }
2461#endif
2462          if (isInPairsetL(i-1,strat->L[j].p1,strat->L[i].p1,&l,strat)
2463          && (pNext(strat->L[l].p) == strat->tail)
2464          && (!pLmEqual(strat->L[i].p,strat->L[l].p))
2465          && pDivisibleBy(p,strat->L[l].lcm))
2466          {
2467            /*
2468            *"NOT equal(...)" because in case of "equal" the element L[l]
2469            *is "older" and has to be from theoretical point of view behind
2470            *L[i], but we do not want to reorder L
2471            */
2472            strat->L[i].p2 = strat->tail;
2473            /*
2474            *L[l] will be canceled, we cannot cancel L[i] later on,
2475            *so we mark it with "tail"
2476            */
2477            deleteInL(strat->L,&strat->Ll,l,strat);
2478            i--;
2479          }
2480          else
2481          {
2482            deleteInL(strat->L,&strat->Ll,i,strat);
2483          }
2484          j--;
2485        }
2486        i--;
2487      }
2488    }
2489    else if (strat->L[j].p2 == strat->tail)
2490    {
2491      /*now L[j] cannot be canceled any more and the tail can be removed*/
2492      strat->L[j].p2 = p;
2493    }
2494    j--;
2495  }
2496}
2497#endif
2498
2499#ifdef HAVE_RINGS
2500long ind2(long arg)
2501{
2502  long ind = 0;
2503  if (arg <= 0) return 0;
2504  while (arg%2 == 0)
2505  {
2506    arg = arg / 2;
2507    ind++;
2508  }
2509  return ind;
2510}
2511
2512long ind_fact_2(long arg)
2513{
2514  long ind = 0;
2515  if (arg <= 0) return 0;
2516  if (arg%2 == 1) { arg--; }
2517  while (arg > 0)
2518  {
2519    ind += ind2(arg);
2520    arg = arg - 2;
2521  }
2522  return ind;
2523}
2524#endif
2525
2526#ifdef HAVE_VANIDEAL
2527long twoPow(long arg)
2528{
2529  return 1L << arg;
2530}
2531
2532/*2
2533* put the pair (p, f) in B and f in T
2534*/
2535void enterOneZeroPairRing (poly f, poly t_p, poly p, int ecart, kStrategy strat, int atR = -1)
2536{
2537  int      l,j,compare,compareCoeff;
2538  LObject  Lp;
2539
2540  if (strat->interred_flag) return;
2541#ifdef KDEBUG
2542  Lp.ecart=0; Lp.length=0;
2543#endif
2544  /*- computes the lcm(s[i],p) -*/
2545  Lp.lcm = pInit();
2546
2547  pLcm(p,f,Lp.lcm);
2548  pSetm(Lp.lcm);
2549  pSetCoeff(Lp.lcm, nLcm(pGetCoeff(p), pGetCoeff(f), currRing));
2550  assume(!strat->sugarCrit);
2551  assume(!strat->fromT);
2552  /*
2553  *the set B collects the pairs of type (S[j],p)
2554  *suppose (r,p) is in B and (s,p) is the new pair and lcm(s,p) != lcm(r,p)
2555  *if the leading term of s devides lcm(r,p) then (r,p) will be canceled
2556  *if the leading term of r devides lcm(s,p) then (s,p) will not enter B
2557  */
2558  for(j = strat->Bl;j>=0;j--)
2559  {
2560    compare=pDivCompRing(strat->B[j].lcm,Lp.lcm);
2561    compareCoeff = nDivComp(pGetCoeff(strat->B[j].lcm), pGetCoeff(Lp.lcm));
2562    if (compareCoeff == 0 || compare == compareCoeff)
2563    {
2564      if (compare == 1)
2565      {
2566        strat->c3++;
2567        pLmDelete(Lp.lcm);
2568        return;
2569      }
2570      else
2571      if (compare == -1)
2572      {
2573        deleteInL(strat->B,&strat->Bl,j,strat);
2574        strat->c3++;
2575      }
2576    }
2577    if (compare == pDivComp_EQUAL)
2578    {
2579      // Add hint for same LM and direction of LC (later) (TODO Oliver)
2580      if (compareCoeff == 1)
2581      {
2582        strat->c3++;
2583        pLmDelete(Lp.lcm);
2584        return;
2585      }
2586      else
2587      if (compareCoeff == -1)
2588      {
2589        deleteInL(strat->B,&strat->Bl,j,strat);
2590        strat->c3++;
2591      }
2592    }
2593  }
2594  /*
2595  *the pair (S[i],p) enters B if the spoly != 0
2596  */
2597  /*-  compute the short s-polynomial -*/
2598  if ((f==NULL) || (p==NULL)) return;
2599  pNorm(p);
2600  {
2601    Lp.p = ksCreateShortSpoly(f, p, strat->tailRing);
2602  }
2603  if (Lp.p == NULL) //deactivated, as we are adding pairs with zeropoly and not from S
2604  {
2605    /*- the case that the s-poly is 0 -*/
2606//    if (strat->pairtest==NULL) initPairtest(strat);
2607//    strat->pairtest[i] = TRUE;/*- hint for spoly(S^[i],p)=0 -*/
2608//    strat->pairtest[strat->sl+1] = TRUE;
2609    /*hint for spoly(S[i],p) == 0 for some i,0 <= i <= sl*/
2610    /*
2611    *suppose we have (s,r),(r,p),(s,p) and spoly(s,p) == 0 and (r,p) is
2612    *still in B (i.e. lcm(r,p) == lcm(s,p) or the leading term of s does not
2613    *devide lcm(r,p)). In the last case (s,r) can be canceled if the leading
2614    *term of p devides the lcm(s,r)
2615    *(this canceling should be done here because
2616    *the case lcm(s,p) == lcm(s,r) is not covered in chainCrit)
2617    *the first case is handeled in chainCrit
2618    */
2619    if (Lp.lcm!=NULL) pLmDelete(Lp.lcm);
2620  }
2621  else
2622  {
2623    /*- the pair (S[i],p) enters B -*/
2624    Lp.p1 = f;
2625    Lp.p2 = p;
2626
2627    pNext(Lp.p) = strat->tail;
2628
2629    LObject tmp_h(f, currRing, strat->tailRing);
2630    tmp_h.SetShortExpVector();
2631    strat->initEcart(&tmp_h);
2632    tmp_h.sev = pGetShortExpVector(tmp_h.p);
2633    tmp_h.t_p = t_p;
2634
2635    enterT(tmp_h, strat, strat->tl + 1);
2636
2637    if (atR >= 0)
2638    {
2639      Lp.i_r2 = atR;
2640      Lp.i_r1 = strat->tl;
2641    }
2642
2643    strat->initEcartPair(&Lp,f,p,0/*strat->ecartS[i]*/,ecart);     // Attention: TODO: break ecart
2644    l = strat->posInL(strat->B,strat->Bl,&Lp,strat);
2645    enterL(&strat->B, &strat->Bl, &strat->Bmax, Lp, l);
2646  }
2647}
2648
2649/* Helper for kCreateZeroPoly
2650 * enumerating the exponents
2651ring r = 2^2, (a, b, c), lp; ideal G0 = system("createG0"); ideal G = interred(G0); ncols(G0); ncols(G);
2652 */
2653
2654int nextZeroSimplexExponent (long exp[], long ind[], long cexp[], long cind[], long* cabsind, long step[], long bound, long N)
2655/* gives the next exponent from the set H_1 */
2656{
2657  long add = ind2(cexp[1] + 2);
2658  if ((*cabsind < bound) && (*cabsind - step[1] + add < bound))
2659  {
2660    cexp[1] += 2;
2661    cind[1] += add;
2662    *cabsind += add;
2663  }
2664  else
2665  {
2666    // cabsind >= habsind
2667    if (N == 1) return 0;
2668    int i = 1;
2669    while (exp[i] == cexp[i] && i <= N) i++;
2670    cexp[i] = exp[i];
2671    *cabsind -= cind[i];
2672    cind[i] = ind[i];
2673    step[i] = 500000;
2674    *cabsind += cind[i];
2675    // Print("in: %d\n", *cabsind);
2676    i += 1;
2677    if (i > N) return 0;
2678    do
2679    {
2680      step[1] = 500000;
2681      for (int j = i + 1; j <= N; j++)
2682      {
2683        if (step[1] > step[j]) step[1] = step[j];
2684      }
2685      add = ind2(cexp[i] + 2);
2686      if (*cabsind - step[1] + add >= bound)
2687      {
2688        cexp[i] = exp[i];
2689        *cabsind -= cind[i];
2690        cind[i] = ind[i];
2691        *cabsind += cind[i];
2692        step[i] = 500000;
2693        i += 1;
2694        if (i > N) return 0;
2695      }
2696      else step[1] = -1;
2697    } while (step[1] != -1);
2698    step[1] = 500000;
2699    cexp[i] += 2;
2700    cind[i] += add;
2701    *cabsind += add;
2702    if (add < step[i]) step[i] = add;
2703    for (i = 2; i <= N; i++)
2704    {
2705      if (step[1] > step[i]) step[1] = step[i];
2706    }
2707  }
2708  return 1;
2709}
2710
2711/*
2712 * Creates the zero Polynomial on position exp
2713 * long exp[] : exponent of leading term
2714 * cabsind    : total 2-ind of exp (if -1 will be computed)
2715 * poly* t_p  : will hold the LT in tailRing
2716 * leadRing   : ring for the LT
2717 * tailRing   : ring for the tail
2718 */
2719
2720poly kCreateZeroPoly(long exp[], long cabsind, poly* t_p, ring leadRing, ring tailRing)
2721{
2722
2723  poly zeroPoly = NULL;
2724
2725  number tmp1;
2726  poly tmp2, tmp3;
2727
2728  if (cabsind == -1)
2729  {
2730    cabsind = 0;
2731    for (int i = 1; i <= leadRing->N; i++)
2732    {
2733      cabsind += ind_fact_2(exp[i]);
2734    }
2735//    Print("cabsind: %d\n", cabsind);
2736  }
2737  if (cabsind < leadRing->ch)
2738  {
2739    zeroPoly = p_ISet(twoPow(leadRing->ch - cabsind), tailRing);
2740  }
2741  else
2742  {
2743    zeroPoly = p_ISet(1, tailRing);
2744  }
2745  for (int i = 1; i <= leadRing->N; i++)
2746  {
2747    for (long j = 1; j <= exp[i]; j++)
2748    {
2749      tmp1 = nInit(j);
2750      tmp2 = p_ISet(1, tailRing);
2751      p_SetExp(tmp2, i, 1, tailRing);
2752      p_Setm(tmp2, tailRing);
2753      if (nIsZero(tmp1))
2754      { // should nowbe obsolet, test ! TODO OLIVER
2755        zeroPoly = p_Mult_q(zeroPoly, tmp2, tailRing);
2756      }
2757      else
2758      {
2759        tmp3 = p_NSet(nCopy(tmp1), tailRing);
2760        zeroPoly = p_Mult_q(zeroPoly, p_Add_q(tmp3, tmp2, tailRing), tailRing);
2761      }
2762    }
2763  }
2764  tmp2 = p_NSet(nCopy(pGetCoeff(zeroPoly)), leadRing);
2765  for (int i = 1; i <= leadRing->N; i++)
2766  {
2767    pSetExp(tmp2, i, p_GetExp(zeroPoly, i, tailRing));
2768  }
2769  p_Setm(tmp2, leadRing);
2770  *t_p = zeroPoly;
2771  zeroPoly = pNext(zeroPoly);
2772  pNext(*t_p) = NULL;
2773  pNext(tmp2) = zeroPoly;
2774  return tmp2;
2775}
2776
2777// #define OLI_DEBUG
2778
2779/*
2780 * Generate the s-polynomial for the virtual set of zero-polynomials
2781 */
2782
2783void initenterzeropairsRing (poly p, int ecart, kStrategy strat, int atR)
2784{
2785  // Initialize
2786  long exp[50];            // The exponent of \hat{X} (basepoint)
2787  long cexp[50];           // The current exponent for iterating over all
2788  long ind[50];            // The power of 2 in the i-th component of exp
2789  long cind[50];           // analog for cexp
2790  long mult[50];           // How to multiply the elements of G
2791  long cabsind = 0;        // The abs. index of cexp, i.e. the sum of cind
2792  long habsind = 0;        // The abs. index of the coefficient of h
2793  long step[50];           // The last increases
2794  for (int i = 1; i <= currRing->N; i++)
2795  {
2796    exp[i] = p_GetExp(p, i, currRing);
2797    if (exp[i] & 1 != 0)
2798    {
2799      exp[i] = exp[i] - 1;
2800      mult[i] = 1;
2801    }
2802    cexp[i] = exp[i];
2803    ind[i] = ind_fact_2(exp[i]);
2804    cabsind += ind[i];
2805    cind[i] = ind[i];
2806    step[i] = 500000;
2807  }
2808  step[1] = 500000;
2809  habsind = ind2((long) p_GetCoeff(p, currRing));
2810  long bound = currRing->ch - habsind;
2811#ifdef OLI_DEBUG
2812  PrintS("-------------\npoly  :");
2813  wrp(p);
2814  Print("\nexp   : (%d, %d)\n", exp[1] + mult[1], exp[2] + mult[1]);
2815  Print("cexp  : (%d, %d)\n", cexp[1], cexp[2]);
2816  Print("cind  : (%d, %d)\n", cind[1], cind[2]);
2817  Print("bound : %d\n", bound);
2818  Print("cind  : %d\n", cabsind);
2819#endif
2820  if (cabsind == 0)
2821  {
2822    if (!(nextZeroSimplexExponent(exp, ind, cexp, cind, &cabsind, step, bound, currRing->N)))
2823    {
2824      return;
2825    }
2826  }
2827  // Now the whole simplex
2828  do
2829  {
2830    // Build s-polynomial
2831    // 2**ind-def * mult * g - exp-def * h
2832    poly t_p;
2833    poly zeroPoly = kCreateZeroPoly(cexp, cabsind, &t_p, currRing, strat->tailRing);
2834#ifdef OLI_DEBUG
2835    Print("%d, (%d, %d), ind = (%d, %d)\n", cabsind, cexp[1], cexp[2], cind[1], cind[2]);
2836    Print("zPoly : ");
2837    wrp(zeroPoly);
2838    Print("\n");
2839#endif
2840    enterOneZeroPairRing(zeroPoly, t_p, p, ecart, strat, atR);
2841  }
2842  while ( nextZeroSimplexExponent(exp, ind, cexp, cind, &cabsind, step, bound, currRing->N) );
2843}
2844
2845/*
2846 * Create the Groebner basis of the vanishing polynomials.
2847 */
2848
2849ideal createG0()
2850{
2851  // Initialize
2852  long exp[50];            // The exponent of \hat{X} (basepoint)
2853  long cexp[50];           // The current exponent for iterating over all
2854  long ind[50];            // The power of 2 in the i-th component of exp
2855  long cind[50];           // analog for cexp
2856  long mult[50];           // How to multiply the elements of G
2857  long cabsind = 0;        // The abs. index of cexp, i.e. the sum of cind
2858  long habsind = 0;        // The abs. index of the coefficient of h
2859  long step[50];           // The last increases
2860  for (int i = 1; i <= currRing->N; i++)
2861  {
2862    exp[i] = 0;
2863    cexp[i] = exp[i];
2864    ind[i] = 0;
2865    step[i] = 500000;
2866    cind[i] = ind[i];
2867  }
2868  long bound = currRing->ch;
2869  step[1] = 500000;
2870#ifdef OLI_DEBUG
2871  PrintS("-------------\npoly  :");
2872//  wrp(p);
2873  Print("\nexp   : (%d, %d)\n", exp[1] + mult[1], exp[2] + mult[1]);
2874  Print("cexp  : (%d, %d)\n", cexp[1], cexp[2]);
2875  Print("cind  : (%d, %d)\n", cind[1], cind[2]);
2876  Print("bound : %d\n", bound);
2877  Print("cind  : %d\n", cabsind);
2878#endif
2879  if (cabsind == 0)
2880  {
2881    if (!(nextZeroSimplexExponent(exp, ind, cexp, cind, &cabsind, step, bound, currRing->N)))
2882    {
2883      return idInit(1, 1);
2884    }
2885  }
2886  ideal G0 = idInit(1, 1);
2887  // Now the whole simplex
2888  do
2889  {
2890    // Build s-polynomial
2891    // 2**ind-def * mult * g - exp-def * h
2892    poly t_p;
2893    poly zeroPoly = kCreateZeroPoly(cexp, cabsind, &t_p, currRing, currRing);
2894#ifdef OLI_DEBUG
2895    Print("%d, (%d, %d), ind = (%d, %d)\n", cabsind, cexp[1], cexp[2], cind[1], cind[2]);
2896    Print("zPoly : ");
2897    wrp(zeroPoly);
2898    Print("\n");
2899#endif
2900    // Add to ideal
2901    pEnlargeSet(&(G0->m), IDELEMS(G0), 1);
2902    IDELEMS(G0) += 1;
2903    G0->m[IDELEMS(G0) - 1] = zeroPoly;
2904  }
2905  while ( nextZeroSimplexExponent(exp, ind, cexp, cind, &cabsind, step, bound, currRing->N) );
2906  idSkipZeroes(G0);
2907  return G0;
2908}
2909#endif
2910
2911#ifdef HAVE_RINGS
2912/*2
2913*(s[0],h),...,(s[k],h) will be put to the pairset L
2914*/
2915void initenterstrongPairs (poly h,int k,int ecart,int isFromQ,kStrategy strat, int atR = -1)
2916{
2917  const int iCompH = pGetComp(h);
2918  if (!nIsOne(pGetCoeff(h)))
2919  {
2920    int j;
2921    BOOLEAN new_pair=FALSE;
2922
2923    for (j=0; j<=k; j++)
2924    {
2925      // Print("j:%d, Ll:%d\n",j,strat->Ll);
2926//      if (((unsigned long) pGetCoeff(h) % (unsigned long) pGetCoeff(strat->S[j]) != 0) &&
2927//         ((unsigned long) pGetCoeff(strat->S[j]) % (unsigned long) pGetCoeff(h) != 0))
2928      if ( iCompH == pGetComp(strat->S[j]) )
2929      {
2930        {
2931          if (enterOneStrongPoly(j,h,ecart,isFromQ,strat, atR))
2932            new_pair=TRUE;
2933        }
2934      }
2935    }
2936  }
2937/*
2938ring r=256,(x,y,z),dp;
2939ideal I=12xz-133y, 2xy-z;
2940*/
2941
2942}
2943
2944/*2
2945* Generates spoly(0, h) if applicable. Assumes ring in Z/2^n.
2946*/
2947void enterExtendedSpoly(poly h,kStrategy strat)
2948{
2949  if (nIsOne(pGetCoeff(h))) return;
2950  number gcd;
2951  bool go = false;
2952  if (nDivBy((number) 0, pGetCoeff(h)))
2953  {
2954    gcd = nIntDiv((number) 0, pGetCoeff(h));
2955    go = true;
2956  }
2957  else
2958    gcd = nGcd((number) 0, pGetCoeff(h), strat->tailRing);
2959  if (go || !nIsOne(gcd))
2960  {
2961    poly p = h->next;
2962    if (!go)
2963    {
2964      number tmp = gcd;
2965      gcd = nIntDiv(0, gcd);
2966      nDelete(&tmp);
2967    }
2968    p_Test(p,strat->tailRing);
2969    p = pp_Mult_nn(p, gcd, strat->tailRing);
2970    nDelete(&gcd);
2971
2972    if (p != NULL)
2973    {
2974      if (TEST_OPT_PROT)
2975      {
2976        PrintS("Z");
2977      }
2978#ifdef KDEBUG
2979      if (TEST_OPT_DEBUG)
2980      {
2981        PrintS("--- create zero spoly: ");
2982        p_wrp(h,currRing,strat->tailRing);
2983        PrintS(" ---> ");
2984      }
2985#endif
2986      poly tmp = pInit();
2987      pSetCoeff0(tmp, pGetCoeff(p));
2988      for (int i = 1; i <= rVar(currRing); i++)
2989      {
2990        pSetExp(tmp, i, p_GetExp(p, i, strat->tailRing));
2991      }
2992      if (rRing_has_Comp(currRing) && rRing_has_Comp(strat->tailRing))
2993      {
2994        p_SetComp(tmp, p_GetComp(p, strat->tailRing), currRing);
2995      }
2996      p_Setm(tmp, currRing);
2997      p = p_LmFreeAndNext(p, strat->tailRing);
2998      pNext(tmp) = p;
2999      LObject h;
3000      h.Init();
3001      h.p = tmp;
3002      h.tailRing = strat->tailRing;
3003      int posx;
3004      if (h.p!=NULL)
3005      {
3006        if (TEST_OPT_INTSTRATEGY)
3007        {
3008          //pContent(h.p);
3009          h.pCleardenom(); // also does a pContent
3010        }
3011        else
3012        {
3013          h.pNorm();
3014        }
3015        strat->initEcart(&h);
3016        if (strat->Ll==-1)
3017          posx =0;
3018        else
3019          posx = strat->posInL(strat->L,strat->Ll,&h,strat);
3020        h.sev = pGetShortExpVector(h.p);
3021        if (strat->tailRing != currRing)
3022        {
3023          h.t_p = k_LmInit_currRing_2_tailRing(h.p, strat->tailRing);
3024        }
3025#ifdef KDEBUG
3026        if (TEST_OPT_DEBUG)
3027        {
3028          p_wrp(tmp,currRing,strat->tailRing);
3029          PrintLn();
3030        }
3031#endif
3032        enterL(&strat->L,&strat->Ll,&strat->Lmax,h,posx);
3033      }
3034    }
3035  }
3036  nDelete(&gcd);
3037}
3038
3039void clearSbatch (poly h,int k,int pos,kStrategy strat)
3040{
3041  int j = pos;
3042  if ( (!strat->fromT)
3043  && (1//(strat->syzComp==0)
3044    //||(pGetComp(h)<=strat->syzComp)))
3045  ))
3046  {
3047    // Print("start clearS k=%d, pos=%d, sl=%d\n",k,pos,strat->sl);
3048    unsigned long h_sev = pGetShortExpVector(h);
3049    loop
3050    {
3051      if (j > k) break;
3052      clearS(h,h_sev, &j,&k,strat);
3053      j++;
3054    }
3055    // Print("end clearS sl=%d\n",strat->sl);
3056  }
3057}
3058
3059/*2
3060* Generates a sufficient set of spolys (maybe just a finite generating
3061* set of the syzygys)
3062*/
3063void superenterpairs (poly h,int k,int ecart,int pos,kStrategy strat, int atR)
3064{
3065    assume (rField_is_Ring(currRing));
3066    // enter also zero divisor * poly, if this is non zero and of smaller degree
3067    if (!(rField_is_Domain(currRing))) enterExtendedSpoly(h, strat);
3068    initenterpairs(h, k, ecart, 0, strat, atR);
3069    initenterstrongPairs(h, k, ecart, 0, strat, atR);
3070    clearSbatch(h, k, pos, strat);
3071}
3072#endif
3073
3074/*2
3075*(s[0],h),...,(s[k],h) will be put to the pairset L(via initenterpairs)
3076*superfluous elements in S will be deleted
3077*/
3078void enterpairs (poly h,int k,int ecart,int pos,kStrategy strat, int atR)
3079{
3080  int j=pos;
3081
3082#ifdef HAVE_RINGS
3083  assume (!rField_is_Ring(currRing));
3084#endif
3085
3086  initenterpairs(h,k,ecart,0,strat, atR);
3087  if ( (!strat->fromT)
3088  && ((strat->syzComp==0)
3089    ||(pGetComp(h)<=strat->syzComp)))
3090  {
3091    //Print("start clearS k=%d, pos=%d, sl=%d\n",k,pos,strat->sl);
3092    unsigned long h_sev = pGetShortExpVector(h);
3093    loop
3094    {
3095      if (j > k) break;
3096      clearS(h,h_sev, &j,&k,strat);
3097      j++;
3098    }
3099    //Print("end clearS sl=%d\n",strat->sl);
3100  }
3101 // PrintS("end enterpairs\n");
3102}
3103
3104/*2
3105*(s[0],h),...,(s[k],h) will be put to the pairset L(via initenterpairs)
3106*superfluous elements in S will be deleted
3107*/
3108void enterpairsSpecial (poly h,int k,int ecart,int pos,kStrategy strat, int atR = -1)
3109{
3110  int j;
3111  const int iCompH = pGetComp(h);
3112
3113  for (j=0; j<=k; j++)
3114  {
3115    const int iCompSj = pGetComp(strat->S[j]);
3116    if ((iCompH==iCompSj)
3117        || (0==iCompH) // TODO: what about this case???
3118        || (0==iCompSj))
3119    {
3120      enterOnePairSpecial(j,h,ecart,strat, atR);
3121    }
3122  }
3123
3124  if (strat->noClearS) return;
3125
3126//   #ifdef HAVE_PLURAL
3127/*
3128  if (rIsPluralRing(currRing))
3129  {
3130    j=pos;
3131    loop
3132    {
3133      if (j > k) break;
3134
3135      if (pLmDivisibleBy(h, strat->S[j]))
3136      {
3137        deleteInS(j, strat);
3138        j--;
3139        k--;
3140      }
3141
3142      j++;
3143    }
3144  }
3145  else
3146*/
3147//   #endif // ??? Why was the following cancelation disabled for non-commutative rings?
3148  {
3149    j=pos;
3150    loop
3151    {
3152      unsigned long h_sev = pGetShortExpVector(h);
3153      if (j > k) break;
3154      clearS(h,h_sev,&j,&k,strat);
3155      j++;
3156    }
3157  }
3158}
3159
3160/*2
3161*reorders  s with respect to posInS,
3162*suc is the first changed index or zero
3163*/
3164
3165void reorderS (int* suc,kStrategy strat)
3166{
3167  int i,j,at,ecart, s2r;
3168  int fq=0;
3169  unsigned long sev;
3170  poly  p;
3171  int new_suc=strat->sl+1;
3172  i= *suc;
3173  if (i<0) i=0;
3174
3175  for (; i<=strat->sl; i++)
3176  {
3177    at = posInS(strat,i-1,strat->S[i],strat->ecartS[i]);
3178    if (at != i)
3179    {
3180      if (new_suc > at) new_suc = at;
3181      p = strat->S[i];
3182      ecart = strat->ecartS[i];
3183      sev = strat->sevS[i];
3184      s2r = strat->S_2_R[i];
3185      if (strat->fromQ!=NULL) fq=strat->fromQ[i];
3186      for (j=i; j>=at+1; j--)
3187      {
3188        strat->S[j] = strat->S[j-1];
3189        strat->ecartS[j] = strat->ecartS[j-1];
3190        strat->sevS[j] = strat->sevS[j-1];
3191        strat->S_2_R[j] = strat->S_2_R[j-1];
3192      }
3193      strat->S[at] = p;
3194      strat->ecartS[at] = ecart;
3195      strat->sevS[at] = sev;
3196      strat->S_2_R[at] = s2r;
3197      if (strat->fromQ!=NULL)
3198      {
3199        for (j=i; j>=at+1; j--)
3200        {
3201          strat->fromQ[j] = strat->fromQ[j-1];
3202        }
3203        strat->fromQ[at]=fq;
3204      }
3205    }
3206  }
3207  if (new_suc <= strat->sl) *suc=new_suc;
3208  else                      *suc=-1;
3209}
3210
3211
3212/*2
3213*looks up the position of p in set
3214*set[0] is the smallest with respect to the ordering-procedure deg/pComp
3215* Assumption: posInS only depends on the leading term
3216*             otherwise, bba has to be changed
3217*/
3218int posInS (const kStrategy strat, const int length,const poly p,
3219            const int ecart_p)
3220{
3221  if(length==-1) return 0;
3222  polyset set=strat->S;
3223  int i;
3224  int an = 0;
3225  int en = length;
3226  int cmp_int = pOrdSgn;
3227  int pc=pGetComp(p);
3228  if ((currRing->MixedOrder)
3229#ifdef HAVE_PLURAL
3230  && (currRing->real_var_start==0)
3231#endif
3232#if 0
3233  || ((strat->ak>0) && ((currRing->order[0]==ringorder_c)||((currRing->order[0]==ringorder_C))))
3234#endif
3235  )
3236  {
3237    int o=pWTotaldegree(p);
3238    int oo=pWTotaldegree(set[length]);
3239
3240    if ((oo<o)
3241    || ((o==oo) && (pLmCmp(set[length],p)!= cmp_int)))
3242      return length+1;
3243
3244    loop
3245    {
3246      if (an >= en-1)
3247      {
3248        if ((pWTotaldegree(set[an])>=o) && (pLmCmp(set[an],p) == cmp_int))
3249        {
3250          return an;
3251        }
3252        return en;
3253      }
3254      i=(an+en) / 2;
3255      if ((pWTotaldegree(set[i])>=o) && (pLmCmp(set[i],p) == cmp_int)) en=i;
3256      else                              an=i;
3257    }
3258  }
3259  else
3260  {
3261#ifdef HAVE_RINGS
3262    if (rField_is_Ring(currRing))
3263    {
3264      if (pLmCmp(set[length],p)== -cmp_int)
3265        return length+1;
3266      int cmp;
3267      loop
3268      {
3269        if (an >= en-1)
3270        {
3271          cmp = pLmCmp(set[an],p);
3272          if (cmp == cmp_int)  return an;
3273          if (cmp == -cmp_int) return en;
3274          if (nDivBy(pGetCoeff(p), pGetCoeff(set[an]))) return en;
3275          return an;
3276        }
3277        i = (an+en) / 2;
3278        cmp = pLmCmp(set[i],p);
3279        if (cmp == cmp_int)         en = i;
3280        else if (cmp == -cmp_int)   an = i;
3281        else
3282        {
3283          if (nDivBy(pGetCoeff(p), pGetCoeff(set[i]))) an = i;
3284          else en = i;
3285        }
3286      }
3287    }
3288    else
3289#endif
3290    if (pLmCmp(set[length],p)== -cmp_int)
3291      return length+1;
3292
3293    loop
3294    {
3295      if (an >= en-1)
3296      {
3297        if (pLmCmp(set[an],p) == cmp_int) return an;
3298        if (pLmCmp(set[an],p) == -cmp_int) return en;
3299        if ((cmp_int!=1)
3300        && ((strat->ecartS[an])>ecart_p))
3301          return an;
3302        return en;
3303      }
3304      i=(an+en) / 2;
3305      if (pLmCmp(set[i],p) == cmp_int) en=i;
3306      else if (pLmCmp(set[i],p) == -cmp_int) an=i;
3307      else
3308      {
3309        if ((cmp_int!=1)
3310        &&((strat->ecartS[i])<ecart_p))
3311          en=i;
3312        else
3313          an=i;
3314      }
3315    }
3316  }
3317}
3318
3319
3320/*2
3321* looks up the position of p in set
3322* the position is the last one
3323*/
3324int posInT0 (const TSet set,const int length,LObject &p)
3325{
3326  return (length+1);
3327}
3328
3329
3330/*2
3331* looks up the position of p in T
3332* set[0] is the smallest with respect to the ordering-procedure
3333* pComp
3334*/
3335int posInT1 (const TSet set,const int length,LObject &p)
3336{
3337  if (length==-1) return 0;
3338
3339  if (pLmCmp(set[length].p,p.p)!= pOrdSgn) return length+1;
3340
3341  int i;
3342  int an = 0;
3343  int en= length;
3344
3345  loop
3346  {
3347    if (an >= en-1)
3348    {
3349      if (pLmCmp(set[an].p,p.p) == pOrdSgn) return an;
3350      return en;
3351    }
3352    i=(an+en) / 2;
3353    if (pLmCmp(set[i].p,p.p) == pOrdSgn) en=i;
3354    else                                 an=i;
3355  }
3356}
3357
3358/*2
3359* looks up the position of p in T
3360* set[0] is the smallest with respect to the ordering-procedure
3361* length
3362*/
3363int posInT2 (const TSet set,const int length,LObject &p)
3364{
3365  p.GetpLength();
3366  if (length==-1)
3367    return 0;
3368  if (set[length].length<p.length)
3369    return length+1;
3370
3371  int i;
3372  int an = 0;
3373  int en= length;
3374
3375  loop
3376  {
3377    if (an >= en-1)
3378    {
3379      if (set[an].length>p.length) return an;
3380      return en;
3381    }
3382    i=(an+en) / 2;
3383    if (set[i].length>p.length) en=i;
3384    else                        an=i;
3385  }
3386}
3387
3388/*2
3389* looks up the position of p in T
3390* set[0] is the smallest with respect to the ordering-procedure
3391* totaldegree,pComp
3392*/
3393int posInT11 (const TSet set,const int length,LObject &p)
3394/*{
3395 * int j=0;
3396 * int o;
3397 *
3398 * o = p.GetpFDeg();
3399 * loop
3400 * {
3401 *   if ((pFDeg(set[j].p) > o)
3402 *   || ((pFDeg(set[j].p) == o) && (pLmCmp(set[j].p,p.p) == pOrdSgn)))
3403 *   {
3404 *     return j;
3405 *   }
3406 *   j++;
3407 *   if (j > length) return j;
3408 * }
3409 *}
3410 */
3411{
3412  if (length==-1) return 0;
3413
3414  int o = p.GetpFDeg();
3415  int op = set[length].GetpFDeg();
3416
3417  if ((op < o)
3418  || ((op == o) && (pLmCmp(set[length].p,p.p) != pOrdSgn)))
3419    return length+1;
3420
3421  int i;
3422  int an = 0;
3423  int en= length;
3424
3425  loop
3426  {
3427    if (an >= en-1)
3428    {
3429      op= set[an].GetpFDeg();
3430      if ((op > o)
3431      || (( op == o) && (pLmCmp(set[an].p,p.p) == pOrdSgn)))
3432        return an;
3433      return en;
3434    }
3435    i=(an+en) / 2;
3436    op = set[i].GetpFDeg();
3437    if (( op > o)
3438    || (( op == o) && (pLmCmp(set[i].p,p.p) == pOrdSgn)))
3439      en=i;
3440    else
3441      an=i;
3442  }
3443}
3444
3445/*2 Pos for rings T: Here I am
3446* looks up the position of p in T
3447* set[0] is the smallest with respect to the ordering-procedure
3448* totaldegree,pComp
3449*/
3450int posInTrg0 (const TSet set,const int length,LObject &p)
3451{
3452  if (length==-1) return 0;
3453  int o = p.GetpFDeg();
3454  int op = set[length].GetpFDeg();
3455  int i;
3456  int an = 0;
3457  int en = length;
3458  int cmp_int = pOrdSgn;
3459  if ((op < o) || (pLmCmp(set[length].p,p.p)== -cmp_int))
3460    return length+1;
3461  int cmp;
3462  loop
3463  {
3464    if (an >= en-1)
3465    {
3466      op = set[an].GetpFDeg();
3467      if (op > o) return an;
3468      if (op < 0) return en;
3469      cmp = pLmCmp(set[an].p,p.p);
3470      if (cmp == cmp_int)  return an;
3471      if (cmp == -cmp_int) return en;
3472      if (nGreater(pGetCoeff(p.p), pGetCoeff(set[an].p))) return en;
3473      return an;
3474    }
3475    i = (an + en) / 2;
3476    op = set[i].GetpFDeg();
3477    if (op > o)       en = i;
3478    else if (op < o)  an = i;
3479    else
3480    {
3481      cmp = pLmCmp(set[i].p,p.p);
3482      if (cmp == cmp_int)                                     en = i;
3483      else if (cmp == -cmp_int)                               an = i;
3484      else if (nGreater(pGetCoeff(p.p), pGetCoeff(set[i].p))) an = i;
3485      else                                                    en = i;
3486    }
3487  }
3488}
3489/*
3490  int o = p.GetpFDeg();
3491  int op = set[length].GetpFDeg();
3492
3493  if ((op < o)
3494  || ((op == o) && (pLmCmp(set[length].p,p.p) != pOrdSgn)))
3495    return length+1;
3496
3497  int i;
3498  int an = 0;
3499  int en= length;
3500
3501  loop
3502  {
3503    if (an >= en-1)
3504    {
3505      op= set[an].GetpFDeg();
3506      if ((op > o)
3507      || (( op == o) && (pLmCmp(set[an].p,p.p) == pOrdSgn)))
3508        return an;
3509      return en;
3510    }
3511    i=(an+en) / 2;
3512    op = set[i].GetpFDeg();
3513    if (( op > o)
3514    || (( op == o) && (pLmCmp(set[i].p,p.p) == pOrdSgn)))
3515      en=i;
3516    else
3517      an=i;
3518  }
3519}
3520  */
3521/*2
3522* looks up the position of p in T
3523* set[0] is the smallest with respect to the ordering-procedure
3524* totaldegree,pComp
3525*/
3526int posInT110 (const TSet set,const int length,LObject &p)
3527{
3528  p.GetpLength();
3529  if (length==-1) return 0;
3530
3531  int o = p.GetpFDeg();
3532  int op = set[length].GetpFDeg();
3533
3534  if (( op < o)
3535  || (( op == o) && (set[length].length<p.length))
3536  || (( op == o) && (set[length].length == p.length)
3537     && (pLmCmp(set[length].p,p.p) != pOrdSgn)))
3538    return length+1;
3539
3540  int i;
3541  int an = 0;
3542  int en= length;
3543  loop
3544  {
3545    if (an >= en-1)
3546    {
3547      op = set[an].GetpFDeg();
3548      if (( op > o)
3549      || (( op == o) && (set[an].length > p.length))
3550      || (( op == o) && (set[an].length == p.length)
3551         && (pLmCmp(set[an].p,p.p) == pOrdSgn)))
3552        return an;
3553      return en;
3554    }
3555    i=(an+en) / 2;
3556    op = set[i].GetpFDeg();
3557    if (( op > o)
3558    || (( op == o) && (set[i].length > p.length))
3559    || (( op == o) && (set[i].length == p.length)
3560       && (pLmCmp(set[i].p,p.p) == pOrdSgn)))
3561      en=i;
3562    else
3563      an=i;
3564  }
3565}
3566
3567/*2
3568* looks up the position of p in set
3569* set[0] is the smallest with respect to the ordering-procedure
3570* pFDeg
3571*/
3572int posInT13 (const TSet set,const int length,LObject &p)
3573{
3574  if (length==-1) return 0;
3575
3576  int o = p.GetpFDeg();
3577
3578  if (set[length].GetpFDeg() <= o)
3579    return length+1;
3580
3581  int i;
3582  int an = 0;
3583  int en= length;
3584  loop
3585  {
3586    if (an >= en-1)
3587    {
3588      if (set[an].GetpFDeg() > o)
3589        return an;
3590      return en;
3591    }
3592    i=(an+en) / 2;
3593    if (set[i].GetpFDeg() > o)
3594      en=i;
3595    else
3596      an=i;
3597  }
3598}
3599
3600// determines the position based on: 1.) Ecart 2.) pLength
3601int posInT_EcartpLength(const TSet set,const int length,LObject &p)
3602{
3603  int ol = p.GetpLength();
3604  if (length==-1) return 0;
3605
3606  int op=p.ecart;
3607
3608  int oo=set[length].ecart;
3609  if ((oo < op) || ((oo==op) && (set[length].length < ol)))
3610    return length+1;
3611
3612  int i;
3613  int an = 0;
3614  int en= length;
3615  loop
3616    {
3617      if (an >= en-1)
3618      {
3619        int oo=set[an].ecart;
3620        if((oo > op)
3621           || ((oo==op) && (set[an].pLength > ol)))
3622          return an;
3623        return en;
3624      }
3625      i=(an+en) / 2;
3626      int oo=set[i].ecart;
3627      if ((oo > op)
3628          || ((oo == op) && (set[i].pLength > ol)))
3629        en=i;
3630      else
3631        an=i;
3632    }
3633}
3634
3635/*2
3636* looks up the position of p in set
3637* set[0] is the smallest with respect to the ordering-procedure
3638* maximaldegree, pComp
3639*/
3640int posInT15 (const TSet set,const int length,LObject &p)
3641/*{
3642 *int j=0;
3643 * int o;
3644 *
3645 * o = p.GetpFDeg()+p.ecart;
3646 * loop
3647 * {
3648 *   if ((set[j].GetpFDeg()+set[j].ecart > o)
3649 *   || ((set[j].GetpFDeg()+set[j].ecart == o)
3650 *     && (pLmCmp(set[j].p,p.p) == pOrdSgn)))
3651 *   {
3652 *     return j;
3653 *   }
3654 *   j++;
3655 *   if (j > length) return j;
3656 * }
3657 *}
3658 */
3659{
3660  if (length==-1) return 0;
3661
3662  int o = p.GetpFDeg() + p.ecart;
3663  int op = set[length].GetpFDeg()+set[length].ecart;
3664
3665  if ((op < o)
3666  || ((op == o)
3667     && (pLmCmp(set[length].p,p.p) != pOrdSgn)))
3668    return length+1;
3669
3670  int i;
3671  int an = 0;
3672  int en= length;
3673  loop
3674  {
3675    if (an >= en-1)
3676    {
3677      op = set[an].GetpFDeg()+set[an].ecart;
3678      if (( op > o)
3679      || (( op  == o) && (pLmCmp(set[an].p,p.p) == pOrdSgn)))
3680        return an;
3681      return en;
3682    }
3683    i=(an+en) / 2;
3684    op = set[i].GetpFDeg()+set[i].ecart;
3685    if (( op > o)
3686    || (( op == o) && (pLmCmp(set[i].p,p.p) == pOrdSgn)))
3687      en=i;
3688    else
3689      an=i;
3690  }
3691}
3692
3693/*2
3694* looks up the position of p in set
3695* set[0] is the smallest with respect to the ordering-procedure
3696* pFDeg+ecart, ecart, pComp
3697*/
3698int posInT17 (const TSet set,const int length,LObject &p)
3699/*
3700*{
3701* int j=0;
3702* int  o;
3703*
3704*  o = p.GetpFDeg()+p.ecart;
3705*  loop
3706*  {
3707*    if ((pFDeg(set[j].p)+set[j].ecart > o)
3708*    || (((pFDeg(set[j].p)+set[j].ecart == o)
3709*      && (set[j].ecart < p.ecart)))
3710*    || ((pFDeg(set[j].p)+set[j].ecart == o)
3711*      && (set[j].ecart==p.ecart)
3712*      && (pLmCmp(set[j].p,p.p)==pOrdSgn)))
3713*      return j;
3714*    j++;
3715*    if (j > length) return j;
3716*  }
3717* }
3718*/
3719{
3720  if (length==-1) return 0;
3721
3722  int o = p.GetpFDeg() + p.ecart;
3723  int op = set[length].GetpFDeg()+set[length].ecart;
3724
3725  if ((op < o)
3726  || (( op == o) && (set[length].ecart > p.ecart))
3727  || (( op == o) && (set[length].ecart==p.ecart)
3728     && (pLmCmp(set[length].p,p.p) != pOrdSgn)))
3729    return length+1;
3730
3731  int i;
3732  int an = 0;
3733  int en= length;
3734  loop
3735  {
3736    if (an >= en-1)
3737    {
3738      op = set[an].GetpFDeg()+set[an].ecart;
3739      if (( op > o)
3740      || (( op == o) && (set[an].ecart < p.ecart))
3741      || (( op  == o) && (set[an].ecart==p.ecart)
3742         && (pLmCmp(set[an].p,p.p) == pOrdSgn)))
3743        return an;
3744      return en;
3745    }
3746    i=(an+en) / 2;
3747    op = set[i].GetpFDeg()+set[i].ecart;
3748    if ((op > o)
3749    || (( op == o) && (set[i].ecart < p.ecart))
3750    || (( op == o) && (set[i].ecart == p.ecart)
3751       && (pLmCmp(set[i].p,p.p) == pOrdSgn)))
3752      en=i;
3753    else
3754      an=i;
3755  }
3756}
3757/*2
3758* looks up the position of p in set
3759* set[0] is the smallest with respect to the ordering-procedure
3760* pGetComp, pFDeg+ecart, ecart, pComp
3761*/
3762int posInT17_c (const TSet set,const int length,LObject &p)
3763{
3764  if (length==-1) return 0;
3765
3766  int cc = (-1+2*currRing->order[0]==ringorder_c);
3767  /* cc==1 for (c,..), cc==-1 for (C,..) */
3768  int o = p.GetpFDeg() + p.ecart;
3769  int c = pGetComp(p.p)*cc;
3770
3771  if (pGetComp(set[length].p)*cc < c)
3772    return length+1;
3773  if (pGetComp(set[length].p)*cc == c)
3774  {
3775    int op = set[length].GetpFDeg()+set[length].ecart;
3776    if ((op < o)
3777    || ((op == o) && (set[length].ecart > p.ecart))
3778    || ((op == o) && (set[length].ecart==p.ecart)
3779       && (pLmCmp(set[length].p,p.p) != pOrdSgn)))
3780      return length+1;
3781  }
3782
3783  int i;
3784  int an = 0;
3785  int en= length;
3786  loop
3787  {
3788    if (an >= en-1)
3789    {
3790      if (pGetComp(set[an].p)*cc < c)
3791        return en;
3792      if (pGetComp(set[an].p)*cc == c)
3793      {
3794        int op = set[an].GetpFDeg()+set[an].ecart;
3795        if ((op > o)
3796        || ((op == o) && (set[an].ecart < p.ecart))
3797        || ((op == o) && (set[an].ecart==p.ecart)
3798           && (pLmCmp(set[an].p,p.p) == pOrdSgn)))
3799          return an;
3800      }
3801      return en;
3802    }
3803    i=(an+en) / 2;
3804    if (pGetComp(set[i].p)*cc > c)
3805      en=i;
3806    else if (pGetComp(set[i].p)*cc == c)
3807    {
3808      int op = set[i].GetpFDeg()+set[i].ecart;
3809      if ((op > o)
3810      || ((op == o) && (set[i].ecart < p.ecart))
3811      || ((op == o) && (set[i].ecart == p.ecart)
3812         && (pLmCmp(set[i].p,p.p) == pOrdSgn)))
3813        en=i;
3814      else
3815        an=i;
3816    }
3817    else
3818      an=i;
3819  }
3820}
3821
3822/*2
3823* looks up the position of p in set
3824* set[0] is the smallest with respect to
3825* ecart, pFDeg, length
3826*/
3827int posInT19 (const TSet set,const int length,LObject &p)
3828{
3829  p.GetpLength();
3830  if (length==-1) return 0;
3831
3832  int o = p.ecart;
3833  int op=p.GetpFDeg();
3834
3835  if (set[length].ecart < o)
3836    return length+1;
3837  if (set[length].ecart == o)
3838  {
3839     int oo=set[length].GetpFDeg();
3840     if ((oo < op) || ((oo==op) && (set[length].length < p.length)))
3841       return length+1;
3842  }
3843
3844  int i;
3845  int an = 0;
3846  int en= length;
3847  loop
3848  {
3849    if (an >= en-1)
3850    {
3851      if (set[an].ecart > o)
3852        return an;
3853      if (set[an].ecart == o)
3854      {
3855         int oo=set[an].GetpFDeg();
3856         if((oo > op)
3857         || ((oo==op) && (set[an].length > p.length)))
3858           return an;
3859      }
3860      return en;
3861    }
3862    i=(an+en) / 2;
3863    if (set[i].ecart > o)
3864      en=i;
3865    else if (set[i].ecart == o)
3866    {
3867       int oo=set[i].GetpFDeg();
3868       if ((oo > op)
3869       || ((oo == op) && (set[i].length > p.length)))
3870         en=i;
3871       else
3872        an=i;
3873    }
3874    else
3875      an=i;
3876  }
3877}
3878
3879/*2
3880*looks up the position of polynomial p in set
3881*set[length] is the smallest element in set with respect
3882*to the ordering-procedure pComp
3883*/
3884int posInLSpecial (const LSet set, const int length,
3885                   LObject *p,const kStrategy strat)
3886{
3887  if (length<0) return 0;
3888
3889  int d=p->GetpFDeg();
3890  int op=set[length].GetpFDeg();
3891
3892  if ((op > d)
3893  || ((op == d) && (p->p1!=NULL)&&(set[length].p1==NULL))
3894  || (pLmCmp(set[length].p,p->p)== pOrdSgn))
3895     return length+1;
3896
3897  int i;
3898  int an = 0;
3899  int en= length;
3900  loop
3901  {
3902    if (an >= en-1)
3903    {
3904      op=set[an].GetpFDeg();
3905      if ((op > d)
3906      || ((op == d) && (p->p1!=NULL) && (set[an].p1==NULL))
3907      || (pLmCmp(set[an].p,p->p)== pOrdSgn))
3908         return en;
3909      return an;
3910    }
3911    i=(an+en) / 2;
3912    op=set[i].GetpFDeg();
3913    if ((op>d)
3914    || ((op==d) && (p->p1!=NULL) && (set[i].p1==NULL))
3915    || (pLmCmp(set[i].p,p->p) == pOrdSgn))
3916      an=i;
3917    else
3918      en=i;
3919  }
3920}
3921
3922/*2
3923*looks up the position of polynomial p in set
3924*set[length] is the smallest element in set with respect
3925*to the ordering-procedure pComp
3926*/
3927int posInL0 (const LSet set, const int length,
3928             LObject* p,const kStrategy strat)
3929{
3930  if (length<0) return 0;
3931
3932  if (pLmCmp(set[length].p,p->p)== pOrdSgn)
3933    return length+1;
3934
3935  int i;
3936  int an = 0;
3937  int en= length;
3938  loop
3939  {
3940    if (an >= en-1)
3941    {
3942      if (pLmCmp(set[an].p,p->p) == pOrdSgn) return en;
3943      return an;
3944    }
3945    i=(an+en) / 2;
3946    if (pLmCmp(set[i].p,p->p) == pOrdSgn) an=i;
3947    else                                 en=i;
3948    /*aend. fuer lazy == in !=- machen */
3949  }
3950}
3951
3952/*2
3953* looks up the position of polynomial p in set
3954* e is the ecart of p
3955* set[length] is the smallest element in set with respect
3956* to the ordering-procedure totaldegree,pComp
3957*/
3958int posInL11 (const LSet set, const int length,
3959              LObject* p,const kStrategy strat)
3960/*{
3961 * int j=0;
3962 * int o;
3963 *
3964 * o = p->GetpFDeg();
3965 * loop
3966 * {
3967 *   if (j > length)            return j;
3968 *   if ((set[j].GetpFDeg() < o)) return j;
3969 *   if ((set[j].GetpFDeg() == o) && (pLmCmp(set[j].p,p->p) == -pOrdSgn))
3970 *   {
3971 *     return j;
3972 *   }
3973 *   j++;
3974 * }
3975 *}
3976 */
3977{
3978  if (length<0) return 0;
3979
3980  int o = p->GetpFDeg();
3981  int op = set[length].GetpFDeg();
3982
3983  if ((op > o)
3984  || ((op == o) && (pLmCmp(set[length].p,p->p) != -pOrdSgn)))
3985    return length+1;
3986  int i;
3987  int an = 0;
3988  int en= length;
3989  loop
3990  {
3991    if (an >= en-1)
3992    {
3993      op = set[an].GetpFDeg();
3994      if ((op > o)
3995      || ((op == o) && (pLmCmp(set[an].p,p->p) != -pOrdSgn)))
3996        return en;
3997      return an;
3998    }
3999    i=(an+en) / 2;
4000    op = set[i].GetpFDeg();
4001    if ((op > o)
4002    || ((op == o) && (pLmCmp(set[i].p,p->p) != -pOrdSgn)))
4003      an=i;
4004    else
4005      en=i;
4006  }
4007}
4008
4009/*2 Position for rings L: Here I am
4010* looks up the position of polynomial p in set
4011* e is the ecart of p
4012* set[length] is the smallest element in set with respect
4013* to the ordering-procedure totaldegree,pComp
4014*/
4015inline int getIndexRng(long coeff)
4016{
4017  if (coeff == 0) return -1;
4018  long tmp = coeff;
4019  int ind = 0;
4020  while (tmp % 2 == 0)
4021  {
4022    tmp = tmp / 2;
4023    ind++;
4024  }
4025  return ind;
4026}
4027
4028int posInLrg0 (const LSet set, const int length,
4029              LObject* p,const kStrategy strat)
4030/*          if (nGreater(pGetCoeff(p), pGetCoeff(set[an]))) return en;
4031        if (pLmCmp(set[i],p) == cmp_int)         en = i;
4032        else if (pLmCmp(set[i],p) == -cmp_int)   an = i;
4033        else
4034        {
4035          if (nGreater(pGetCoeff(p), pGetCoeff(set[i]))) an = i;
4036          else en = i;
4037        }*/
4038{
4039  if (length < 0) return 0;
4040
4041  int o = p->GetpFDeg();
4042  int op = set[length].GetpFDeg();
4043
4044  if ((op > o) || ((op == o) && (pLmCmp(set[length].p,p->p) != -pOrdSgn)))
4045    return length + 1;
4046  int i;
4047  int an = 0;
4048  int en = length;
4049  loop
4050  {
4051    if (an >= en - 1)
4052    {
4053      op = set[an].GetpFDeg();
4054      if ((op > o) || ((op == o) && (pLmCmp(set[an].p,p->p) != -pOrdSgn)))
4055        return en;
4056      return an;
4057    }
4058    i = (an+en) / 2;
4059    op = set[i].GetpFDeg();
4060    if ((op > o) || ((op == o) && (pLmCmp(set[i].p,p->p) != -pOrdSgn)))
4061      an = i;
4062    else
4063      en = i;
4064  }
4065}
4066
4067/*{
4068  if (length < 0) return 0;
4069
4070  int o = p->GetpFDeg();
4071  int op = set[length].GetpFDeg();
4072
4073  int inde = getIndexRng((unsigned long) pGetCoeff(set[length].p));
4074  int indp = getIndexRng((unsigned long) pGetCoeff(p->p));
4075  int inda;
4076  int indi;
4077
4078  if ((inda > indp) || ((inda == inde) && ((op > o) || ((op == o) && (pLmCmp(set[length].p,p->p) != -pOrdSgn)))))
4079    return length + 1;
4080  int i;
4081  int an = 0;
4082  inda = getIndexRng((unsigned long) pGetCoeff(set[an].p));
4083  int en = length;
4084  loop
4085  {
4086    if (an >= en-1)
4087    {
4088      op = set[an].GetpFDeg();
4089      if ((indp > inda) || ((indp == inda) && ((op > o) || ((op == o) && (pLmCmp(set[an].p,p->p) != -pOrdSgn)))))
4090        return en;
4091      return an;
4092    }
4093    i = (an + en) / 2;
4094    indi = getIndexRng((unsigned long) pGetCoeff(set[i].p));
4095    op = set[i].GetpFDeg();
4096    if ((indi > indp) || ((indi == indp) && ((op > o) || ((op == o) && (pLmCmp(set[i].p,p->p) != -pOrdSgn)))))
4097    // if ((op > o) || ((op == o) && (pLmCmp(set[i].p,p->p) != -pOrdSgn)))
4098    {
4099      an = i;
4100      inda = getIndexRng((unsigned long) pGetCoeff(set[an].p));
4101    }
4102    else
4103      en = i;
4104  }
4105} */
4106
4107/*2
4108* looks up the position of polynomial p in set
4109* set[length] is the smallest element in set with respect
4110* to the ordering-procedure totaldegree,pLength0
4111*/
4112int posInL110 (const LSet set, const int length,
4113               LObject* p,const kStrategy strat)
4114{
4115  if (length<0) return 0;
4116
4117  int o = p->GetpFDeg();
4118  int op = set[length].GetpFDeg();
4119
4120  if ((op > o)
4121  || ((op == o) && (set[length].length >p->length))
4122  || ((op == o) && (set[length].length <= p->length)
4123     && (pLmCmp(set[length].p,p->p) != -pOrdSgn)))
4124    return length+1;
4125  int i;
4126  int an = 0;
4127  int en= length;
4128  loop
4129  {
4130    if (an >= en-1)
4131    {
4132      op = set[an].GetpFDeg();
4133      if ((op > o)
4134      || ((op == o) && (set[an].length >p->length))
4135      || ((op == o) && (set[an].length <=p->length)
4136         && (pLmCmp(set[an].p,p->p) != -pOrdSgn)))
4137        return en;
4138      return an;
4139    }
4140    i=(an+en) / 2;
4141    op = set[i].GetpFDeg();
4142    if ((op > o)
4143    || ((op == o) && (set[i].length > p->length))
4144    || ((op == o) && (set[i].length <= p->length)
4145       && (pLmCmp(set[i].p,p->p) != -pOrdSgn)))
4146      an=i;
4147    else
4148      en=i;
4149  }
4150}
4151
4152/*2
4153* looks up the position of polynomial p in set
4154* e is the ecart of p
4155* set[length] is the smallest element in set with respect
4156* to the ordering-procedure totaldegree
4157*/
4158int posInL13 (const LSet set, const int length,
4159              LObject* p,const kStrategy strat)
4160{
4161  if (length<0) return 0;
4162
4163  int o = p->GetpFDeg();
4164
4165  if (set[length].GetpFDeg() > o)
4166    return length+1;
4167
4168  int i;
4169  int an = 0;
4170  int en= length;
4171  loop
4172  {
4173    if (an >= en-1)
4174    {
4175      if (set[an].GetpFDeg() >= o)
4176        return en;
4177      return an;
4178    }
4179    i=(an+en) / 2;
4180    if (set[i].GetpFDeg() >= o)
4181      an=i;
4182    else
4183      en=i;
4184  }
4185}
4186
4187/*2
4188* looks up the position of polynomial p in set
4189* e is the ecart of p
4190* set[length] is the smallest element in set with respect
4191* to the ordering-procedure maximaldegree,pComp
4192*/
4193int posInL15 (const LSet set, const int length,
4194              LObject* p,const kStrategy strat)
4195/*{
4196 * int j=0;
4197 * int o;
4198 *
4199 * o = p->ecart+p->GetpFDeg();
4200 * loop
4201 * {
4202 *   if (j > length)                       return j;
4203 *   if (set[j].GetpFDeg()+set[j].ecart < o) return j;
4204 *   if ((set[j].GetpFDeg()+set[j].ecart == o)
4205 *   && (pLmCmp(set[j].p,p->p) == -pOrdSgn))
4206 *   {
4207 *     return j;
4208 *   }
4209 *   j++;
4210 * }
4211 *}
4212 */
4213{
4214  if (length<0) return 0;
4215
4216  int o = p->GetpFDeg() + p->ecart;
4217  int op = set[length].GetpFDeg() + set[length].ecart;
4218
4219  if ((op > o)
4220  || ((op == o) && (pLmCmp(set[length].p,p->p) != -pOrdSgn)))
4221    return length+1;
4222  int i;
4223  int an = 0;
4224  int en= length;
4225  loop
4226  {
4227    if (an >= en-1)
4228    {
4229      op = set[an].GetpFDeg() + set[an].ecart;
4230      if ((op > o)
4231      || ((op == o) && (pLmCmp(set[an].p,p->p) != -pOrdSgn)))
4232        return en;
4233      return an;
4234    }
4235    i=(an+en) / 2;
4236    op = set[i].GetpFDeg() + set[i].ecart;
4237    if ((op > o)
4238    || ((op == o) && (pLmCmp(set[i].p,p->p) != -pOrdSgn)))
4239      an=i;
4240    else
4241      en=i;
4242  }
4243}
4244
4245/*2
4246* looks up the position of polynomial p in set
4247* e is the ecart of p
4248* set[length] is the smallest element in set with respect
4249* to the ordering-procedure totaldegree
4250*/
4251int posInL17 (const LSet set, const int length,
4252              LObject* p,const kStrategy strat)
4253{
4254  if (length<0) return 0;
4255
4256  int o = p->GetpFDeg() + p->ecart;
4257
4258  if ((set[length].GetpFDeg() + set[length].ecart > o)
4259  || ((set[length].GetpFDeg() + set[length].ecart == o)
4260     && (set[length].ecart > p->ecart))
4261  || ((set[length].GetpFDeg() + set[length].ecart == o)
4262     && (set[length].ecart == p->ecart)
4263     && (pLmCmp(set[length].p,p->p) != -pOrdSgn)))
4264    return length+1;
4265  int i;
4266  int an = 0;
4267  int en= length;
4268  loop
4269  {
4270    if (an >= en-1)
4271    {
4272      if ((set[an].GetpFDeg() + set[an].ecart > o)
4273      || ((set[an].GetpFDeg() + set[an].ecart == o)
4274         && (set[an].ecart > p->ecart))
4275      || ((set[an].GetpFDeg() + set[an].ecart == o)
4276         && (set[an].ecart == p->ecart)
4277         && (pLmCmp(set[an].p,p->p) != -pOrdSgn)))
4278        return en;
4279      return an;
4280    }
4281    i=(an+en) / 2;
4282    if ((set[i].GetpFDeg() + set[i].ecart > o)
4283    || ((set[i].GetpFDeg() + set[i].ecart == o)
4284       && (set[i].ecart > p->ecart))
4285    || ((set[i].GetpFDeg() +set[i].ecart == o)
4286       && (set[i].ecart == p->ecart)
4287       && (pLmCmp(set[i].p,p->p) != -pOrdSgn)))
4288      an=i;
4289    else
4290      en=i;
4291  }
4292}
4293/*2
4294* looks up the position of polynomial p in set
4295* e is the ecart of p
4296* set[length] is the smallest element in set with respect
4297* to the ordering-procedure pComp
4298*/
4299int posInL17_c (const LSet set, const int length,
4300                LObject* p,const kStrategy strat)
4301{
4302  if (length<0) return 0;
4303
4304  int cc = (-1+2*currRing->order[0]==ringorder_c);
4305  /* cc==1 for (c,..), cc==-1 for (C,..) */
4306  int c = pGetComp(p->p)*cc;
4307  int o = p->GetpFDeg() + p->ecart;
4308
4309  if (pGetComp(set[length].p)*cc > c)
4310    return length+1;
4311  if (pGetComp(set[length].p)*cc == c)
4312  {
4313    if ((set[length].GetpFDeg() + set[length].ecart > o)
4314    || ((set[length].GetpFDeg() + set[length].ecart == o)
4315       && (set[length].ecart > p->ecart))
4316    || ((set[length].GetpFDeg() + set[length].ecart == o)
4317       && (set[length].ecart == p->ecart)
4318       && (pLmCmp(set[length].p,p->p) != -pOrdSgn)))
4319      return length+1;
4320  }
4321  int i;
4322  int an = 0;
4323  int en= length;
4324  loop
4325  {
4326    if (an >= en-1)
4327    {
4328      if (pGetComp(set[an].p)*cc > c)
4329        return en;
4330      if (pGetComp(set[an].p)*cc == c)
4331      {
4332        if ((set[an].GetpFDeg() + set[an].ecart > o)
4333        || ((set[an].GetpFDeg() + set[an].ecart == o)
4334           && (set[an].ecart > p->ecart))
4335        || ((set[an].GetpFDeg() + set[an].ecart == o)
4336           && (set[an].ecart == p->ecart)
4337           && (pLmCmp(set[an].p,p->p) != -pOrdSgn)))
4338          return en;
4339      }
4340      return an;
4341    }
4342    i=(an+en) / 2;
4343    if (pGetComp(set[i].p)*cc > c)
4344      an=i;
4345    else if (pGetComp(set[i].p)*cc == c)
4346    {
4347      if ((set[i].GetpFDeg() + set[i].ecart > o)
4348      || ((set[i].GetpFDeg() + set[i].ecart == o)
4349         && (set[i].ecart > p->ecart))
4350      || ((set[i].GetpFDeg() +set[i].ecart == o)
4351         && (set[i].ecart == p->ecart)
4352         && (pLmCmp(set[i].p,p->p) != -pOrdSgn)))
4353        an=i;
4354      else
4355        en=i;
4356    }
4357    else
4358      en=i;
4359  }
4360}
4361
4362/***************************************************************
4363 *
4364 * Tail reductions
4365 *
4366 ***************************************************************/
4367TObject*
4368kFindDivisibleByInS(kStrategy strat, int pos, LObject* L, TObject *T,
4369                    long ecart)
4370{
4371  int j = 0;
4372  const unsigned long not_sev = ~L->sev;
4373  const unsigned long* sev = strat->sevS;
4374  poly p;
4375  ring r;
4376  L->GetLm(p, r);
4377
4378  assume(~not_sev == p_GetShortExpVector(p, r));
4379
4380  if (r == currRing)
4381  {
4382    loop
4383    {
4384      if (j > pos) return NULL;
4385#if defined(PDEBUG) || defined(PDIV_DEBUG)
4386      if (p_LmShortDivisibleBy(strat->S[j], sev[j], p, not_sev, r) &&
4387          (ecart== LONG_MAX || ecart>= strat->ecartS[j]))
4388        break;
4389#else
4390      if (!(sev[j] & not_sev) &&
4391          (ecart== LONG_MAX || ecart>= strat->ecartS[j]) &&
4392          p_LmDivisibleBy(strat->S[j], p, r))
4393        break;
4394
4395#endif
4396      j++;
4397    }
4398    // if called from NF, T objects do not exist:
4399    if (strat->tl < 0 || strat->S_2_R[j] == -1)
4400    {
4401      T->Set(strat->S[j], r, strat->tailRing);
4402      return T;
4403    }
4404    else
4405    {
4406/////      assume (j >= 0 && j <= strat->tl && strat->S_2_T(j) != NULL
4407/////      && strat->S_2_T(j)->p == strat->S[j]); // wrong?
4408//      assume (j >= 0 && j <= strat->sl && strat->S_2_T(j) != NULL && strat->S_2_T(j)->p == strat->S[j]);
4409      return strat->S_2_T(j);
4410    }
4411  }
4412  else
4413  {
4414    TObject* t;
4415    loop
4416    {
4417      if (j > pos) return NULL;
4418      assume(strat->S_2_R[j] != -1);
4419#if defined(PDEBUG) || defined(PDIV_DEBUG)
4420      t = strat->S_2_T(j);
4421      assume(t != NULL && t->t_p != NULL && t->tailRing == r);
4422      if (p_LmShortDivisibleBy(t->t_p, sev[j], p, not_sev, r) &&
4423          (ecart== LONG_MAX || ecart>= strat->ecartS[j]))
4424        return t;
4425#else
4426      if (! (sev[j] & not_sev) && (ecart== LONG_MAX || ecart>= strat->ecartS[j]))
4427      {
4428        t = strat->S_2_T(j);
4429        assume(t != NULL && t->t_p != NULL && t->tailRing == r && t->p == strat->S[j]);
4430        if (p_LmDivisibleBy(t->t_p, p, r)) return t;
4431      }
4432#endif
4433      j++;
4434    }
4435  }
4436}
4437
4438poly redtail (LObject* L, int pos, kStrategy strat)
4439{
4440  poly h, hn;
4441  int j;
4442  unsigned long not_sev;
4443  strat->redTailChange=FALSE;
4444
4445  poly p = L->p;
4446  if (strat->noTailReduction || pNext(p) == NULL)
4447    return p;
4448
4449  LObject Ln(strat->tailRing);
4450  TObject* With;
4451  // placeholder in case strat->tl < 0
4452  TObject  With_s(strat->tailRing);
4453  h = p;
4454  hn = pNext(h);
4455  long op = strat->tailRing->pFDeg(hn, strat->tailRing);
4456  long e;
4457  int l;
4458  BOOLEAN save_HE=strat->kHEdgeFound;
4459  strat->kHEdgeFound |=
4460    ((Kstd1_deg>0) && (op<=Kstd1_deg)) || TEST_OPT_INFREDTAIL;
4461
4462  while(hn != NULL)
4463  {
4464    op = strat->tailRing->pFDeg(hn, strat->tailRing);
4465    if ((Kstd1_deg>0)&&(op>Kstd1_deg)) goto all_done;
4466    e = strat->tailRing->pLDeg(hn, &l, strat->tailRing) - op;
4467    loop
4468    {
4469      Ln.Set(hn, strat->tailRing);
4470      Ln.sev = p_GetShortExpVector(hn, strat->tailRing);
4471      if (strat->kHEdgeFound)
4472        With = kFindDivisibleByInS(strat, pos, &Ln, &With_s);
4473      else
4474        With = kFindDivisibleByInS(strat, pos, &Ln, &With_s, e);
4475      if (With == NULL) break;
4476      With->length=0;
4477      With->pLength=0;
4478      strat->redTailChange=TRUE;
4479      if (ksReducePolyTail(L, With, h, strat->kNoetherTail()))
4480      {
4481        // reducing the tail would violate the exp bound
4482        if (kStratChangeTailRing(strat, L))
4483        {
4484          strat->kHEdgeFound = save_HE;
4485          return redtail(L, pos, strat);
4486        }
4487        else
4488          return NULL;
4489      }
4490      hn = pNext(h);
4491      if (hn == NULL) goto all_done;
4492      op = strat->tailRing->pFDeg(hn, strat->tailRing);
4493      if ((Kstd1_deg>0)&&(op>Kstd1_deg)) goto all_done;
4494      e = strat->tailRing->pLDeg(hn, &l, strat->tailRing) - op;
4495    }
4496    h = hn;
4497    hn = pNext(h);
4498  }
4499
4500  all_done:
4501  if (strat->redTailChange)
4502  {
4503    L->last = NULL;
4504    L->pLength = 0;
4505  }
4506  strat->kHEdgeFound = save_HE;
4507  return p;
4508}
4509
4510poly redtail (poly p, int pos, kStrategy strat)
4511{
4512  LObject L(p, currRing);
4513  return redtail(&L, pos, strat);
4514}
4515
4516poly redtailBba (LObject* L, int pos, kStrategy strat, BOOLEAN withT, BOOLEAN normalize)
4517{
4518#define REDTAIL_CANONICALIZE 100
4519  strat->redTailChange=FALSE;
4520  if (strat->noTailReduction) return L->GetLmCurrRing();
4521  poly h, p;
4522  p = h = L->GetLmTailRing();
4523  if ((h==NULL) || (pNext(h)==NULL))
4524    return L->GetLmCurrRing();
4525
4526  TObject* With;
4527  // placeholder in case strat->tl < 0
4528  TObject  With_s(strat->tailRing);
4529
4530  LObject Ln(pNext(h), strat->tailRing);
4531  Ln.pLength = L->GetpLength() - 1;
4532
4533  pNext(h) = NULL;
4534  if (L->p != NULL) pNext(L->p) = NULL;
4535  L->pLength = 1;
4536
4537  Ln.PrepareRed(strat->use_buckets);
4538
4539  int cnt=REDTAIL_CANONICALIZE;
4540  while(!Ln.IsNull())
4541  {
4542    loop
4543    {
4544      Ln.SetShortExpVector();
4545      if (withT)
4546      {
4547        int j;
4548        j = kFindDivisibleByInT(strat->T, strat->sevT, strat->tl, &Ln);
4549        if (j < 0) break;
4550        With = &(strat->T[j]);
4551      }
4552      else
4553      {
4554        With = kFindDivisibleByInS(strat, pos, &Ln, &With_s);
4555        if (With == NULL) break;
4556      }
4557      cnt--;
4558      if (cnt==0)
4559      {
4560        cnt=REDTAIL_CANONICALIZE;
4561        poly tmp=Ln.CanonicalizeP();
4562        if (normalize)
4563        {
4564          Ln.Normalize();
4565          //pNormalize(tmp);
4566          //if (TEST_OPT_PROT) { PrintS("n"); mflush(); }
4567        }
4568      }
4569      if (normalize && (!TEST_OPT_INTSTRATEGY) && (!nIsOne(pGetCoeff(With->p))))
4570      {
4571        With->pNorm();
4572      }
4573      strat->redTailChange=TRUE;
4574      if (ksReducePolyTail(L, With, &Ln))
4575      {
4576        // reducing the tail would violate the exp bound
4577        //  set a flag and hope for a retry (in bba)
4578        strat->completeReduce_retry=TRUE;
4579        if ((Ln.p != NULL) && (Ln.t_p != NULL)) Ln.p=NULL;
4580        do
4581        {
4582          pNext(h) = Ln.LmExtractAndIter();
4583          pIter(h);
4584          L->pLength++;
4585        } while (!Ln.IsNull());
4586        goto all_done;
4587      }
4588      if (Ln.IsNull()) goto all_done;
4589      if (! withT) With_s.Init(currRing);
4590    }
4591    pNext(h) = Ln.LmExtractAndIter();
4592    pIter(h);
4593    pNormalize(h);
4594    L->pLength++;
4595  }
4596
4597  all_done:
4598  Ln.Delete();
4599  if (L->p != NULL) pNext(L->p) = pNext(p);
4600
4601  if (strat->redTailChange)
4602  {
4603    L->last = NULL;
4604    L->length = 0;
4605  }
4606
4607  //if (TEST_OPT_PROT) { PrintS("N"); mflush(); }
4608  //L->Normalize(); // HANNES: should have a test
4609  kTest_L(L);
4610  return L->GetLmCurrRing();
4611}
4612
4613#ifdef HAVE_RINGS
4614poly redtailBba_Z (LObject* L, int pos, kStrategy strat )
4615// normalize=FALSE, withT=FALSE, coeff=Z
4616{
4617  strat->redTailChange=FALSE;
4618  if (strat->noTailReduction) return L->GetLmCurrRing();
4619  poly h, p;
4620  p = h = L->GetLmTailRing();
4621  if ((h==NULL) || (pNext(h)==NULL))
4622    return L->GetLmCurrRing();
4623
4624  TObject* With;
4625  // placeholder in case strat->tl < 0
4626  TObject  With_s(strat->tailRing);
4627
4628  LObject Ln(pNext(h), strat->tailRing);
4629  Ln.pLength = L->GetpLength() - 1;
4630
4631  pNext(h) = NULL;
4632  if (L->p != NULL) pNext(L->p) = NULL;
4633  L->pLength = 1;
4634
4635  Ln.PrepareRed(strat->use_buckets);
4636
4637  int cnt=REDTAIL_CANONICALIZE;
4638  while(!Ln.IsNull())
4639  {
4640    loop
4641    {
4642      Ln.SetShortExpVector();
4643      With = kFindDivisibleByInS(strat, pos, &Ln, &With_s);
4644      if (With == NULL) break;
4645      cnt--;
4646      if (cnt==0)
4647      {
4648        cnt=REDTAIL_CANONICALIZE;
4649        poly tmp=Ln.CanonicalizeP();
4650      }
4651      // we are in Z, do not Ccall pNorm
4652      strat->redTailChange=TRUE;
4653      // test divisibility of coefs:
4654      poly p_Ln=Ln.GetLmCurrRing();
4655      poly p_With=With->GetLmCurrRing();
4656      number z=nIntMod(pGetCoeff(p_Ln),pGetCoeff(p_With));
4657      if (!nIsZero(z))
4658      {
4659        // subtract z*Ln, add z.Ln to L
4660        poly m=pHead(p_Ln);
4661        pSetCoeff(m,z);
4662        poly mm=pHead(m);
4663        pNext(h) = m;
4664        pIter(h);
4665        L->pLength++;
4666        mm=pNeg(mm);
4667        if (Ln.bucket!=NULL)
4668        {
4669          int dummy=1;
4670          kBucket_Add_q(Ln.bucket,mm,&dummy);
4671        }
4672        else
4673        {
4674          if (Ln.p!=NULL) Ln.p=pAdd(Ln.p,mm);
4675          else if (Ln.t_p!=NULL)  Ln.t_p=p_Add_q(Ln.t_p,mm,strat->tailRing);
4676        }
4677      }
4678      else
4679        nDelete(&z);
4680
4681      if (ksReducePolyTail(L, With, &Ln))
4682      {
4683        // reducing the tail would violate the exp bound
4684        //  set a flag and hope for a retry (in bba)
4685        strat->completeReduce_retry=TRUE;
4686        if ((Ln.p != NULL) && (Ln.t_p != NULL)) Ln.p=NULL;
4687        do
4688        {
4689          pNext(h) = Ln.LmExtractAndIter();
4690          pIter(h);
4691          L->pLength++;
4692        } while (!Ln.IsNull());
4693        goto all_done;
4694      }
4695      if (Ln.IsNull()) goto all_done;
4696      With_s.Init(currRing);
4697    }
4698    pNext(h) = Ln.LmExtractAndIter();
4699    pIter(h);
4700    pNormalize(h);
4701    L->pLength++;
4702  }
4703
4704  all_done:
4705  Ln.Delete();
4706  if (L->p != NULL) pNext(L->p) = pNext(p);
4707
4708  if (strat->redTailChange)
4709  {
4710    L->last = NULL;
4711    L->length = 0;
4712  }
4713
4714  //if (TEST_OPT_PROT) { PrintS("N"); mflush(); }
4715  //L->Normalize(); // HANNES: should have a test
4716  kTest_L(L);
4717  return L->GetLmCurrRing();
4718}
4719#endif
4720
4721/*2
4722*checks the change degree and write progress report
4723*/
4724void message (int i,int* reduc,int* olddeg,kStrategy strat, int red_result)
4725{
4726  if (i != *olddeg)
4727  {
4728    Print("%d",i);
4729    *olddeg = i;
4730  }
4731  if (TEST_OPT_OLDSTD)
4732  {
4733    if (strat->Ll != *reduc)
4734    {
4735      if (strat->Ll != *reduc-1)
4736        Print("(%d)",strat->Ll+1);
4737      else
4738        PrintS("-");
4739      *reduc = strat->Ll;
4740    }
4741    else
4742      PrintS(".");
4743    mflush();
4744  }
4745  else
4746  {
4747    if (red_result == 0)
4748      PrintS("-");
4749    else if (red_result < 0)
4750      PrintS(".");
4751    if ((red_result > 0) || ((strat->Ll % 100)==99))
4752    {
4753      if (strat->Ll != *reduc && strat->Ll > 0)
4754      {
4755        Print("(%d)",strat->Ll+1);
4756        *reduc = strat->Ll;
4757      }
4758    }
4759  }
4760}
4761
4762/*2
4763*statistics
4764*/
4765void messageStat (int srmax,int lrmax,int hilbcount,kStrategy strat)
4766{
4767  //PrintS("\nUsage/Allocation of temporary storage:\n");
4768  //Print("%d/%d polynomials in standard base\n",srmax,IDELEMS(Shdl));
4769  //Print("%d/%d polynomials in set L (for lazy alg.)",lrmax+1,strat->Lmax);
4770  Print("product criterion:%d chain criterion:%d\n",strat->cp,strat->c3);
4771  if (hilbcount!=0) Print("hilbert series criterion:%d\n",hilbcount);
4772  /* in usual case strat->cv is 0, it gets changed only in shift routines */
4773  if (strat->cv!=0) Print("shift V criterion:%d\n",strat->cv);
4774  /*mflush();*/
4775}
4776
4777#ifdef KDEBUG
4778/*2
4779*debugging output: all internal sets, if changed
4780*for testing purpuse only/has to be changed for later use
4781*/
4782void messageSets (kStrategy strat)
4783{
4784  int i;
4785  if (strat->news)
4786  {
4787    PrintS("set S");
4788    for (i=0; i<=strat->sl; i++)
4789    {
4790      Print("\n  %d:",i);
4791      p_wrp(strat->S[i], currRing, strat->tailRing);
4792    }
4793    strat->news = FALSE;
4794  }
4795  if (strat->newt)
4796  {
4797    PrintS("\nset T");
4798    for (i=0; i<=strat->tl; i++)
4799    {
4800      Print("\n  %d:",i);
4801      strat->T[i].wrp();
4802      Print(" o:%ld e:%d l:%d",
4803        strat->T[i].pFDeg(),strat->T[i].ecart,strat->T[i].length);
4804    }
4805    strat->newt = FALSE;
4806  }
4807  PrintS("\nset L");
4808  for (i=strat->Ll; i>=0; i--)
4809  {
4810    Print("\n%d:",i);
4811    p_wrp(strat->L[i].p1, currRing, strat->tailRing);
4812    PrintS("  ");
4813    p_wrp(strat->L[i].p2, currRing, strat->tailRing);
4814    PrintS(" lcm: ");p_wrp(strat->L[i].lcm, currRing);
4815    PrintS("\n  p : ");
4816    strat->L[i].wrp();
4817    Print("  o:%ld e:%d l:%d",
4818          strat->L[i].pFDeg(),strat->L[i].ecart,strat->L[i].length);
4819  }
4820  PrintLn();
4821}
4822
4823#endif
4824
4825
4826/*2
4827*construct the set s from F
4828*/
4829void initS (ideal F, ideal Q, kStrategy strat)
4830{
4831  int   i,pos;
4832
4833  if (Q!=NULL) i=((IDELEMS(F)+IDELEMS(Q)+(setmaxTinc-1))/setmaxTinc)*setmaxTinc;
4834  else         i=((IDELEMS(F)+(setmaxTinc-1))/setmaxTinc)*setmaxTinc;
4835  strat->ecartS=initec(i);
4836  strat->sevS=initsevS(i);
4837  strat->S_2_R=initS_2_R(i);
4838  strat->fromQ=NULL;
4839  strat->Shdl=idInit(i,F->rank);
4840  strat->S=strat->Shdl->m;
4841  /*- put polys into S -*/
4842  if (Q!=NULL)
4843  {
4844    strat->fromQ=initec(i);
4845    memset(strat->fromQ,0,i*sizeof(int));
4846    for (i=0; i<IDELEMS(Q); i++)
4847    {
4848      if (Q->m[i]!=NULL)
4849      {
4850        LObject h;
4851        h.p = pCopy(Q->m[i]);
4852        if (TEST_OPT_INTSTRATEGY)
4853        {
4854          //pContent(h.p);
4855          h.pCleardenom(); // also does a pContent
4856        }
4857        else
4858        {
4859          h.pNorm();
4860        }
4861        if (pOrdSgn==-1)
4862        {
4863          deleteHC(&h, strat);
4864        }
4865        if (h.p!=NULL)
4866        {
4867          strat->initEcart(&h);
4868          if (strat->sl==-1)
4869            pos =0;
4870          else
4871          {
4872            pos = posInS(strat,strat->sl,h.p,h.ecart);
4873          }
4874          h.sev = pGetShortExpVector(h.p);
4875          strat->enterS(h,pos,strat,-1);
4876          strat->fromQ[pos]=1;
4877        }
4878      }
4879    }
4880  }
4881  for (i=0; i<IDELEMS(F); i++)
4882  {
4883    if (F->m[i]!=NULL)
4884    {
4885      LObject h;
4886      h.p = pCopy(F->m[i]);
4887      if (pOrdSgn==-1)
4888      {
4889        cancelunit(&h);  /*- tries to cancel a unit -*/
4890        deleteHC(&h, strat);
4891      }
4892      if (h.p!=NULL)
4893      // do not rely on the input being a SB!
4894      {
4895        if (TEST_OPT_INTSTRATEGY)
4896        {
4897          //pContent(h.p);
4898          h.pCleardenom(); // also does a pContent
4899        }
4900        else
4901        {
4902          h.pNorm();
4903        }
4904        strat->initEcart(&h);
4905        if (strat->sl==-1)
4906          pos =0;
4907        else
4908          pos = posInS(strat,strat->sl,h.p,h.ecart);
4909        h.sev = pGetShortExpVector(h.p);
4910        strat->enterS(h,pos,strat,-1);
4911      }
4912    }
4913  }
4914  /*- test, if a unit is in F -*/
4915  if ((strat->sl>=0)
4916#ifdef HAVE_RINGS
4917       && nIsUnit(pGetCoeff(strat->S[0]))
4918#endif
4919       && pIsConstant(strat->S[0]))
4920  {
4921    while (strat->sl>0) deleteInS(strat->sl,strat);
4922  }
4923}
4924
4925void initSL (ideal F, ideal Q,kStrategy strat)
4926{
4927  int   i,pos;
4928
4929  if (Q!=NULL) i=((IDELEMS(Q)+(setmaxTinc-1))/setmaxTinc)*setmaxTinc;
4930  else i=setmaxT;
4931  strat->ecartS=initec(i);
4932  strat->sevS=initsevS(i);
4933  strat->S_2_R=initS_2_R(i);
4934  strat->fromQ=NULL;
4935  strat->Shdl=idInit(i,F->rank);
4936  strat->S=strat->Shdl->m;
4937  /*- put polys into S -*/
4938  if (Q!=NULL)
4939  {
4940    strat->fromQ=initec(i);
4941    memset(strat->fromQ,0,i*sizeof(int));
4942    for (i=0; i<IDELEMS(Q); i++)
4943    {
4944      if (Q->m[i]!=NULL)
4945      {
4946        LObject h;
4947        h.p = pCopy(Q->m[i]);
4948        if (pOrdSgn==-1)
4949        {
4950          deleteHC(&h,strat);
4951        }
4952        if (TEST_OPT_INTSTRATEGY)
4953        {
4954          //pContent(h.p);
4955          h.pCleardenom(); // also does a pContent
4956        }
4957        else
4958        {
4959          h.pNorm();
4960        }
4961        if (h.p!=NULL)
4962        {
4963          strat->initEcart(&h);
4964          if (strat->sl==-1)
4965            pos =0;
4966          else
4967          {
4968            pos = posInS(strat,strat->sl,h.p,h.ecart);
4969          }
4970          h.sev = pGetShortExpVector(h.p);
4971          strat->enterS(h,pos,strat,-1);
4972          strat->fromQ[pos]=1;
4973        }
4974      }
4975    }
4976  }
4977  for (i=0; i<IDELEMS(F); i++)
4978  {
4979    if (F->m[i]!=NULL)
4980    {
4981      LObject h;
4982      h.p = pCopy(F->m[i]);
4983      if (h.p!=NULL)
4984      {
4985        if (pOrdSgn==-1)
4986        {
4987          cancelunit(&h);  /*- tries to cancel a unit -*/
4988          deleteHC(&h, strat);
4989        }
4990        if (h.p!=NULL)
4991        {
4992          if (TEST_OPT_INTSTRATEGY)
4993          {
4994            //pContent(h.p);
4995            h.pCleardenom(); // also does a pContent
4996          }
4997          else
4998          {
4999            h.pNorm();
5000          }
5001          strat->initEcart(&h);
5002          if (strat->Ll==-1)
5003            pos =0;
5004          else
5005            pos = strat->posInL(strat->L,strat->Ll,&h,strat);
5006          h.sev = pGetShortExpVector(h.p);
5007          enterL(&strat->L,&strat->Ll,&strat->Lmax,h,pos);
5008        }
5009      }
5010    }
5011  }
5012  /*- test, if a unit is in F -*/
5013
5014  if ((strat->Ll>=0)
5015#ifdef HAVE_RINGS
5016       && nIsUnit(pGetCoeff(strat->L[strat->Ll].p))
5017#endif
5018       && pIsConstant(strat->L[strat->Ll].p))
5019  {
5020    while (strat->Ll>0) deleteInL(strat->L,&strat->Ll,strat->Ll-1,strat);
5021  }
5022}
5023
5024
5025/*2
5026*construct the set s from F and {P}
5027*/
5028void initSSpecial (ideal F, ideal Q, ideal P,kStrategy strat)
5029{
5030  int   i,pos;
5031
5032  if (Q!=NULL) i=((IDELEMS(Q)+(setmaxTinc-1))/setmaxTinc)*setmaxTinc;
5033  else i=setmaxT;
5034  i=((i+IDELEMS(F)+IDELEMS(P)+15)/16)*16;
5035  strat->ecartS=initec(i);
5036  strat->sevS=initsevS(i);
5037  strat->S_2_R=initS_2_R(i);
5038  strat->fromQ=NULL;
5039  strat->Shdl=idInit(i,F->rank);
5040  strat->S=strat->Shdl->m;
5041
5042  /*- put polys into S -*/
5043  if (Q!=NULL)
5044  {
5045    strat->fromQ=initec(i);
5046    memset(strat->fromQ,0,i*sizeof(int));
5047    for (i=0; i<IDELEMS(Q); i++)
5048    {
5049      if (Q->m[i]!=NULL)
5050      {
5051        LObject h;
5052        h.p = pCopy(Q->m[i]);
5053        //if (TEST_OPT_INTSTRATEGY)
5054        //{
5055        //  //pContent(h.p);
5056        //  h.pCleardenom(); // also does a pContent
5057        //}
5058        //else
5059        //{
5060        //  h.pNorm();
5061        //}
5062        if (pOrdSgn==-1)
5063        {
5064          deleteHC(&h,strat);
5065        }
5066        if (h.p!=NULL)
5067        {
5068          strat->initEcart(&h);
5069          if (strat->sl==-1)
5070            pos =0;
5071          else
5072          {
5073            pos = posInS(strat,strat->sl,h.p,h.ecart);
5074          }
5075          h.sev = pGetShortExpVector(h.p);
5076          strat->enterS(h,pos,strat, strat->tl+1);
5077          enterT(h, strat);
5078          strat->fromQ[pos]=1;
5079        }
5080      }
5081    }
5082  }
5083  /*- put polys into S -*/
5084  for (i=0; i<IDELEMS(F); i++)
5085  {
5086    if (F->m[i]!=NULL)
5087    {
5088      LObject h;
5089      h.p = pCopy(F->m[i]);
5090      if (pOrdSgn==-1)
5091      {
5092        deleteHC(&h,strat);
5093      }
5094      else
5095      {
5096        h.p=redtailBba(h.p,strat->sl,strat);
5097      }
5098      if (h.p!=NULL)
5099      {
5100        strat->initEcart(&h);
5101        if (strat->sl==-1)
5102          pos =0;
5103        else
5104          pos = posInS(strat,strat->sl,h.p,h.ecart);
5105        h.sev = pGetShortExpVector(h.p);
5106        strat->enterS(h,pos,strat, strat->tl+1);
5107        enterT(h,strat);
5108      }
5109    }
5110  }
5111  for (i=0; i<IDELEMS(P); i++)
5112  {
5113    if (P->m[i]!=NULL)
5114    {
5115      LObject h;
5116      h.p=pCopy(P->m[i]);
5117      if (TEST_OPT_INTSTRATEGY)
5118      {
5119        h.pCleardenom();
5120      }
5121      else
5122      {
5123        h.pNorm();
5124      }
5125      if(strat->sl>=0)
5126      {
5127        if (pOrdSgn==1)
5128        {
5129          h.p=redBba(h.p,strat->sl,strat);
5130          if (h.p!=NULL)
5131          {
5132            h.p=redtailBba(h.p,strat->sl,strat);
5133          }
5134        }
5135        else
5136        {
5137          h.p=redMora(h.p,strat->sl,strat);
5138        }
5139        if(h.p!=NULL)
5140        {
5141          strat->initEcart(&h);
5142          if (TEST_OPT_INTSTRATEGY)
5143          {
5144            h.pCleardenom();
5145          }
5146          else
5147          {
5148            h.is_normalized = 0;
5149            h.pNorm();
5150          }
5151          h.sev = pGetShortExpVector(h.p);
5152          h.SetpFDeg();
5153          pos = posInS(strat,strat->sl,h.p,h.ecart);
5154          enterpairsSpecial(h.p,strat->sl,h.ecart,pos,strat,strat->tl+1);
5155          strat->enterS(h,pos,strat, strat->tl+1);
5156          enterT(h,strat);
5157        }
5158      }
5159      else
5160      {
5161        h.sev = pGetShortExpVector(h.p);
5162        strat->initEcart(&h);
5163        strat->enterS(h,0,strat, strat->tl+1);
5164        enterT(h,strat);
5165      }
5166    }
5167  }
5168}
5169/*2
5170* reduces h using the set S
5171* procedure used in cancelunit1
5172*/
5173static poly redBba1 (poly h,int maxIndex,kStrategy strat)
5174{
5175  int j = 0;
5176  unsigned long not_sev = ~ pGetShortExpVector(h);
5177
5178  while (j <= maxIndex)
5179  {
5180    if (pLmShortDivisibleBy(strat->S[j],strat->sevS[j],h, not_sev))
5181       return ksOldSpolyRedNew(strat->S[j],h,strat->kNoetherTail());
5182    else j++;
5183  }
5184  return h;
5185}
5186
5187/*2
5188*tests if p.p=monomial*unit and cancels the unit
5189*/
5190void cancelunit1 (LObject* p,int *suc, int index,kStrategy strat )
5191{
5192  int k;
5193  poly r,h,h1,q;
5194
5195  if (!pIsVector((*p).p) && ((*p).ecart != 0))
5196  {
5197#ifdef HAVE_RINGS_LOC
5198    // Leading coef have to be a unit
5199    if ( !(nIsUnit(p_GetCoeff((*p).p, r))) ) return;
5200#endif
5201    k = 0;
5202    h1 = r = pCopy((*p).p);
5203    h =pNext(r);
5204    loop
5205    {
5206      if (h==NULL)
5207      {
5208        pDelete(&r);
5209        pDelete(&(pNext((*p).p)));
5210        (*p).ecart = 0;
5211        (*p).length = 1;
5212#ifdef HAVE_RINGS_LOC
5213        (*p).pLength = 1;  // Why wasn't this set already?
5214#endif
5215        (*suc)=0;
5216        return;
5217      }
5218      if (!pDivisibleBy(r,h))
5219      {
5220        q=redBba1(h,index ,strat);
5221        if (q != h)
5222        {
5223          k++;
5224          pDelete(&h);
5225          pNext(h1) = h = q;
5226        }
5227        else
5228        {
5229          pDelete(&r);
5230          return;
5231        }
5232      }
5233      else
5234      {
5235        h1 = h;
5236        pIter(h);
5237      }
5238      if (k > 10)
5239      {
5240        pDelete(&r);
5241        return;
5242      }
5243    }
5244  }
5245}
5246
5247#if 0
5248/*2
5249* reduces h using the elements from Q in the set S
5250* procedure used in updateS
5251* must not be used for elements of Q or elements of an ideal !
5252*/
5253static poly redQ (poly h, int j, kStrategy strat)
5254{
5255  int start;
5256  unsigned long not_sev = ~ pGetShortExpVector(h);
5257  while ((j <= strat->sl) && (pGetComp(strat->S[j])!=0)) j++;
5258  start=j;
5259  while (j<=strat->sl)
5260  {
5261    if (pLmShortDivisibleBy(strat->S[j],strat->sevS[j], h, not_sev))
5262    {
5263      h = ksOldSpolyRed(strat->S[j],h,strat->kNoetherTail());
5264      if (h==NULL) return NULL;
5265      j = start;
5266      not_sev = ~ pGetShortExpVector(h);
5267    }
5268    else j++;
5269  }
5270  return h;
5271}
5272#endif
5273
5274/*2
5275* reduces h using the set S
5276* procedure used in updateS
5277*/
5278static poly redBba (poly h,int maxIndex,kStrategy strat)
5279{
5280  int j = 0;
5281  unsigned long not_sev = ~ pGetShortExpVector(h);
5282
5283  while (j <= maxIndex)
5284  {
5285    if (pLmShortDivisibleBy(strat->S[j],strat->sevS[j], h, not_sev))
5286    {
5287      h = ksOldSpolyRed(strat->S[j],h,strat->kNoetherTail());
5288      if (h==NULL) return NULL;
5289      j = 0;
5290      not_sev = ~ pGetShortExpVector(h);    }
5291    else j++;
5292  }
5293  return h;
5294}
5295
5296/*2
5297* reduces h using the set S
5298*e is the ecart of h
5299*procedure used in updateS
5300*/
5301static poly redMora (poly h,int maxIndex,kStrategy strat)
5302{
5303  int  j=0;
5304  int  e,l;
5305  unsigned long not_sev = ~ pGetShortExpVector(h);
5306
5307  if (maxIndex >= 0)
5308  {
5309    e = pLDeg(h,&l,currRing)-pFDeg(h,currRing);
5310    do
5311    {
5312      if (pLmShortDivisibleBy(strat->S[j],strat->sevS[j], h, not_sev)
5313      && ((e >= strat->ecartS[j]) || strat->kHEdgeFound))
5314      {
5315#ifdef KDEBUG
5316        if (TEST_OPT_DEBUG)
5317          {PrintS("reduce ");wrp(h);Print(" with S[%d] (",j);wrp(strat->S[j]);}
5318#endif
5319        h = ksOldSpolyRed(strat->S[j],h,strat->kNoetherTail());
5320#ifdef KDEBUG
5321        if(TEST_OPT_DEBUG)
5322          {PrintS(")\nto "); wrp(h); PrintLn();}
5323#endif
5324        // pDelete(&h);
5325        if (h == NULL) return NULL;
5326        e = pLDeg(h,&l,currRing)-pFDeg(h,currRing);
5327        j = 0;
5328        not_sev = ~ pGetShortExpVector(h);
5329      }
5330      else j++;
5331    }
5332    while (j <= maxIndex);
5333  }
5334  return h;
5335}
5336
5337/*2
5338*updates S:
5339*the result is a set of polynomials which are in
5340*normalform with respect to S
5341*/
5342void updateS(BOOLEAN toT,kStrategy strat)
5343{
5344  LObject h;
5345  int i, suc=0;
5346  poly redSi=NULL;
5347  BOOLEAN change,any_change;
5348//  Print("nach initS: updateS start mit sl=%d\n",(strat->sl));
5349//  for (i=0; i<=(strat->sl); i++)
5350//  {
5351//    Print("s%d:",i);
5352//    if (strat->fromQ!=NULL) Print("(Q:%d) ",strat->fromQ[i]);
5353//    pWrite(strat->S[i]);
5354//  }
5355//  Print("pOrdSgn=%d\n", pOrdSgn);
5356  any_change=FALSE;
5357  if (pOrdSgn==1)
5358  {
5359    while (suc != -1)
5360    {
5361      i=suc+1;
5362      while (i<=strat->sl)
5363      {
5364        change=FALSE;
5365        if (((strat->fromQ==NULL) || (strat->fromQ[i]==0)) && (i>0))
5366        {
5367          redSi = pHead(strat->S[i]);
5368          strat->S[i] = redBba(strat->S[i],i-1,strat);
5369          //if ((strat->ak!=0)&&(strat->S[i]!=NULL))
5370          //  strat->S[i]=redQ(strat->S[i],i+1,strat); /*reduce S[i] mod Q*/
5371          if (pCmp(redSi,strat->S[i])!=0)
5372          {
5373            change=TRUE;
5374            any_change=TRUE;
5375            #ifdef KDEBUG
5376            if (TEST_OPT_DEBUG)
5377            {
5378              PrintS("reduce:");
5379              wrp(redSi);PrintS(" to ");p_wrp(strat->S[i], currRing, strat->tailRing);PrintLn();
5380            }
5381            #endif
5382            if (TEST_OPT_PROT)
5383            {
5384              if (strat->S[i]==NULL)
5385                PrintS("V");
5386              else
5387                PrintS("v");
5388              mflush();
5389            }
5390          }
5391          pLmDelete(&redSi);
5392          if (strat->S[i]==NULL)
5393          {
5394            deleteInS(i,strat);
5395            i--;
5396          }
5397          else if (change)
5398          {
5399            if (TEST_OPT_INTSTRATEGY)
5400            {
5401              //pContent(strat->S[i]);
5402              strat->S[i]=p_Cleardenom(strat->S[i], currRing);// also does a pContent
5403            }
5404            else
5405            {
5406              pNorm(strat->S[i]);
5407            }
5408            strat->sevS[i] = pGetShortExpVector(strat->S[i]);
5409          }
5410        }
5411        i++;
5412      }
5413      if (any_change) reorderS(&suc,strat);
5414      else break;
5415    }
5416    if (toT)
5417    {
5418      for (i=0; i<=strat->sl; i++)
5419      {
5420        if ((strat->fromQ==NULL) || (strat->fromQ[i]==0))
5421        {
5422          h.p = redtailBba(strat->S[i],i-1,strat);
5423          if (TEST_OPT_INTSTRATEGY)
5424          {
5425            h.pCleardenom();// also does a pContent
5426          }
5427        }
5428        else
5429        {
5430          h.p = strat->S[i];
5431        }
5432        strat->initEcart(&h);
5433        if (strat->honey)
5434        {
5435          strat->ecartS[i] = h.ecart;
5436        }
5437        if (strat->sevS[i] == 0) {strat->sevS[i] = pGetShortExpVector(h.p);}
5438        else assume(strat->sevS[i] == pGetShortExpVector(h.p));
5439        h.sev = strat->sevS[i];
5440        /*puts the elements of S also to T*/
5441        strat->initEcart(&h);
5442        enterT(h,strat);
5443        strat->S_2_R[i] = strat->tl;
5444      }
5445    }
5446  }
5447  else
5448  {
5449    while (suc != -1)
5450    {
5451      i=suc;
5452      while (i<=strat->sl)
5453      {
5454        change=FALSE;
5455        if (((strat->fromQ==NULL) || (strat->fromQ[i]==0)) && (i>0))
5456        {
5457          redSi=pHead((strat->S)[i]);
5458          (strat->S)[i] = redMora((strat->S)[i],i-1,strat);
5459          if ((strat->S)[i]==NULL)
5460          {
5461            deleteInS(i,strat);
5462            i--;
5463          }
5464          else if (pCmp((strat->S)[i],redSi)!=0)
5465          {
5466            any_change=TRUE;
5467            h.p = strat->S[i];
5468            strat->initEcart(&h);
5469            strat->ecartS[i] = h.ecart;
5470            if (TEST_OPT_INTSTRATEGY)
5471            {
5472              strat->S[i]=p_Cleardenom(strat->S[i], currRing);// also does a pContent
5473            }
5474            else
5475            {
5476              pNorm(strat->S[i]); // == h.p
5477            }
5478            h.sev =  pGetShortExpVector(h.p);
5479            strat->sevS[i] = h.sev;
5480          }
5481          pLmDelete(&redSi);
5482          kTest(strat);
5483        }
5484        i++;
5485      }
5486#ifdef KDEBUG
5487      kTest(strat);
5488#endif
5489      if (any_change) reorderS(&suc,strat);
5490      else { suc=-1; break; }
5491      if (h.p!=NULL)
5492      {
5493        if (!strat->kHEdgeFound)
5494        {
5495          /*strat->kHEdgeFound =*/ HEckeTest(h.p,strat);
5496        }
5497        if (strat->kHEdgeFound)
5498          newHEdge(strat->S,strat);
5499      }
5500    }
5501    for (i=0; i<=strat->sl; i++)
5502    {
5503      if ((strat->fromQ==NULL) || (strat->fromQ[i]==0))
5504      {
5505        strat->S[i] = h.p = redtail(strat->S[i],strat->sl,strat);
5506        strat->initEcart(&h);
5507        strat->ecartS[i] = h.ecart;
5508        h.sev = pGetShortExpVector(h.p);
5509        strat->sevS[i] = h.sev;
5510      }
5511      else
5512      {
5513        h.p = strat->S[i];
5514        h.ecart=strat->ecartS[i];
5515        h.sev = strat->sevS[i];
5516        h.length = h.pLength = pLength(h.p);
5517      }
5518      if ((strat->fromQ==NULL) || (strat->fromQ[i]==0))
5519        cancelunit1(&h,&suc,strat->sl,strat);
5520      h.SetpFDeg();
5521      /*puts the elements of S also to T*/
5522      enterT(h,strat);
5523      strat->S_2_R[i] = strat->tl;
5524    }
5525    if (suc!= -1) updateS(toT,strat);
5526  }
5527#ifdef KDEBUG
5528  kTest(strat);
5529#endif
5530}
5531
5532
5533/*2
5534* -puts p to the standardbasis s at position at
5535* -saves the result in S
5536*/
5537void enterSBba (LObject p,int atS,kStrategy strat, int atR)
5538{
5539  int i;
5540  strat->news = TRUE;
5541  /*- puts p to the standardbasis s at position at -*/
5542  if (strat->sl == IDELEMS(strat->Shdl)-1)
5543  {
5544    strat->sevS = (unsigned long*) omRealloc0Size(strat->sevS,
5545                                    IDELEMS(strat->Shdl)*sizeof(unsigned long),
5546                                    (IDELEMS(strat->Shdl)+setmaxTinc)
5547                                                  *sizeof(unsigned long));
5548    strat->ecartS = (intset)omReallocSize(strat->ecartS,
5549                                          IDELEMS(strat->Shdl)*sizeof(int),
5550                                          (IDELEMS(strat->Shdl)+setmaxTinc)
5551                                                  *sizeof(int));
5552    strat->S_2_R = (int*) omRealloc0Size(strat->S_2_R,
5553                                         IDELEMS(strat->Shdl)*sizeof(int),
5554                                         (IDELEMS(strat->Shdl)+setmaxTinc)
5555                                                  *sizeof(int));
5556    if (strat->lenS!=NULL)
5557      strat->lenS=(int*)omRealloc0Size(strat->lenS,
5558                                       IDELEMS(strat->Shdl)*sizeof(int),
5559                                       (IDELEMS(strat->Shdl)+setmaxTinc)
5560                                                 *sizeof(int));
5561    if (strat->lenSw!=NULL)
5562      strat->lenSw=(wlen_type*)omRealloc0Size(strat->lenSw,
5563                                       IDELEMS(strat->Shdl)*sizeof(wlen_type),
5564                                       (IDELEMS(strat->Shdl)+setmaxTinc)
5565                                                 *sizeof(wlen_type));
5566    if (strat->fromQ!=NULL)
5567    {
5568      strat->fromQ = (intset)omReallocSize(strat->fromQ,
5569                                    IDELEMS(strat->Shdl)*sizeof(int),
5570                                    (IDELEMS(strat->Shdl)+setmaxTinc)*sizeof(int));
5571    }
5572    pEnlargeSet(&strat->S,IDELEMS(strat->Shdl),setmaxTinc);
5573    IDELEMS(strat->Shdl)+=setmaxTinc;
5574    strat->Shdl->m=strat->S;
5575  }
5576  if (atS <= strat->sl)
5577  {
5578#ifdef ENTER_USE_MEMMOVE
5579// #if 0
5580    memmove(&(strat->S[atS+1]), &(strat->S[atS]),
5581            (strat->sl - atS + 1)*sizeof(poly));
5582    memmove(&(strat->ecartS[atS+1]), &(strat->ecartS[atS]),
5583            (strat->sl - atS + 1)*sizeof(int));
5584    memmove(&(strat->sevS[atS+1]), &(strat->sevS[atS]),
5585            (strat->sl - atS + 1)*sizeof(unsigned long));
5586    memmove(&(strat->S_2_R[atS+1]), &(strat->S_2_R[atS]),
5587            (strat->sl - atS + 1)*sizeof(int));
5588    if (strat->lenS!=NULL)
5589    memmove(&(strat->lenS[atS+1]), &(strat->lenS[atS]),
5590            (strat->sl - atS + 1)*sizeof(int));
5591    if (strat->lenSw!=NULL)
5592    memmove(&(strat->lenSw[atS+1]), &(strat->lenSw[atS]),
5593            (strat->sl - atS + 1)*sizeof(wlen_type));
5594#else
5595    for (i=strat->sl+1; i>=atS+1; i--)
5596    {
5597      strat->S[i] = strat->S[i-1];
5598      strat->ecartS[i] = strat->ecartS[i-1];
5599      strat->sevS[i] = strat->sevS[i-1];
5600      strat->S_2_R[i] = strat->S_2_R[i-1];
5601    }
5602    if (strat->lenS!=NULL)
5603    for (i=strat->sl+1; i>=atS+1; i--)
5604      strat->lenS[i] = strat->lenS[i-1];
5605    if (strat->lenSw!=NULL)
5606    for (i=strat->sl+1; i>=atS+1; i--)
5607      strat->lenSw[i] = strat->lenSw[i-1];
5608#endif
5609  }
5610  if (strat->fromQ!=NULL)
5611  {
5612#ifdef ENTER_USE_MEMMOVE
5613    memmove(&(strat->fromQ[atS+1]), &(strat->fromQ[atS]),
5614                  (strat->sl - atS + 1)*sizeof(int));
5615#else
5616    for (i=strat->sl+1; i>=atS+1; i--)
5617    {
5618      strat->fromQ[i] = strat->fromQ[i-1];
5619    }
5620#endif
5621    strat->fromQ[atS]=0;
5622  }
5623
5624  /*- save result -*/
5625  strat->S[atS] = p.p;
5626  if (strat->honey) strat->ecartS[atS] = p.ecart;
5627  if (p.sev == 0)
5628    p.sev = pGetShortExpVector(p.p);
5629  else
5630    assume(p.sev == pGetShortExpVector(p.p));
5631  strat->sevS[atS] = p.sev;
5632  strat->ecartS[atS] = p.ecart;
5633  strat->S_2_R[atS] = atR;
5634  strat->sl++;
5635}
5636
5637/*2
5638* puts p to the set T at position atT
5639*/
5640void enterT(LObject p, kStrategy strat, int atT)
5641{
5642  int i;
5643
5644  pp_Test(p.p, currRing, p.tailRing);
5645  assume(strat->tailRing == p.tailRing);
5646  // redMoraNF complains about this -- but, we don't really
5647  // neeed this so far
5648  assume(p.pLength == 0 || pLength(p.p) == p.pLength || rIsSyzIndexRing(currRing)); // modulo syzring
5649  assume(p.FDeg == p.pFDeg());
5650  assume(!p.is_normalized || nIsOne(pGetCoeff(p.p)));
5651
5652#ifdef KDEBUG
5653  // do not put an LObject twice into T:
5654  for(i=strat->tl;i>=0;i--)
5655  {
5656    if (p.p==strat->T[i].p)
5657    {
5658      printf("already in T at pos %d of %d, atT=%d\n",i,strat->tl,atT);
5659      return;
5660    }
5661  }
5662#endif
5663  strat->newt = TRUE;
5664  if (atT < 0)
5665    atT = strat->posInT(strat->T, strat->tl, p);
5666  if (strat->tl == strat->tmax-1)
5667    enlargeT(strat->T,strat->R,strat->sevT,strat->tmax,setmaxTinc);
5668  if (atT <= strat->tl)
5669  {
5670#ifdef ENTER_USE_MEMMOVE
5671    memmove(&(strat->T[atT+1]), &(strat->T[atT]),
5672            (strat->tl-atT+1)*sizeof(TObject));
5673    memmove(&(strat->sevT[atT+1]), &(strat->sevT[atT]),
5674            (strat->tl-atT+1)*sizeof(unsigned long));
5675#endif
5676    for (i=strat->tl+1; i>=atT+1; i--)
5677    {
5678#ifndef ENTER_USE_MEMMOVE
5679      strat->T[i] = strat->T[i-1];
5680      strat->sevT[i] = strat->sevT[i-1];
5681#endif
5682      strat->R[strat->T[i].i_r] = &(strat->T[i]);
5683    }
5684  }
5685
5686  if (strat->tailBin != NULL && (pNext(p.p) != NULL))
5687  {
5688    pNext(p.p)=p_ShallowCopyDelete(pNext(p.p),
5689                                   (strat->tailRing != NULL ?
5690                                    strat->tailRing : currRing),
5691                                   strat->tailBin);
5692    if (p.t_p != NULL) pNext(p.t_p) = pNext(p.p);
5693  }
5694  strat->T[atT] = (TObject) p;
5695
5696  if (strat->tailRing != currRing && pNext(p.p) != NULL)
5697    strat->T[atT].max = p_GetMaxExpP(pNext(p.p), strat->tailRing);
5698  else
5699    strat->T[atT].max = NULL;
5700
5701  strat->tl++;
5702  strat->R[strat->tl] = &(strat->T[atT]);
5703  strat->T[atT].i_r = strat->tl;
5704  assume(p.sev == 0 || pGetShortExpVector(p.p) == p.sev);
5705  strat->sevT[atT] = (p.sev == 0 ? pGetShortExpVector(p.p) : p.sev);
5706  kTest_T(&(strat->T[atT]));
5707}
5708
5709void initHilbCrit(ideal F, ideal Q, intvec **hilb,kStrategy strat)
5710{
5711  if (strat->homog!=isHomog)
5712  {
5713    *hilb=NULL;
5714  }
5715}
5716
5717void initBuchMoraCrit(kStrategy strat)
5718{
5719  strat->enterOnePair=enterOnePairNormal;
5720  strat->chainCrit=chainCritNormal;
5721#ifdef HAVE_RINGS
5722  if (rField_is_Ring(currRing))
5723  {
5724    strat->enterOnePair=enterOnePairRing;
5725    strat->chainCrit=chainCritRing;
5726  }
5727#endif
5728#ifdef HAVE_RATGRING
5729  if (rIsRatGRing(currRing))
5730  {
5731     strat->chainCrit=chainCritPart;
5732     /* enterOnePairNormal get rational part in it */
5733  }
5734#endif
5735
5736  strat->sugarCrit =        TEST_OPT_SUGARCRIT;
5737  strat->Gebauer =          strat->homog || strat->sugarCrit;
5738  strat->honey =            !strat->homog || strat->sugarCrit || TEST_OPT_WEIGHTM;
5739  if (TEST_OPT_NOT_SUGAR) strat->honey = FALSE;
5740  strat->pairtest = NULL;
5741  /* alway use tailreduction, except:
5742  * - in local rings, - in lex order case, -in ring over extensions */
5743  strat->noTailReduction = !TEST_OPT_REDTAIL;
5744
5745#ifdef HAVE_PLURAL
5746  // and r is plural_ring
5747  //  hence this holds for r a rational_plural_ring
5748  if( rIsPluralRing(currRing) || (rIsSCA(currRing) && !strat->z2homog) )
5749  {    //or it has non-quasi-comm type... later
5750    strat->sugarCrit = FALSE;
5751    strat->Gebauer = FALSE;
5752    strat->honey = FALSE;
5753  }
5754#endif
5755
5756#ifdef HAVE_RINGS
5757  // Coefficient ring?
5758  if (rField_is_Ring(currRing))
5759  {
5760    strat->sugarCrit = FALSE;
5761    strat->Gebauer = FALSE ;
5762    strat->honey = FALSE;
5763  }
5764#endif
5765  #ifdef KDEBUG
5766  if (TEST_OPT_DEBUG)
5767  {
5768    if (strat->homog) PrintS("ideal/module is homogeneous\n");
5769    else              PrintS("ideal/module is not homogeneous\n");
5770  }
5771  #endif
5772}
5773
5774BOOLEAN kPosInLDependsOnLength(int (*pos_in_l)
5775                               (const LSet set, const int length,
5776                                LObject* L,const kStrategy strat))
5777{
5778  if (pos_in_l == posInL110 ||
5779      pos_in_l == posInL10)
5780    return TRUE;
5781
5782  return FALSE;
5783}
5784
5785void initBuchMoraPos (kStrategy strat)
5786{
5787  if (pOrdSgn==1)
5788  {
5789    if (strat->honey)
5790    {
5791      strat->posInL = posInL15;
5792      // ok -- here is the deal: from my experiments for Singular-2-0
5793      // I conclude that that posInT_EcartpLength is the best of
5794      // posInT15, posInT_EcartFDegpLength, posInT_FDegLength, posInT_pLength
5795      // see the table at the end of this file
5796      if (TEST_OPT_OLDSTD)
5797        strat->posInT = posInT15;
5798      else
5799        strat->posInT = posInT_EcartpLength;
5800    }
5801    else if (pLexOrder && !TEST_OPT_INTSTRATEGY)
5802    {
5803      strat->posInL = posInL11;
5804      strat->posInT = posInT11;
5805    }
5806    else if (TEST_OPT_INTSTRATEGY)
5807    {
5808      strat->posInL = posInL11;
5809      strat->posInT = posInT11;
5810    }
5811    else
5812    {
5813      strat->posInL = posInL0;
5814      strat->posInT = posInT0;
5815    }
5816    //if (strat->minim>0) strat->posInL =posInLSpecial;
5817    if (strat->homog)
5818    {
5819       strat->posInL = posInL110;
5820       strat->posInT = posInT110;
5821    }
5822  }
5823  else
5824  {
5825    if (strat->homog)
5826    {
5827      strat->posInL = posInL11;
5828      strat->posInT = posInT11;
5829    }
5830    else
5831    {
5832      if ((currRing->order[0]==ringorder_c)
5833      ||(currRing->order[0]==ringorder_C))
5834      {
5835        strat->posInL = posInL17_c;
5836        strat->posInT = posInT17_c;
5837      }
5838      else
5839      {
5840        strat->posInL = posInL17;
5841        strat->posInT = posInT17;
5842      }
5843    }
5844  }
5845  if (strat->minim>0) strat->posInL =posInLSpecial;
5846  // for further tests only
5847  if ((BTEST1(11)) || (BTEST1(12)))
5848    strat->posInL = posInL11;
5849  else if ((BTEST1(13)) || (BTEST1(14)))
5850    strat->posInL = posInL13;
5851  else if ((BTEST1(15)) || (BTEST1(16)))
5852    strat->posInL = posInL15;
5853  else if ((BTEST1(17)) || (BTEST1(18)))
5854    strat->posInL = posInL17;
5855  if (BTEST1(11))
5856    strat->posInT = posInT11;
5857  else if (BTEST1(13))
5858    strat->posInT = posInT13;
5859  else if (BTEST1(15))
5860    strat->posInT = posInT15;
5861  else if ((BTEST1(17)))
5862    strat->posInT = posInT17;
5863  else if ((BTEST1(19)))
5864    strat->posInT = posInT19;
5865  else if (BTEST1(12) || BTEST1(14) || BTEST1(16) || BTEST1(18))
5866    strat->posInT = posInT1;
5867#ifdef HAVE_RINGS
5868  if (rField_is_Ring(currRing))
5869  {
5870    strat->posInL = posInL11;
5871    strat->posInT = posInT11;
5872  }
5873#endif
5874  strat->posInLDependsOnLength = kPosInLDependsOnLength(strat->posInL);
5875}
5876
5877void initBuchMora (ideal F,ideal Q,kStrategy strat)
5878{
5879  strat->interpt = BTEST1(OPT_INTERRUPT);
5880  strat->kHEdge=NULL;
5881  if (pOrdSgn==1) strat->kHEdgeFound=FALSE;
5882  /*- creating temp data structures------------------- -*/
5883  strat->cp = 0;
5884  strat->c3 = 0;
5885  strat->tail = pInit();
5886  /*- set s -*/
5887  strat->sl = -1;
5888  /*- set L -*/
5889  strat->Lmax = ((IDELEMS(F)+setmaxLinc-1)/setmaxLinc)*setmaxLinc;
5890  strat->Ll = -1;
5891  strat->L = initL(((IDELEMS(F)+setmaxLinc-1)/setmaxLinc)*setmaxLinc);
5892  /*- set B -*/
5893  strat->Bmax = setmaxL;
5894  strat->Bl = -1;
5895  strat->B = initL();
5896  /*- set T -*/
5897  strat->tl = -1;
5898  strat->tmax = setmaxT;
5899  strat->T = initT();
5900  strat->R = initR();
5901  strat->sevT = initsevT();
5902  /*- init local data struct.---------------------------------------- -*/
5903  strat->P.ecart=0;
5904  strat->P.length=0;
5905  if (pOrdSgn==-1)
5906  {
5907    if (strat->kHEdge!=NULL) pSetComp(strat->kHEdge, strat->ak);
5908    if (strat->kNoether!=NULL) pSetComp(strat->kNoetherTail(), strat->ak);
5909  }
5910  if(TEST_OPT_SB_1)
5911  {
5912    int i;
5913    ideal P=idInit(IDELEMS(F)-strat->newIdeal,F->rank);
5914    for (i=strat->newIdeal;i<IDELEMS(F);i++)
5915    {
5916      P->m[i-strat->newIdeal] = F->m[i];
5917      F->m[i] = NULL;
5918    }
5919    initSSpecial(F,Q,P,strat);
5920    for (i=strat->newIdeal;i<IDELEMS(F);i++)
5921    {
5922      F->m[i] = P->m[i-strat->newIdeal];
5923      P->m[i-strat->newIdeal] = NULL;
5924    }
5925    idDelete(&P);
5926  }
5927  else
5928  {
5929    /*Shdl=*/initSL(F, Q,strat); /*sets also S, ecartS, fromQ */
5930    // /*Shdl=*/initS(F, Q,strat); /*sets also S, ecartS, fromQ */
5931  }
5932  strat->kIdeal = NULL;
5933  strat->fromT = FALSE;
5934  strat->noTailReduction = !TEST_OPT_REDTAIL;
5935  if (!TEST_OPT_SB_1)
5936  {
5937    updateS(TRUE,strat);
5938  }
5939  if (strat->fromQ!=NULL) omFreeSize(strat->fromQ,IDELEMS(strat->Shdl)*sizeof(int));
5940  strat->fromQ=NULL;
5941}
5942
5943void exitBuchMora (kStrategy strat)
5944{
5945  /*- release temp data -*/
5946  cleanT(strat);
5947  omFreeSize(strat->T,(strat->tmax)*sizeof(TObject));
5948  omFreeSize(strat->R,(strat->tmax)*sizeof(TObject*));
5949  omFreeSize(strat->sevT, (strat->tmax)*sizeof(unsigned long));
5950  omFreeSize(strat->ecartS,IDELEMS(strat->Shdl)*sizeof(int));
5951  omFreeSize(strat->sevS,IDELEMS(strat->Shdl)*sizeof(int));
5952  omFreeSize(strat->S_2_R,IDELEMS(strat->Shdl)*sizeof(int));
5953  /*- set L: should be empty -*/
5954  omFreeSize(strat->L,(strat->Lmax)*sizeof(LObject));
5955  /*- set B: should be empty -*/
5956  omFreeSize(strat->B,(strat->Bmax)*sizeof(LObject));
5957  pLmDelete(&strat->tail);
5958  strat->syzComp=0;
5959  if (strat->kIdeal!=NULL)
5960  {
5961    omFreeBin(strat->kIdeal, sleftv_bin);
5962    strat->kIdeal=NULL;
5963  }
5964}
5965
5966/*2
5967* in the case of a standardbase of a module over a qring:
5968* replace polynomials in i by ak vectors,
5969* (the polynomial * unit vectors gen(1)..gen(ak)
5970* in every case (also for ideals:)
5971* deletes divisible vectors/polynomials
5972*/
5973void updateResult(ideal r,ideal Q, kStrategy strat)
5974{
5975  int l;
5976  if (strat->ak>0)
5977  {
5978    for (l=IDELEMS(r)-1;l>=0;l--)
5979    {
5980      if ((r->m[l]!=NULL) && (pGetComp(r->m[l])==0))
5981      {
5982        pDelete(&r->m[l]); // and set it to NULL
5983      }
5984    }
5985    int q;
5986    poly p;
5987    for (l=IDELEMS(r)-1;l>=0;l--)
5988    {
5989      if ((r->m[l]!=NULL)
5990      //&& (strat->syzComp>0)
5991      //&& (pGetComp(r->m[l])<=strat->syzComp)
5992      )
5993      {
5994        for(q=IDELEMS(Q)-1; q>=0;q--)
5995        {
5996          if ((Q->m[q]!=NULL)
5997          &&(pLmDivisibleBy(Q->m[q],r->m[l])))
5998          {
5999            if (TEST_OPT_REDSB)
6000            {
6001              p=r->m[l];
6002              r->m[l]=kNF(Q,NULL,p);
6003              pDelete(&p);
6004            }
6005            else
6006            {
6007              pDelete(&r->m[l]); // and set it to NULL
6008            }
6009            break;
6010          }
6011        }
6012      }
6013    }
6014  }
6015  else
6016  {
6017    int q;
6018    poly p;
6019    BOOLEAN reduction_found=FALSE;
6020    for (l=IDELEMS(r)-1;l>=0;l--)
6021    {
6022      if (r->m[l]!=NULL)
6023      {
6024        for(q=IDELEMS(Q)-1; q>=0;q--)
6025        {
6026          if ((Q->m[q]!=NULL)
6027          &&(pLmEqual(r->m[l],Q->m[q])))
6028          {
6029            if (TEST_OPT_REDSB)
6030            {
6031              p=r->m[l];
6032              r->m[l]=kNF(Q,NULL,p);
6033              pDelete(&p);
6034              reduction_found=TRUE;
6035            }
6036            else
6037            {
6038              pDelete(&r->m[l]); // and set it to NULL
6039            }
6040            break;
6041          }
6042        }
6043      }
6044    }
6045    if (/*TEST_OPT_REDSB &&*/ reduction_found)
6046    {
6047      for (l=IDELEMS(r)-1;l>=0;l--)
6048      {
6049        if (r->m[l]!=NULL)
6050        {
6051          for(q=IDELEMS(r)-1;q>=0;q--)
6052          {
6053            if ((l!=q)
6054            && (r->m[q]!=NULL)
6055            &&(pLmDivisibleBy(r->m[l],r->m[q])))
6056            {
6057