source: git/kernel/tgb.cc @ b16a613

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