source: git/kernel/tgb.cc @ ebe987

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