source: git/kernel/tgb.cc @ 19370c

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