source: git/Singular/tgb.cc @ 1142b2

fieker-DuValspielwiese
Last change on this file since 1142b2 was 1142b2, checked in by Michael Brickenstein <bricken@…>, 21 years ago
*bricken: res freed git-svn-id: file:///usr/local/Singular/svn/trunk@6682 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 63.5 KB
Line 
1
2// #define OM_CHECK 3
3// #define OM_TRACK 5
4// #define OM_KEEP  1
5// TODO:
6//       deg -> poly_crit
7//       multiple rings
8//       shorten_tails und dessen Aufrufe pruefen wlength!!!
9//       calculating with formal sums
10//       try to create spolys as formal sums
11
12#include "tgb.h"
13
14#define LEN_VAR1
15#define degbound(p) assume(pTotaldegree(p)<10)
16#define inDebug(p) assume((debug_Ideal==NULL)||(kNF(debug_Ideal,NULL,p,0,0)==0))
17#ifdef LEN_VAR1
18// erste Variante: Laenge: Anzahl der Monome
19int pSLength(poly p, int l) {
20  return l; }
21int kSBucketLength(kBucket* bucket) {return bucket_guess(bucket);}
22#endif
23
24#ifdef LEN_VAR2
25// 2. Variante: Laenge: Platz fuer die Koeff.
26int pSLength(poly p,int l)
27{
28  int s=0;
29  while (p!=NULL) { s+=nSize(pGetCoeff(p));pIter(p); }
30  return s;
31}
32int kSBucketLength(kBucket* b)
33{
34  int s=0;
35  int i;
36  for (i=MAX_BUCKET;i>=0;i--)
37  {
38    s+=pSLength(b->buckets[i],0);
39  }
40  return s;
41}
42#endif
43
44#ifdef LEN_VAR3
45// 3.Variante: Laenge: Platz fuer Leitk * Monomanzahl
46int pSLength(poly p,int l)
47{
48  int c=nSize(pGetCoeff(p));
49  return c*l /*pLength(p)*/;
50}
51int kSBucketLength(kBucket* b)
52{
53  int s=0;
54  int c=nSize(pGetCoeff(kBucketGetLm(b)))+1;
55  int i;
56  for (i=MAX_BUCKET;i>0;i--)
57  {
58    s+=b->buckets_length[i] /*pLength(b->buckets[i])*/;
59  }
60  return s*c;
61}
62#endif
63
64#ifdef LEN_VAR4
65// 4.Variante: Laenge: Platz fuer Leitk * (1+Platz fuer andere Koeff.)
66int pSLength(poly p, int l)
67{
68  int s=1;
69  int c=nSize(pGetCoeff(p));
70  pIter(p);
71  while (p!=NULL) { s+=nSize(pGetCoeff(p));pIter(p); }
72  return s*c;
73}
74int kSBucketLength(kBucket* b)
75{
76  int s=1;
77  int c=nSize(pGetCoeff(kBucketGetLm(b)));
78  int i;
79  for (i=MAX_BUCKET;i>0;i--)
80  {
81    if(b->buckets[i]==NULL) continue;
82    s+=pSLength(b->buckets[i],0);
83  }
84  return s*c;
85}
86#endif
87
88static int LObject_better_gen(const void* ap, const void* bp)
89{
90  LObject* a=*(LObject**)ap;
91  LObject* b=*(LObject**)bp;
92  return(pLmCmp(a->p,b->p));
93}
94static int red_object_better_gen(const void* ap, const void* bp)
95{
96
97
98  return(pLmCmp(((red_object*) ap)->p,((red_object*) bp)->p));
99}
100static int pair_better_gen2(const void* ap,const void* bp){
101  return(-pair_better_gen(ap,bp));
102}
103static int kFindDivisibleByInS_easy(kStrategy strat,const red_object & obj){
104  int i;
105  long not_sev=~obj.sev;
106  poly p=obj.p;
107  for(i=0;i<=strat->sl;i++){
108    if (pLmShortDivisibleBy(strat->S[i],strat->sevS[i],p,not_sev))
109      return i;
110  }
111  return -1;
112}
113static int posInPairs (sorted_pair_node**  p, int pn, sorted_pair_node* qe,calc_dat* c,int an=0)
114{
115  if(pn==0) return 0;
116
117  int length=pn-1;
118  int i;
119  //int an = 0;
120  int en= length;
121
122  if (pair_better(qe,p[en],c))
123    return length+1;
124
125  while(1)
126    {
127      //if (an >= en-1)
128      if(en-1<=an)
129      {
130        if (pair_better(p[an],qe,c)) return an;
131        return en;
132      }
133      i=(an+en) / 2;
134        if (pair_better(p[i],qe,c))
135          en=i;
136      else an=i;
137    }
138}
139static BOOLEAN  ascending(int* i,int top){
140  if(top<1) return TRUE;
141  if(i[top]<i[top-1]) return FALSE;
142  return ascending(i,top-1);
143}
144
145sorted_pair_node**  merge(sorted_pair_node** p, int pn,sorted_pair_node **q, int qn,calc_dat* c){
146  int i;
147  int* a= (int*) omalloc(qn*sizeof(int));
148//   int mc;
149//   PrintS("Debug\n");
150//   for(mc=0;mc<qn;mc++)
151// {
152
153//     wrp(q[mc]->lcm_of_lm);
154//     PrintS("\n");
155// }
156//    PrintS("Debug they are in\n");
157//   for(mc=0;mc<pn;mc++)
158// {
159
160//     wrp(p[mc]->lcm_of_lm);
161//     PrintS("\n");
162// }
163  int lastpos=0;
164  for(i=0;i<qn;i++){
165    lastpos=posInPairs(p,pn,q[i],c, max(lastpos-1,0));
166    //   cout<<lastpos<<"\n";
167    a[i]=lastpos;
168
169  }
170  if((pn+qn)>c->max_pairs){
171    p=(sorted_pair_node**) omrealloc(p,2*(pn+qn)*sizeof(sorted_pair_node*));
172    c->max_pairs=2*(pn+qn);
173  }
174  for(i=qn-1;i>=0;i--){
175    size_t size;
176    if(qn-1>i)
177      size=(a[i+1]-a[i])*sizeof(sorted_pair_node*);
178    else
179      size=(pn-a[i])*sizeof(sorted_pair_node*); //as indices begin with 0
180    memmove (p+a[i]+(1+i), p+a[i], size);
181    p[a[i]+i]=q[i];
182  }
183  omfree(a);
184  return p;
185}
186
187
188static BOOLEAN trivial_syzygie(int pos1,int pos2,poly bound,calc_dat* c){
189
190
191  poly p1=c->S->m[pos1];
192  poly p2=c->S->m[pos2];
193  ring r=c->r;
194 
195
196  if (pGetComp(p1) > 0 || pGetComp(p2) > 0)
197    return FALSE;
198  int i = 1;
199  poly m=NULL;
200  poly gcd1=c->gcd_of_terms[pos1];
201  poly gcd2=c->gcd_of_terms[pos2];
202 
203  if((gcd1!=NULL) && (gcd2!=NULL)) 
204    {
205      gcd1->next=gcd2; //may ordered incorrect
206      poly m=gcd_of_terms(gcd1,c->r);
207      gcd1->next=NULL;
208     
209    } 
210
211  if (m==NULL) 
212  {
213     loop
214      {
215        if (pGetExp(p1, i)+ pGetExp(p2, i) > pGetExp(bound,i))   return FALSE;
216        if (i == pVariables){
217          PrintS("trivial");
218          return TRUE;
219        }
220        i++;
221      }
222  }
223  else 
224  {
225    loop
226      {
227        if (pGetExp(p1, i)-pGetExp(m,i) + pGetExp(p2, i) > pGetExp(bound,i))   return FALSE;
228        if (i == pVariables){
229          pDelete(&m);
230          PrintS("trivial");
231          return TRUE;
232        }
233        i++;
234      }
235  }
236
237 
238
239 
240}
241
242///returns position sets w as weight
243int find_best(red_object* r,int l, int u, int &w, calc_dat* c){
244
245  int sz=u-l+1;
246  int n=sz/10+1;
247  int filled=0;
248  int* indizes=(int*) omalloc(n*sizeof(int));
249  int* weight=(int*) omalloc(n*sizeof(int));
250  int worst=-1;
251  int i; 
252  for(i=l;i<=u;i++){
253    int q=r[i].guess_quality(c);
254    if ((filled<n)||(q<worst)){
255      if(filled<n){
256        worst=max(q,worst);
257        indizes[filled]=i;
258        weight[filled]=q;
259        filled++;
260      }
261    }
262    else{
263      int j;
264      for(j=0;j<filled;j++){
265        if (worst==weight[j]){
266          weight[j]=q;
267          indizes[j]=i;
268        }
269      }
270      worst=-1;
271      for(j=0;j<filled;j++){
272        if (worst<weight[j]){
273          worst=weight[j];
274        }
275      }
276    }
277  }
278  assume(filled==n);
279  int pos=0;
280
281  for(i=0;i<filled;i++){ 
282    r[indizes[i]].canonicalize();
283    weight[i]=r[indizes[i]].guess_quality(c);
284    if(weight[i]<weight[pos]) pos=i;
285  }
286  w=weight[pos];
287  pos=indizes[pos];
288
289  omfree(indizes);
290  omfree(weight);
291
292  assume(w==r[pos].guess_quality(c));
293  assume(l<=pos);
294  assume(u>=pos);
295  return pos;
296 
297}
298void red_object::canonicalize(){
299  kBucketCanonicalize(bucket);
300  if(sum)
301    kBucketCanonicalize(sum->ac->bucket);
302 
303}
304BOOLEAN good_has_t_rep(int i, int j,calc_dat* c){
305  assume(i>=0);
306    assume(j>=0);
307  if (has_t_rep(i,j,c)) return TRUE;
308  poly lm=pOne();
309
310  pLcm(c->S->m[i], c->S->m[j], lm);
311  pSetm(lm);
312  assume(lm!=NULL);
313  int deciding_deg= pTotaldegree(lm);
314  int* i_con =make_connections(i,j,lm,c);
315  p_Delete(&lm,c->r);
316
317
318  for (int n=0;((n<c->n) && (i_con[n]>=0));n++){
319    if (i_con[n]==j){
320      now_t_rep(i,j,c);
321      omfree(i_con);
322
323      return TRUE;
324    }
325  }
326  omfree(i_con);
327
328  return FALSE;
329}
330BOOLEAN lenS_correct(kStrategy strat){
331  int i;
332  for(i=0;i<=strat->sl;i++){
333    if (strat->lenS[i]!=pLength(strat->S[i]))
334      return FALSE;
335  }
336  return TRUE;
337}
338
339static void notice_miss(int i, int j, calc_dat* c){
340  PrintS("-");
341 
342}
343
344static void cleanS(kStrategy strat){
345  int i=0;
346  LObject P;
347  while(i<=strat->sl){
348    P.p=strat->S[i];
349    P.sev=strat->sevS[i];
350    if(kFindDivisibleByInS(strat->S,strat->sevS,strat->sl,&P)!=i){
351      deleteInS(i,strat);
352      //remember destroying poly
353    }
354    else i++;
355  }
356}
357static int bucket_guess(kBucket* bucket){
358  int sum=0;
359  int i;
360  for (i=MAX_BUCKET;i>=0;i--){
361    sum+=bucket->buckets_length[i];
362  }
363  return sum;
364}
365
366
367
368
369
370
371static int add_to_reductors(calc_dat* c, poly h, int len){
372  inDebug(h);
373  assume(lenS_correct(c->strat));
374 
375  int i;
376  if (c->is_char0)
377       i=simple_posInS(c->strat,h,pSLength(h,len),c->is_char0);
378  else
379    i=simple_posInS(c->strat,h,len,c->is_char0);
380 
381  LObject P; memset(&P,0,sizeof(P));
382  P.tailRing=c->r;
383  P.p=h; /*p_Copy(h,c->r);*/
384  P.FDeg=pFDeg(P.p,c->r);
385  if (!rField_is_Zp(c->r)){ 
386    pCleardenom(P.p);
387    pContent(P.p); //is a duplicate call, but belongs here
388  }
389 
390  else                     
391    pNorm(P.p);
392 
393 
394  c->strat->enterS(P,i,c->strat);
395 
396 
397
398  c->strat->lenS[i]=len;
399 
400  if(c->strat->lenSw)
401    c->strat->lenSw[i]=pSLength(P.p,len);
402  return i;
403 
404}
405static void length_one_crit(calc_dat* c, int pos, int len)
406{
407  if (len==1)
408  {
409    int i;
410    for ( i=0;i<pos;i++)
411    {
412      if (c->lengths[i]==1)
413        c->states[pos][i]=HASTREP;
414    }
415    for ( i=pos+1;i<c->n;i++){
416      if (c->lengths[i]==1)
417        c->states[i][pos]=HASTREP;
418    }
419    shorten_tails(c,c->S->m[pos]);
420  }
421}
422static sorted_pair_node* find_next_pair2(calc_dat* c, BOOLEAN go_higher){
423  clean_top_of_pair_list(c);
424  sorted_pair_node* s=pop_pair(c);
425
426
427  return s;
428}
429
430static void move_forward_in_S(int old_pos, int new_pos,kStrategy strat, BOOLEAN is_char0)
431{
432  assume(old_pos>=new_pos);
433  poly p=strat->S[old_pos];
434  int ecart=strat->ecartS[old_pos];
435  long sev=strat->sevS[old_pos];
436  int s_2_r=strat->S_2_R[old_pos];
437  int length=strat->lenS[old_pos];
438  int length_w;
439  if(is_char0)
440    length_w=strat->lenSw[old_pos];
441  int i;
442  for (i=old_pos; i>new_pos; i--)
443  {
444    strat->S[i] = strat->S[i-1];
445    strat->ecartS[i] = strat->ecartS[i-1];
446    strat->sevS[i] = strat->sevS[i-1];
447    strat->S_2_R[i] = strat->S_2_R[i-1];
448  }
449  if (strat->lenS!=NULL)
450    for (i=old_pos; i>new_pos; i--)
451      strat->lenS[i] = strat->lenS[i-1];
452  if (strat->lenSw!=NULL)
453    for (i=old_pos; i>new_pos; i--)
454      strat->lenSw[i] = strat->lenSw[i-1];
455
456  strat->S[new_pos]=p;
457  strat->ecartS[new_pos]=ecart;
458  strat->sevS[new_pos]=sev;
459  strat->S_2_R[new_pos]=s_2_r;
460  strat->lenS[new_pos]=length;
461  if(is_char0)
462    strat->lenSw[new_pos]=length_w;
463  //assume(lenS_correct(strat));
464}
465static void replace_pair(int & i, int & j, calc_dat* c)
466{
467  c->soon_free=NULL;
468  int curr_deg;
469  poly lm=pOne();
470
471  pLcm(c->S->m[i], c->S->m[j], lm);
472  pSetm(lm);
473  int deciding_deg= pTotaldegree(lm);
474  int* i_con =make_connections(i,j,lm,c);
475  int z=0;
476
477  for (int n=0;((n<c->n) && (i_con[n]>=0));n++){
478    if (i_con[n]==j){
479      now_t_rep(i,j,c);
480      omfree(i_con);
481      p_Delete(&lm,c->r);
482      return;
483    }
484  }
485
486  int* j_con =make_connections(j,lm,c);
487  i= i_con[0];
488  j=j_con[0];
489  if(c->n>1){
490    if (i_con[1]>=0)
491      i=i_con[1];
492    else {
493      if (j_con[1]>=0)
494        j=j_con[1];
495    }
496  }
497  pLcm(c->S->m[i], c->S->m[j], lm);
498  pSetm(lm);
499  poly short_s;
500  curr_deg=pTotaldegree(lm);
501  int_pair_node* last=NULL;
502
503  for (int n=0;((n<c->n) && (j_con[n]>=0));n++){
504    for (int m=0;((m<c->n) && (i_con[m]>=0));m++){
505      pLcm(c->S->m[i_con[m]], c->S->m[j_con[n]], lm);
506      pSetm(lm);
507      if (pTotaldegree(lm)>=deciding_deg)
508      {
509        soon_t_rep(i_con[m],j_con[n],c);
510        int_pair_node* h= (int_pair_node*)omalloc(sizeof(int_pair_node));
511        if (last!=NULL)
512          last->next=h;
513        else
514          c->soon_free=h;
515        h->next=NULL;
516        h->a=i_con[m];
517        h->b=j_con[n];
518        last=h;
519      }
520      //      if ((comp_deg<curr_deg)
521      //  ||
522      //  ((comp_deg==curr_deg) &&
523      short_s=ksCreateShortSpoly(c->S->m[i_con[m]],c->S->m[j_con[n]],c->r);
524      if (short_s==NULL) {
525        i=i_con[m];
526        j=j_con[n];
527        now_t_rep(i_con[m],j_con[n],c);
528        p_Delete(&lm,c->r);
529        omfree(i_con);
530        omfree(j_con);
531
532        return;
533      }
534#ifdef QUICK_SPOLY_TEST
535      for (int dz=0;dz<=c->n;dz++){
536        if (dz==c->n) {
537          //have found not head reducing pair
538          i=i_con[m];
539          j=j_con[n];
540          p_Delete(&short_s,c->r);
541          p_Delete(&lm,c->r);
542          omfree(i_con);
543          omfree(j_con);
544
545          return;
546        }
547        if (p_LmDivisibleBy(c->S->m[dz],short_s,c->r)) break;
548      }
549#endif
550      int comp_deg(pTotaldegree(short_s));
551      p_Delete(&short_s,c->r);
552      if ((comp_deg<curr_deg))
553         
554      {
555        curr_deg=comp_deg;
556        i=i_con[m];
557        j=j_con[n];
558      }
559    }
560  }
561  p_Delete(&lm,c->r);
562  omfree(i_con);
563  omfree(j_con);
564  return;
565}
566
567
568static int* make_connections(int from, poly bound, calc_dat* c)
569{
570  ideal I=c->S;
571  int s=pTotaldegree(bound);
572  int* cans=(int*) omalloc(c->n*sizeof(int));
573  int* connected=(int*) omalloc(c->n*sizeof(int));
574  int cans_length=0;
575  connected[0]=from;
576  int connected_length=1;
577  long neg_bounds_short= ~p_GetShortExpVector(bound,c->r);
578  for (int i=0;i<c->n;i++){
579    if (c->T_deg[i]>s) continue;
580    if (i!=from){
581      if(p_LmShortDivisibleBy(I->m[i],c->short_Exps[i],bound,neg_bounds_short,c->r)){
582        cans[cans_length]=i;
583        cans_length++;
584      }
585    }
586  }
587  int not_yet_found=cans_length;
588  int con_checked=0;
589  int pos;
590  while((not_yet_found>0) && (con_checked<connected_length)){
591    pos=connected[con_checked];
592    for(int i=0;i<cans_length;i++){
593      if (cans[i]<0) continue;
594      if (has_t_rep(pos,cans[i],c))
595      {
596        connected[connected_length]=cans[i];
597        connected_length++;
598        cans[i]=-1;
599        --not_yet_found;
600      }
601    }
602    con_checked++;
603  }
604  if (connected_length<c->n){
605    connected[connected_length]=-1;
606  }
607  omfree(cans);
608  return connected;
609}
610static int* make_connections(int from, int to, poly bound, calc_dat* c)
611{
612  ideal I=c->S;
613  int s=pTotaldegree(bound);
614  int* cans=(int*) omalloc(c->n*sizeof(int));
615  int* connected=(int*) omalloc(c->n*sizeof(int));
616  cans[0]=to;
617  int cans_length=1;
618  connected[0]=from;
619  int last_cans_pos=-1;
620  int connected_length=1;
621  long neg_bounds_short= ~p_GetShortExpVector(bound,c->r);
622
623  int not_yet_found=cans_length;
624  int con_checked=0;
625  int pos;
626  BOOLEAN can_find_more=TRUE;
627  while(((not_yet_found>0) && (con_checked<connected_length))||can_find_more){
628    if ((con_checked<connected_length)&& (not_yet_found>0)){
629      pos=connected[con_checked];
630      for(int i=0;i<cans_length;i++){
631        if (cans[i]<0) continue;
632        if (has_t_rep(pos,cans[i],c))//||(trivial_syzygie(pos,cans[i],bound,c))
633{
634
635          connected[connected_length]=cans[i];
636          connected_length++;
637          cans[i]=-1;
638          --not_yet_found;
639
640          if (connected[connected_length-1]==to){
641            if (connected_length<c->n){
642              connected[connected_length]=-1;
643            }
644            omfree(cans);
645            return connected;
646          }
647        }
648      }
649      con_checked++;
650    }
651    else
652    {
653      for(last_cans_pos++;last_cans_pos<=c->n;last_cans_pos++){
654        if (last_cans_pos==c->n){
655          if (connected_length<c->n){
656            connected[connected_length]=-1;
657          }
658          omfree(cans);
659          return connected;
660        }
661        if ((last_cans_pos==from)||(last_cans_pos==to))
662          continue;
663        if(p_LmShortDivisibleBy(I->m[last_cans_pos],c->short_Exps[last_cans_pos],bound,neg_bounds_short,c->r)){
664          cans[cans_length]=last_cans_pos;
665          cans_length++;
666          break;
667        }
668      }
669      not_yet_found++;
670      for (int i=0;i<con_checked;i++){
671        if (has_t_rep(connected[i],last_cans_pos,c)){
672
673          connected[connected_length]=last_cans_pos;
674          connected_length++;
675          cans[cans_length-1]=-1;
676
677          --not_yet_found;
678          if (connected[connected_length-1]==to){
679            if (connected_length<c->n){
680              connected[connected_length]=-1;
681            }
682
683            omfree(cans);
684            return connected;
685          }
686          break;
687        }
688      }
689    }
690  }
691  if (connected_length<c->n){
692    connected[connected_length]=-1;
693  }
694
695  omfree(cans);
696  return connected;
697}
698#ifdef HEAD_BIN
699static inline poly p_MoveHead(poly p, omBin b)
700{
701  poly np;
702  omTypeAllocBin(poly, np, b);
703  memmove(np, p, b->sizeW*sizeof(long));
704  omFreeBinAddr(p);
705  return np;
706}
707#endif
708
709
710
711//len should be weighted length in char 0
712static int simple_posInS (kStrategy strat, poly p,int len, BOOLEAN is_char0)
713{
714
715
716  if(strat->sl==-1) return 0;
717  polyset set=strat->S;
718  intset setL=strat->lenS;
719  if (is_char0) setL=strat->lenSw;
720  int length=strat->sl;
721  int i;
722  int an = 0;
723  int en= length;
724
725  if ((len>setL[length])
726      || ((len==setL[length]) && (pLmCmp(set[length],p)== -1)))
727    return length+1;
728
729  loop
730  {
731    if (an >= en-1)
732    {
733      if ((len<setL[an])
734          || ((len==setL[an]) && (pLmCmp(set[an],p) == 1))) return an;
735      return en;
736    }
737    i=(an+en) / 2;
738    if ((len<setL[i])
739        || ((len==setL[i]) && (pLmCmp(set[i],p) == 1))) en=i;
740    //else if ((len>setL[i])
741    //|| ((len==setL[i]) && (pLmCmp(set[i],p) == -1))) an=i;
742    else an=i;
743  }
744}
745/*2
746 *if the leading term of p
747 *divides the leading term of some S[i] it will be canceled
748 */
749static inline void clearS (poly p, unsigned long p_sev,int l, int* at, int* k,
750                           kStrategy strat)
751{
752  assume(p_sev == pGetShortExpVector(p));
753  if (!pLmShortDivisibleBy(p,p_sev, strat->S[*at], ~ strat->sevS[*at])) return;
754  if (l>=strat->lenS[*at]) return;
755  PrintS("!");mflush();
756  //pDelete(&strat->S[*at]);
757  deleteInS((*at),strat);
758  (*at)--;
759  (*k)--;
760//  assume(lenS_correct(strat));
761}
762static sorted_pair_node** add_to_basis(poly h, int i_pos, int j_pos,calc_dat* c, int* ip)
763{
764
765  assume(h!=NULL);
766//  BOOLEAN corr=lenS_correct(c->strat);
767  BOOLEAN R_found=FALSE;
768  void* hp;
769
770  ++(c->n);
771  ++(c->S->ncols);
772  int i,j;
773  i=c->n-1;
774  sorted_pair_node** nodes=(sorted_pair_node**) omalloc(sizeof(sorted_pair_node*)*i);
775  int spc=0;
776  c->T_deg=(int*) omrealloc(c->T_deg,c->n*sizeof(int));
777  c->T_deg[i]=pTotaldegree(h);
778  hp=omrealloc(c->rep, c->n *sizeof(int));
779  if (hp!=NULL){
780    c->rep=(int*) hp;
781  } else {
782    exit(1);
783  }
784  c->short_Exps=(long *) omrealloc(c->short_Exps ,c->n*sizeof(long));
785
786  hp=omrealloc(c->lengths, c->n *sizeof(int));
787  if (hp!=NULL){
788    c->lengths=(int*) hp;
789  } else {
790    exit(1);
791  }
792  c->lengths[i]=pLength(h);
793  hp=omrealloc(c->states, c->n * sizeof(char*));
794 
795    c->states=(char**) hp;
796  c->gcd_of_terms=(poly*) omrealloc(c->gcd_of_terms, c->n *sizeof(poly));
797  c->gcd_of_terms[i]=gcd_of_terms(h,c->r);
798  c->rep[i]=i;
799  hp=omalloc(i*sizeof(char));
800  if (hp!=NULL){
801    c->states[i]=(char*) hp;
802  } else {
803    exit(1);
804  }
805  hp=omrealloc(c->S->m,c->n*sizeof(poly));
806  if (hp!=NULL){
807    c->S->m=(poly*) hp;
808  } else {
809    exit(1);
810  }
811  c->S->m[i]=h;
812  c->short_Exps[i]=p_GetShortExpVector(h,c->r);
813  for (j=0;j<i;j++){
814    if (c->rep[j]==j){
815      //check product criterion
816
817      c->states[i][j]=UNCALCULATED;
818
819      //lies I[i] under I[j] ?
820      if(p_LmShortDivisibleBy(c->S->m[i],c->short_Exps[i],c->S->m[j],~(c->short_Exps[j]),c->r)){
821        c->rep[j]=i;
822       
823        PrintS("R"); R_found=TRUE;
824
825        c->Rcounter++;
826        if((i_pos>=0) && (j_pos>=0)){
827       
828        }
829        for(int z=0;z<j;z++){
830          if(c->rep[z]!=z) continue;
831          if (c->states[j][z]==UNCALCULATED){
832            c->states[j][z]=UNIMPORTANT;
833          }
834        }
835        for(int z=j+1;z<i;z++){
836          if(c->rep[z]!=z) continue;
837          if (c->states[z][j]==UNCALCULATED){
838            c->states[z][j]=UNIMPORTANT;
839          }
840        }
841      }
842    }
843    else {
844      c->states[i][j]=UNIMPORTANT;
845    }
846    if ((c->lengths[i]==1) && (c->lengths[j]==1))
847      c->states[i][j]=HASTREP;
848    else if (pHasNotCF(c->S->m[i],c->S->m[j])){
849      c->easy_product_crit++;
850      c->states[i][j]=HASTREP;
851    }
852                        else if(extended_product_criterion(c->S->m[i],c->gcd_of_terms[i],c->S->m[j],c->gcd_of_terms[j],c)){
853                                        c->states[i][j]=HASTREP;
854                                        c->extended_product_crit++;
855                                        //PrintS("E");
856                        }
857    if (c->states[i][j]==UNCALCULATED){
858
859     
860      poly short_s=ksCreateShortSpoly(c->S->m[i],c->S->m[j],c->r);
861      if (short_s)
862      {
863        sorted_pair_node* s=(sorted_pair_node*) omalloc(sizeof(sorted_pair_node));
864        s->i=max(i,j);
865        s->j=min(i,j);
866        s->expected_length=c->lengths[i]+c->lengths[j]-2;
867        s->deg=pTotaldegree(short_s);
868        poly lm=pOne();
869
870        pLcm(c->S->m[i], c->S->m[j], lm);
871        pSetm(lm);
872        s->lcm_of_lm=lm;
873          pDelete(&short_s);
874        //assume(lm!=NULL);
875        nodes[spc]=s;
876        spc++;
877      }
878      else
879      {
880        c->states[i][j]=HASTREP;
881      }
882    }
883  }
884
885  add_to_reductors(c, h, c->lengths[c->n-1]);
886  //i=posInS(c->strat,c->strat->sl,h,0 /*ecart*/);
887
888  if (c->lengths[c->n-1]==1)
889    shorten_tails(c,c->S->m[c->n-1]);
890  //you should really update c->lengths, c->strat->lenS, and the oder of polys in strat if you sort after lengths
891
892  //for(i=c->strat->sl; i>0;i--)
893  //  if(c->strat->lenS[i]<c->strat->lenS[i-1]) printf("fehler bei %d\n",i);
894  if (c->Rcounter>50) {
895    c->Rcounter=0;
896    cleanS(c->strat);
897  }
898  if(!ip){
899    qsort(nodes,spc,sizeof(sorted_pair_node*),pair_better_gen2);
900 
901   
902    c->apairs=merge(c->apairs,c->pair_top+1,nodes,spc,c);
903    c->pair_top+=spc;
904    clean_top_of_pair_list(c);
905    omfree(nodes);
906    return NULL;
907  }
908  {
909    *ip=spc;
910    return nodes;
911  }
912
913 
914
915}
916#if 0
917static poly redNF (poly h,kStrategy strat)
918{
919  int j = 0;
920  int z = 3;
921  unsigned long not_sev;
922
923  if (0 > strat->sl)
924  {
925    return h;
926  }
927  not_sev = ~ pGetShortExpVector(h);
928  loop
929    {
930      if (pLmShortDivisibleBy(strat->S[j], strat->sevS[j], h, not_sev))
931      {
932        //if (strat->interpt) test_int_std(strat->kIdeal);
933        /*- compute the s-polynomial -*/
934#ifdef KDEBUG
935        if (TEST_OPT_DEBUG)
936        {
937          PrintS("red:");
938          wrp(h);
939          PrintS(" with ");
940          wrp(strat->S[j]);
941        }
942#endif
943        h = ksOldSpolyRed(strat->S[j],h,strat->kNoether);
944#ifdef KDEBUG
945        if (TEST_OPT_DEBUG)
946        {
947          PrintS("\nto:");
948          wrp(h);
949          PrintLn();
950        }
951#endif
952        if (h == NULL) return NULL;
953        z++;
954        if (z>=10)
955        {
956          z=0;
957          pNormalize(h);
958        }
959        /*- try to reduce the s-polynomial -*/
960        j = 0;
961        not_sev = ~ pGetShortExpVector(h);
962      }
963      else
964      {
965        if (j >= strat->sl) return h;
966        j++;
967      }
968    }
969}
970#else
971
972static poly redNF2 (poly h,calc_dat* c , int &len, number&  m,int n)
973{
974  m=nInit(1);
975  if (h==NULL) return NULL;
976
977  assume(len==pLength(h));
978  kStrategy strat=c->strat;
979  if (0 > strat->sl)
980  {
981    return h;
982  }
983  int j;
984 
985  LObject P(h);
986  P.SetShortExpVector();
987  P.bucket = kBucketCreate(currRing);
988  // BOOLEAN corr=lenS_correct(strat);
989  kBucketInit(P.bucket,P.p,len /*pLength(P.p)*/);
990  //int max_pos=simple_posInS(strat,P.p);
991  loop
992    {
993
994      j=kFindDivisibleByInS(strat->S,strat->sevS,strat->sl,&P);
995      if ((j>=0) && ((!n)||(strat->lenS[j]<=n)))
996      {
997
998       
999        nNormalize(pGetCoeff(P.p));
1000#ifdef KDEBUG
1001        if (TEST_OPT_DEBUG)
1002        {
1003          PrintS("red:");
1004          wrp(h);
1005          PrintS(" with ");
1006          wrp(strat->S[j]);
1007        }
1008#endif
1009       
1010        number coef=kBucketPolyRed(P.bucket,strat->S[j],
1011                                   strat->lenS[j]/*pLength(strat->S[j])*/,
1012                                   strat->kNoether);
1013        number m2=nMult(m,coef);
1014        nDelete(&m);
1015        m=m2;
1016        nDelete(&coef);
1017        h = kBucketGetLm(P.bucket);
1018       
1019        if (h==NULL) {
1020          len=0;
1021          return 
1022          NULL;}
1023        P.p=h;
1024        P.t_p=NULL;
1025        P.SetShortExpVector();
1026#ifdef KDEBUG
1027        if (TEST_OPT_DEBUG)
1028        {
1029          PrintS("\nto:");
1030          wrp(h);
1031          PrintLn();
1032        }
1033#endif
1034      }
1035      else
1036      {
1037        kBucketClear(P.bucket,&(P.p),&len);
1038        kBucketDestroy(&P.bucket);
1039        pNormalize(P.p);
1040        assume(len==(pLength(P.p)));
1041        return P.p;
1042      }
1043    }
1044}
1045
1046
1047static poly redTailShort(poly h, kStrategy strat){
1048
1049  int sl=strat->sl;
1050  int i;
1051  int len=pLength(h);
1052  for(i=0;i<=strat->sl;i++){
1053    if(strat->lenS[i]>2)
1054      break;
1055  }
1056  return(redNFTail(h,i-1,strat, len));
1057}
1058
1059static void line_of_extended_prod(int fixpos,calc_dat* c){
1060    if (c->gcd_of_terms[fixpos]==NULL)
1061  {
1062    c->gcd_of_terms[fixpos]=gcd_of_terms(c->S->m[fixpos],c->r);
1063    if (c->gcd_of_terms[fixpos])
1064    {
1065      int i;
1066      for(i=0;i<fixpos;i++)
1067        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)))
1068{
1069          c->states[fixpos][i]=HASTREP;
1070          c->extended_product_crit++;
1071}     
1072      for(i=fixpos+1;i<c->n;i++)
1073        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)))
1074        {        c->states[i][fixpos]=HASTREP;
1075        c->extended_product_crit++;
1076        }
1077    }
1078  }
1079}
1080static void c_S_element_changed_hook(int pos, calc_dat* c){
1081  length_one_crit(c,pos, c->lengths[pos]);
1082  line_of_extended_prod(pos,c);
1083}
1084class poly_tree_node {
1085public:
1086  poly p;
1087  poly_tree_node* l;
1088  poly_tree_node* r;
1089  int n;
1090  poly_tree_node(int sn):l(NULL),r(NULL),n(sn){}
1091};
1092class exp_number_builder{
1093public:
1094  poly_tree_node* top_level;
1095  int n;
1096  int get_n(poly p);
1097  exp_number_builder():top_level(0),n(0){}
1098};
1099int exp_number_builder::get_n(poly p){
1100  poly_tree_node** node=&top_level;
1101  while(*node!=NULL){
1102    int c=pLmCmp(p,(*node)->p);
1103    if (c==0) return (*node)->n;
1104    if (c==-1) node=&((*node)->r);
1105    else
1106      node=&((*node)->l);
1107  }
1108  (*node)= new poly_tree_node(n);
1109  n++;
1110  (*node)->p=pLmInit(p);
1111  return (*node)->n;
1112}
1113class mac_poly_r{
1114public:
1115  number coef;
1116  mac_poly_r* next;
1117  int exp;
1118  mac_poly_r():next(NULL){}
1119};
1120//mac_polys exp are smaller iff they are greater by monomial ordering
1121//corresponding to solving linear equations notation
1122struct int_poly_pair{
1123  poly p;
1124  int n;
1125};
1126typedef mac_poly_r* mac_poly;
1127  void t2ippa_rec(poly* ip,int* ia, poly_tree_node* k, int &offset){
1128    if(!k) return;
1129    t2ippa_rec(ip,ia,k->l,offset);
1130    ip[offset]=k->p;
1131    ia[k->n]=offset;
1132    ++offset;
1133
1134    t2ippa_rec(ip,ia,k->r,offset);
1135    delete k;
1136  }
1137void t2ippa(poly* ip,int* ia,exp_number_builder & e){
1138
1139  int o=0;
1140  t2ippa_rec(ip,ia,e.top_level,o);
1141}
1142int anti_poly_order(const void* a, const void* b){
1143  return -pLmCmp(((int_poly_pair*) a)->p,((int_poly_pair*) b)->p );
1144}
1145mac_poly mac_p_add_ff_qq(mac_poly a, number f,mac_poly b){
1146  mac_poly erg;
1147  mac_poly* set_this;
1148  set_this=&erg;
1149  while((a!=NULL) &&(b!=NULL)){
1150    if (a->exp<b->exp){
1151      (*set_this)=a;
1152      a=a->next;
1153      set_this= &((*set_this)->next);
1154    } 
1155    else{
1156      if (a->exp>b->exp){
1157        mac_poly in =new mac_poly_r();
1158        in->exp=b->exp;
1159        in->coef=nMult(b->coef,f);
1160        (*set_this)=in;
1161        b=b->next;
1162        set_this= &((*set_this)->next);
1163      }
1164      else {
1165        //a->exp==b->ecp
1166        number n=nMult(b->coef,f);
1167        number n2=nAdd(a->coef,n);
1168        nDelete(&n);
1169        nDelete(&(a->coef));
1170        if (nIsZero(n2)){
1171          nDelete(&n2);
1172          mac_poly ao=a;
1173          a=a->next;
1174          delete ao;
1175          b=b->next;
1176         
1177        } else {
1178          a->coef=n2;
1179          b=b->next;
1180          (*set_this)=a;
1181          a=a->next;
1182          set_this= &((*set_this)->next);
1183        }
1184 
1185      }
1186   
1187    }
1188  }
1189  if((a==NULL)&&(b==NULL)){
1190    (*set_this)=NULL;
1191    return erg;
1192  }
1193  if (b==NULL) {
1194    (*set_this=a);
1195    return erg;
1196  }
1197 
1198  //a==NULL
1199  while(b!=NULL){
1200    mac_poly mp= new mac_poly_r();
1201    mp->exp=b->exp;
1202    mp->coef=nMult(f,b->coef);
1203    (*set_this)=mp;
1204    set_this=&(mp->next);
1205    b=b->next;
1206  }
1207  (*set_this)=NULL;
1208  return erg;
1209 
1210}
1211void mac_mult_cons(mac_poly p,number c){
1212  while(p){
1213    number m=nMult(p->coef,c);
1214    nDelete(&(p->coef));
1215    p->coef=m;
1216    p=p->next;
1217  }
1218 
1219}
1220int mac_length(mac_poly p){
1221  int l=0;
1222  while(p){
1223    l++;
1224    p=p->next;
1225  }
1226  return l;
1227}
1228void pre_comp(poly* p,int & pn,calc_dat* c){
1229  if(!(pn))
1230    return;
1231  mac_poly* q=(mac_poly*) omalloc(pn*sizeof(mac_poly)); 
1232  int i;
1233  exp_number_builder e;
1234  for(i=0;i<pn;i++){
1235    assume(p[i]!=NULL);
1236    q[i]=new mac_poly_r();
1237    q[i]->exp=e.get_n(p[i]);
1238    q[i]->coef=nCopy(p[i]->coef);
1239   
1240    mac_poly qa=q[i];
1241    poly pa=p[i];
1242    while(pa->next!=NULL){
1243      qa->next = new mac_poly_r();
1244      qa=qa->next;
1245      pa=pa->next,
1246      qa->exp=e.get_n(pa);
1247      qa->coef=nCopy(pa->coef);
1248     
1249     
1250    }
1251    qa->next=NULL;
1252    pDelete(&p[i]);
1253  }
1254  poly* ip= (poly*)omalloc (e.n*sizeof(poly));
1255  int* ia=(int*) omalloc (e.n*sizeof(int));
1256  t2ippa(ip,ia,e);
1257  for(i=0;i<pn;i++){
1258    mac_poly mp=q[i];
1259    while(mp!=NULL){
1260      mp->exp=ia[mp->exp];
1261      mp=mp->next;
1262    }
1263   
1264  }
1265//gaus reduction start
1266  int col, row;
1267  col=0;
1268  row=0;
1269  assume(pn>0);
1270  while(row<pn-1){
1271    //row is the row where pivot should be
1272    // row== pn-1 means we have only to act on one row so no red nec.
1273    //we assume further all rows till the pn-1 row are non-zero
1274   
1275    //select column
1276    int i;
1277    col=q[row]->exp;
1278    int found_in_row=row;
1279    for(i=row;i<pn;i++){
1280      if(q[i]->exp<col){
1281        col=q[i]->exp;
1282        found_in_row=i;
1283      }
1284     
1285    }
1286    //select pivot
1287    int act_l=mac_length(q[found_in_row]);
1288    for(i=row+1;i<pn;i++){
1289      if((q[i]->exp==col)&&(mac_length(q[i])<act_l)){
1290        found_in_row=i;
1291        act_l=mac_length(q[i]);//should be optimized here
1292      }
1293    }
1294    mac_poly h=q[row];
1295    q[row]=q[found_in_row];
1296    q[found_in_row]=h;
1297
1298    //reduction
1299    for(i=row+1;i<pn;i++){
1300      if(q[i]->exp==q[row]->exp){
1301       
1302        number c1=nNeg(q[i]->coef);
1303        number c2=q[row]->coef;
1304        //use checkcoeff later
1305        mac_mult_cons(q[i],c2);
1306        q[i]=mac_p_add_ff_qq(q[i],c1,q[row]);
1307      }
1308         
1309       
1310       
1311    }
1312    for(i=row+1;i<pn;i++){
1313      if(q[i]==NULL){
1314        q[i]=q[pn-1];
1315        pn--;
1316        if(i!=pn){i--;}
1317      }
1318    }
1319 
1320
1321
1322
1323
1324
1325
1326
1327    row++;
1328  }
1329
1330
1331//gaus reduction end 
1332
1333
1334
1335
1336
1337  for(i=0;i<pn;i++){
1338    poly pa;
1339    mac_poly qa;
1340    p[i]=pLmInit(ip[q[i]->exp]);
1341    pSetCoeff(p[i],q[i]->coef);
1342    pa=p[i];
1343    qa=q[i];
1344    while(qa->next!=NULL){
1345      qa=qa->next;
1346      pa->next=pLmInit(ip[qa->exp]);
1347      pa=pa->next;
1348      pa->coef=qa->coef;
1349    }
1350    pa->next=NULL;
1351  }
1352  for(i=0;i<e.n;i++){
1353    pDelete(&ip[i]); 
1354  }
1355  omfree(ip);
1356  omfree(ia);
1357}
1358
1359static void go_on (calc_dat* c){
1360  //set limit of 1000 for multireductions, at the moment for
1361  //programming reasons
1362  int i=0;
1363
1364  poly* p=(poly*) omalloc((PAR_N+1)*sizeof(poly));//nullterminated
1365
1366  int curr_deg=-1;
1367  while(i<PAR_N){
1368    sorted_pair_node* s=top_pair(c);
1369    if (!s) break;
1370    if(curr_deg>=0){
1371      if (s->deg >curr_deg) break;
1372    }
1373
1374    else curr_deg=s->deg;
1375    quick_pop_pair(c);
1376    if(s->i>=0){
1377      //replace_pair(s->i,s->j,c);
1378    if(s->i==s->j) {
1379      free_sorted_pair_node(s,c->r);
1380      continue;
1381        }
1382    }
1383    poly h;
1384    if(s->i>=0)
1385      h=ksOldCreateSpoly(c->S->m[s->i], c->S->m[s->j], NULL, c->r);
1386    else
1387      h=s->lcm_of_lm;
1388    if(s->i>=0)
1389      now_t_rep(s->j,s->i,c);
1390    number coef;
1391    int mlen=pLength(h);
1392    h=redNF2(h,c,mlen,coef,2);
1393    redTailShort(h,c->strat);
1394    nDelete(&coef);
1395
1396    free_sorted_pair_node(s,c->r);
1397    if(!h) continue;
1398    int len=pLength(h);
1399    p[i]=h;
1400   
1401    i++;
1402  }
1403  p[i]=NULL;
1404  pre_comp(p,i,c);
1405  red_object* buf=(red_object*) omalloc(i*sizeof(red_object));
1406  c->normal_forms+=i;
1407  int j;
1408  for(j=0;j<i;j++){
1409    buf[j].p=p[j];
1410    buf[j].sev=pGetShortExpVector(p[j]);
1411    buf[j].sum=NULL;
1412    buf[j].bucket = kBucketCreate(currRing);
1413    int len=pLength(p[j]);
1414    kBucketInit(buf[j].bucket,buf[j].p,len);
1415  }
1416  omfree(p);
1417  qsort(buf,i,sizeof(red_object),red_object_better_gen);
1418//    Print("\ncurr_deg:%i\n",curr_deg);
1419  Print("M[%i, ",i); 
1420  multi_reduction(buf, i, c);
1421  Print("%i]",i); 
1422 //  for(j=0;j<i;j++){
1423//     if(buf[j].p==NULL) PrintS("\n ZERO ALERT \n");
1424//     int z;
1425//      for(z=0;z<j;z++){
1426//       if (pLmEqual(buf[z].p, buf[j].p))
1427//      PrintS("\n Critical Warning!!!! \n");
1428     
1429//     }
1430//   }
1431  int* ibuf=(int*) omalloc(i*sizeof(int));
1432  sorted_pair_node*** sbuf=(sorted_pair_node***) omalloc(i*sizeof(sorted_pair_node**));
1433  for(j=0;j<i;j++){
1434 
1435 
1436    int len;
1437    poly p;
1438    buf[j].flatten();
1439    kBucketClear(buf[j].bucket,&p, &len);
1440    kBucketDestroy(&buf[j].bucket);
1441    // delete buf[j];
1442    //remember to free res here
1443    p=redTailShort(p, c->strat);
1444    sbuf[j]=add_to_basis(p,-1,-1,c,ibuf+j);
1445   
1446  }
1447  int sum=0;
1448  for(j=0;j<i;j++){
1449    sum+=ibuf[j];
1450  }
1451  sorted_pair_node** big_sbuf=(sorted_pair_node**) omalloc(sum*sizeof(sorted_pair_node*));
1452  int partsum=0;
1453  for(j=0;j<i;j++){
1454    memmove(big_sbuf+partsum, sbuf[j],ibuf[j]*sizeof(sorted_pair_node*));
1455    omfree(sbuf[j]);
1456    partsum+=ibuf[j];
1457  }
1458
1459  qsort(big_sbuf,sum,sizeof(sorted_pair_node*),pair_better_gen2);
1460  c->apairs=merge(c->apairs,c->pair_top+1,big_sbuf,sum,c);
1461  c->pair_top+=sum;
1462  clean_top_of_pair_list(c);
1463  omfree(big_sbuf);
1464  omfree(sbuf);
1465  omfree(ibuf);
1466  omfree(buf);
1467  return;
1468}
1469
1470static poly redNF (poly h,kStrategy strat, int &len)
1471{
1472  len=0;
1473  if (h==NULL) return NULL;
1474  int j;
1475
1476  len=pLength(h);
1477  if (0 > strat->sl)
1478  {
1479    return h;
1480  }
1481  LObject P(h);
1482  P.SetShortExpVector();
1483  P.bucket = kBucketCreate(currRing);
1484  kBucketInit(P.bucket,P.p,len /*pLength(P.p)*/);
1485  //int max_pos=simple_posInS(strat,P.p);
1486  loop
1487    {
1488      j=kFindDivisibleByInS(strat->S,strat->sevS,strat->sl,&P);
1489      if (j>=0)
1490      {
1491        nNormalize(pGetCoeff(P.p));
1492#ifdef KDEBUG
1493        if (TEST_OPT_DEBUG)
1494        {
1495          PrintS("red:");
1496          wrp(h);
1497          PrintS(" with ");
1498          wrp(strat->S[j]);
1499        }
1500#endif
1501        number coef=kBucketPolyRed(P.bucket,strat->S[j],
1502                                   strat->lenS[j]/*pLength(strat->S[j])*/,
1503                                   strat->kNoether);
1504        nDelete(&coef);
1505        h = kBucketGetLm(P.bucket);
1506        if (h==NULL) return NULL;
1507        P.p=h;
1508        P.t_p=NULL;
1509        P.SetShortExpVector();
1510#ifdef KDEBUG
1511        if (TEST_OPT_DEBUG)
1512        {
1513          PrintS("\nto:");
1514          wrp(h);
1515          PrintLn();
1516        }
1517#endif
1518      }
1519      else
1520      {
1521        kBucketClear(P.bucket,&(P.p),&len);
1522        kBucketDestroy(&P.bucket);
1523        pNormalize(P.p);
1524        return P.p;
1525      }
1526    }
1527}
1528#endif
1529#ifdef REDTAIL_S
1530
1531static poly redNFTail (poly h,const int sl,kStrategy strat, int len)
1532{
1533  if (h==NULL) return NULL;
1534  pTest(h);
1535  if (0 > sl)
1536    return h;
1537  if (pNext(h)==NULL) return h;
1538
1539  int j;
1540  poly res=h;
1541  poly act=res;
1542  LObject P(pNext(h));
1543  pNext(res)=NULL;
1544  P.bucket = kBucketCreate(currRing);
1545  len--;
1546  h=P.p;
1547  if (len <=0) len=pLength(h);
1548  kBucketInit(P.bucket,h /*P.p*/,len /*pLength(P.p)*/);
1549  pTest(h);
1550  loop
1551    {
1552      P.p=h;
1553      P.t_p=NULL;
1554      P.SetShortExpVector();
1555      loop
1556        {
1557          j=kFindDivisibleByInS(strat->S,strat->sevS,sl,&P);
1558          if (j>=0)
1559          {
1560#ifdef REDTAIL_PROT
1561            PrintS("r");
1562#endif
1563            nNormalize(pGetCoeff(P.p));
1564#ifdef KDEBUG
1565            if (TEST_OPT_DEBUG)
1566            {
1567              PrintS("red tail:");
1568              wrp(h);
1569              PrintS(" with ");
1570              wrp(strat->S[j]);
1571            }
1572#endif
1573            number coef;
1574            pTest(strat->S[j]);
1575            coef=kBucketPolyRed(P.bucket,strat->S[j],
1576                                strat->lenS[j]/*pLength(strat->S[j])*/,strat->kNoether);
1577            pMult_nn(res,coef);
1578            nDelete(&coef);
1579            h = kBucketGetLm(P.bucket);
1580            pTest(h);
1581            if (h==NULL)
1582            {
1583#ifdef REDTAIL_PROT
1584              PrintS(" ");
1585#endif
1586              return res;
1587            }
1588            pTest(h);
1589            P.p=h;
1590            P.t_p=NULL;
1591            P.SetShortExpVector();
1592#ifdef KDEBUG
1593            if (TEST_OPT_DEBUG)
1594            {
1595              PrintS("\nto tail:");
1596              wrp(h);
1597              PrintLn();
1598            }
1599#endif
1600          }
1601          else
1602          {
1603#ifdef REDTAIL_PROT
1604            PrintS("n");
1605#endif
1606            break;
1607          }
1608        } /* end loop current mon */
1609      poly tmp=pHead(h /*kBucketGetLm(P.bucket)*/);
1610      act->next=tmp;pIter(act);
1611      poly tmp2=pHead(h);
1612      pNeg(tmp2);
1613      int ltmp2=1;
1614      pTest(tmp2);
1615      kBucket_Add_q(P.bucket, tmp2, &ltmp2);
1616
1617      h = kBucketGetLm(P.bucket);
1618      if (h==NULL)
1619      {
1620#ifdef REDTAIL_PROT
1621        PrintS(" ");
1622#endif
1623        return res;
1624      }
1625      pTest(h);
1626    }
1627}
1628#endif
1629
1630static void do_this_spoly_stuff(int i,int j,calc_dat* c){
1631  poly f=c->S->m[i];
1632  poly g=c->S->m[j];
1633  poly h=ksOldCreateSpoly(f, g, NULL, c->r);
1634  poly hr=NULL;
1635#ifdef FULLREDUCTIONS
1636  if (h!=NULL)
1637  {
1638    int len;
1639
1640//    hr=redNF2(h,c,len);
1641//      hr=redNF(h,c->strat,len);
1642
1643    if (hr!=NULL)
1644#ifdef REDTAIL_S
1645      hr = redNFTail(hr,c->strat->sl,c->strat,len);
1646#else
1647    hr = redtailBba(hr,c->strat->sl,c->strat);
1648#endif
1649
1650  }
1651#else
1652  if (h!=NULL)
1653  {
1654    int len;
1655//    hr=redNF2(h,c,len);
1656  }
1657#endif
1658  c->normal_forms++;
1659  if (hr==NULL)
1660  {
1661    notice_miss(i, j, c);
1662
1663  }
1664  else
1665  {
1666
1667#ifdef HEAD_BIN
1668    hr=p_MoveHead(hr,c->HeadBin);
1669#endif
1670    add_to_basis(hr, i,j,c);
1671  }
1672}
1673//try to fill, return FALSE iff queue is empty
1674
1675static int poly_crit(const void* ap1, const void* ap2){
1676  poly p1,p2;
1677  p1=*((poly*) ap1);
1678  p2=*((poly*)ap2);
1679
1680  int c=pLmCmp(p1,p2);
1681  if (c !=0) return c;
1682  int l1=pLength(p1);
1683  int l2=pLength(p2);
1684  if (l1<l2) return -1;
1685  if (l1>l2) return 1;
1686  return 0;
1687}
1688ideal t_rep_gb(ring r,ideal arg_I, ideal arg_debug_Ideal){
1689  debug_Ideal=arg_debug_Ideal;
1690  if (debug_Ideal) PrintS("DebugIdeal received\n");
1691   Print("Idelems %i \n----------\n",IDELEMS(arg_I));
1692  ideal I=idCompactify(arg_I);
1693  qsort(I->m,IDELEMS(I),sizeof(poly),poly_crit);
1694  Print("Idelems %i \n----------\n",IDELEMS(I));
1695  calc_dat* c=(calc_dat*) omalloc(sizeof(calc_dat));
1696  c->r=currRing;
1697  void* h;
1698  poly hp;
1699  int i,j;
1700  c->easy_product_crit=0;
1701  c->extended_product_crit=0;
1702  c->is_char0=(rChar()==0);
1703  c->reduction_steps=0;
1704  c->last_index=-1;
1705
1706
1707
1708  c->Rcounter=0;
1709
1710  c->soon_free=NULL;
1711
1712
1713  c->normal_forms=0;
1714  c->current_degree=1;
1715 
1716  c->max_pairs=5*I->idelems();
1717 
1718  c->apairs=(sorted_pair_node**) omalloc(sizeof(sorted_pair_node*)*c->max_pairs);
1719  c->pair_top=-1;
1720  int n=I->idelems();
1721  for (i=0;i<n;i++){
1722    wrp(I->m[i]);
1723    PrintS("\n");
1724  }
1725    i=0;
1726  c->n=0;
1727  c->T_deg=(int*) omalloc(n*sizeof(int));
1728 
1729#ifdef HEAD_BIN
1730  c->HeadBin=omGetSpecBin(POLYSIZE + (currRing->ExpL_Size)*sizeof(long));
1731#endif
1732  /* omUnGetSpecBin(&(c->HeadBin)); */
1733  h=omalloc(n*sizeof(char*));
1734  c->states=(char**) h;
1735  h=omalloc(n*sizeof(int));
1736  c->lengths=(int*) h;
1737  h=omalloc(n*sizeof(int));
1738        c->gcd_of_terms=(poly*) omalloc(n*sizeof(poly));
1739  c->rep=(int*) h;
1740  c->short_Exps=(long*) omalloc(n*sizeof(long));
1741  c->S=idInit(n,1);
1742  c->strat=new skStrategy;
1743  c->strat->syzComp = 0;
1744  initBuchMoraCrit(c->strat);
1745  initBuchMoraPos(c->strat);
1746  c->strat->initEcart = initEcartBBA;
1747  c->strat->enterS = enterSBba;
1748  c->strat->sl = -1;
1749  i=n;
1750  /* initS(c->S,NULL,c->strat); */
1751/* intS start: */
1752  i=((i+IDELEMS(c->S)+15)/16)*16;
1753  c->strat->ecartS=(intset)omAlloc(i*sizeof(int)); /*initec(i);*/
1754  c->strat->sevS=(unsigned long*)omAlloc0(i*sizeof(unsigned long));
1755  /*initsevS(i);*/
1756  c->strat->S_2_R=(int*)omAlloc0(i*sizeof(int));/*initS_2_R(i);*/
1757  c->strat->fromQ=NULL;
1758  c->strat->Shdl=idInit(1,1);
1759  c->strat->S=c->strat->Shdl->m;
1760  c->strat->lenS=(int*)omAlloc0(i*sizeof(int));
1761  if(c->is_char0)
1762    c->strat->lenSw=(int*)omAlloc0(i*sizeof(int));
1763  else
1764    c->strat->lenSw=NULL;
1765  sorted_pair_node* si;
1766  assume(n>0);
1767  add_to_basis(I->m[0],-1,-1,c);
1768
1769  assume(c->strat->sl==c->strat->Shdl->idelems()-1);
1770
1771  for (i=1;i<n;i++)//the 1 is wanted, because first element is added to basis
1772   {
1773//     add_to_basis(I->m[i],-1,-1,c);
1774     si=(sorted_pair_node*) omalloc(sizeof(sorted_pair_node));
1775      si->i=-1;
1776      si->j=-1;
1777      si->expected_length=pLength(I->m[i]);
1778      si->deg=pTotaldegree(I->m[i]);
1779      si->lcm_of_lm=I->m[i];
1780
1781//      c->apairs[n-1-i]=si;
1782      c->apairs[n-i-1]=si;
1783      ++(c->pair_top);
1784   }
1785 
1786
1787  while(c->pair_top>=0)
1788    go_on(c);
1789
1790  for(int z=0;z<c->n;z++){
1791    omfree(c->states[z]);
1792  }
1793  omfree(c->states);
1794  omfree(c->lengths);
1795
1796
1797  omfree(c->short_Exps);
1798  omfree(c->T_deg);
1799
1800     omFree(c->strat->ecartS);
1801     omFree(c->strat->sevS);
1802//   /*initsevS(i);*/
1803   omFree(c->strat->S_2_R);
1804   
1805
1806  omFree(c->strat->lenS);
1807
1808   if(c->strat->lenSw)  omFree(c->strat->lenSw);
1809
1810
1811
1812
1813  for(i=0;i<c->n;i++){
1814    if(c->gcd_of_terms[i])
1815      pDelete(&(c->gcd_of_terms[i]));
1816  }
1817  omfree(c->gcd_of_terms);
1818
1819  omfree(c->apairs);
1820  printf("calculated %d NFs\n",c->normal_forms);
1821  printf("applied %i product crit, %i extended_product crit \n", c->easy_product_crit, c->extended_product_crit);
1822  int deleted_form_c_s=0;
1823
1824  for(i=0;i<c->n;i++){
1825    if (c->rep[i]!=i){
1826      for(j=0;j<=c->strat->sl;j++){
1827        if(c->strat->S[j]==c->S->m[i]){
1828          c->strat->S[j]=NULL;
1829          break;
1830        }
1831      }
1832      PrintS("R_delete");
1833      pDelete(&c->S->m[i]);
1834    }
1835  }
1836  for(i=0;i<=c->strat->sl;i++){
1837    if (!c->strat->S[i]) continue;
1838    BOOLEAN found=FALSE;
1839    for(j=0;j<c->n;j++){
1840      if (c->S->m[j]==c->strat->S[i]){
1841        found=TRUE;
1842        break;
1843      }
1844    }
1845    if(!found) pDelete(&c->strat->S[i]);
1846  }
1847  omfree(c->rep);
1848  I=c->S;
1849 
1850  IDELEMS(I)=c->n;
1851
1852  idSkipZeroes(c->S);
1853  for(i=0;i<=c->strat->sl;i++)
1854    c->strat->S[i]=NULL;
1855  id_Delete(&c->strat->Shdl,c->r);
1856  delete c->strat;
1857  omfree(c);
1858
1859  return(I);
1860}
1861static void now_t_rep(const int & arg_i, const int & arg_j, calc_dat* c){
1862  int i,j;
1863  if (arg_i==arg_j){
1864    return;
1865  }
1866  if (arg_i>arg_j){
1867    i=arg_j;
1868    j=arg_i;
1869  } else {
1870    i=arg_i;
1871    j=arg_j;
1872  }
1873  c->states[j][i]=HASTREP;
1874}
1875static void soon_t_rep(const int& arg_i, const int& arg_j, calc_dat* c)
1876{
1877  assume(0<=arg_i);
1878  assume(0<=arg_j);
1879  assume(arg_i<c->n);
1880  assume(arg_j<c->n);
1881  int i,j;
1882  if (arg_i==arg_j){
1883    return;
1884  }
1885  if (arg_i>arg_j){
1886    i=arg_j;
1887    j=arg_i;
1888  } else {
1889    i=arg_i;
1890    j=arg_j;
1891  }
1892  if (!
1893      (c->states[j][i]==HASTREP))
1894    c->states[j][i]=SOONTREP;
1895}
1896static BOOLEAN has_t_rep(const int & arg_i, const  int & arg_j, calc_dat* state){
1897  assume(0<=arg_i);
1898  assume(0<=arg_j);
1899  assume(arg_i<state->n);
1900  assume(arg_j<state->n);
1901  if (arg_i==arg_j)
1902  {
1903    return (TRUE);
1904  }
1905  if (arg_i>arg_j)
1906  {
1907    return (state->states[arg_i][arg_j]==HASTREP);
1908  } else
1909  {
1910    return (state->states[arg_j][arg_i]==HASTREP);
1911  }
1912}
1913static int pLcmDeg(poly a, poly b)
1914{
1915  int i;
1916  int n=0;
1917  for (i=pVariables; i; i--)
1918  {
1919    n+=max( pGetExp(a,i), pGetExp(b,i));
1920  }
1921  return n;
1922
1923}
1924
1925
1926
1927static void shorten_tails(calc_dat* c, poly monom)
1928{
1929  return;
1930// BOOLEAN corr=lenS_correct(c->strat);
1931  for(int i=0;i<c->n;i++)
1932  {
1933    //enter tail
1934    if (c->rep[i]!=i) continue;
1935    if (c->S->m[i]==NULL) continue;
1936    poly tail=c->S->m[i]->next;
1937    poly prev=c->S->m[i];
1938    BOOLEAN did_something=FALSE;
1939    while((tail!=NULL)&& (pLmCmp(tail, monom)>=0))
1940    {
1941      if (p_LmDivisibleBy(monom,tail,c->r))
1942      {
1943        did_something=TRUE;
1944        prev->next=tail->next;
1945        tail->next=NULL;
1946        p_Delete(& tail,c->r);
1947        tail=prev;
1948        //PrintS("Shortened");
1949        c->lengths[i]--;
1950      }
1951      prev=tail;
1952      tail=tail->next;
1953    }
1954    if (did_something)
1955    {
1956      int new_pos;
1957      int q;
1958      q=quality(c->S->m[i],c->lengths[i],c);
1959      new_pos=simple_posInS(c->strat,c->S->m[i],q,c->is_char0);
1960
1961      int old_pos=-1;
1962      //assume new_pos<old_pos
1963      for (int z=0;z<=c->strat->sl;z++)
1964      {
1965        if (c->strat->S[z]==c->S->m[i])
1966        {
1967          old_pos=z;
1968          break;
1969        }
1970      }
1971      if (old_pos== -1)
1972        for (int z=new_pos-1;z>=0;z--)
1973        {
1974          if (c->strat->S[z]==c->S->m[i])
1975          {
1976            old_pos=z;
1977            break;
1978          }
1979        }
1980      assume(old_pos>=0);
1981      assume(new_pos<=old_pos);
1982      assume(pLength(c->strat->S[old_pos])==c->lengths[i]);
1983      c->strat->lenS[old_pos]=c->lengths[i];
1984      if (c->strat->lenSw)
1985        c->strat->lenSw[old_pos]=q;
1986
1987      if (new_pos<old_pos)
1988        move_forward_in_S(old_pos,new_pos,c->strat, c->is_char0);
1989
1990      length_one_crit(c,i,c->lengths[i]);
1991    }
1992  }
1993}
1994static sorted_pair_node* pop_pair(calc_dat* c){
1995  clean_top_of_pair_list(c);
1996
1997  if(c->pair_top<0) return NULL;
1998  else return (c->apairs[c->pair_top--]);
1999}
2000static sorted_pair_node* top_pair(calc_dat* c){
2001  super_clean_top_of_pair_list(c);//yeah, I know, it's odd that I use a different proc here
2002
2003  if(c->pair_top<0) return NULL;
2004  else return (c->apairs[c->pair_top]);
2005}
2006static sorted_pair_node* quick_pop_pair(calc_dat* c){
2007  if(c->pair_top<0) return NULL;
2008  else return (c->apairs[c->pair_top--]);
2009}
2010static BOOLEAN no_pairs(calc_dat* c){
2011  clean_top_of_pair_list(c);
2012  return (c->pair_top==-1);
2013}
2014
2015
2016static void super_clean_top_of_pair_list(calc_dat* c){
2017  while((c->pair_top>=0) && (c->apairs[c->pair_top]->i>=0) && (good_has_t_rep(c->apairs[c->pair_top]->j, c->apairs[c->pair_top]->i,c))){
2018
2019    free_sorted_pair_node(c->apairs[c->pair_top],c->r);
2020    c->pair_top--;
2021
2022  }
2023}
2024static void clean_top_of_pair_list(calc_dat* c){
2025  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))){
2026
2027    free_sorted_pair_node(c->apairs[c->pair_top],c->r);
2028    c->pair_top--;
2029
2030  }
2031}
2032static BOOLEAN state_is(calc_state state, const int & arg_i, const  int & arg_j, calc_dat* c){
2033  assume(0<=arg_i);
2034  assume(0<=arg_j);
2035  assume(arg_i<c->n);
2036  assume(arg_j<c->n);
2037  if (arg_i==arg_j)
2038  {
2039    return (TRUE);
2040  }
2041  if (arg_i>arg_j)
2042  {
2043    return (c->states[arg_i][arg_j]==state);
2044  }
2045  else return(c->states[arg_j][arg_i]==state);
2046}
2047static void free_sorted_pair_node(sorted_pair_node* s, ring r){
2048  if (s->i>=0)
2049    p_Delete(&s->lcm_of_lm,r);
2050  omfree(s);
2051}
2052static BOOLEAN pair_better(sorted_pair_node* a,sorted_pair_node* b, calc_dat* c){
2053  if (a->deg<b->deg) return TRUE;
2054  if (a->deg>b->deg) return FALSE;
2055
2056  if (a->expected_length<b->expected_length) return TRUE;
2057  if (a->expected_length>b->expected_length) return FALSE;
2058  int comp=pLmCmp(a->lcm_of_lm, b->lcm_of_lm);
2059  if (comp==1) return FALSE;
2060  if (-1==comp) return TRUE;
2061  if (a->i<b->i) return TRUE;
2062  if (a->j<b->j) return TRUE;
2063  return FALSE;
2064}
2065
2066static int pair_better_gen(const void* ap,const void* bp){
2067
2068  sorted_pair_node* a=*((sorted_pair_node**)ap);
2069  sorted_pair_node* b=*((sorted_pair_node**)bp);
2070  assume(a->i>a->j);
2071  assume(b->i>b->j);
2072  if (a->deg<b->deg) return -1;
2073  if (a->deg>b->deg) return 1;
2074
2075
2076  if (a->expected_length<b->expected_length) return -1;
2077  if (a->expected_length>b->expected_length) return 1;
2078 int comp=pLmCmp(a->lcm_of_lm, b->lcm_of_lm);
2079 
2080  if (comp==1) return 1;
2081  if (-1==comp) return -1;
2082  if (a->i<b->i) return -1;
2083  if (a->j<b->j) return -1;
2084  return 0;
2085}
2086
2087
2088static poly gcd_of_terms(poly p, ring r){
2089  int max_g_0=0;
2090  assume(p!=NULL);
2091  int i;
2092  poly m=pOne();
2093  poly t;
2094  for (i=pVariables; i; i--)
2095  {
2096      pSetExp(m,i, pGetExp(p,i));
2097      if (max_g_0==0)
2098        if (pGetExp(m,i)>0)
2099          max_g_0=i;
2100  }
2101 
2102  t=p->next;
2103  while (t!=NULL){
2104   
2105    if (max_g_0==0) break;
2106    for (i=pVariables; i; i--)
2107    {
2108      pSetExp(m,i, min(pGetExp(t,i),pGetExp(m,i)));
2109      if (max_g_0==i)
2110        if (pGetExp(m,i)==0)
2111          max_g_0=0;
2112      if ((max_g_0==0) && (pGetExp(m,i)>0)){
2113        max_g_0=i;
2114      }
2115    }
2116                t=t->next;
2117  }
2118        for (i=pVariables;i;i--)
2119        {
2120                if(pGetExp(m,i)>0)
2121                        return m;
2122  }
2123        pDelete(&m);
2124  return NULL;
2125}
2126static BOOLEAN pHasNotCFExtended(poly p1, poly p2, poly m)
2127{
2128
2129  if (pGetComp(p1) > 0 || pGetComp(p2) > 0)
2130    return FALSE;
2131  int i = 1;
2132  loop
2133  {
2134    if ((pGetExp(p1, i)-pGetExp(m,i) >0) && (pGetExp(p2, i) -pGetExp(m,i)> 0))   return FALSE;
2135    if (i == pVariables)                                return TRUE;
2136    i++;
2137  }
2138}
2139
2140
2141//for impl reasons may return false if the the normal product criterion matches
2142static BOOLEAN extended_product_criterion(poly p1, poly gcd1, poly p2, poly gcd2, calc_dat* c){
2143  if(gcd1==NULL) return FALSE;
2144        if(gcd2==NULL) return FALSE;
2145        gcd1->next=gcd2; //may ordered incorrect
2146        poly m=gcd_of_terms(gcd1,c->r);
2147        gcd1->next=NULL;
2148        if (m==NULL) return FALSE;
2149
2150        BOOLEAN erg=pHasNotCFExtended(p1,p2,m);
2151        pDelete(&m);
2152        return erg;
2153}
2154static poly kBucketGcd(kBucket* b, ring r)
2155{
2156  int s=0;
2157  int i;
2158  poly m, n;
2159  BOOLEAN initialized=FALSE;
2160  for (i=MAX_BUCKET-1;i>=0;i--)
2161  { 
2162    if (b->buckets[i]!=NULL){
2163      if (!initialized){
2164        m=gcd_of_terms(b->buckets[i],r);
2165        initialized=TRUE;
2166        if (m==NULL) return NULL;
2167      }
2168      else
2169        {
2170          n=gcd_of_terms(b->buckets[i],r);
2171          if (n==NULL) {
2172            pDelete(&m);
2173            return NULL;   
2174          }
2175          n->next=m;
2176          poly t=gcd_of_terms(n,r);
2177          n->next=NULL;
2178          pDelete(&m);
2179          pDelete(&n);
2180          m=t;
2181          if (m==NULL) return NULL;
2182         
2183        }
2184    }
2185  }
2186  return m;
2187}
2188
2189
2190
2191static int guess_quality(const red_object & p, calc_dat* c){
2192  //looks only on bucket
2193  if (c->is_char0) return kSBucketLength(p.bucket);
2194  return (bucket_guess(p.bucket));
2195}
2196static int quality_of_pos_in_strat_S(int pos, calc_dat* c){
2197  if (c->is_char0) return c->strat->lenSw[pos];
2198  return c->strat->lenS[pos];
2199}
2200static int quality(poly p, int len, calc_dat* c){
2201  if (c->is_char0) return pSLength(p,len);
2202  return pLength(p);
2203}
2204static void multi_reduction_lls_trick(red_object* los, int losl,calc_dat* c,find_erg & erg){
2205  erg.expand=NULL;
2206  BOOLEAN swap_roles; //from reduce_by, to_reduce_u if fromS
2207  if(erg.fromS){
2208    if(pLmEqual(c->strat->S[erg.reduce_by],los[erg.to_reduce_u].p))
2209    {
2210      int i;
2211      int quality_a=quality_of_pos_in_strat_S(erg.reduce_by,c);
2212      int best=erg.to_reduce_u+1;
2213/*
2214      for (i=erg.to_reduce_u;i>=erg.to_reduce_l;i--){
2215        int qc=los[i].guess_quality(c);
2216        if (qc<quality_a){
2217          best=i;
2218          quality_a=qc;
2219        }
2220      }
2221      if(best!=erg.to_reduce_u+1){*/
2222      int qc;
2223      best=find_best(los,erg.to_reduce_l,erg.to_reduce_u,qc,c);
2224      assume(qc==los[best].guess_quality(c));
2225      if(qc<quality_a){
2226        red_object h=los[erg.to_reduce_u];
2227        los[erg.to_reduce_u]=los[best];
2228        los[best]=h;
2229        swap_roles=TRUE;
2230      }
2231      else{
2232       
2233        swap_roles=FALSE;
2234      }
2235 
2236    }
2237      else
2238    {
2239      if (erg.to_reduce_u>erg.to_reduce_l){
2240
2241        int i;
2242        int quality_a=quality_of_pos_in_strat_S(erg.reduce_by,c);
2243        int best=erg.to_reduce_u+1;
2244        int qc;
2245        best=find_best(los,erg.to_reduce_l,erg.to_reduce_u,qc,c);
2246        assume(qc==los[best].guess_quality(c));
2247        if(qc<quality_a){
2248          //(best!=erg.to_reduce_u+1){
2249          red_object h=los[erg.to_reduce_l];
2250          los[erg.to_reduce_l]=los[best];
2251          los[best]=h;
2252          erg.reduce_by=erg.to_reduce_l;
2253          erg.fromS=FALSE;
2254          erg.to_reduce_l++;
2255         
2256        }
2257      }
2258      else 
2259      {
2260        assume(erg.to_reduce_u==erg.to_reduce_l);
2261        int quality_a=quality_of_pos_in_strat_S(erg.reduce_by,c);
2262        int qc=los[erg.to_reduce_u].guess_quality(c);
2263        if(qc<quality_a){
2264          BOOLEAN exp=FALSE;
2265          if(qc<=2)
2266            exp=TRUE;
2267          else {
2268            if (qc<quality_a/2)
2269              exp=TRUE;
2270            else
2271              if(erg.reduce_by<c->n/4)
2272                exp=TRUE;
2273          }
2274          if (exp){
2275            poly clear_into;
2276            los[erg.to_reduce_u].flatten();
2277            kBucketClear(los[erg.to_reduce_u].bucket,&clear_into,&erg.expand_length);
2278            erg.expand=pCopy(clear_into);
2279            kBucketInit(los[erg.to_reduce_u].bucket,clear_into,erg.expand_length);
2280            PrintS("e");
2281           
2282          }
2283        }
2284
2285       
2286      }
2287     
2288      swap_roles=FALSE;
2289      return;
2290      }
2291   
2292  }
2293  else{
2294    if(erg.reduce_by>erg.to_reduce_u){
2295      //then lm(rb)>= lm(tru) so =
2296      assume(erg.reduce_by==erg.to_reduce_u+1);
2297      int best=erg.reduce_by;
2298      int quality_a=los[erg.reduce_by].guess_quality(c);
2299      int qc;
2300      best=find_best(los,erg.to_reduce_l,erg.to_reduce_u,qc,c);
2301     
2302      int i;
2303      if(qc<quality_a){
2304          red_object h=los[erg.reduce_by];
2305          los[erg.reduce_by]=los[best];
2306          los[best]=h;
2307        }
2308        swap_roles=FALSE;
2309        return;
2310       
2311         
2312    }
2313    else
2314    {
2315      assume(!pLmEqual(los[erg.reduce_by].p,los[erg.to_reduce_l].p));
2316      //further assume, that reduce_by is the above all other polys
2317      //with same leading term
2318      int il=erg.reduce_by;
2319      int quality_a =los[erg.reduce_by].guess_quality(c);
2320      int qc;
2321      while((il>0) && pLmEqual(los[il-1].p,los[il].p)){
2322        il--;
2323        qc=los[il].guess_quality(c);
2324        if (qc<quality_a){
2325          quality_a=qc;
2326          erg.reduce_by=il;
2327        }
2328      }
2329      swap_roles=FALSE;
2330    }
2331 
2332  }
2333  if(swap_roles){
2334    PrintS("b");
2335    poly clear_into;
2336    int dummy_len;
2337    int new_length;
2338    int bp=erg.to_reduce_u;//bucket_positon
2339    //kBucketClear(los[bp].bucket,&clear_into,&new_length);
2340    new_length=los[bp].clear_to_poly();
2341    clear_into=los[bp].p;
2342    poly p=c->strat->S[erg.reduce_by];
2343    int j=erg.reduce_by;
2344    int old_length=c->strat->lenS[j];// in view of S
2345    los[bp].p=p;
2346    kBucketInit(los[bp].bucket,p,old_length);
2347    int qal=quality(clear_into,new_length,c);
2348    int pos_in_c=-1;   
2349    int z;
2350    int new_pos;
2351    new_pos=simple_posInS(c->strat,clear_into,qal,c->is_char0);
2352    assume(new_pos<=j);
2353    for (z=c->n;z;z--)
2354    {
2355      if(p==c->S->m[z-1])
2356      {
2357        pos_in_c=z-1;
2358        break;
2359      }
2360    }
2361    if(pos_in_c>=0)
2362    {
2363      c->S->m[pos_in_c]=clear_into;
2364      c->lengths[pos_in_c]=new_length;
2365      c_S_element_changed_hook(pos_in_c,c);
2366    }
2367    c->strat->S[j]=clear_into;
2368    c->strat->lenS[j]=new_length;
2369    if(c->strat->lenSw)
2370      c->strat->lenS[j]=qal;
2371    if(c->is_char0)
2372    {
2373      pContent(clear_into);
2374      pCleardenom(clear_into);
2375    }
2376  else                     
2377    pNorm(clear_into);
2378    if (new_pos<j){
2379      move_forward_in_S(j,new_pos,c->strat,c->is_char0);
2380      erg.reduce_by=new_pos;
2381    }
2382  }
2383}
2384static void multi_reduction_find(red_object* los, int losl,calc_dat* c,int startf,find_erg & erg){
2385  kStrategy strat=c->strat;
2386
2387  assume(startf<=losl);
2388  assume((startf==losl-1)||(pLmCmp(los[startf].p,los[startf+1].p)==-1));
2389  int i=startf;
2390 
2391  int j;
2392  while(i>=0){
2393    assume((i==losl-1)||(pLmCmp(los[i].p,los[i+1].p)<=0));
2394    j=kFindDivisibleByInS_easy(strat,los[i]);
2395    if(j>=0){
2396     
2397      erg.to_reduce_u=i;
2398      erg.reduce_by=j;
2399      erg.fromS=TRUE;
2400      int i2;
2401      for(i2=i-1;i2>=0;i2--){
2402        if(!pLmEqual(los[i].p,los[i2].p))
2403          break;
2404      }
2405      erg.to_reduce_l=i2+1;
2406      assume((i==losl-1)||(pLmCmp(los[i].p,los[i+1].p)==-1));
2407      return;
2408    }
2409    if (j<0){
2410     
2411      //not reduceable, try to use this for reducing higher terms
2412      int i2;
2413      i2=i;
2414      while((i2>0)&&(pLmEqual(los[i].p,los[i2-1].p)))
2415        i2--;
2416      if(i2!=i){
2417       
2418        erg.to_reduce_u=i-1;
2419        erg.to_reduce_l=i2;
2420        erg.reduce_by=i;
2421        erg.fromS=FALSE;
2422        assume((i==losl-1)||(pLmCmp(los[i].p,los[i+1].p)==-1));
2423        return;
2424      }
2425 
2426      for (i2=i+1;i2<losl;i2++){
2427        if (p_LmShortDivisibleBy(los[i].p,los[i].sev,los[i2].p,~los[i2].sev,
2428                                c->r)){
2429          int i3=i2;
2430          while((i3+1<losl) && (pLmEqual(los[i2].p, los[i3+1].p)))
2431            i3++;
2432          erg.to_reduce_u=i3;
2433          erg.to_reduce_l=i2;
2434          erg.reduce_by=i;
2435          erg.fromS=FALSE;
2436          assume((i==losl-1)||(pLmCmp(los[i].p,los[i+1].p)==-1));
2437          return;
2438        }
2439//      else {assume(!p_LmDivisibleBy(los[i].p, los[i2].p,c->r));}
2440      }
2441
2442      i--;
2443    }
2444  }
2445  erg.reduce_by=-1;//error code
2446  return;
2447}
2448
2449 //  nicht reduzierbare eintraege in ergebnisliste schreiben
2450//   nullen loeschen
2451//   while(finde_groessten leitterm reduzierbar(c,erg)){
2452 
2453static int multi_reduction_clear_zeroes(red_object* los, int  losl, int l, int u)
2454{
2455
2456
2457  int deleted=0;
2458  int  i=l;
2459  int last=-1;
2460  while(i<=u)
2461  {
2462   
2463    if(los[i].p==NULL){
2464      kBucketDestroy(&los[i].bucket);
2465//      delete los[i];//here we assume los are constructed with new
2466      //destroy resources, must be added here   
2467     if (last>=0)
2468     {
2469       memmove(los+(int)(last+1-deleted),los+(last+1),sizeof(red_object)*(i-1-last));
2470     }
2471     last=i;
2472     deleted++;
2473    }
2474    i++;
2475  }
2476  if((last>=0)&&(last!=losl-1))
2477      memmove(los+(int)(last+1-deleted),los+last+1,sizeof(red_object)*(losl-1-last));
2478  return deleted;
2479 
2480}
2481
2482static void sort_region_down(red_object* los, int l, int u, calc_dat* c)
2483{
2484  qsort(los+l,u-l+1,sizeof(red_object),red_object_better_gen);
2485  int i;
2486
2487  for(i=l;i<=u;i++)
2488  {
2489    BOOLEAN moved=FALSE;
2490    int j;
2491    for(j=i;j;j--)
2492    {
2493      if(pLmCmp(los[j].p,los[j-1].p)==-1){
2494        red_object h=los[j];
2495        los[j]=los[j-1];
2496        los[j-1]=h;
2497        moved=TRUE;
2498      }
2499      else break;
2500    }
2501    if(!moved) return;
2502  }
2503}
2504
2505//assume that los is ordered ascending by leading term, all non zero
2506static void multi_reduction(red_object* los, int & losl, calc_dat* c)
2507{
2508 
2509  //initialize;
2510  assume(c->strat->sl>=0);
2511  assume(losl>0);
2512  int i;
2513  for(i=0;i<losl;i++){
2514    los[i].sev=pGetShortExpVector(los[i].p);
2515//SetShortExpVector();
2516    los[i].p=kBucketGetLm(los[i].bucket);
2517  }
2518
2519  kStrategy strat=c->strat;
2520  int curr_pos=losl-1;
2521
2522
2523//  nicht reduzierbare einträge in ergebnisliste schreiben
2524  // nullen loeschen
2525  while(curr_pos>=0){
2526    find_erg erg;
2527    multi_reduction_find(los, losl,c,curr_pos,erg);//last argument should be curr_pos
2528   //  PrintS("\n erg:\n");
2529//     Print("upper:%i\n",erg.to_reduce_u);
2530//     Print("lower:%i\n",erg.to_reduce_l);
2531//     Print("reduce_by:%i\n",erg.reduce_by);
2532//     Print("fromS:%i\n",erg.fromS);
2533    if(erg.reduce_by<0) break;
2534    multi_reduction_lls_trick(los,losl,c,erg);
2535    int sum=0;
2536
2537   
2538    int i;
2539    int len;
2540    multi_reduce_step(erg,los,c);
2541//     reduction_step *rs=create_reduction_step(erg, los, c);
2542//     rs->reduce(los,erg.to_reduce_l,erg.to_reduce_u);
2543//     finalize_reduction_step(rs);
2544                 
2545    int deleted=multi_reduction_clear_zeroes(los, losl, erg.to_reduce_l, erg.to_reduce_u);
2546    curr_pos=erg.to_reduce_u;
2547    losl -= deleted;
2548    curr_pos -= deleted;
2549
2550    //Print("deleted %i \n",deleted);
2551    sort_region_down(los, erg.to_reduce_l, erg.to_reduce_u-deleted, c);
2552//   sort_region_down(los, 0, losl-1, c);
2553    //  qsort(los,losl,sizeof(red_object),red_object_better_gen);
2554    if(erg.expand)
2555      add_to_reductors(c,erg.expand,erg.expand_length);
2556  }
2557  return;
2558}
2559void red_object::flatten(){
2560  if (sum!=NULL)
2561  {
2562
2563 
2564    if(kBucketGetLm(sum->ac->bucket)!=NULL){
2565      number mult_my=n_Mult(sum->c_my,sum->ac->multiplied,currRing);
2566      poly add_this;
2567      if(!nIsOne(mult_my))
2568        kBucket_Mult_n(bucket,mult_my);
2569      int len;
2570      poly clear_into;
2571      kBucketClear(sum->ac->bucket,&clear_into,&len);
2572      if(sum->ac->counter>1){
2573        add_this=pCopy(clear_into);
2574        kBucketInit(sum->ac->bucket,clear_into,len);
2575      }
2576      else
2577        add_this=clear_into;
2578      pMult_nn(add_this, sum->c_ac);
2579     
2580      kBucket_Add_q(bucket,add_this, &len);
2581       nDelete(&mult_my);
2582     
2583    }
2584    nDelete(&sum->c_ac);
2585    nDelete(&sum->c_my);
2586   
2587    sum->ac->decrease_counter();
2588    delete sum;
2589    p=kBucketGetLm(bucket);
2590    sum=NULL;
2591  }
2592  assume(p==kBucketGetLm(bucket));
2593  assume(sum==NULL);
2594}
2595void red_object::validate(){
2596  BOOLEAN flattened=FALSE;
2597  if(sum!=NULL)
2598  {
2599    poly lm=kBucketGetLm(bucket);
2600    poly lm_ac=kBucketGetLm(sum->ac->bucket);
2601    if ((lm_ac==NULL)||((lm!=NULL) && (pLmCmp(lm,lm_ac)!=-1))){
2602      flatten();
2603      assume(sum==NULL);
2604      flattened=TRUE;
2605      p=kBucketGetLm(bucket);
2606      if (p!=NULL)
2607        sev=pGetShortExpVector(p);
2608    } 
2609    else
2610    {
2611 
2612      p=lm_ac;
2613      assume(sum->ac->sev==pGetShortExpVector(p));
2614      sev=sum->ac->sev;
2615    }
2616   
2617  }
2618  else{
2619    p=kBucketGetLm(bucket);
2620    if(p)
2621    sev=pGetShortExpVector(p);
2622  }
2623  assume((sum==NULL)||(kBucketGetLm(sum->ac->bucket)!=NULL));
2624}
2625int red_object::clear_to_poly(){
2626  flatten();
2627  int l;
2628  kBucketClear(bucket,&p,&l);
2629  return l;
2630}
2631
2632 
2633
2634void red_object::adjust_coefs(number c_r, number c_ac_r){
2635  assume(this->sum!=NULL);
2636  number n1=nMult(sum->c_my, c_ac_r);
2637  number n2=nMult(sum->c_ac,c_r);
2638  nDelete(&sum->c_my);
2639  nDelete(&sum->c_ac);
2640 
2641  int ct = ksCheckCoeff(&n1, &n2);
2642  sum->c_my=n1;
2643  sum->c_ac=nNeg(n2);
2644  nDelete(&n2);
2645 
2646
2647}
2648int red_object::guess_quality(calc_dat* c){
2649    //works at the moment only for lenvar 1, because in different
2650    //case, you have to look on coefs
2651    int s=0;
2652    if (c->is_char0)
2653      s=kSBucketLength(bucket);
2654    else 
2655      s=bucket_guess(bucket);
2656    if (sum!=NULL){
2657      if (c->is_char0)
2658      s+=kSBucketLength(sum->ac->bucket);
2659    else 
2660      s+=bucket_guess(sum->ac->bucket);
2661    }
2662    return s;
2663}
2664void reduction_step::reduce(red_object* r, int l, int u){}
2665void simple_reducer::target_is_no_sum_reduce(red_object & ro){
2666  kBucketPolyRed(ro.bucket,p,
2667                 p_len,
2668                 c->strat->kNoether);
2669}
2670
2671void simple_reducer::target_is_a_sum_reduce(red_object & ro){
2672  pTest(p);
2673  kbTest(ro.bucket);
2674  kbTest(ro.sum->ac->bucket);
2675  assume(ro.sum!=NULL);
2676  assume(ro.sum->ac!=NULL);
2677  if(ro.sum->ac->last_reduction_id!=reduction_id){
2678    number n1=kBucketPolyRed(ro.sum->ac->bucket,p, p_len, c->strat->kNoether);
2679    number n2=nMult(n1,ro.sum->ac->multiplied);
2680    nDelete(&ro.sum->ac->multiplied);
2681    nDelete(&n1);
2682    ro.sum->ac->multiplied=n2;
2683    poly lm=kBucketGetLm(ro.sum->ac->bucket);
2684    if (lm)
2685      ro.sum->ac->sev=pGetShortExpVector(lm);
2686  }
2687  ro.sev=ro.sum->ac->sev;
2688}
2689void simple_reducer::reduce(red_object* r, int l, int u){
2690  this->pre_reduce(r,l,u);
2691  int i;
2692  for(i=l;i<=u;i++){
2693    if(r[i].sum==NULL)
2694      this->target_is_no_sum_reduce(r[i]);
2695
2696    else 
2697    {
2698      this->target_is_a_sum_reduce(r[i]);
2699      //reduce and adjust multiplied
2700      r[i].sum->ac->last_reduction_id=reduction_id;
2701     
2702    }
2703    //most elegant would be multimethods at this point and subclassing
2704    //red_object for sum
2705 
2706  }
2707  for(i=l;i<=u;i++){
2708    r[i].validate();
2709      }
2710}
2711reduction_step::~reduction_step(){}
2712simple_reducer::~simple_reducer(){
2713  if(fill_back!=NULL)
2714  {
2715    kBucketInit(fill_back,p,p_len);
2716  }
2717   
2718}
2719
2720void multi_reduce_step(find_erg & erg, red_object* r, calc_dat* c){
2721  static int id=0;
2722  id++;
2723//   if(id==0)
2724//     PrintS("warning: integer overflow");
2725//   BOOLEAN join=FALSE;
2726//   int i;
2727//   int sum;
2728//   for(i=erg.to_reduce_l;i<=erg.to_reduce_u;i++){
2729//     if(!r[i].sum) sum++;
2730//     if (sum>=AC_NEW_MIN) {join=TRUE;break;}
2731//   }
2732//   simple_reducer* pointer;
2733//   poly p;
2734
2735//   if (erg.fromS){
2736   
2737//     p=c->strat->S[erg.reduce_by];
2738//     if ((!join)||(p->next==NULL)){
2739     
2740//       pointer=new simple_reducer(p,c->strat->lenS[erg.reduce_by]);
2741   
2742//     }
2743//     else{
2744//       join_simple_reducer* jp;
2745//       pointer= jp=new join_simple_reducer(p,c->strat->lenS[erg.reduce_by],r[erg.to_reduce_l].p);
2746
2747
2748     
2749//     }
2750   
2751//     pointer->fill_back=NULL;
2752//   }
2753//   else
2754//   {
2755//     if(r[erg.reduce_by].sum)
2756//       kbTest(r[erg.reduce_by].sum->ac->bucket);
2757//     r[erg.reduce_by].flatten();
2758//     int len;
2759//     kBucket_pt bucket=r[erg.reduce_by].bucket;
2760//     kBucketClear(bucket,&p,&len);
2761//     if(c->is_char0)
2762//       pContent(p);
2763//     pTest(p);
2764//     if ((!join)||(p->next==NULL))
2765//       pointer=new simple_reducer(p,len);
2766//     else
2767//     {
2768//       join_simple_reducer* jp;
2769//       pointer=jp=new join_simple_reducer(p,len,r[erg.to_reduce_l].p);
2770
2771
2772
2773//     }
2774   
2775//     pointer->p=p;
2776//     pointer->p_len=len;
2777//     pointer->fill_back=bucket;
2778 
2779//   }
2780
2781//   pointer->reduction_id=id;
2782//   pointer->c=c;
2783 
2784//   return pointer;
2785  int rn=erg.reduce_by;
2786  poly red;
2787  int red_len;
2788  simple_reducer* pointer;
2789  BOOLEAN woc=FALSE;
2790  if(erg.fromS){
2791    red=c->strat->S[rn];
2792    red_len=c->strat->lenS[rn];
2793  }
2794  else
2795  {
2796    kBucketClear(r[rn].bucket,&red,&red_len);
2797  }
2798  #ifdef HANS_IDEA
2799  if(erg.to_reduce_u-erg.to_reduce_l>5){
2800    poly my=pOne();
2801    int ol=red_len;
2802    for(int i=1;i<=pVariables;i++)
2803      pSetExp(my,i,(pGetExp(r[erg.to_reduce_u].p, i)-pGetExp(red,i)));
2804   
2805    pSetm(my);
2806    number m;
2807    poly a=ppMult_mm(red,my);
2808    woc=TRUE;
2809    red_len--;
2810
2811//    a->next=redNF2(a->next,c,red_len, m,2);
2812    redTailShort(a,c->strat);
2813//    pSetCoeff(a,nMult(a->coef,m));
2814    red_len++;
2815    red_len=pLength(a);
2816    assume(red_len==pLength(a));
2817    if(!erg.fromS) kBucketInit(r[rn].bucket,red, ol);
2818    red=a;
2819    red_len=pLength(a);
2820  }
2821  #endif
2822  pointer=new simple_reducer(red,red_len,c);
2823
2824  if ((!woc) && (!erg.fromS))
2825    pointer->fill_back=r[rn].bucket;
2826  else
2827    pointer->fill_back=NULL;
2828  pointer->reduce(r,erg.to_reduce_l, erg.to_reduce_u);
2829  if(woc) pDelete(&pointer->p);
2830  delete pointer;
2831 
2832};
2833
2834void join_simple_reducer::target_is_no_sum_reduce(red_object & ro){
2835 
2836  ro.sum=new formal_sum_descriptor();
2837  ro.sum->ac=ac;
2838  ac->counter++;
2839  kBucket_pt bucket=ro.bucket;
2840  poly a1 = pNext(p), lm = kBucketExtractLm(bucket);
2841  BOOLEAN reset_vec=FALSE;
2842  number rn;
2843  assume(a1!=NULL);
2844  number an = pGetCoeff(p), bn = pGetCoeff(lm);
2845  lm->next=NULL;
2846  int ct = ksCheckCoeff(&an, &bn);
2847  ro.sum->c_ac=nNeg(bn);
2848  ro.sum->c_my=an;
2849  assume(nIsZero(nAdd(nMult(ro.sum->c_my,lm->coef),nMult(p->coef,ro.sum->c_ac) )));
2850  if (p_GetComp(p, bucket->bucket_ring) != p_GetComp(lm, bucket->bucket_ring))
2851  {
2852    p_SetCompP(a1, p_GetComp(lm, bucket->bucket_ring), bucket->bucket_ring);
2853    reset_vec = TRUE;
2854    p_SetComp(lm, p_GetComp(p, bucket->bucket_ring), bucket->bucket_ring);
2855    p_Setm(lm, bucket->bucket_ring);
2856  }
2857
2858
2859 
2860
2861  p_DeleteLm(&lm, bucket->bucket_ring);
2862  if (reset_vec) p_SetCompP(a1, 0, bucket->bucket_ring);
2863  kbTest(bucket);
2864
2865}
2866
2867  reduction_accumulator::reduction_accumulator(poly p, int p_len, poly high_to){
2868    //sev needs to be removed from interfaces,makes no sense
2869   
2870    degbound(p);
2871    degbound(p->next);
2872    degbound(high_to);
2873    poly my=pOne();
2874    counter=0;
2875   
2876    for(int i=1;i<=pVariables;i++)
2877      pSetExp(my,i,(pGetExp(high_to, i)-pGetExp(p,i)));
2878    degbound(my);
2879   
2880   
2881    pSetm(my);
2882    last_reduction_id=-1;
2883    multiplied=nInit(1);
2884    bucket=kBucketCreate(currRing);
2885    poly a=pMult_mm(pCopy(p->next),my);
2886    degbound(a);
2887    this->sev=pGetShortExpVector(a);
2888    kBucketInit(bucket, a,p_len-1);
2889    pDelete(&my);
2890  }
2891void simple_reducer:: pre_reduce(red_object* r, int l, int u){}
2892void join_simple_reducer:: pre_reduce(red_object* r, int l, int u){
2893  for(int i=l;i<=u;i++)
2894    {
2895      if (r[i].sum){
2896        if(r[i].sum->ac->counter<=15) r[i].flatten();
2897       
2898      }
2899    }
2900}
Note: See TracBrowser for help on using the repository browser.