source: git/kernel/tgb.cc @ 47e836b

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