source: git/kernel/tgb.cc @ 2fd387d

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