source: git/kernel/tgb.cc @ 391323

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