source: git/kernel/tgb.cc @ e65be8e

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