source: git/kernel/tgb.cc @ 06662e

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