source: git/kernel/tgb.cc @ e9ade0f

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