source: git/kernel/tgb.cc @ 9cf06bf

spielwiese
Last change on this file since 9cf06bf was 9cf06bf, checked in by Michael Brickenstein <bricken@…>, 18 years ago
*bricken: more careful for int overflows in multipl. of length and coef size git-svn-id: file:///usr/local/Singular/svn/trunk@9187 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 78.7 KB
Line 
1//! \file tgb.cc
2//       multiple rings
3//       shorten_tails und dessen Aufrufe pruefen wlength!!!
4/****************************************
5*  Computer Algebra System SINGULAR     *
6****************************************/
7/* $Id: tgb.cc,v 1.95 2006-06-08 08:01:51 bricken Exp $ */
8/*
9* ABSTRACT: slimgb and F4 implementation
10*/
11//#include <vector>
12//using namespace std;
13#include "mod2.h"
14#include "tgb.h"
15#include "tgb_internal.h"
16#include "tgbgauss.h"
17#include "F4.h"
18#include "digitech.h"
19#include "gring.h"
20
21#include "longrat.h"
22#define SR_HDL(A) ((long)(A))
23static const int bundle_size=100;
24static const int delay_factor=3;
25int QlogSize(number n);
26#define ADD_LATER_SIZE 500
27#if 1
28static omBin lm_bin=NULL;
29
30static void simplify_poly(poly p, ring r) {
31     assume(r==currRing);
32     if (!rField_is_Zp(r))
33     {
34        pCleardenom(p);
35        pContent(p); //is a duplicate call, but belongs here
36     }
37     else
38       pNorm(p);
39}
40//static const BOOLEAN up_to_radical=TRUE;
41
42int slim_nsize(number n, ring r)
43{
44  if (rField_is_Zp(r))
45  {
46    return 1;
47  }
48  if (rField_is_Q(r))
49  {
50    return QlogSize(n);
51  }
52  else
53  {
54    return n_Size(n,r);
55  }
56}
57static BOOLEAN monomial_root(poly m, ring r){
58    BOOLEAN changed=FALSE;
59    int i;
60    for(i=1;i<=rVar(r);i++){
61        int e=p_GetExp(m,i,r);
62        if (e>1){
63            p_SetExp(m,i,1,r);
64            changed=TRUE;
65        }
66    }
67    if (changed) {
68        p_Setm(m,r);
69    }
70    return changed;
71}
72static BOOLEAN polynomial_root(poly h, ring r){
73  poly got=gcd_of_terms(h,r);
74  BOOLEAN changed=FALSE;
75  if((got!=NULL) &&(TEST_V_UPTORADICAL)) {
76    poly copy=p_Copy(got,r);
77    //p_wrp(got,c->r);
78    changed=monomial_root(got,r);
79    if (changed)
80    {
81         poly div_by=pDivide(copy, got);
82         poly iter=h;
83         while(iter){
84            pExpVectorSub(iter,div_by);
85            pIter(iter);
86         }
87         p_Delete(&div_by, r);
88         if (TEST_OPT_PROT)
89             PrintS("U");
90    }
91    p_Delete(&copy,r);
92  }
93  p_Delete(&got,r);
94  return changed;
95}
96static inline poly p_Init_Special(const ring r)
97{
98  return p_Init(r,lm_bin);
99}
100static inline poly pOne_Special(const ring r=currRing)
101{
102  poly rc = p_Init_Special(r);
103  pSetCoeff0(rc,r->cf->nInit(1));
104  return rc;
105}
106// zum Initialiseren: in t_rep_gb plazieren:
107
108#endif
109#define LEN_VAR3
110#define degbound(p) assume(pTotaldegree(p)<10)
111//#define inDebug(p) assume((debug_Ideal==NULL)||(kNF(debug_Ideal,NULL,p,0,0)==0))
112
113//die meisten Varianten stossen sich an coef_buckets
114
115
116
117#ifdef LEN_VAR1
118// erste Variante: Laenge: Anzahl der Monome
119inline int pSLength(poly p, int l) {
120  return l; }
121inline int kSBucketLength(kBucket* bucket, poly lm) {return bucket_guess(bucket);}
122#endif
123
124#ifdef LEN_VAR2
125// 2. Variante: Laenge: Platz fuer die Koeff.
126int pSLength(poly p,int l)
127{
128  int s=0;
129  while (p!=NULL) { s+=nSize(pGetCoeff(p));pIter(p); }
130  return s;
131}
132int kSBucketLength(kBucket* b, poly lm)
133{
134  int s=0;
135  int i;
136  for (i=MAX_BUCKET;i>=0;i--)
137  {
138    s+=pSLength(b->buckets[i],0);
139  }
140  return s;
141}
142#endif
143
144
145// 3.Variante: Laenge: Platz fuer Leitk * Monomanzahl
146
147
148
149
150int QlogSize(number n){
151
152    if (SR_HDL(n) & SR_INT){
153       long i=SR_TO_INT(n);
154       if (i==0) return 0;
155
156       unsigned long v;
157       v=(i>=0)?i:-i;
158       int r = 0;
159
160       while (v >>= 1)
161       {
162        r++;
163       }
164       return r+1;
165    }
166    //assume denominator is 0
167    return mpz_sizeinbase(&n->z,2);
168}
169
170
171#ifdef LEN_VAR3
172inline wlen_type pSLength(poly p,int l)
173{
174  int c;
175  number coef=pGetCoeff(p);
176  if (rField_is_Q(currRing)){
177    c=QlogSize(coef);
178  }
179  else
180    c=nSize(coef);
181  if (!(TEST_V_COEFSTRAT))
182      return (wlen_type)c*(wlen_type)l /*pLength(p)*/;
183  else {
184    wlen_type res=l;
185    res*=c;
186    res*=c;
187    return res;
188  }
189}
190//! TODO CoefBuckets berücksichtigen
191wlen_type kSBucketLength(kBucket* b, poly lm=NULL)
192{
193  int s=0;
194  int c;
195  number coef;
196  if(lm==NULL)
197    coef=pGetCoeff(kBucketGetLm(b));
198    //c=nSize(pGetCoeff(kBucketGetLm(b)));
199  else
200    coef=pGetCoeff(lm);
201    //c=nSize(pGetCoeff(lm));
202  if (rField_is_Q(currRing)){
203    c=QlogSize(coef);
204  }
205  else
206    c=nSize(coef);
207
208  int i;
209  for (i=b->buckets_used;i>=0;i--)
210  {
211    assume((b->buckets_length[i]==0)||(b->buckets[i]!=NULL));
212    s+=b->buckets_length[i] /*pLength(b->buckets[i])*/;
213  }
214  #ifdef HAVE_COEF_BUCKETS
215  assume(b->buckets[0]==kBucketGetLm(b));
216  if (b->coef[0]!=NULL){
217
218    if (rField_is_Q(currRing)){
219      int modifier=QlogSize(pGetCoeff(b->coef[0]));
220      c+=modifier;
221  }
222    else{
223      int modifier=nSize(pGetCoeff(b->coef[0]));
224      c*=modifier;}
225    }
226  #endif
227  if (!(TEST_V_COEFSTRAT)){
228  return s*c;
229  } else
230  {
231    wlen_type res=s;
232    res*=c;
233    res*=c;
234    return res;
235  }
236}
237#endif
238#ifdef LEN_VAR5
239inline wlen_type pSLength(poly p,int l)
240{
241  int c;
242  number coef=pGetCoeff(p);
243  if (rField_is_Q(currRing)){
244    c=QlogSize(coef);
245  }
246  else
247    c=nSize(coef);
248  wlen_type erg=l;
249  erg*=c;
250  erg*=c;
251  //PrintS("lenvar 5");
252  assume(erg>=0);
253  return erg; /*pLength(p)*/;
254}
255//! TODO CoefBuckets berücksichtigen
256wlen_type kSBucketLength(kBucket* b, poly lm=NULL)
257{
258  wlen_type s=0;
259  int c;
260  number coef;
261  if(lm==NULL)
262    coef=pGetCoeff(kBucketGetLm(b));
263    //c=nSize(pGetCoeff(kBucketGetLm(b)));
264  else
265    coef=pGetCoeff(lm);
266    //c=nSize(pGetCoeff(lm));
267  if (rField_is_Q(currRing)){
268    c=QlogSize(coef);
269  }
270  else
271    c=nSize(coef);
272
273  int i;
274  for (i=b->buckets_used;i>=0;i--)
275  {
276    assume((b->buckets_length[i]==0)||(b->buckets[i]!=NULL));
277    s+=b->buckets_length[i] /*pLength(b->buckets[i])*/;
278  }
279  #ifdef HAVE_COEF_BUCKETS
280  assume(b->buckets[0]==kBucketGetLm(b));
281  if (b->coef[0]!=NULL){
282
283    if (rField_is_Q(currRing)){
284      int modifier=QlogSize(pGetCoeff(b->coef[0]));
285      c+=modifier;
286  }
287    else{
288      int modifier=nSize(pGetCoeff(b->coef[0]));
289      c*=modifier;}
290    }
291  #endif
292  wlen_type erg=s;
293  erg*=c;
294  erg*=c;
295  return erg;
296}
297#endif
298
299#ifdef LEN_VAR4
300// 4.Variante: Laenge: Platz fuer Leitk * (1+Platz fuer andere Koeff.)
301int pSLength(poly p, int l)
302{
303  int s=1;
304  int c=nSize(pGetCoeff(p));
305  pIter(p);
306  while (p!=NULL) { s+=nSize(pGetCoeff(p));pIter(p); }
307  return s*c;
308}
309int kSBucketLength(kBucket* b)
310{
311  int s=1;
312  int c=nSize(pGetCoeff(kBucketGetLm(b)));
313  int i;
314  for (i=MAX_BUCKET;i>0;i--)
315  {
316    if(b->buckets[i]==NULL) continue;
317    s+=pSLength(b->buckets[i],0);
318  }
319  return s*c;
320}
321#endif
322static int do_pELength(poly p, slimgb_alg* c, int dlm=-1){
323  if(p==NULL) return 0;
324  int s=0;
325  poly pi=p;
326  if(dlm<0){
327    dlm=pTotaldegree(p,c->r);
328    s=1;
329    pi=p->next;
330  }
331
332  while(pi){
333    int d=pTotaldegree(pi,c->r);
334    if(d>dlm)
335      s+=1+d-dlm;
336    else
337      ++s;
338    pi=pi->next;
339  }
340  return s;
341}
342
343wlen_type pELength(poly p, ring r){
344  if(p==NULL) return 0;
345  int s=0;
346  poly pi=p;
347  int dlm;
348    dlm=pTotaldegree(p,r);
349    s=1;
350    pi=p->next;
351
352
353  while(pi){
354    int d=pTotaldegree(pi,r);
355    if(d>dlm)
356      s+=1+d-dlm;
357    else
358      ++s;
359    pi=pi->next;
360  }
361  return s;
362}
363
364int kEBucketLength(kBucket* b, poly lm,slimgb_alg* ca)
365{
366  int s=0;
367  if(lm==NULL){
368    lm=kBucketGetLm(b);
369  }
370  if(lm==NULL) return 0;
371  int d=pTotaldegree(lm,ca->r);
372  int i;
373  for (i=b->buckets_used;i>=0;i--)
374  {
375    if(b->buckets[i]==NULL) continue;
376    s+=do_pELength(b->buckets[i],ca,d);
377  }
378  return s;
379}
380
381static inline int pELength(poly p, slimgb_alg* c,int l){
382  if (p==NULL) return 0;
383  return do_pELength(p,c);
384}
385
386
387
388
389static inline wlen_type pQuality(poly p, slimgb_alg* c, int l=-1){
390
391  if(l<0)
392    l=pLength(p);
393  if(c->isDifficultField) {
394    if(c->eliminationProblem){
395      int cs;
396      number coef=pGetCoeff(p);
397      if (rField_is_Q(currRing)){
398         cs=QlogSize(coef);
399      }
400      else
401     cs=nSize(coef);
402     wlen_type erg=cs;
403     //erg*=cs;//for quadratic
404     erg*=pELength(p,c,l);
405    //FIXME: not quadratic coeff size
406      //return cs*pELength(p,c,l);
407      return erg;
408    }
409    //Print("I am here");
410    wlen_type r=pSLength(p,l);
411    assume(r>=0);
412    return r;
413  }
414  if(c->eliminationProblem) return pELength(p,c,l);
415  return l;
416}
417
418static inline int pTotaldegree_full(poly p){
419  int r=0;
420  while(p){
421    int d=pTotaldegree(p);
422    r=si_max(r,d);
423    pIter(p);
424  }
425  return r;
426}
427
428wlen_type red_object::guess_quality(slimgb_alg* c){
429    //works at the moment only for lenvar 1, because in different
430    //case, you have to look on coefs
431    wlen_type s=0;
432    if (c->isDifficultField){
433      //s=kSBucketLength(bucket,this->p);
434      if(c->eliminationProblem){
435    int cs;
436    number coef;
437
438    coef=pGetCoeff(kBucketGetLm(bucket));
439    //c=nSize(pGetCoeff(kBucketGetLm(b)));
440
441    //c=nSize(pGetCoeff(lm));
442    if (rField_is_Q(currRing)){
443      cs=QlogSize(coef);
444    }
445    else
446      cs=nSize(coef);
447    #ifdef HAVE_COEF_BUCKETS
448    if (bucket->coef[0]!=NULL){
449      if (rField_is_Q(currRing)){
450        int modifier=QlogSize(pGetCoeff(bucket->coef[0]));
451        cs+=modifier;
452      }
453      else{
454        int modifier=nSize(pGetCoeff(bucket->coef[0]));
455        cs*=modifier;}
456    }
457    #endif
458    //FIXME:not quadratic
459    wlen_type erg=kEBucketLength(this->bucket,this->p,c);
460    //erg*=cs;//for quadratic
461    erg*=cs;
462    //return cs*kEBucketLength(this->bucket,this->p,c);
463    return erg;
464      }
465      s=kSBucketLength(bucket,NULL);
466    }
467    else
468    {
469      if(c->eliminationProblem)
470  //if (false)
471  s=kEBucketLength(this->bucket,this->p,c);
472      else s=bucket_guess(bucket);
473    }
474
475    return s;
476}
477
478
479
480static void finalize_reduction_step(reduction_step* r){
481  delete r;
482}
483static int LObject_better_gen(const void* ap, const void* bp)
484{
485  LObject* a=*(LObject**)ap;
486  LObject* b=*(LObject**)bp;
487  return(pLmCmp(a->p,b->p));
488}
489static int red_object_better_gen(const void* ap, const void* bp)
490{
491
492
493  return(pLmCmp(((red_object*) ap)->p,((red_object*) bp)->p));
494}
495
496
497static int pLmCmp_func_inverted(const void* ap1, const void* ap2){
498    poly p1,p2;
499  p1=*((poly*) ap1);
500  p2=*((poly*)ap2);
501
502  return -pLmCmp(p1,p2);
503}
504
505int tgb_pair_better_gen2(const void* ap,const void* bp){
506  return(-tgb_pair_better_gen(ap,bp));
507}
508int kFindDivisibleByInS_easy(kStrategy strat,const red_object & obj){
509  int i;
510  long not_sev=~obj.sev;
511  poly p=obj.p;
512  for(i=0;i<=strat->sl;i++){
513    if (pLmShortDivisibleBy(strat->S[i],strat->sevS[i],p,not_sev))
514      return i;
515  }
516  return -1;
517}
518int kFindDivisibleByInS_easy(kStrategy strat,poly p, long sev){
519  int i;
520  long not_sev=~sev;
521  for(i=0;i<=strat->sl;i++){
522    if (pLmShortDivisibleBy(strat->S[i],strat->sevS[i],p,not_sev))
523      return i;
524  }
525  return -1;
526}
527static int posInPairs (sorted_pair_node**  p, int pn, sorted_pair_node* qe,slimgb_alg* c,int an=0)
528{
529  if(pn==0) return 0;
530
531  int length=pn-1;
532  int i;
533  //int an = 0;
534  int en= length;
535
536  if (pair_better(qe,p[en],c))
537    return length+1;
538
539  while(1)
540    {
541      //if (an >= en-1)
542      if(en-1<=an)
543      {
544        if (pair_better(p[an],qe,c)) return an;
545        return en;
546      }
547      i=(an+en) / 2;
548        if (pair_better(p[i],qe,c))
549          en=i;
550      else an=i;
551    }
552}
553
554static BOOLEAN  ascending(int* i,int top){
555  if(top<1) return TRUE;
556  if(i[top]<i[top-1]) return FALSE;
557  return ascending(i,top-1);
558}
559
560sorted_pair_node**  spn_merge(sorted_pair_node** p, int pn,sorted_pair_node **q, int qn,slimgb_alg* c){
561  int i;
562  int* a= (int*) omalloc(qn*sizeof(int));
563//   int mc;
564//   PrintS("Debug\n");
565//   for(mc=0;mc<qn;mc++)
566// {
567
568//     wrp(q[mc]->lcm_of_lm);
569//     PrintS("\n");
570// }
571//    PrintS("Debug they are in\n");
572//   for(mc=0;mc<pn;mc++)
573// {
574
575//     wrp(p[mc]->lcm_of_lm);
576//     PrintS("\n");
577// }
578  int lastpos=0;
579  for(i=0;i<qn;i++){
580    lastpos=posInPairs(p,pn,q[i],c, si_max(lastpos-1,0));
581    //   cout<<lastpos<<"\n";
582    a[i]=lastpos;
583
584  }
585  if((pn+qn)>c->max_pairs){
586    p=(sorted_pair_node**) omrealloc(p,2*(pn+qn)*sizeof(sorted_pair_node*));
587    c->max_pairs=2*(pn+qn);
588  }
589  for(i=qn-1;i>=0;i--){
590    size_t size;
591    if(qn-1>i)
592      size=(a[i+1]-a[i])*sizeof(sorted_pair_node*);
593    else
594      size=(pn-a[i])*sizeof(sorted_pair_node*); //as indices begin with 0
595    memmove (p+a[i]+(1+i), p+a[i], size);
596    p[a[i]+i]=q[i];
597  }
598  omfree(a);
599  return p;
600}
601
602
603static BOOLEAN trivial_syzygie(int pos1,int pos2,poly bound,slimgb_alg* c){
604
605
606  poly p1=c->S->m[pos1];
607  poly p2=c->S->m[pos2];
608
609
610  if (pGetComp(p1) > 0 || pGetComp(p2) > 0)
611    return FALSE;
612  int i = 1;
613  poly m=NULL;
614  poly gcd1=c->gcd_of_terms[pos1];
615  poly gcd2=c->gcd_of_terms[pos2];
616
617  if((gcd1!=NULL) && (gcd2!=NULL))
618    {
619      gcd1->next=gcd2; //may ordered incorrect
620      m=gcd_of_terms(gcd1,c->r);
621      gcd1->next=NULL;
622
623    }
624
625  if (m==NULL)
626  {
627     loop
628      {
629  if (pGetExp(p1, i)+ pGetExp(p2, i) > pGetExp(bound,i))   return FALSE;
630  if (i == pVariables){
631    //PrintS("trivial");
632    return TRUE;
633  }
634  i++;
635      }
636  }
637  else
638  {
639    loop
640      {
641  if (pGetExp(p1, i)-pGetExp(m,i) + pGetExp(p2, i) > pGetExp(bound,i))   return FALSE;
642  if (i == pVariables){
643    pDelete(&m);
644    //PrintS("trivial");
645    return TRUE;
646  }
647  i++;
648      }
649  }
650
651
652
653
654}
655
656//! returns position sets w as weight
657int find_best(red_object* r,int l, int u, int &w, slimgb_alg* c){
658  int best=l;
659  int i;
660  w=r[l].guess_quality(c);
661  for(i=l+1;i<=u;i++){
662    int w2=r[i].guess_quality(c);
663    if(w2<w){
664      w=w2;
665      best=i;
666    }
667
668  }
669 return best;
670}
671
672
673void red_object::canonicalize(){
674  kBucketCanonicalize(bucket);
675
676
677}
678BOOLEAN good_has_t_rep(int i, int j,slimgb_alg* c){
679  assume(i>=0);
680    assume(j>=0);
681  if (has_t_rep(i,j,c)) return TRUE;
682  //poly lm=pOne();
683  assume (c->tmp_lm!=NULL);
684  poly lm=c->tmp_lm;
685
686  pLcm(c->S->m[i], c->S->m[j], lm);
687  pSetm(lm);
688  assume(lm!=NULL);
689  //int deciding_deg= pTotaldegree(lm);
690  int* i_con =make_connections(i,j,lm,c);
691  //p_Delete(&lm,c->r);
692
693
694  for (int n=0;((n<c->n) && (i_con[n]>=0));n++){
695    if (i_con[n]==j){
696      now_t_rep(i,j,c);
697      omfree(i_con);
698
699      return TRUE;
700    }
701  }
702  omfree(i_con);
703
704  return FALSE;
705}
706BOOLEAN lenS_correct(kStrategy strat){
707  int i;
708  for(i=0;i<=strat->sl;i++){
709    if (strat->lenS[i]!=pLength(strat->S[i]))
710      return FALSE;
711  }
712  return TRUE;
713}
714
715
716static void cleanS(kStrategy strat, slimgb_alg* c){
717  int i=0;
718  LObject P;
719  while(i<=strat->sl){
720    P.p=strat->S[i];
721    P.sev=strat->sevS[i];
722    if(kFindDivisibleByInS(strat->S,strat->sevS,strat->sl,&P)!=i)
723    {
724      deleteInS(i,strat);
725      //remember destroying poly
726      BOOLEAN found=FALSE;
727      int j;
728      for(j=0;j<c->n;j++)
729  if(c->S->m[j]==P.p)
730  {
731    found=TRUE;
732    break;
733  }
734      if (!found)
735  pDelete(&P.p);
736      //remember additional reductors
737    }
738    else i++;
739  }
740}
741static int bucket_guess(kBucket* bucket){
742  int sum=0;
743  int i;
744  for (i=bucket->buckets_used;i>=0;i--){
745    if(bucket->buckets[i])
746       sum+=bucket->buckets_length[i];
747  }
748  return sum;
749}
750
751
752
753
754
755
756static int add_to_reductors(slimgb_alg* c, poly h, int len, BOOLEAN simplified){
757  //inDebug(h);
758  assume(lenS_correct(c->strat));
759  assume(len==pLength(h));
760  int i;
761//   if (c->isDifficultField)
762//        i=simple_posInS(c->strat,h,pSLength(h,len),c->isDifficultField);
763//   else
764//     i=simple_posInS(c->strat,h,len,c->isDifficultField);
765
766  LObject P; memset(&P,0,sizeof(P));
767  P.tailRing=c->r;
768  P.p=h; /*p_Copy(h,c->r);*/
769  P.FDeg=pFDeg(P.p,c->r);
770  if (!(simplified)){
771      if (!rField_is_Zp(c->r)){
772        pCleardenom(P.p);
773        pContent(P.p); //is a duplicate call, but belongs here
774
775      }
776      else
777        pNorm(P.p);
778    pNormalize(P.p);
779  }
780  wlen_type pq=pQuality(h,c,len);
781  i=simple_posInS(c->strat,h,len,pq);
782  c->strat->enterS(P,i,c->strat,-1);
783
784
785
786  c->strat->lenS[i]=len;
787  assume(pLength(c->strat->S[i])==c->strat->lenS[i]);
788  if(c->strat->lenSw!=NULL)
789    c->strat->lenSw[i]=pq;
790
791  return i;
792
793}
794static void length_one_crit(slimgb_alg* c, int pos, int len)
795{
796  if (c->nc)
797    return;
798  if (len==1)
799  {
800    int i;
801    for ( i=0;i<pos;i++)
802    {
803      if (c->lengths[i]==1)
804        c->states[pos][i]=HASTREP;
805    }
806    for ( i=pos+1;i<c->n;i++){
807      if (c->lengths[i]==1)
808        c->states[i][pos]=HASTREP;
809    }
810    if (!c->nc)
811      shorten_tails(c,c->S->m[pos]);
812  }
813}
814
815
816static void move_forward_in_S(int old_pos, int new_pos,kStrategy strat)
817{
818  assume(old_pos>=new_pos);
819  poly p=strat->S[old_pos];
820  int ecart=strat->ecartS[old_pos];
821  long sev=strat->sevS[old_pos];
822  int s_2_r=strat->S_2_R[old_pos];
823  int length=strat->lenS[old_pos];
824  assume(length==pLength(strat->S[old_pos]));
825  wlen_type length_w;
826  if(strat->lenSw!=NULL)
827    length_w=strat->lenSw[old_pos];
828  int i;
829  for (i=old_pos; i>new_pos; i--)
830  {
831    strat->S[i] = strat->S[i-1];
832    strat->ecartS[i] = strat->ecartS[i-1];
833    strat->sevS[i] = strat->sevS[i-1];
834    strat->S_2_R[i] = strat->S_2_R[i-1];
835  }
836  if (strat->lenS!=NULL)
837    for (i=old_pos; i>new_pos; i--)
838      strat->lenS[i] = strat->lenS[i-1];
839  if (strat->lenSw!=NULL)
840    for (i=old_pos; i>new_pos; i--)
841      strat->lenSw[i] = strat->lenSw[i-1];
842
843  strat->S[new_pos]=p;
844  strat->ecartS[new_pos]=ecart;
845  strat->sevS[new_pos]=sev;
846  strat->S_2_R[new_pos]=s_2_r;
847  strat->lenS[new_pos]=length;
848  if(strat->lenSw!=NULL)
849    strat->lenSw[new_pos]=length_w;
850  //assume(lenS_correct(strat));
851}
852
853static int* make_connections(int from, int to, poly bound, slimgb_alg* c)
854{
855  ideal I=c->S;
856  int* cans=(int*) omalloc(c->n*sizeof(int));
857  int* connected=(int*) omalloc(c->n*sizeof(int));
858  cans[0]=to;
859  int cans_length=1;
860  connected[0]=from;
861  int last_cans_pos=-1;
862  int connected_length=1;
863  long neg_bounds_short= ~p_GetShortExpVector(bound,c->r);
864
865  int not_yet_found=cans_length;
866  int con_checked=0;
867  int pos;
868
869  while(TRUE){
870    if ((con_checked<connected_length)&& (not_yet_found>0)){
871      pos=connected[con_checked];
872      for(int i=0;i<cans_length;i++){
873        if (cans[i]<0) continue;
874        //FIXME: triv. syz. does not hold on noncommutative, check it for modules
875        if ((has_t_rep(pos,cans[i],c)) ||((!rIsPluralRing(c->r))&&(trivial_syzygie(pos,cans[i],bound,c))))
876{
877
878          connected[connected_length]=cans[i];
879          connected_length++;
880          cans[i]=-1;
881          --not_yet_found;
882
883          if (connected[connected_length-1]==to){
884            if (connected_length<c->n){
885              connected[connected_length]=-1;
886            }
887            omfree(cans);
888            return connected;
889          }
890        }
891      }
892      con_checked++;
893    }
894    else
895    {
896      for(last_cans_pos++;last_cans_pos<=c->n;last_cans_pos++){
897        if (last_cans_pos==c->n){
898          if (connected_length<c->n){
899            connected[connected_length]=-1;
900          }
901          omfree(cans);
902          return connected;
903        }
904        if ((last_cans_pos==from)||(last_cans_pos==to))
905          continue;
906        if(p_LmShortDivisibleBy(I->m[last_cans_pos],c->short_Exps[last_cans_pos],bound,neg_bounds_short,c->r)){
907          cans[cans_length]=last_cans_pos;
908          cans_length++;
909          break;
910        }
911      }
912      not_yet_found++;
913      for (int i=0;i<con_checked;i++){
914        if (has_t_rep(connected[i],last_cans_pos,c)){
915
916          connected[connected_length]=last_cans_pos;
917          connected_length++;
918          cans[cans_length-1]=-1;
919
920          --not_yet_found;
921          if (connected[connected_length-1]==to){
922            if (connected_length<c->n){
923              connected[connected_length]=-1;
924            }
925
926            omfree(cans);
927            return connected;
928          }
929          break;
930        }
931      }
932    }
933  }
934  if (connected_length<c->n){
935    connected[connected_length]=-1;
936  }
937
938  omfree(cans);
939  return connected;
940}
941#ifdef HEAD_BIN
942static inline poly p_MoveHead(poly p, omBin b)
943{
944  poly np;
945  omTypeAllocBin(poly, np, b);
946  memmove(np, p, b->sizeW*sizeof(long));
947  omFreeBinAddr(p);
948  return np;
949}
950#endif
951
952static void replace_pair(int & i, int & j,slimgb_alg* c)
953{
954  if (i<0) return;
955  c->soon_free=NULL;
956  int syz_deg;
957  poly lm=pOne();
958
959  pLcm(c->S->m[i], c->S->m[j], lm);
960  pSetm(lm);
961
962  int* i_con =make_connections(i,j,lm,c);
963
964  for (int n=0;((n<c->n) && (i_con[n]>=0));n++){
965    if (i_con[n]==j){
966      now_t_rep(i,j,c);
967      omfree(i_con);
968      p_Delete(&lm,c->r);
969      return;
970    }
971  }
972
973  int* j_con =make_connections(j,i,lm,c);
974
975//   if(c->n>1){
976//     if (i_con[1]>=0)
977//       i=i_con[1];
978//     else {
979//       if (j_con[1]>=0)
980//         j=j_con[1];
981//     }
982 // }
983
984  int sugar=pTotaldegree(lm);
985  p_Delete(&lm, c->r);
986    if(c->T_deg_full)//Sugar
987    {
988      int t_i=c->T_deg_full[i]-c->T_deg[i];
989      int t_j=c->T_deg_full[j]-c->T_deg[j];
990      sugar+=si_max(t_i,t_j);
991      //Print("\n max: %d\n",max(t_i,t_j));
992    }
993
994
995
996
997
998  for (int m=0;((m<c->n) && (i_con[m]>=0));m++){
999    if(c->T_deg_full!=NULL){
1000        int s1=c->T_deg_full[i_con[m]]+syz_deg-c->T_deg[i_con[m]];
1001        if (s1>sugar) continue;
1002    }
1003    if (c->weighted_lengths[i_con[m]]<c->weighted_lengths[i])
1004        i=i_con[m];
1005    }
1006    for (int m=0;((m<c->n) && (j_con[m]>=0));m++){
1007        if (c->T_deg_full!=NULL){
1008        int s1=c->T_deg_full[j_con[m]]+syz_deg-c->T_deg[j_con[m]];
1009        if (s1>sugar) continue;}
1010        if (c->weighted_lengths[j_con[m]]<c->weighted_lengths[j])
1011            j=j_con[m];
1012    }
1013
1014     //can also try dependend search
1015  omfree(i_con);
1016  omfree(j_con);
1017  return;
1018}
1019
1020
1021static void add_later(poly p, char* prot, slimgb_alg* c){
1022    int i=0;
1023    //check, if it is already in the queue
1024
1025
1026    while(c->add_later->m[i]!=NULL){
1027        if (p_LmEqual(c->add_later->m[i],p,c->r))
1028            return;
1029        i++;
1030    }
1031    if (TEST_OPT_PROT)
1032        PrintS(prot);
1033    c->add_later->m[i]=p;
1034}
1035static int simple_posInS (kStrategy strat, poly p,int len, wlen_type wlen)
1036{
1037
1038
1039  if(strat->sl==-1) return 0;
1040  if (strat->lenSw) return pos_helper(strat,p,(wlen_type) wlen,(wlen_set) strat->lenSw,strat->S);
1041  return pos_helper(strat,p,len,strat->lenS,strat->S);
1042
1043}
1044
1045/*2
1046 *if the leading term of p
1047 *divides the leading term of some S[i] it will be canceled
1048 */
1049static inline void clearS (poly p, unsigned long p_sev,int l, int* at, int* k,
1050                           kStrategy strat)
1051{
1052  assume(p_sev == pGetShortExpVector(p));
1053  if (!pLmShortDivisibleBy(p,p_sev, strat->S[*at], ~ strat->sevS[*at])) return;
1054  if (l>=strat->lenS[*at]) return;
1055  if (TEST_OPT_PROT)
1056    PrintS("!");
1057  mflush();
1058  //pDelete(&strat->S[*at]);
1059  deleteInS((*at),strat);
1060  (*at)--;
1061  (*k)--;
1062//  assume(lenS_correct(strat));
1063}
1064
1065
1066
1067static int iq_crit(const void* ap,const void* bp){
1068
1069  sorted_pair_node* a=*((sorted_pair_node**)ap);
1070  sorted_pair_node* b=*((sorted_pair_node**)bp);
1071  assume(a->i>a->j);
1072  assume(b->i>b->j);
1073
1074
1075  if (a->deg<b->deg) return -1;
1076  if (a->deg>b->deg) return 1;
1077  int comp=pLmCmp(a->lcm_of_lm, b->lcm_of_lm);
1078  if(comp!=0)
1079    return comp;
1080  if (a->expected_length<b->expected_length) return -1;
1081  if (a->expected_length>b->expected_length) return 1;
1082  if (a->j>b->j) return 1;
1083  if (a->j<b->j) return -1;
1084  return 0;
1085}
1086static wlen_type coeff_mult_size_estimate(int s1, int s2, ring r){
1087    if (rField_is_Q(r)) return s1+s2;
1088    else return s1*s2;
1089}
1090static wlen_type pair_weighted_length(int i, int j, slimgb_alg* c){
1091    if ((c->isDifficultField) && (c->eliminationProblem))  {
1092        int c1=slim_nsize(p_GetCoeff(c->S->m[i],c->r),c->r);
1093        int c2=slim_nsize(p_GetCoeff(c->S->m[j],c->r),c->r);
1094        wlen_type el1=c->weighted_lengths[i]/c1;
1095        assume(el1!=0);
1096        assume(c->weighted_lengths[i] %c1==0);
1097        wlen_type el2=c->weighted_lengths[j]/c2;
1098        assume(el2!=0);
1099        assume(c->weighted_lengths[j] %c2==0);
1100        //should be * for function fields
1101        //return (c1+c2) * (el1+el2-2);
1102        wlen_type res=coeff_mult_size_estimate(c1,c2,c->r);
1103        res*=el1+el2-2;
1104        return res;
1105
1106    }
1107    if (c->isDifficultField) {
1108        //int cs=slim_nsize(p_GetCoeff(c->S->m[i],c->r),c->r)+
1109        //    slim_nsize(p_GetCoeff(c->S->m[j],c->r),c->r);
1110        if(!(TEST_V_COEFSTRAT)){
1111        wlen_type cs=
1112            coeff_mult_size_estimate(
1113                slim_nsize(p_GetCoeff(c->S->m[i],c->r),c->r),
1114                slim_nsize(p_GetCoeff(c->S->m[j],c->r),c->r),c->r);
1115        return (wlen_type)(c->lengths[i]+c->lengths[j]-2)*
1116            (wlen_type)cs;}
1117            else {
1118
1119            wlen_type cs=
1120            coeff_mult_size_estimate(
1121                slim_nsize(p_GetCoeff(c->S->m[i],c->r),c->r),
1122                slim_nsize(p_GetCoeff(c->S->m[j],c->r),c->r),c->r);
1123            cs*=cs;
1124        return (wlen_type)(c->lengths[i]+c->lengths[j]-2)*
1125            (wlen_type)cs;
1126            }
1127    }
1128    if (c->eliminationProblem) {
1129
1130        return (c->weighted_lengths[i]+c->weighted_lengths[j]-2);
1131    }
1132    return c->lengths[i]+c->lengths[j]-2;
1133
1134}
1135sorted_pair_node** add_to_basis_ideal_quotient(poly h, slimgb_alg* c, int* ip)
1136{
1137
1138  assume(h!=NULL);
1139  poly got=gcd_of_terms(h,c->r);
1140  if((got!=NULL) &&(TEST_V_UPTORADICAL)) {
1141    poly copy=p_Copy(got,c->r);
1142    //p_wrp(got,c->r);
1143    BOOLEAN changed=monomial_root(got,c->r);
1144    if (changed)
1145    {
1146         poly div_by=pDivide(copy, got);
1147         poly iter=h;
1148         while(iter){
1149            pExpVectorSub(iter,div_by);
1150            pIter(iter);
1151         }
1152         p_Delete(&div_by, c->r);
1153         PrintS("U");
1154    }
1155    p_Delete(&copy,c->r);
1156  }
1157
1158#define ENLARGE(pointer, type) pointer=(type*) omrealloc(pointer, c->array_lengths*sizeof(type))
1159//  BOOLEAN corr=lenS_correct(c->strat);
1160  BOOLEAN R_found=FALSE;
1161  void* hp;
1162
1163  ++(c->n);
1164  ++(c->S->ncols);
1165  int i,j;
1166  i=c->n-1;
1167  sorted_pair_node** nodes=(sorted_pair_node**) omalloc(sizeof(sorted_pair_node*)*i);
1168  int spc=0;
1169  if(c->n>c->array_lengths){
1170    c->array_lengths=c->array_lengths*2;
1171    assume(c->array_lengths>=c->n);
1172    ENLARGE(c->T_deg, int);
1173    ENLARGE(c->tmp_pair_lm,poly);
1174    ENLARGE(c->tmp_spn,sorted_pair_node*);
1175
1176    ENLARGE(c->short_Exps,long);
1177    ENLARGE(c->lengths,int);
1178    #ifndef HAVE_BOOST
1179    #ifndef USE_STDVECBOOL
1180   
1181    ENLARGE(c->states, char*);
1182    #endif
1183    #endif
1184    ENLARGE(c->gcd_of_terms,poly);
1185    //if (c->weighted_lengths!=NULL) {
1186    ENLARGE(c->weighted_lengths,wlen_type);
1187    //}
1188    //ENLARGE(c->S->m,poly);
1189
1190  }
1191  pEnlargeSet(&c->S->m,c->n-1,1);
1192  if (c->T_deg_full)
1193    ENLARGE(c->T_deg_full,int);
1194  c->T_deg[i]=pTotaldegree(h);
1195  if(c->T_deg_full){
1196    c->T_deg_full[i]=pTotaldegree_full(h);
1197  }
1198
1199
1200  c->tmp_pair_lm[i]=pOne_Special(c->r);
1201
1202
1203  c->tmp_spn[i]=(sorted_pair_node*) omalloc(sizeof(sorted_pair_node));
1204
1205
1206  c->lengths[i]=pLength(h);
1207
1208  //necessary for correct weighted length
1209
1210   if (!rField_is_Zp(c->r)){
1211    pCleardenom(h);
1212    pContent(h); //is a duplicate call, but belongs here
1213
1214  }
1215  else
1216    pNorm(h);
1217  pNormalize(h);
1218
1219  c->weighted_lengths[i]=pQuality(h, c, c->lengths[i]);
1220  c->gcd_of_terms[i]=got;
1221  #ifdef HAVE_BOOST
1222    c->states.push_back(dynamic_bitset<>(i));
1223
1224  #else
1225  #ifdef USE_STDVECBOOL
1226   
1227    c->states.push_back(vector<bool>(i));
1228   
1229 
1230  #else
1231  if (i>0)
1232    c->states[i]=(char*)  omalloc(i*sizeof(char));
1233  else
1234    c->states[i]=NULL;
1235  #endif
1236  #endif
1237
1238  c->S->m[i]=h;
1239  c->short_Exps[i]=p_GetShortExpVector(h,c->r);
1240
1241#undef ENLARGE
1242  if (p_GetComp(h,currRing)<=c->syz_comp){
1243  for (j=0;j<i;j++){
1244
1245
1246    #ifndef HAVE_BOOST
1247    c->states[i][j]=UNCALCULATED;
1248    #endif
1249    assume(p_LmDivisibleBy(c->S->m[i],c->S->m[j],c->r)==
1250     p_LmShortDivisibleBy(c->S->m[i],c->short_Exps[i],c->S->m[j],~(c->short_Exps[j]),c->r));
1251
1252
1253    if (_p_GetComp(c->S->m[i],c->r)!=_p_GetComp(c->S->m[j],c->r)){
1254      //c->states[i][j]=UNCALCULATED;
1255      //WARNUNG: be careful
1256      continue;
1257    } else
1258    if ((!c->nc) && (c->lengths[i]==1) && (c->lengths[j]==1)){
1259      c->states[i][j]=HASTREP;
1260
1261      }
1262    else if ((!(c->nc)) &&  (pHasNotCF(c->S->m[i],c->S->m[j])))
1263    {
1264      c->easy_product_crit++;
1265      c->states[i][j]=HASTREP;
1266      continue;
1267    }
1268    else if(extended_product_criterion(c->S->m[i],c->gcd_of_terms[i],c->S->m[j],c->gcd_of_terms[j],c))
1269    {
1270      c->states[i][j]=HASTREP;
1271      c->extended_product_crit++;
1272
1273      //PrintS("E");
1274    }
1275      //  if (c->states[i][j]==UNCALCULATED){
1276
1277    if ((TEST_V_FINDMONOM) &&(!c->nc)) {
1278        //PrintS("COMMU");
1279       //  if (c->lengths[i]==c->lengths[j]){
1280//             poly short_s=ksCreateShortSpoly(c->S->m[i],c->S->m[j],c->r);
1281//             if (short_s==NULL){
1282//                 c->states[i][j]=HASTREP;
1283//             } else
1284//             {
1285//                 p_Delete(&short_s, currRing);
1286//             }
1287//         }
1288        if (c->lengths[i]+c->lengths[j]==3){
1289
1290
1291             poly short_s=ksCreateShortSpoly(c->S->m[i],c->S->m[j],c->r);
1292            if (short_s==NULL){
1293                c->states[i][j]=HASTREP;
1294            } else
1295            {
1296                assume(pLength(short_s)==1);
1297                if (TEST_V_UPTORADICAL)
1298                   monomial_root(short_s,c->r);
1299                int iS=
1300                   kFindDivisibleByInS_easy(c->strat,short_s, p_GetShortExpVector(short_s,c->r));
1301                if (iS<0){
1302                    //PrintS("N");
1303                    if (TRUE) {
1304                    c->states[i][j]=HASTREP;
1305                    add_later(short_s,"N",c);
1306                    } else p_Delete(&short_s,currRing);
1307                }
1308                else {
1309                    if (c->strat->lenS[iS]>1){
1310                        //PrintS("O");
1311                        if (TRUE) {
1312                        c->states[i][j]=HASTREP;
1313                        add_later(short_s,"O",c);
1314                        } else p_Delete(&short_s,currRing);
1315                    }
1316                    else
1317                     p_Delete(&short_s, currRing);
1318                     c->states[i][j]=HASTREP;
1319                }
1320
1321
1322            }
1323        }
1324    }
1325      //    if (short_s)
1326      //    {
1327    assume(spc<=j);
1328    sorted_pair_node* s=c->tmp_spn[spc];//(sorted_pair_node*) omalloc(sizeof(sorted_pair_node));
1329    s->i=si_max(i,j);
1330    s->j=si_min(i,j);
1331    assume(s->j==j);
1332    s->expected_length=pair_weighted_length(i,j,c);//c->lengths[i]+c->lengths[j]-2;
1333
1334    poly lm=c->tmp_pair_lm[spc];//=pOne_Special();
1335
1336    pLcm(c->S->m[i], c->S->m[j], lm);
1337    pSetm(lm);
1338    s->deg=pTotaldegree(lm);
1339
1340    if(c->T_deg_full)//Sugar
1341    {
1342      int t_i=c->T_deg_full[s->i]-c->T_deg[s->i];
1343      int t_j=c->T_deg_full[s->j]-c->T_deg[s->j];
1344      s->deg+=si_max(t_i,t_j);
1345      //Print("\n max: %d\n",max(t_i,t_j));
1346    }
1347    s->lcm_of_lm=lm;
1348    //          pDelete(&short_s);
1349    //assume(lm!=NULL);
1350    nodes[spc]=s;
1351    spc++;
1352
1353  // }
1354  //else
1355  //{
1356        //c->states[i][j]=HASTREP;
1357  //}
1358  }
1359  }//if syz_comp end
1360 
1361 
1362 
1363 
1364  assume(spc<=i);
1365  //now ideal quotient crit
1366  qsort(nodes,spc,sizeof(sorted_pair_node*),iq_crit);
1367
1368    sorted_pair_node** nodes_final=(sorted_pair_node**) omalloc(sizeof(sorted_pair_node*)*i);
1369  int spc_final=0;
1370  j=0;
1371  while(j<spc)
1372  {
1373    int lower=j;
1374    int upper;
1375    BOOLEAN has=FALSE;
1376    for(upper=lower+1;upper<spc;upper++)
1377    {
1378
1379      if(!pLmEqual(nodes[lower]->lcm_of_lm,nodes[upper]->lcm_of_lm))
1380      {
1381  break;
1382      }
1383      if (has_t_rep(nodes[upper]->i,nodes[upper]->j,c))
1384  has=TRUE;
1385
1386    }
1387    upper=upper-1;
1388    int z;
1389    assume(spc_final<=j);
1390    for(z=0;z<spc_final;z++)
1391    {
1392      if(p_LmDivisibleBy(nodes_final[z]->lcm_of_lm,nodes[lower]->lcm_of_lm,c->r))
1393      {
1394  has=TRUE;
1395  break;
1396      }
1397    }
1398
1399    if(has)
1400    {
1401      for(;lower<=upper;lower++)
1402      {
1403  //free_sorted_pair_node(nodes[lower],c->r);
1404  //omfree(nodes[lower]);
1405  nodes[lower]=NULL;
1406      }
1407      j=upper+1;
1408      continue;
1409    }
1410    else
1411    {
1412      nodes[lower]->lcm_of_lm=pCopy(nodes[lower]->lcm_of_lm);
1413      assume(_p_GetComp(c->S->m[nodes[lower]->i],c->r)==_p_GetComp(c->S->m[nodes[lower]->j],c->r));
1414      nodes_final[spc_final]=(sorted_pair_node*) omalloc(sizeof(sorted_pair_node));
1415
1416      *(nodes_final[spc_final++])=*(nodes[lower]);
1417      //c->tmp_spn[nodes[lower]->j]=(sorted_pair_node*) omalloc(sizeof(sorted_pair_node));
1418      nodes[lower]=NULL;
1419      for(lower=lower+1;lower<=upper;lower++)
1420      {
1421  //  free_sorted_pair_node(nodes[lower],c->r);
1422  //omfree(nodes[lower]);
1423  nodes[lower]=NULL;
1424      }
1425      j=upper+1;
1426      continue;
1427    }
1428  }
1429
1430  //  Print("i:%d,spc_final:%d",i,spc_final);
1431
1432
1433
1434
1435  assume(spc_final<=spc);
1436  omfree(nodes);
1437  nodes=NULL;
1438
1439  add_to_reductors(c, h, c->lengths[c->n-1], TRUE);
1440  //i=posInS(c->strat,c->strat->sl,h,0 ecart);
1441  if (!(c->nc)){
1442    if (c->lengths[c->n-1]==1)
1443      shorten_tails(c,c->S->m[c->n-1]);
1444  }
1445  //you should really update c->lengths, c->strat->lenS, and the oder of polys in strat if you sort after lengths
1446
1447  //for(i=c->strat->sl; i>0;i--)
1448  //  if(c->strat->lenS[i]<c->strat->lenS[i-1]) printf("fehler bei %d\n",i);
1449  if (c->Rcounter>50) {
1450    c->Rcounter=0;
1451    cleanS(c->strat,c);
1452  }
1453  if(!ip){
1454    qsort(nodes_final,spc_final,sizeof(sorted_pair_node*),tgb_pair_better_gen2);
1455
1456
1457    c->apairs=spn_merge(c->apairs,c->pair_top+1,nodes_final,spc_final,c);
1458    c->pair_top+=spc_final;
1459    clean_top_of_pair_list(c);
1460    omfree(nodes_final);
1461    return NULL;
1462  }
1463  {
1464    *ip=spc_final;
1465    return nodes_final;
1466  }
1467
1468
1469
1470}
1471
1472
1473static poly redNF2 (poly h,slimgb_alg* c , int &len, number&  m,int n)
1474{
1475  m=nInit(1);
1476  if (h==NULL) return NULL;
1477
1478  assume(len==pLength(h));
1479  kStrategy strat=c->strat;
1480  if (0 > strat->sl)
1481  {
1482    return h;
1483  }
1484  int j;
1485
1486  LObject P(h);
1487  P.SetShortExpVector();
1488  P.bucket = kBucketCreate(currRing);
1489  // BOOLEAN corr=lenS_correct(strat);
1490  kBucketInit(P.bucket,P.p,len /*pLength(P.p)*/);
1491  //wlen_set lenSw=(wlen_set) c->strat->lenS;
1492  //FIXME: plainly wrong
1493  //strat->lenS;
1494  //if (strat->lenSw!=NULL)
1495  //  lenSw=strat->lenSw;
1496  //int max_pos=simple_posInS(strat,P.p);
1497  loop
1498    {
1499
1500      j=kFindDivisibleByInS(strat->S,strat->sevS,strat->sl,&P);
1501      if ((j>=0) && ((!n)||
1502        ((strat->lenS[j]<=n) &&
1503         ((strat->lenSw==NULL)||
1504         (strat->lenSw[j]<=n)))))
1505      {
1506
1507
1508        nNormalize(pGetCoeff(P.p));
1509#ifdef KDEBUG
1510        if (TEST_OPT_DEBUG)
1511        {
1512          PrintS("red:");
1513          wrp(h);
1514          PrintS(" with ");
1515          wrp(strat->S[j]);
1516        }
1517#endif
1518
1519        number coef=kBucketPolyRed(P.bucket,strat->S[j],
1520                                   strat->lenS[j]/*pLength(strat->S[j])*/,
1521                                   strat->kNoether);
1522  number m2=nMult(m,coef);
1523  nDelete(&m);
1524  m=m2;
1525        nDelete(&coef);
1526        h = kBucketGetLm(P.bucket);
1527
1528  if (h==NULL) {
1529    len=0;
1530    kBucketDestroy(&P.bucket);
1531    return
1532    NULL;}
1533        P.p=h;
1534        P.t_p=NULL;
1535        P.SetShortExpVector();
1536#ifdef KDEBUG
1537        if (TEST_OPT_DEBUG)
1538        {
1539          PrintS("\nto:");
1540          wrp(h);
1541          PrintLn();
1542        }
1543#endif
1544      }
1545      else
1546      {
1547        kBucketClear(P.bucket,&(P.p),&len);
1548        kBucketDestroy(&P.bucket);
1549        pNormalize(P.p);
1550  assume(len==(pLength(P.p)));
1551        return P.p;
1552      }
1553    }
1554}
1555
1556
1557
1558static poly redTailShort(poly h, kStrategy strat){
1559  if (h==NULL) return NULL;//n_Init(1,currRing);
1560  if (TEST_V_MODPSOLVSB){
1561    bit_reduce(pNext(h), strat->tailRing);
1562  }
1563  int sl=strat->sl;
1564  int i;
1565  int len=pLength(h);
1566  for(i=0;i<=strat->sl;i++){
1567    if((strat->lenS[i]>2) || ((strat->lenSw!=NULL) && (strat->lenSw[i]>2)))
1568      break;
1569  }
1570  return(redNFTail(h,i-1,strat, len));
1571}
1572
1573static void line_of_extended_prod(int fixpos,slimgb_alg* c){
1574    if (c->gcd_of_terms[fixpos]==NULL)
1575  {
1576    c->gcd_of_terms[fixpos]=gcd_of_terms(c->S->m[fixpos],c->r);
1577    if (c->gcd_of_terms[fixpos])
1578    {
1579      int i;
1580      for(i=0;i<fixpos;i++)
1581        if((c->states[fixpos][i]!=HASTREP)&& (extended_product_criterion(c->S->m[fixpos],c->gcd_of_terms[fixpos], c->S->m[i],c->gcd_of_terms[i],c)))
1582{
1583          c->states[fixpos][i]=HASTREP;
1584    c->extended_product_crit++;
1585}
1586      for(i=fixpos+1;i<c->n;i++)
1587        if((c->states[i][fixpos]!=HASTREP)&& (extended_product_criterion(c->S->m[fixpos],c->gcd_of_terms[fixpos], c->S->m[i],c->gcd_of_terms[i],c)))
1588  {        c->states[i][fixpos]=HASTREP;
1589  c->extended_product_crit++;
1590  }
1591    }
1592  }
1593}
1594static void c_S_element_changed_hook(int pos, slimgb_alg* c){
1595  length_one_crit(c,pos, c->lengths[pos]);
1596  if (!c->nc)
1597    line_of_extended_prod(pos,c);
1598}
1599class poly_tree_node {
1600public:
1601  poly p;
1602  poly_tree_node* l;
1603  poly_tree_node* r;
1604  int n;
1605  poly_tree_node(int sn):l(NULL),r(NULL),n(sn){}
1606};
1607class exp_number_builder{
1608public:
1609  poly_tree_node* top_level;
1610  int n;
1611  int get_n(poly p);
1612  exp_number_builder():top_level(0),n(0){}
1613};
1614int exp_number_builder::get_n(poly p){
1615  poly_tree_node** node=&top_level;
1616  while(*node!=NULL){
1617    int c=pLmCmp(p,(*node)->p);
1618    if (c==0) return (*node)->n;
1619    if (c==-1) node=&((*node)->r);
1620    else
1621      node=&((*node)->l);
1622  }
1623  (*node)= new poly_tree_node(n);
1624  n++;
1625  (*node)->p=pLmInit(p);
1626  return (*node)->n;
1627}
1628
1629//mac_polys exp are smaller iff they are greater by monomial ordering
1630//corresponding to solving linear equations notation
1631
1632//! obsolete
1633struct int_poly_pair{
1634  poly p;
1635  int n;
1636};
1637
1638
1639//! obsolete
1640void t2ippa_rec(poly* ip,int* ia, poly_tree_node* k, int &offset){
1641    if(!k) return;
1642    t2ippa_rec(ip,ia,k->l,offset);
1643    ip[offset]=k->p;
1644    ia[k->n]=offset;
1645    ++offset;
1646
1647    t2ippa_rec(ip,ia,k->r,offset);
1648    delete k;
1649  }
1650
1651//! obsolete
1652void t2ippa(poly* ip,int* ia,exp_number_builder & e){
1653
1654  int o=0;
1655  t2ippa_rec(ip,ia,e.top_level,o);
1656}
1657int anti_poly_order(const void* a, const void* b){
1658  return -pLmCmp(((int_poly_pair*) a)->p,((int_poly_pair*) b)->p );
1659}
1660
1661BOOLEAN is_valid_ro(red_object & ro){
1662  red_object r2=ro;
1663  ro.validate();
1664  if ((r2.p!=ro.p)||(r2.sev!=ro.sev)) return FALSE;
1665  return TRUE;
1666}
1667
1668
1669
1670static void go_on (slimgb_alg* c){
1671  //set limit of 1000 for multireductions, at the moment for
1672  //programming reasons
1673  int i=0;
1674  c->average_length=0;
1675  for(i=0;i<c->n;i++){
1676    c->average_length+=c->lengths[i];
1677  }
1678  c->average_length=c->average_length/c->n;
1679  i=0;
1680  poly* p=(poly*) omalloc((PAR_N+1)*sizeof(poly));//nullterminated
1681
1682  int curr_deg=-1;
1683  while(i<PAR_N){
1684    sorted_pair_node* s=top_pair(c);//here is actually chain criterium done
1685
1686
1687    if (!s) break;
1688
1689    if(curr_deg>=0){
1690      if (s->deg >curr_deg) break;
1691    }
1692
1693    else curr_deg=s->deg;
1694    quick_pop_pair(c);
1695    if(s->i>=0){
1696      //be careful replace_pair use createShortSpoly which is not noncommutative
1697      now_t_rep(s->i,s->j,c);
1698    replace_pair(s->i,s->j,c);
1699
1700    if(s->i==s->j) {
1701      free_sorted_pair_node(s,c->r);
1702      continue;
1703    }
1704    now_t_rep(s->i,s->j,c);
1705    }
1706    poly h;
1707    if(s->i>=0){
1708      if (!c->nc)
1709  h=ksOldCreateSpoly(c->S->m[s->i], c->S->m[s->j], NULL, c->r);
1710      else
1711  h= nc_CreateSpoly(c->S->m[s->i], c->S->m[s->j], NULL, c->r);
1712    }
1713    else
1714      h=s->lcm_of_lm;
1715    // if(s->i>=0)
1716//       now_t_rep(s->j,s->i,c);
1717    number coef;
1718    int mlen=pLength(h);
1719    if (!c->nc){
1720      h=redNF2(h,c,mlen,coef,2);
1721      redTailShort(h,c->strat);
1722      nDelete(&coef);
1723    }
1724    free_sorted_pair_node(s,c->r);
1725    if(!h) continue;
1726    int len=pLength(h);
1727    p[i]=h;
1728
1729    i++;
1730  }
1731  p[i]=NULL;
1732//  pre_comp(p,i,c);
1733  if(i==0){
1734    omfree(p);
1735    return;
1736  }
1737  #ifdef TGB_RESORT_PAIRS
1738  c->replaced=new bool[c->n];
1739  c->used_b=FALSE;
1740  #endif
1741  red_object* buf=(red_object*) omalloc(i*sizeof(red_object));
1742  c->normal_forms+=i;
1743  int j;
1744  for(j=0;j<i;j++){
1745    buf[j].p=p[j];
1746    buf[j].sev=pGetShortExpVector(p[j]);
1747    buf[j].bucket = kBucketCreate(currRing);
1748
1749    int len=pLength(p[j]);
1750    kBucketInit(buf[j].bucket,buf[j].p,len);
1751    buf[j].initial_quality=buf[j].guess_quality(c);
1752    assume(buf[j].initial_quality>=0);
1753  }
1754  omfree(p);
1755  qsort(buf,i,sizeof(red_object),red_object_better_gen);
1756//    Print("\ncurr_deg:%i\n",curr_deg);
1757  if (TEST_OPT_PROT)
1758  {
1759    Print("%dM[%d,",curr_deg,i);
1760  }
1761#ifdef FIND_DETERMINISTIC
1762  c->modifiedS=(BOOLEAN*) omalloc((c->strat->sl+1)*sizeof(BOOLEAN));
1763  c->expandS=(poly*) omalloc((1)*sizeof(poly));
1764  c->expandS[0]=NULL;
1765  int z2;
1766  for(z2=0;z2<=c->strat->sl;z2++)
1767    c->modifiedS[z2]=FALSE;
1768#endif
1769  multi_reduction(buf, i, c);
1770  #ifdef TGB_RESORT_PAIRS
1771  if (c->used_b) {
1772    if (TEST_OPT_PROT)
1773        PrintS("B");
1774    int e;
1775    for(e=0;e<=c->pair_top;e++){
1776        if(c->apairs[e]->i<0) continue;
1777        assume(c->apairs[e]->j>=0);
1778        if ((c->replaced[c->apairs[e]->i])||(c->replaced[c->apairs[e]->j])) {
1779            sorted_pair_node* s=c->apairs[e];
1780            s->expected_length=pair_weighted_length(s->i,s->j,c);
1781        }
1782    }
1783    qsort(c->apairs,c->pair_top+1,sizeof(sorted_pair_node*),tgb_pair_better_gen2);
1784  }
1785  #endif
1786#ifdef TGB_DEBUG
1787 {
1788   int k;
1789   for(k=0;k<i;k++)
1790   {
1791     assume(kFindDivisibleByInS_easy(c->strat,buf[k])<0);
1792     int k2;
1793     for(k2=0;k2<i;k2++)
1794     {
1795       if(k==k2) continue;
1796       assume((!(p_LmDivisibleBy(buf[k].p,buf[k2].p,c->r)))||(wrp(buf[k].p),Print(" k %d k2 %d ",k,k2),wrp(buf[k2].p),FALSE));
1797     }
1798   }
1799 }
1800#endif
1801  //resort S
1802#ifdef FIND_DETERMINISTIC
1803  for(z2=0;z2<=c->strat->sl;z2++)
1804  {
1805    if (c->modifiedS[z2])
1806    {
1807      wlen_type qual;
1808      int new_pos;
1809      if (c->strat->lenSw!=NULL)
1810          new_pos=simple_posInS(c->strat,c->strat->S[z2],strat->lenS[z2],strat->Sw[z2]);
1811      else
1812          new_pos=simple_posInS(c->strat,c->strat->S[z2],strat->lenS[z2],lenS[z2]);
1813
1814      if (new_pos<z2)
1815      {
1816         move_forward_in_S(z2,new_pos,c->strat);
1817      }
1818
1819      assume(new_pos<=z2);
1820    }
1821  }
1822  for(z2=0;c->expandS[z2]!=NULL;z2++)
1823  {
1824    add_to_reductors(c,c->expandS[z2],pLength(c->expandS[z2]));
1825    // PrintS("E");
1826  }
1827  omfree(c->modifiedS);
1828  c->modifiedS=NULL;
1829  omfree(c->expandS);
1830  c->expandS=NULL;
1831#endif
1832  if (TEST_OPT_PROT)
1833      Print("%i]",i);
1834
1835  int* ibuf=(int*) omalloc(i*sizeof(int));
1836  sorted_pair_node*** sbuf=(sorted_pair_node***) omalloc(i*sizeof(sorted_pair_node**));
1837  for(j=0;j<i;j++)
1838  {
1839    int len;
1840    poly p;
1841    buf[j].flatten();
1842    kBucketClear(buf[j].bucket,&p, &len);
1843    kBucketDestroy(&buf[j].bucket);
1844
1845    if (!c->nc) {
1846      if ((TEST_OPT_REDTAIL)&&(c->S->rank<=1)){
1847      p=redNFTail(p,c->strat->sl,c->strat, 0);
1848      } else {
1849      p=redTailShort(p, c->strat);
1850      }
1851      }
1852    sbuf[j]=add_to_basis_ideal_quotient(p,c,ibuf+j);
1853    //sbuf[j]=add_to_basis(p,-1,-1,c,ibuf+j);
1854  }
1855  int sum=0;
1856  for(j=0;j<i;j++){
1857    sum+=ibuf[j];
1858  }
1859  sorted_pair_node** big_sbuf=(sorted_pair_node**) omalloc(sum*sizeof(sorted_pair_node*));
1860  int partsum=0;
1861  for(j=0;j<i;j++)
1862  {
1863    memmove(big_sbuf+partsum, sbuf[j],ibuf[j]*sizeof(sorted_pair_node*));
1864    omfree(sbuf[j]);
1865    partsum+=ibuf[j];
1866  }
1867
1868  qsort(big_sbuf,sum,sizeof(sorted_pair_node*),tgb_pair_better_gen2);
1869  c->apairs=spn_merge(c->apairs,c->pair_top+1,big_sbuf,sum,c);
1870  c->pair_top+=sum;
1871  clean_top_of_pair_list(c);
1872  omfree(big_sbuf);
1873  omfree(sbuf);
1874  omfree(ibuf);
1875  omfree(buf);
1876#ifdef TGB_DEBUG
1877  int z;
1878  for(z=1;z<=c->pair_top;z++)
1879  {
1880    assume(pair_better(c->apairs[z],c->apairs[z-1],c));
1881  }
1882#endif
1883  if (TEST_OPT_PROT)
1884      Print("(%d)",c->pair_top+1);
1885  while(!(idIs0(c->add_later))){
1886    ideal add=c->add_later;
1887    idSkipZeroes(add);
1888    for(j=0;j<add->idelems();j++){
1889        assume(pLength(add->m[j])==1);
1890        p_SetCoeff(add->m[j],n_Init(1,currRing),currRing);
1891
1892    }
1893    ideal add2=kInterRed(add,NULL);
1894    id_Delete(&add,currRing);
1895    idSkipZeroes(add2);
1896    c->add_later=idInit(ADD_LATER_SIZE,c->S->rank);
1897    memset(c->add_later->m,0,ADD_LATER_SIZE*sizeof(poly));
1898//     for(i=0;i<add2->idelems();i++){
1899//       if (add2->m[i]!=NULL)
1900//           add_to_basis_ideal_quotient(add2->m[i],-1,-1,c,NULL);
1901//       add2->m[i]=NULL;
1902//     }
1903    int i=add2->idelems();
1904    int* ibuf=(int*) omalloc(i*sizeof(int));
1905  sorted_pair_node*** sbuf=(sorted_pair_node***) omalloc(i*sizeof(sorted_pair_node**));
1906
1907  for(j=0;j<i;j++)
1908  {
1909    int len;
1910    poly p;
1911    //buf[j].flatten();
1912    //kBucketClear(buf[j].bucket,&p, &len);
1913    //kBucketDestroy(&buf[j].bucket);
1914    p=add2->m[j];
1915    add2->m[j]=NULL;
1916
1917    sbuf[j]=add_to_basis_ideal_quotient(p,c,ibuf+j);
1918    //sbuf[j]=add_to_basis(p,-1,-1,c,ibuf+j);
1919  }
1920  int sum=0;
1921  for(j=0;j<i;j++){
1922    sum+=ibuf[j];
1923  }
1924  sorted_pair_node** big_sbuf=(sorted_pair_node**) omalloc(sum*sizeof(sorted_pair_node*));
1925  int partsum=0;
1926  for(j=0;j<i;j++)
1927  {
1928    memmove(big_sbuf+partsum, sbuf[j],ibuf[j]*sizeof(sorted_pair_node*));
1929    omfree(sbuf[j]);
1930    partsum+=ibuf[j];
1931  }
1932
1933  qsort(big_sbuf,sum,sizeof(sorted_pair_node*),tgb_pair_better_gen2);
1934  c->apairs=spn_merge(c->apairs,c->pair_top+1,big_sbuf,sum,c);
1935  c->pair_top+=sum;
1936  clean_top_of_pair_list(c);
1937  omfree(big_sbuf);
1938  omfree(sbuf);
1939  omfree(ibuf);
1940  //omfree(buf);
1941    id_Delete(&add2, c->r);
1942  }
1943  #ifdef TGB_RESORT_PAIRS
1944  delete c->replaced;
1945  c->replaced=NULL;
1946  c->used_b=FALSE;
1947  #endif
1948  return;
1949}
1950
1951
1952
1953#ifdef REDTAIL_S
1954
1955static poly redNFTail (poly h,const int sl,kStrategy strat, int len)
1956{
1957  if (h==NULL) return NULL;
1958  pTest(h);
1959  if (0 > sl)
1960    return h;
1961  if (pNext(h)==NULL) return h;
1962
1963  int j;
1964  poly res=h;
1965  poly act=res;
1966  LObject P(pNext(h));
1967  pNext(res)=NULL;
1968  P.bucket = kBucketCreate(currRing);
1969  len--;
1970  h=P.p;
1971  if (len <=0) len=pLength(h);
1972  kBucketInit(P.bucket,h /*P.p*/,len /*pLength(P.p)*/);
1973  pTest(h);
1974  loop
1975  {
1976      P.p=h;
1977      P.t_p=NULL;
1978      P.SetShortExpVector();
1979      loop
1980      {
1981          j=kFindDivisibleByInS(strat->S,strat->sevS,sl,&P);
1982          if (j>=0)
1983          {
1984#ifdef REDTAIL_PROT
1985            PrintS("r");
1986#endif
1987            nNormalize(pGetCoeff(P.p));
1988#ifdef KDEBUG
1989            if (TEST_OPT_DEBUG)
1990            {
1991              PrintS("red tail:");
1992              wrp(h);
1993              PrintS(" with ");
1994              wrp(strat->S[j]);
1995            }
1996#endif
1997            number coef;
1998            pTest(strat->S[j]);
1999            coef=kBucketPolyRed(P.bucket,strat->S[j],
2000                                strat->lenS[j]/*pLength(strat->S[j])*/,strat->kNoether);
2001            pMult_nn(res,coef);
2002            nDelete(&coef);
2003            h = kBucketGetLm(P.bucket);
2004            pTest(h);
2005            if (h==NULL)
2006            {
2007#ifdef REDTAIL_PROT
2008              PrintS(" ");
2009#endif
2010        kBucketDestroy(&P.bucket);
2011              return res;
2012            }
2013            pTest(h);
2014            P.p=h;
2015            P.t_p=NULL;
2016            P.SetShortExpVector();
2017#ifdef KDEBUG
2018            if (TEST_OPT_DEBUG)
2019            {
2020              PrintS("\nto tail:");
2021              wrp(h);
2022              PrintLn();
2023            }
2024#endif
2025          }
2026          else
2027          {
2028#ifdef REDTAIL_PROT
2029            PrintS("n");
2030#endif
2031            break;
2032          }
2033      } /* end loop current mon */
2034      //   poly tmp=pHead(h /*kBucketGetLm(P.bucket)*/);
2035      //act->next=tmp;pIter(act);
2036      act->next=kBucketExtractLm(P.bucket);pIter(act);
2037      h = kBucketGetLm(P.bucket);
2038      if (h==NULL)
2039      {
2040#ifdef REDTAIL_PROT
2041        PrintS(" ");
2042#endif
2043  kBucketDestroy(&P.bucket);
2044        return res;
2045      }
2046      pTest(h);
2047  }
2048}
2049#endif
2050
2051
2052//try to fill, return FALSE iff queue is empty
2053
2054//transfers ownership of m to mat
2055void init_with_mac_poly(tgb_sparse_matrix* mat, int row, mac_poly m){
2056  assume(mat->mp[row]==NULL);
2057  mat->mp[row]=m;
2058#ifdef TGB_DEBUG
2059  mac_poly r=m;
2060  while(r){
2061    assume(r->exp<mat->columns);
2062    r=r->next;
2063  }
2064#endif
2065}
2066poly free_row_to_poly(tgb_sparse_matrix* mat, int row, poly* monoms, int monom_index){
2067  poly p=NULL;
2068  poly* set_this=&p;
2069  mac_poly r=mat->mp[row];
2070  mat->mp[row]=NULL;
2071  while(r)
2072  {
2073    (*set_this)=pLmInit(monoms[monom_index-1-r->exp]);
2074    pSetCoeff((*set_this),r->coef);
2075    set_this=&((*set_this)->next);
2076    mac_poly old=r;
2077    r=r->next;
2078    delete old;
2079
2080  }
2081  return p;
2082
2083}
2084
2085
2086
2087
2088
2089
2090
2091
2092
2093
2094
2095static int poly_crit(const void* ap1, const void* ap2){
2096  poly p1,p2;
2097  p1=*((poly*) ap1);
2098  p2=*((poly*)ap2);
2099
2100  int c=pLmCmp(p1,p2);
2101  if (c !=0) return c;
2102  int l1=pLength(p1);
2103  int l2=pLength(p2);
2104  if (l1<l2) return -1;
2105  if (l1>l2) return 1;
2106  return 0;
2107}
2108slimgb_alg::slimgb_alg(ideal I, int syz_comp,BOOLEAN F4){
2109  completed=FALSE;
2110  this->syz_comp=syz_comp;
2111  r=currRing;
2112  nc=rIsPluralRing(r);
2113
2114  is_homog=TRUE;
2115  {
2116    int hz;
2117    for(hz=0;hz<IDELEMS(I);hz++){
2118      assume(I->m[hz]!=NULL);
2119      int d=pTotaldegree(I->m[hz]);
2120      poly t=I->m[hz]->next;
2121      while(t)
2122      {
2123        if (d!=pTotaldegree(t,r))
2124        {
2125          is_homog=FALSE;
2126          break;
2127        }
2128        t=t->next;
2129      }
2130      if(!(is_homog)) break;
2131    }
2132  }
2133  eliminationProblem=((!(is_homog))&&((pLexOrder)));
2134  //  Print("is homog:%d",c->is_homog);
2135  void* h;
2136  poly hp;
2137  int i,j;
2138  to_destroy=NULL;
2139  easy_product_crit=0;
2140  extended_product_crit=0;
2141  if (rField_is_Zp(r))
2142    isDifficultField=FALSE;
2143  else
2144    isDifficultField=TRUE;
2145  //not fully correct
2146  //(rChar()==0);
2147  F4_mode=F4;
2148
2149  reduction_steps=0;
2150  last_index=-1;
2151
2152  F=NULL;
2153  F_minus=NULL;
2154
2155  Rcounter=0;
2156
2157  soon_free=NULL;
2158
2159  tmp_lm=pOne();
2160
2161  normal_forms=0;
2162  current_degree=1;
2163
2164  max_pairs=5*I->idelems();
2165
2166  apairs=(sorted_pair_node**) omalloc(sizeof(sorted_pair_node*)*max_pairs);
2167  pair_top=-1;
2168
2169  int n=I->idelems();
2170  array_lengths=n;
2171
2172
2173  i=0;
2174  this->n=0;
2175  T_deg=(int*) omalloc(n*sizeof(int));
2176  if(eliminationProblem)
2177    T_deg_full=(int*) omalloc(n*sizeof(int));
2178  else
2179    T_deg_full=NULL;
2180  tmp_pair_lm=(poly*) omalloc(n*sizeof(poly));
2181  tmp_spn=(sorted_pair_node**) omalloc(n*sizeof(sorted_pair_node*));
2182  lm_bin=omGetSpecBin(POLYSIZE + (r->ExpL_Size)*sizeof(long));
2183#ifdef HEAD_BIN
2184  HeadBin=omGetSpecBin(POLYSIZE + (currRing->ExpL_Size)*sizeof(long));
2185#endif
2186  /* omUnGetSpecBin(&(c->HeadBin)); */
2187  #ifndef HAVE_BOOST
2188  #ifdef USE_STDVECBOOL
2189  #else
2190  h=omalloc(n*sizeof(char*));
2191
2192  states=(char**) h;
2193  #endif
2194  #endif
2195  h=omalloc(n*sizeof(int));
2196  lengths=(int*) h;
2197  weighted_lengths=(wlen_type*)omalloc(n*sizeof(wlen_type));
2198  gcd_of_terms=(poly*) omalloc(n*sizeof(poly));
2199
2200  short_Exps=(long*) omalloc(n*sizeof(long));
2201  if (F4_mode)
2202    S=idInit(n,I->rank);
2203  else
2204    S=idInit(1,I->rank);
2205  strat=new skStrategy;
2206  strat->syzComp = 0;
2207  initBuchMoraCrit(strat);
2208  initBuchMoraPos(strat);
2209  strat->initEcart = initEcartBBA;
2210  strat->tailRing=r;
2211  strat->enterS = enterSBba;
2212  strat->sl = -1;
2213  i=n;
2214  i=1;//some strange bug else
2215  /* initS(c->S,NULL,c->strat); */
2216  /* intS start: */
2217  // i=((i+IDELEMS(c->S)+15)/16)*16;
2218  strat->ecartS=(intset)omAlloc(i*sizeof(int)); /*initec(i);*/
2219  strat->sevS=(unsigned long*)omAlloc0(i*sizeof(unsigned long));
2220  /*initsevS(i);*/
2221  strat->S_2_R=(int*)omAlloc0(i*sizeof(int));/*initS_2_R(i);*/
2222  strat->fromQ=NULL;
2223  strat->Shdl=idInit(1,1);
2224  strat->S=strat->Shdl->m;
2225  strat->lenS=(int*)omAlloc0(i*sizeof(int));
2226  if((isDifficultField)||(eliminationProblem))
2227    strat->lenSw=(wlen_type*)omAlloc0(i*sizeof(wlen_type));
2228  else
2229    strat->lenSw=NULL;
2230  sorted_pair_node* si;
2231  assume(n>0);
2232  add_to_basis_ideal_quotient(I->m[0],this,NULL);
2233
2234  assume(strat->sl==strat->Shdl->idelems()-1);
2235  if(!(F4_mode))
2236  {
2237    for (i=1;i<n;i++)//the 1 is wanted, because first element is added to basis
2238    {
2239      //     add_to_basis(I->m[i],-1,-1,c);
2240      si=(sorted_pair_node*) omalloc(sizeof(sorted_pair_node));
2241      si->i=-1;
2242      si->j=-2;
2243      si->expected_length=pQuality(I->m[i],this,pLength(I->m[i]));
2244      si->deg=pTotaldegree(I->m[i]);
2245      if (!rField_is_Zp(r)){
2246        pCleardenom(I->m[i]);
2247      }
2248      si->lcm_of_lm=I->m[i];
2249
2250      //      c->apairs[n-1-i]=si;
2251      apairs[n-i-1]=si;
2252      ++(pair_top);
2253    }
2254  }
2255  else
2256  {
2257    for (i=1;i<n;i++)//the 1 is wanted, because first element is added to basis
2258      add_to_basis_ideal_quotient(I->m[i],this,NULL);
2259  }
2260  for(i=0;i<I->idelems();i++)
2261  {
2262    I->m[i]=NULL;
2263  }
2264  idDelete(&I);
2265  add_later=idInit(ADD_LATER_SIZE,S->rank);
2266  memset(add_later->m,0,ADD_LATER_SIZE*sizeof(poly));
2267}
2268slimgb_alg::~slimgb_alg(){
2269 
2270  if (!(completed)){
2271      poly* add=(poly*) omalloc((pair_top+2)*sizeof(poly));
2272      int piter;
2273      int pos=0;
2274      for(piter=0;piter<=pair_top;piter++){
2275        sorted_pair_node* s=apairs[piter];
2276        if (s->i<0){
2277            //delayed element
2278            if (s->lcm_of_lm!=NULL){
2279                add[pos]=s->lcm_of_lm;
2280                pos++;
2281
2282            }
2283        }
2284   
2285        free_sorted_pair_node(s,r);
2286        apairs[piter]=NULL;
2287      }
2288      pair_top=-1;
2289      add[pos]=NULL;
2290      pos=0;
2291      while(add[pos]!=NULL){
2292        add_to_basis_ideal_quotient(add[pos],this,NULL);
2293        pos++;
2294      }
2295      for(piter=0;piter<=pair_top;piter++){
2296        sorted_pair_node* s=apairs[piter];
2297        assume(s->i>=0);
2298        free_sorted_pair_node(s,r);
2299        apairs[piter]=NULL;
2300      }
2301      pair_top=-1;
2302  }
2303  id_Delete(&add_later,r);
2304  int i,j;
2305  slimgb_alg* c=this;
2306  while(c->to_destroy)
2307  {
2308    pDelete(&(c->to_destroy->p));
2309    poly_list_node* old=c->to_destroy;
2310    c->to_destroy=c->to_destroy->next;
2311    omfree(old);
2312  }
2313  while(c->F)
2314  {
2315    for(i=0;i<c->F->size;i++){
2316      pDelete(&(c->F->mp[i].m));
2317    }
2318    omfree(c->F->mp);
2319    c->F->mp=NULL;
2320    mp_array_list* old=c->F;
2321    c->F=c->F->next;
2322    omfree(old);
2323  }
2324  while(c->F_minus)
2325  {
2326    for(i=0;i<c->F_minus->size;i++){
2327      pDelete(&(c->F_minus->p[i]));
2328    }
2329    omfree(c->F_minus->p);
2330    c->F_minus->p=NULL;
2331    poly_array_list* old=c->F_minus;
2332    c->F_minus=c->F_minus->next;
2333    omfree(old);
2334  }
2335  #ifndef HAVE_BOOST
2336  #ifndef USE_STDVECBOOL
2337  for(int z=1 /* zero length at 0 */;z<c->n;z++){
2338    omfree(c->states[z]);
2339  }
2340  omfree(c->states);
2341  #endif
2342  #endif
2343
2344  omfree(c->lengths);
2345  omfree(c->weighted_lengths);
2346  for(int z=0;z<c->n;z++)
2347  {
2348    pDelete(&c->tmp_pair_lm[z]);
2349    omfree(c->tmp_spn[z]);
2350  }
2351  omfree(c->tmp_pair_lm);
2352  omfree(c->tmp_spn);
2353
2354  omfree(c->T_deg);
2355  if(c->T_deg_full)
2356    omfree(c->T_deg_full);
2357
2358  omFree(c->strat->ecartS);
2359  omFree(c->strat->sevS);
2360//   initsevS(i);
2361  omFree(c->strat->S_2_R);
2362
2363
2364  omFree(c->strat->lenS);
2365
2366  if(c->strat->lenSw)  omFree(c->strat->lenSw);
2367
2368
2369
2370
2371  for(i=0;i<c->n;i++){
2372    if(c->gcd_of_terms[i])
2373      pDelete(&(c->gcd_of_terms[i]));
2374  }
2375  omfree(c->gcd_of_terms);
2376
2377  omfree(c->apairs);
2378  if (TEST_OPT_PROT)
2379  {
2380    Print("calculated %d NFs\n",c->normal_forms);
2381    Print("applied %i product crit, %i extended_product crit \n", c->easy_product_crit, c->extended_product_crit);
2382  }
2383  int deleted_form_c_s=0;
2384
2385  for(i=0;i<=c->strat->sl;i++){
2386    if (!c->strat->S[i]) continue;
2387    BOOLEAN found=FALSE;
2388    for(j=0;j<c->n;j++){
2389      if (c->S->m[j]==c->strat->S[i]){
2390        found=TRUE;
2391        break;
2392      }
2393    }
2394    if(!found) pDelete(&c->strat->S[i]);
2395  }
2396//   for(i=0;i<c->n;i++){
2397//     if (c->rep[i]!=i){
2398// //       for(j=0;j<=c->strat->sl;j++){
2399// //   if(c->strat->S[j]==c->S->m[i]){
2400// //     c->strat->S[j]=NULL;
2401// //     break;
2402// //   }
2403// //       }
2404// //      PrintS("R_delete");
2405//       pDelete(&c->S->m[i]);
2406//     }
2407//   }
2408
2409  if (completed){
2410  for(i=0;i<c->n;i++)
2411  {
2412    assume(c->S->m[i]!=NULL);
2413    if (p_GetComp(c->S->m[i],currRing)>this->syz_comp) continue;
2414    for(j=0;j<c->n;j++)
2415    {
2416      if((c->S->m[j]==NULL)||(i==j))
2417        continue;
2418      assume(p_LmShortDivisibleBy(c->S->m[j],c->short_Exps[j],
2419             c->S->m[i],~c->short_Exps[i],
2420             c->r)==p_LmDivisibleBy(c->S->m[j],
2421             c->S->m[i],
2422             c->r));
2423      if (p_LmShortDivisibleBy(c->S->m[j],c->short_Exps[j],
2424          c->S->m[i],~c->short_Exps[i],
2425          c->r))
2426      {
2427        pDelete(&c->S->m[i]);
2428        break;
2429      }
2430    }
2431  }
2432  }
2433  omfree(c->short_Exps);
2434
2435
2436  ideal I=c->S;
2437
2438  IDELEMS(I)=c->n;
2439
2440  idSkipZeroes(I);
2441  for(i=0;i<=c->strat->sl;i++)
2442    c->strat->S[i]=NULL;
2443  id_Delete(&c->strat->Shdl,c->r);
2444  pDelete(&c->tmp_lm);
2445  omUnGetSpecBin(&lm_bin);
2446  delete c->strat;
2447}
2448ideal t_rep_gb(ring r,ideal arg_I, int syz_comp, BOOLEAN F4_mode){
2449
2450  //  Print("QlogSize(0) %d, QlogSize(1) %d,QlogSize(-2) %d, QlogSize(5) %d\n", QlogSize(nlInit(0)),QlogSize(nlInit(1)),QlogSize(nlInit(-2)),QlogSize(nlInit(5)));
2451 
2452  if (TEST_OPT_PROT)
2453    if (F4_mode)
2454      PrintS("F4 Modus \n");
2455
2456  //debug_Ideal=arg_debug_Ideal;
2457  //if (debug_Ideal) PrintS("DebugIdeal received\n");
2458  // Print("Idelems %i \n----------\n",IDELEMS(arg_I));
2459  ideal I=idCompactify(arg_I);
2460   int i;
2461  if (idIs0(I)) return I;
2462  for(i=0;i<IDELEMS(I);i++)
2463  {
2464    assume(I->m[i]!=NULL);
2465    simplify_poly(I->m[i],currRing);
2466  }
2467 
2468
2469  qsort(I->m,IDELEMS(I),sizeof(poly),poly_crit);
2470  //Print("Idelems %i \n----------\n",IDELEMS(I));
2471  //slimgb_alg* c=(slimgb_alg*) omalloc(sizeof(slimgb_alg));
2472  //int syz_comp=arg_I->rank;
2473  slimgb_alg* c=new slimgb_alg(I, syz_comp,F4_mode);
2474
2475
2476  while ((c->pair_top>=0) && ((!(TEST_OPT_DEGBOUND)) || (c->apairs[c->pair_top]->deg<=Kstd1_deg)))
2477  {
2478    if(F4_mode)
2479      go_on_F4(c);
2480    else
2481      go_on(c);
2482  }
2483  if (c->pair_top<0)
2484    c->completed=TRUE;
2485  I=c->S;
2486  delete c;
2487  if (TEST_OPT_REDSB){
2488    ideal erg=kInterRed(I,NULL);
2489    assume(I!=erg);
2490    id_Delete(&I, currRing);
2491    return erg;
2492  }
2493  //qsort(I->m, IDELEMS(I),sizeof(poly),pLmCmp_func);
2494  assume(I->rank>=idRankFreeModule(I));
2495  return(I);
2496}
2497void now_t_rep(const int & arg_i, const int & arg_j, slimgb_alg* c){
2498  int i,j;
2499  if (arg_i==arg_j){
2500    return;
2501  }
2502  if (arg_i>arg_j){
2503    i=arg_j;
2504    j=arg_i;
2505  } else {
2506    i=arg_i;
2507    j=arg_j;
2508  }
2509  c->states[j][i]=HASTREP;
2510}
2511
2512static BOOLEAN has_t_rep(const int & arg_i, const  int & arg_j, slimgb_alg* state){
2513  assume(0<=arg_i);
2514  assume(0<=arg_j);
2515  assume(arg_i<state->n);
2516  assume(arg_j<state->n);
2517  if (arg_i==arg_j)
2518  {
2519    return (TRUE);
2520  }
2521  if (arg_i>arg_j)
2522  {
2523    return (state->states[arg_i][arg_j]==HASTREP);
2524  } else
2525  {
2526    return (state->states[arg_j][arg_i]==HASTREP);
2527  }
2528}
2529static int pLcmDeg(poly a, poly b)
2530{
2531  int i;
2532  int n=0;
2533  for (i=pVariables; i; i--)
2534  {
2535    n+=si_max( pGetExp(a,i), pGetExp(b,i));
2536  }
2537  return n;
2538
2539}
2540
2541
2542
2543static void shorten_tails(slimgb_alg* c, poly monom)
2544{
2545  return;
2546// BOOLEAN corr=lenS_correct(c->strat);
2547  for(int i=0;i<c->n;i++)
2548  {
2549    //enter tail
2550
2551    if (c->S->m[i]==NULL) continue;
2552    poly tail=c->S->m[i]->next;
2553    poly prev=c->S->m[i];
2554    BOOLEAN did_something=FALSE;
2555    while((tail!=NULL)&& (pLmCmp(tail, monom)>=0))
2556    {
2557      if (p_LmDivisibleBy(monom,tail,c->r))
2558      {
2559        did_something=TRUE;
2560        prev->next=tail->next;
2561        tail->next=NULL;
2562        p_Delete(& tail,c->r);
2563        tail=prev;
2564        //PrintS("Shortened");
2565        c->lengths[i]--;
2566      }
2567      prev=tail;
2568      tail=tail->next;
2569    }
2570    if (did_something)
2571    {
2572      int new_pos;
2573      wlen_type q;
2574      q=pQuality(c->S->m[i],c,c->lengths[i]);
2575      new_pos=simple_posInS(c->strat,c->S->m[i],c->lengths[i],q);
2576
2577      int old_pos=-1;
2578      //assume new_pos<old_pos
2579      for (int z=0;z<=c->strat->sl;z++)
2580      {
2581        if (c->strat->S[z]==c->S->m[i])
2582        {
2583          old_pos=z;
2584          break;
2585        }
2586      }
2587      if (old_pos== -1)
2588        for (int z=new_pos-1;z>=0;z--)
2589        {
2590          if (c->strat->S[z]==c->S->m[i])
2591          {
2592            old_pos=z;
2593            break;
2594          }
2595        }
2596      assume(old_pos>=0);
2597      assume(new_pos<=old_pos);
2598      assume(pLength(c->strat->S[old_pos])==c->lengths[i]);
2599      c->strat->lenS[old_pos]=c->lengths[i];
2600      if (c->strat->lenSw)
2601        c->strat->lenSw[old_pos]=q;
2602
2603      if (new_pos<old_pos)
2604        move_forward_in_S(old_pos,new_pos,c->strat);
2605
2606      length_one_crit(c,i,c->lengths[i]);
2607    }
2608  }
2609}
2610static sorted_pair_node* pop_pair(slimgb_alg* c){
2611  clean_top_of_pair_list(c);
2612
2613  if(c->pair_top<0) return NULL;
2614  else return (c->apairs[c->pair_top--]);
2615}
2616sorted_pair_node* top_pair(slimgb_alg* c){
2617  super_clean_top_of_pair_list(c);//yeah, I know, it's odd that I use a different proc here
2618
2619  if(c->pair_top<0) return NULL;
2620  else return (c->apairs[c->pair_top]);
2621}
2622sorted_pair_node* quick_pop_pair(slimgb_alg* c){
2623  if(c->pair_top<0) return NULL;
2624  else return (c->apairs[c->pair_top--]);
2625}
2626
2627
2628
2629static void super_clean_top_of_pair_list(slimgb_alg* c){
2630  while((c->pair_top>=0)
2631  && (c->apairs[c->pair_top]->i>=0)
2632  && (good_has_t_rep(c->apairs[c->pair_top]->j, c->apairs[c->pair_top]->i,c)))
2633  {
2634
2635    free_sorted_pair_node(c->apairs[c->pair_top],c->r);
2636    c->pair_top--;
2637
2638  }
2639}
2640void clean_top_of_pair_list(slimgb_alg* c){
2641  while((c->pair_top>=0) && (c->apairs[c->pair_top]->i>=0) && (!state_is(UNCALCULATED,c->apairs[c->pair_top]->j, c->apairs[c->pair_top]->i,c))){
2642
2643    free_sorted_pair_node(c->apairs[c->pair_top],c->r);
2644    c->pair_top--;
2645
2646  }
2647}
2648static BOOLEAN state_is(calc_state state, const int & arg_i, const  int & arg_j, slimgb_alg* c){
2649  assume(0<=arg_i);
2650  assume(0<=arg_j);
2651  assume(arg_i<c->n);
2652  assume(arg_j<c->n);
2653  if (arg_i==arg_j)
2654  {
2655    return (TRUE);
2656  }
2657  if (arg_i>arg_j)
2658  {
2659    return (c->states[arg_i][arg_j]==state);
2660  }
2661  else return(c->states[arg_j][arg_i]==state);
2662}
2663
2664
2665void free_sorted_pair_node(sorted_pair_node* s, ring r){
2666  if (s->i>=0)
2667    p_Delete(&s->lcm_of_lm,r);
2668  omfree(s);
2669}
2670static BOOLEAN pair_better(sorted_pair_node* a,sorted_pair_node* b, slimgb_alg* c){
2671  if (a->deg<b->deg) return TRUE;
2672  if (a->deg>b->deg) return FALSE;
2673
2674
2675  int comp=pLmCmp(a->lcm_of_lm, b->lcm_of_lm);
2676  if (comp==1) return FALSE;
2677  if (-1==comp) return TRUE;
2678  if (a->expected_length<b->expected_length) return TRUE;
2679  if (a->expected_length>b->expected_length) return FALSE;
2680  if (a->i+a->j<b->i+b->j) return TRUE;
2681   if (a->i+a->j>b->i+b->j) return FALSE;
2682  if (a->i<b->i) return TRUE;
2683  if (a->i>b->i) return FALSE;
2684  return TRUE;
2685}
2686
2687static int tgb_pair_better_gen(const void* ap,const void* bp){
2688
2689  sorted_pair_node* a=*((sorted_pair_node**)ap);
2690  sorted_pair_node* b=*((sorted_pair_node**)bp);
2691  assume((a->i>a->j) || (a->i < 0));
2692  assume((b->i>b->j) || (b->i < 0));
2693  if (a->deg<b->deg) return -1;
2694  if (a->deg>b->deg) return 1;
2695
2696
2697 
2698 int comp=pLmCmp(a->lcm_of_lm, b->lcm_of_lm);
2699
2700  if (comp==1) return 1;
2701  if (-1==comp) return -1;
2702   if (a->expected_length<b->expected_length) return -1;
2703  if (a->expected_length>b->expected_length) return 1;
2704  if (a->i+a->j<b->i+b->j) return -1;
2705   if (a->i+a->j>b->i+b->j) return 1;
2706  if (a->i<b->i) return -1;
2707   if (a->i>b->i) return 1;
2708  return 0;
2709}
2710
2711
2712static poly gcd_of_terms(poly p, ring r){
2713  int max_g_0=0;
2714  assume(p!=NULL);
2715  int i;
2716  poly m=pOne();
2717  poly t;
2718  for (i=pVariables; i; i--)
2719  {
2720      pSetExp(m,i, pGetExp(p,i));
2721      if (max_g_0==0)
2722  if (pGetExp(m,i)>0)
2723    max_g_0=i;
2724  }
2725
2726  t=p->next;
2727  while (t!=NULL){
2728
2729    if (max_g_0==0) break;
2730    for (i=max_g_0; i; i--)
2731    {
2732      pSetExp(m,i, si_min(pGetExp(t,i),pGetExp(m,i)));
2733      if (max_g_0==i)
2734  if (pGetExp(m,i)==0)
2735    max_g_0=0;
2736      if ((max_g_0==0) && (pGetExp(m,i)>0)){
2737  max_g_0=i;
2738      }
2739    }
2740    t=t->next;
2741  }
2742
2743  if (max_g_0>0)
2744    return m;
2745  pDelete(&m);
2746  return NULL;
2747}
2748static inline BOOLEAN pHasNotCFExtended(poly p1, poly p2, poly m)
2749{
2750
2751  if (pGetComp(p1) > 0 || pGetComp(p2) > 0)
2752    return FALSE;
2753  int i = 1;
2754  loop
2755  {
2756    if ((pGetExp(p1, i)-pGetExp(m,i) >0) && (pGetExp(p2, i) -pGetExp(m,i)> 0))   return FALSE;
2757    if (i == pVariables)                                return TRUE;
2758    i++;
2759  }
2760}
2761
2762
2763//for impl reasons may return false if the the normal product criterion matches
2764static inline BOOLEAN extended_product_criterion(poly p1, poly gcd1, poly p2, poly gcd2, slimgb_alg* c){
2765  if (c->nc)
2766    return FALSE;
2767  if(gcd1==NULL) return FALSE;
2768        if(gcd2==NULL) return FALSE;
2769        gcd1->next=gcd2; //may ordered incorrect
2770        poly m=gcd_of_terms(gcd1,c->r);
2771        gcd1->next=NULL;
2772        if (m==NULL) return FALSE;
2773
2774        BOOLEAN erg=pHasNotCFExtended(p1,p2,m);
2775        pDelete(&m);
2776        return erg;
2777}
2778static poly kBucketGcd(kBucket* b, ring r)
2779{
2780  int s=0;
2781  int i;
2782  poly m, n;
2783  BOOLEAN initialized=FALSE;
2784  for (i=MAX_BUCKET-1;i>=0;i--)
2785  {
2786    if (b->buckets[i]!=NULL){
2787      if (!initialized){
2788  m=gcd_of_terms(b->buckets[i],r);
2789  initialized=TRUE;
2790  if (m==NULL) return NULL;
2791      }
2792      else
2793  {
2794    n=gcd_of_terms(b->buckets[i],r);
2795    if (n==NULL) {
2796      pDelete(&m);
2797      return NULL;
2798    }
2799    n->next=m;
2800    poly t=gcd_of_terms(n,r);
2801    n->next=NULL;
2802    pDelete(&m);
2803    pDelete(&n);
2804    m=t;
2805    if (m==NULL) return NULL;
2806
2807  }
2808    }
2809  }
2810  return m;
2811}
2812
2813
2814
2815
2816static inline int quality_of_pos_in_strat_S(int pos, slimgb_alg* c){
2817  if (c->strat->lenSw!=NULL) return c->strat->lenSw[pos];
2818  return c->strat->lenS[pos];
2819}
2820static inline int quality_of_pos_in_strat_S_mult_high(int pos, poly high, slimgb_alg* c)
2821  //meant only for nc
2822{
2823  poly m=pOne();
2824  pExpVectorDiff(m,high ,c->strat->S[pos]);
2825  poly product=nc_mm_Mult_p(m, pCopy(c->strat->S[pos]), c->r);
2826  int erg=pQuality(product,c);
2827  pDelete(&m);
2828  pDelete(&product);
2829  return erg;
2830}
2831
2832static void multi_reduction_lls_trick(red_object* los, int losl,slimgb_alg* c,find_erg & erg){
2833  erg.expand=NULL;
2834  BOOLEAN swap_roles; //from reduce_by, to_reduce_u if fromS
2835  if(erg.fromS){
2836    if(pLmEqual(c->strat->S[erg.reduce_by],los[erg.to_reduce_u].p))
2837    {
2838      int i;
2839      int quality_a=quality_of_pos_in_strat_S(erg.reduce_by,c);
2840      int best=erg.to_reduce_u+1;
2841/*
2842      for (i=erg.to_reduce_u;i>=erg.to_reduce_l;i--){
2843  int qc=los[i].guess_quality(c);
2844  if (qc<quality_a){
2845    best=i;
2846    quality_a=qc;
2847  }
2848      }
2849      if(best!=erg.to_reduce_u+1){*/
2850      int qc;
2851      best=find_best(los,erg.to_reduce_l,erg.to_reduce_u,qc,c);
2852      if(qc<quality_a){
2853  los[best].flatten();
2854  int b_pos=kBucketCanonicalize(los[best].bucket);
2855  los[best].p=los[best].bucket->buckets[b_pos];
2856  qc=pQuality(los[best].bucket->buckets[b_pos],c);
2857  if(qc<quality_a){
2858    red_object h=los[erg.to_reduce_u];
2859    los[erg.to_reduce_u]=los[best];
2860    los[best]=h;
2861    swap_roles=TRUE;
2862  }
2863  else
2864    swap_roles=FALSE;
2865      }
2866      else{
2867
2868  swap_roles=FALSE;
2869      }
2870
2871    }
2872      else
2873    {
2874      if (erg.to_reduce_u>erg.to_reduce_l){
2875
2876  int i;
2877  int quality_a=quality_of_pos_in_strat_S(erg.reduce_by,c);
2878  if (c->nc)
2879    quality_a=quality_of_pos_in_strat_S_mult_high(erg.reduce_by, los[erg.to_reduce_u].p, c);
2880  int best=erg.to_reduce_u+1;
2881  int qc;
2882  best=find_best(los,erg.to_reduce_l,erg.to_reduce_u,qc,c);
2883  assume(qc==los[best].guess_quality(c));
2884  if(qc<quality_a){
2885    los[best].flatten();
2886    int b_pos=kBucketCanonicalize(los[best].bucket);
2887    los[best].p=los[best].bucket->buckets[b_pos];
2888    qc==pQuality(los[best].bucket->buckets[b_pos],c);
2889    //(best!=erg.to_reduce_u+1)
2890    if(qc<quality_a){
2891    red_object h=los[erg.to_reduce_u];
2892    los[erg.to_reduce_u]=los[best];
2893    los[best]=h;
2894    erg.reduce_by=erg.to_reduce_u;
2895    erg.fromS=FALSE;
2896    erg.to_reduce_u--;
2897    }
2898  }
2899      }
2900      else
2901      {
2902  assume(erg.to_reduce_u==erg.to_reduce_l);
2903  wlen_type quality_a=
2904        quality_of_pos_in_strat_S(erg.reduce_by,c);
2905  wlen_type qc=los[erg.to_reduce_u].guess_quality(c);
2906  if (qc<0) PrintS("Wrong wlen_type");
2907  if(qc<quality_a){
2908    int best=erg.to_reduce_u;
2909    los[best].flatten();
2910    int b_pos=kBucketCanonicalize(los[best].bucket);
2911    los[best].p=los[best].bucket->buckets[b_pos];
2912    qc=pQuality(los[best].bucket->buckets[b_pos],c);
2913    assume(qc>=0);
2914    if(qc<quality_a){
2915      BOOLEAN exp=FALSE;
2916      if(qc<=2){
2917         //Print("\n qc is %lld \n",qc);
2918         exp=TRUE;
2919      }
2920
2921      else {
2922         if (qc<quality_a/2)
2923          exp=TRUE;
2924         else
2925       if(erg.reduce_by<c->n/4)
2926          exp=TRUE;
2927      }
2928      if (exp){
2929        poly clear_into;
2930        los[erg.to_reduce_u].flatten();
2931        kBucketClear(los[erg.to_reduce_u].bucket,&clear_into,&erg.expand_length);
2932        erg.expand=pCopy(clear_into);
2933        kBucketInit(los[erg.to_reduce_u].bucket,clear_into,erg.expand_length);
2934        if (TEST_OPT_PROT)
2935    PrintS("e");
2936
2937      }
2938    }
2939  }
2940
2941
2942      }
2943
2944      swap_roles=FALSE;
2945      return;
2946      }
2947
2948  }
2949  else{
2950    if(erg.reduce_by>erg.to_reduce_u){
2951      //then lm(rb)>= lm(tru) so =
2952      assume(erg.reduce_by==erg.to_reduce_u+1);
2953      int best=erg.reduce_by;
2954      wlen_type quality_a=los[erg.reduce_by].guess_quality(c);
2955      int qc;
2956      best=find_best(los,erg.to_reduce_l,erg.to_reduce_u,qc,c);
2957
2958      int i;
2959      if(qc<quality_a){
2960    red_object h=los[erg.reduce_by];
2961    los[erg.reduce_by]=los[best];
2962    los[best]=h;
2963  }
2964  swap_roles=FALSE;
2965  return;
2966
2967
2968    }
2969    else
2970    {
2971      assume(!pLmEqual(los[erg.reduce_by].p,los[erg.to_reduce_l].p));
2972      assume(erg.to_reduce_u==erg.to_reduce_l);
2973      //further assume, that reduce_by is the above all other polys
2974      //with same leading term
2975      int il=erg.reduce_by;
2976      int quality_a =los[erg.reduce_by].guess_quality(c);
2977      int qc;
2978      while((il>0) && pLmEqual(los[il-1].p,los[il].p)){
2979  il--;
2980  qc=los[il].guess_quality(c);
2981  if (qc<quality_a){
2982    quality_a=qc;
2983    erg.reduce_by=il;
2984  }
2985      }
2986      swap_roles=FALSE;
2987    }
2988
2989  }
2990  if(swap_roles)
2991  {
2992    if (TEST_OPT_PROT)
2993      PrintS("b");
2994    poly clear_into;
2995    int dummy_len;
2996    int new_length;
2997    int bp=erg.to_reduce_u;//bucket_positon
2998    //kBucketClear(los[bp].bucket,&clear_into,&new_length);
2999    new_length=los[bp].clear_to_poly();
3000    clear_into=los[bp].p;
3001    poly p=c->strat->S[erg.reduce_by];
3002    int j=erg.reduce_by;
3003    int old_length=c->strat->lenS[j];// in view of S
3004    los[bp].p=p;
3005    kBucketInit(los[bp].bucket,p,old_length);
3006    wlen_type qal=pQuality(clear_into,c,new_length);
3007    int pos_in_c=-1;
3008    int z;
3009    int new_pos;
3010    new_pos=simple_posInS(c->strat,clear_into,new_length, qal);
3011    assume(new_pos<=j);
3012    for (z=c->n;z;z--)
3013    {
3014      if(p==c->S->m[z-1])
3015      {
3016  pos_in_c=z-1;
3017  break;
3018      }
3019    }
3020    if(pos_in_c>=0)
3021    {
3022      #ifdef TGB_RESORT_PAIRS
3023      c->used_b=TRUE;
3024      c->replaced[pos_in_c]=TRUE;
3025      #endif
3026      c->S->m[pos_in_c]=clear_into;
3027      c->lengths[pos_in_c]=new_length;
3028      c->weighted_lengths[pos_in_c]=qal;
3029      if (c->gcd_of_terms[pos_in_c]==NULL)
3030        c->gcd_of_terms[pos_in_c]=gcd_of_terms(clear_into,c->r);
3031      if (c->T_deg_full)
3032  c->T_deg_full[pos_in_c]=pTotaldegree_full(clear_into);
3033      c_S_element_changed_hook(pos_in_c,c);
3034    }
3035    c->strat->S[j]=clear_into;
3036    c->strat->lenS[j]=new_length;
3037
3038    assume(pLength(clear_into)==new_length);
3039    if(c->strat->lenSw!=NULL)
3040      c->strat->lenSw[j]=qal;
3041    if (!rField_is_Zp(c->r))
3042    {
3043      pContent(clear_into);
3044      pCleardenom(clear_into);//should be unnecessary
3045    }
3046    else
3047      pNorm(clear_into);
3048#ifdef FIND_DETERMINISTIC
3049    erg.reduce_by=j;
3050    //resort later see diploma thesis, find_in_S must be deterministic
3051    //during multireduction if spolys are only in the span of the
3052    //input polys
3053#else
3054
3055    if (new_pos<j)
3056    {
3057      move_forward_in_S(j,new_pos,c->strat);
3058      erg.reduce_by=new_pos;
3059    }
3060#endif
3061  }
3062}
3063static int fwbw(red_object* los, int i){
3064   int i2=i;
3065   int step=1;
3066
3067   BOOLEAN bw=FALSE;
3068   BOOLEAN incr=TRUE;
3069
3070   while(1)
3071   {
3072     if(!bw)
3073     {
3074       step=si_min(i2,step);
3075       if (step==0) break;
3076       i2-=step;
3077
3078       if(!pLmEqual(los[i].p,los[i2].p))
3079       {
3080   bw=TRUE;
3081   incr=FALSE;
3082       }
3083       else
3084       {
3085   if ((!incr) &&(step==1)) break;
3086       }
3087
3088
3089     }
3090     else
3091     {
3092
3093       step=si_min(i-i2,step);
3094       if (step==0) break;
3095       i2+=step;
3096       if(pLmEqual(los[i].p,los[i2].p)){
3097   if(step==1) break;
3098   else
3099   {
3100     bw=FALSE;
3101   }
3102       }
3103
3104     }
3105     if (incr)
3106       step*=2;
3107     else
3108     {
3109       if (step%2==1)
3110   step=(step+1)/2;
3111       else
3112   step/=2;
3113
3114     }
3115   }
3116   return i2;
3117}
3118static void canonicalize_region(red_object* los, int l, int u,slimgb_alg* c){
3119    assume(l<=u+1);
3120    int i;
3121    for(i=l;i<=u;i++){
3122        kBucketCanonicalize(los[i].bucket);
3123    }
3124
3125}
3126static void multi_reduction_find(red_object* los, int losl,slimgb_alg* c,int startf,find_erg & erg){
3127  kStrategy strat=c->strat;
3128
3129  assume(startf<=losl);
3130  assume((startf==losl-1)||(pLmCmp(los[startf].p,los[startf+1].p)==-1));
3131  int i=startf;
3132
3133  int j;
3134  while(i>=0){
3135    assume((i==losl-1)||(pLmCmp(los[i].p,los[i+1].p)<=0));
3136    assume(is_valid_ro(los[i]));
3137    j=kFindDivisibleByInS_easy(strat,los[i]);
3138    if(j>=0){
3139
3140      erg.to_reduce_u=i;
3141      erg.reduce_by=j;
3142      erg.fromS=TRUE;
3143      int i2=fwbw(los,i);
3144      assume(pLmEqual(los[i].p,los[i2].p));
3145      assume((i2==0)||(!pLmEqual(los[i2].p,los[i2-1].p)));
3146      assume(i>=i2);
3147
3148
3149      erg.to_reduce_l=i2;
3150      assume((i==losl-1)||(pLmCmp(los[i].p,los[i+1].p)==-1));
3151      canonicalize_region(los,erg.to_reduce_u+1,startf,c);
3152      return;
3153    }
3154    if (j<0){
3155
3156      //not reduceable, try to use this for reducing higher terms
3157      int i2=fwbw(los,i);
3158      assume(pLmEqual(los[i].p,los[i2].p));
3159      assume((i2==0)||(!pLmEqual(los[i2].p,los[i2-1].p)));
3160      assume(i>=i2);
3161      if(i2!=i){
3162
3163
3164  erg.to_reduce_u=i-1;
3165  erg.to_reduce_l=i2;
3166  erg.reduce_by=i;
3167  erg.fromS=FALSE;
3168  assume((i==losl-1)||(pLmCmp(los[i].p,los[i+1].p)==-1));
3169  canonicalize_region(los,erg.to_reduce_u+1,startf,c);
3170  return;
3171      }
3172
3173      i--;
3174    }
3175  }
3176  erg.reduce_by=-1;//error code
3177  return;
3178}
3179
3180 //  nicht reduzierbare eintraege in ergebnisliste schreiben
3181//   nullen loeschen
3182//   while(finde_groessten leitterm reduzierbar(c,erg)){
3183
3184static int multi_reduction_clear_zeroes(red_object* los, int  losl, int l, int u)
3185{
3186
3187
3188  int deleted=0;
3189  int  i=l;
3190  int last=-1;
3191  while(i<=u)
3192  {
3193
3194    if(los[i].p==NULL){
3195      kBucketDestroy(&los[i].bucket);
3196//      delete los[i];//here we assume los are constructed with new
3197      //destroy resources, must be added here
3198     if (last>=0)
3199     {
3200       memmove(los+(int)(last+1-deleted),los+(last+1),sizeof(red_object)*(i-1-last));
3201     }
3202     last=i;
3203     deleted++;
3204    }
3205    i++;
3206  }
3207  if((last>=0)&&(last!=losl-1))
3208      memmove(los+(int)(last+1-deleted),los+last+1,sizeof(red_object)*(losl-1-last));
3209  return deleted;
3210
3211}
3212
3213static void sort_region_down(red_object* los, int l, int u, slimgb_alg* c)
3214{
3215  qsort(los+l,u-l+1,sizeof(red_object),red_object_better_gen);
3216  int i;
3217
3218  for(i=l;i<=u;i++)
3219  {
3220    BOOLEAN moved=FALSE;
3221    int j;
3222    for(j=i;j;j--)
3223    {
3224      if(pLmCmp(los[j].p,los[j-1].p)==-1){
3225  red_object h=los[j];
3226  los[j]=los[j-1];
3227  los[j-1]=h;
3228  moved=TRUE;
3229      }
3230      else break;
3231    }
3232    if(!moved) return;
3233  }
3234}
3235
3236//assume that los is ordered ascending by leading term, all non zero
3237static void multi_reduction(red_object* los, int & losl, slimgb_alg* c)
3238{
3239  poly* delay=(poly*) omalloc(losl*sizeof(poly));
3240  int delay_s=0;
3241  //initialize;
3242  assume(c->strat->sl>=0);
3243  assume(losl>0);
3244  int i;
3245  wlen_type max_initial_quality=0;
3246
3247  for(i=0;i<losl;i++){
3248    los[i].sev=pGetShortExpVector(los[i].p);
3249//SetShortExpVector();
3250    los[i].p=kBucketGetLm(los[i].bucket);
3251    if (los[i].initial_quality>max_initial_quality)
3252        max_initial_quality=los[i].initial_quality;
3253    // else
3254//         Print("init2_qal=%lld;", los[i].initial_quality);
3255//     Print("initial_quality=%lld;",max_initial_quality);
3256  }
3257
3258  kStrategy strat=c->strat;
3259  int curr_pos=losl-1;
3260
3261
3262//  nicht reduzierbare einträge in ergebnisliste schreiben
3263  // nullen loeschen
3264  while(curr_pos>=0){
3265
3266    find_erg erg;
3267    multi_reduction_find(los, losl,c,curr_pos,erg);//last argument should be curr_pos
3268    if(erg.reduce_by<0) break;
3269
3270
3271
3272    erg.expand=NULL;
3273    int d=erg.to_reduce_u-erg.to_reduce_l+1;
3274
3275
3276    multi_reduction_lls_trick(los,losl,c,erg);
3277
3278
3279    int i;
3280    int len;
3281    //    wrp(los[erg.to_reduce_u].p);
3282    //Print("\n");
3283    multi_reduce_step(erg,los,c);
3284
3285
3286    if(!K_TEST_OPT_REDTHROUGH){
3287  for(i=erg.to_reduce_l;i<=erg.to_reduce_u;i++){
3288     if  (los[i].p!=NULL)  //the check (los[i].p!=NULL) might be invalid
3289     {
3290         //
3291         assume(los[i].initial_quality>0);
3292
3293               if(los[i].guess_quality(c)
3294                  >1.5*delay_factor*max_initial_quality){
3295                       if (TEST_OPT_PROT)
3296                           PrintS("v");
3297                       los[i].canonicalize();
3298                       if(los[i].guess_quality(c)
3299                           >delay_factor*max_initial_quality){
3300                               if (TEST_OPT_PROT)
3301                                   PrintS(".");
3302                               los[i].clear_to_poly();
3303                               //delay.push_back(los[i].p);
3304                               delay[delay_s]=los[i].p;
3305                               delay_s++;
3306
3307                               los[i].p=NULL;
3308
3309                      }
3310                  }
3311
3312            }
3313     }
3314  }
3315    int deleted=multi_reduction_clear_zeroes(los, losl, erg.to_reduce_l, erg.to_reduce_u);
3316    if(erg.fromS==FALSE)
3317      curr_pos=si_max(erg.to_reduce_u,erg.reduce_by);
3318    else
3319      curr_pos=erg.to_reduce_u;
3320    losl -= deleted;
3321    curr_pos -= deleted;
3322
3323    //Print("deleted %i \n",deleted);
3324    if ((TEST_V_UPTORADICAL) &&(!(erg.fromS)))
3325        sort_region_down(los,si_min(erg.to_reduce_l,erg.reduce_by),(si_max(erg.to_reduce_u,erg.reduce_by))-deleted,c);
3326    else
3327    sort_region_down(los, erg.to_reduce_l, erg.to_reduce_u-deleted, c);
3328
3329
3330    if(erg.expand)
3331    {
3332#ifdef FIND_DETERMINISTIC
3333      int i;
3334      for(i=0;c->expandS[i];i++);
3335      c->expandS=(poly*) omrealloc(c->expandS,(i+2)*sizeof(poly));
3336      c->expandS[i]=erg.expand;
3337      c->expandS[i+1]=NULL;
3338#else
3339      add_to_reductors(c,erg.expand,erg.expand_length);
3340#endif
3341    }
3342
3343  }
3344
3345
3346  sorted_pair_node** pairs=(sorted_pair_node**)
3347    omalloc(delay_s*sizeof(sorted_pair_node*));
3348  for(i=0;i<delay_s;i++){
3349
3350      poly p=delay[i];
3351      //if (rPar(c->r)==0)
3352      simplify_poly(p,c->r);
3353      sorted_pair_node* si=(sorted_pair_node*) omalloc(sizeof(sorted_pair_node));
3354      si->i=-1;
3355      si->j=-1;
3356       if (!rField_is_Zp(c->r)){
3357        if (!c->nc)
3358            p=redTailShort(p, c->strat);
3359        pCleardenom(p);
3360        pContent(p);
3361      }
3362      si->expected_length=pQuality(p,c,pLength(p));
3363      si->deg=pTotaldegree(p);
3364
3365      si->lcm_of_lm=p;
3366      pairs[i]=si;
3367  }
3368  qsort(pairs,delay_s,sizeof(sorted_pair_node*),tgb_pair_better_gen2);
3369  c->apairs=spn_merge(c->apairs,c->pair_top+1,pairs,delay_s,c);
3370  c->pair_top+=delay_s;
3371  omfree(delay);
3372  return;
3373}
3374void red_object::flatten(){
3375  assume(p==kBucketGetLm(bucket));
3376}
3377void red_object::validate(){
3378  p=kBucketGetLm(bucket);
3379  if(p)
3380    sev=pGetShortExpVector(p);
3381}
3382int red_object::clear_to_poly(){
3383  flatten();
3384  int l;
3385  kBucketClear(bucket,&p,&l);
3386  return l;
3387}
3388
3389
3390
3391
3392
3393void reduction_step::reduce(red_object* r, int l, int u){}
3394void simple_reducer::do_reduce(red_object & ro){
3395  number coef;
3396  if (!c->nc)
3397    coef=kBucketPolyRed(ro.bucket,p,
3398       p_len,
3399       c->strat->kNoether);
3400  else
3401    nc_kBucketPolyRed_Z(ro.bucket, p, &coef);
3402  nDelete(&coef);
3403}
3404
3405
3406void simple_reducer::reduce(red_object* r, int l, int u){
3407  this->pre_reduce(r,l,u);
3408  int i;
3409//debug start
3410  int im;
3411
3412
3413  for(i=l;i<=u;i++){
3414
3415
3416
3417    this->do_reduce(r[i]);
3418
3419  }
3420  for(i=l;i<=u;i++){
3421
3422    kBucketSimpleContent(r[i].bucket);
3423    r[i].validate();
3424    #ifdef TGB_DEBUG
3425    #endif
3426  }
3427}
3428reduction_step::~reduction_step(){}
3429simple_reducer::~simple_reducer(){
3430  if(fill_back!=NULL)
3431  {
3432    kBucketInit(fill_back,p,p_len);
3433  }
3434  fill_back=NULL;
3435
3436}
3437
3438void multi_reduce_step(find_erg & erg, red_object* r, slimgb_alg* c){
3439  static int id=0;
3440  id++;
3441  unsigned long sev;
3442    BOOLEAN lt_changed=FALSE;
3443  int rn=erg.reduce_by;
3444  poly red;
3445  int red_len;
3446  simple_reducer* pointer;
3447  BOOLEAN work_on_copy=FALSE;
3448  if(erg.fromS){
3449    red=c->strat->S[rn];
3450    red_len=c->strat->lenS[rn];
3451    assume(red_len==pLength(red));
3452  }
3453  else
3454  {
3455    r[rn].flatten();
3456    kBucketClear(r[rn].bucket,&red,&red_len);
3457
3458    if (!rField_is_Zp(c->r))
3459    {
3460      pContent(red);
3461      pCleardenom(red);//should be unnecessary
3462
3463    }
3464    pNormalize(red);
3465
3466
3467    if ((!(erg.fromS))&&(TEST_V_UPTORADICAL)){
3468
3469         if (polynomial_root(red,c->r))
3470            lt_changed=TRUE;
3471            sev=p_GetShortExpVector(red,c->r);}
3472    red_len=pLength(red);
3473  }
3474  if (((TEST_V_MODPSOLVSB)&&(red_len>1))||((c->nc)||(erg.to_reduce_u-erg.to_reduce_l>5))){
3475    work_on_copy=TRUE;
3476    // poly m=pOne();
3477    poly m=c->tmp_lm;
3478    pSetCoeff(m,nInit(1));
3479    for(int i=1;i<=pVariables;i++)
3480      pSetExp(m,i,(pGetExp(r[erg.to_reduce_l].p, i)-pGetExp(red,i)));
3481    pSetm(m);
3482    poly red_cp;
3483    if (!c->nc)
3484      red_cp=ppMult_mm(red,m);
3485    else
3486      red_cp=nc_mm_Mult_p(m, pCopy(red), c->r);
3487    if(!erg.fromS){
3488      kBucketInit(r[rn].bucket,red,red_len);
3489    }
3490    //now reduce the copy
3491    //static poly redNF2 (poly h,slimgb_alg* c , int &len, number&  m,int n)
3492
3493    if (!c->nc)
3494      redTailShort(red_cp,c->strat);
3495    //number mul;
3496    // red_len--;
3497//     red_cp->next=redNF2(red_cp->next,c,red_len,mul,c->average_length);
3498//     pSetCoeff(red_cp,nMult(red_cp->coef,mul));
3499//     nDelete(&mul);
3500//     red_len++;
3501    red=red_cp;
3502    red_len=pLength(red);
3503    // pDelete(&m);
3504
3505  }
3506  int i;
3507
3508
3509
3510  assume(red_len==pLength(red));
3511
3512  pointer=new simple_reducer(red,red_len,c);
3513
3514  if ((!work_on_copy) && (!erg.fromS))
3515    pointer->fill_back=r[rn].bucket;
3516  else
3517    pointer->fill_back=NULL;
3518  pointer->reduction_id=id;
3519  pointer->c=c;
3520
3521  pointer->reduce(r,erg.to_reduce_l, erg.to_reduce_u);
3522  if(work_on_copy) pDelete(&pointer->p);
3523  delete pointer;
3524  if (lt_changed){
3525    assume(!erg.fromS);
3526    r[erg.reduce_by].sev=sev;
3527  }
3528
3529};
3530
3531
3532
3533
3534void simple_reducer:: pre_reduce(red_object* r, int l, int u){}
3535
Note: See TracBrowser for help on using the repository browser.