source: git/kernel/tgb.cc @ 89b59f

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