source: git/kernel/tgb.cc @ 9108d3

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