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