source: git/kernel/tgb.cc @ 2fc974

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