source: git/kernel/tgb.cc @ f478f5b

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