1 /****************************************
2 *  Computer Algebra System SINGULAR     *
3 ****************************************/
4 /*
5 *  ABSTRACT -  dimension, multiplicity, HC, kbase
6 */
7 
8 #include "kernel/mod2.h"
9 
10 #include "misc/intvec.h"
11 #include "coeffs/numbers.h"
12 
13 #include "kernel/structs.h"
14 #include "kernel/ideals.h"
15 #include "kernel/polys.h"
16 
17 #include "kernel/combinatorics/hutil.h"
18 #include "kernel/combinatorics/hilb.h"
19 #include "kernel/combinatorics/stairc.h"
20 #include "reporter/reporter.h"
21 
22 #ifdef HAVE_SHIFTBBA
23 #include <vector>
24 #include "Singular/libsingular.h"
25 #endif
26 
27 VAR int  hCo, hMu, hMu2;
28 VAR omBin indlist_bin = omGetSpecBin(sizeof(indlist));
29 
30 /*0 implementation*/
31 
32 // dimension
33 
hDimSolve(scmon pure,int Npure,scfmon rad,int Nrad,varset var,int Nvar)34 void hDimSolve(scmon pure, int Npure, scfmon rad, int Nrad,
35  varset var, int Nvar)
36 {
37   int  dn, iv, rad0, b, c, x;
38   scmon pn;
39   scfmon rn;
40   if (Nrad < 2)
41   {
42     dn = Npure + Nrad;
43     if (dn < hCo)
44       hCo = dn;
45     return;
46   }
47   if (Npure+1 >= hCo)
48     return;
49   iv = Nvar;
50   while(pure[var[iv]]) iv--;
51   hStepR(rad, Nrad, var, iv, &rad0);
52   if (rad0!=0)
53   {
54     iv--;
55     if (rad0 < Nrad)
56     {
57       pn = hGetpure(pure);
58       rn = hGetmem(Nrad, rad, radmem[iv]);
59       hDimSolve(pn, Npure + 1, rn, rad0, var, iv);
60       b = rad0;
61       c = Nrad;
62       hElimR(rn, &rad0, b, c, var, iv);
63       hPure(rn, b, &c, var, iv, pn, &x);
64       hLex2R(rn, rad0, b, c, var, iv, hwork);
65       rad0 += (c - b);
66       hDimSolve(pn, Npure + x, rn, rad0, var, iv);
67     }
68     else
69     {
70       hDimSolve(pure, Npure, rad, Nrad, var, iv);
71     }
72   }
73   else
74     hCo = Npure + 1;
75 }
76 
scDimInt(ideal S,ideal Q)77 int  scDimInt(ideal S, ideal Q)
78 {
79   id_Test(S, currRing);
80   if( Q!=NULL ) id_Test(Q, currRing);
81 
82   int  mc;
83   hexist = hInit(S, Q, &hNexist, currRing);
84   if (!hNexist)
85     return (currRing->N);
86   hwork = (scfmon)omAlloc(hNexist * sizeof(scmon));
87   hvar = (varset)omAlloc(((currRing->N) + 1) * sizeof(int));
88   hpure = (scmon)omAlloc((1 + ((currRing->N) * (currRing->N))) * sizeof(int));
89   mc = hisModule;
90   if (!mc)
91   {
92     hrad = hexist;
93     hNrad = hNexist;
94   }
95   else
96     hrad = (scfmon)omAlloc(hNexist * sizeof(scmon));
97   radmem = hCreate((currRing->N) - 1);
98   hCo = (currRing->N) + 1;
99   loop
100   {
101     if (mc)
102       hComp(hexist, hNexist, mc, hrad, &hNrad);
103     if (hNrad)
104     {
105       hNvar = (currRing->N);
106       hRadical(hrad, &hNrad, hNvar);
107       hSupp(hrad, hNrad, hvar, &hNvar);
108       if (hNvar)
109       {
110         memset(hpure, 0, ((currRing->N) + 1) * sizeof(int));
111         hPure(hrad, 0, &hNrad, hvar, hNvar, hpure, &hNpure);
112         hLexR(hrad, hNrad, hvar, hNvar);
113         hDimSolve(hpure, hNpure, hrad, hNrad, hvar, hNvar);
114       }
115     }
116     else
117     {
118       hCo = 0;
119       break;
120     }
121     mc--;
122     if (mc <= 0)
123       break;
124   }
125   hKill(radmem, (currRing->N) - 1);
126   omFreeSize((ADDRESS)hpure, (1 + ((currRing->N) * (currRing->N))) * sizeof(int));
127   omFreeSize((ADDRESS)hvar, ((currRing->N) + 1) * sizeof(int));
128   omFreeSize((ADDRESS)hwork, hNexist * sizeof(scmon));
129   hDelete(hexist, hNexist);
130   if (hisModule)
131     omFreeSize((ADDRESS)hrad, hNexist * sizeof(scmon));
132   return (currRing->N) - hCo;
133 }
134 
scDimIntRing(ideal vid,ideal Q)135 int  scDimIntRing(ideal vid, ideal Q)
136 {
137 #ifdef HAVE_RINGS
138   if (rField_is_Ring(currRing))
139   {
140     int i = idPosConstant(vid);
141     if ((i != -1) && (n_IsUnit(pGetCoeff(vid->m[i]),currRing->cf)))
142     { /* ideal v contains unit; dim = -1 */
143       return(-1);
144     }
145     ideal vv = id_Head(vid,currRing);
146     idSkipZeroes(vv);
147     i = idPosConstant(vid);
148     int d;
149     if(i == -1)
150     {
151       d = scDimInt(vv, Q);
152       if(rField_is_Z(currRing))
153         d++;
154     }
155     else
156     {
157       if(n_IsUnit(pGetCoeff(vv->m[i]),currRing->cf))
158         d = -1;
159       else
160         d = scDimInt(vv, Q);
161     }
162     //Anne's Idea for std(4,2x) = 0 bug
163     int dcurr = d;
164     for(unsigned ii=0;ii<(unsigned)IDELEMS(vv);ii++)
165     {
166       if(vv->m[ii] != NULL && !n_IsUnit(pGetCoeff(vv->m[ii]),currRing->cf))
167       {
168         ideal vc = idCopy(vv);
169         poly c = pInit();
170         pSetCoeff0(c,nCopy(pGetCoeff(vv->m[ii])));
171         idInsertPoly(vc,c);
172         idSkipZeroes(vc);
173         for(unsigned jj = 0;jj<(unsigned)IDELEMS(vc)-1;jj++)
174         {
175           if((vc->m[jj]!=NULL)
176           && (n_DivBy(pGetCoeff(vc->m[jj]),pGetCoeff(c),currRing->cf)))
177           {
178             pDelete(&vc->m[jj]);
179           }
180         }
181         idSkipZeroes(vc);
182         i = idPosConstant(vc);
183         if (i != -1) pDelete(&vc->m[i]);
184         dcurr = scDimInt(vc, Q);
185         // the following assumes the ground rings to be either zero- or one-dimensional
186         if((i==-1) && rField_is_Z(currRing))
187         {
188           // should also be activated for other euclidean domains as groundfield
189           dcurr++;
190         }
191         idDelete(&vc);
192       }
193       if(dcurr > d)
194           d = dcurr;
195     }
196     idDelete(&vv);
197     return d;
198   }
199 #endif
200   return scDimInt(vid,Q);
201 }
202 
203 // independent set
204 STATIC_VAR scmon hInd;
205 
hIndSolve(scmon pure,int Npure,scfmon rad,int Nrad,varset var,int Nvar)206 static void hIndSolve(scmon pure, int Npure, scfmon rad, int Nrad,
207  varset var, int Nvar)
208 {
209   int  dn, iv, rad0, b, c, x;
210   scmon pn;
211   scfmon rn;
212   if (Nrad < 2)
213   {
214     dn = Npure + Nrad;
215     if (dn < hCo)
216     {
217       hCo = dn;
218       for (iv=(currRing->N); iv; iv--)
219       {
220         if (pure[iv])
221           hInd[iv] = 0;
222         else
223           hInd[iv] = 1;
224       }
225       if (Nrad)
226       {
227         pn = *rad;
228         iv = Nvar;
229         loop
230         {
231           x = var[iv];
232           if (pn[x])
233           {
234             hInd[x] = 0;
235             break;
236           }
237           iv--;
238         }
239       }
240     }
241     return;
242   }
243   if (Npure+1 >= hCo)
244     return;
245   iv = Nvar;
246   while(pure[var[iv]]) iv--;
247   hStepR(rad, Nrad, var, iv, &rad0);
248   if (rad0)
249   {
250     iv--;
251     if (rad0 < Nrad)
252     {
253       pn = hGetpure(pure);
254       rn = hGetmem(Nrad, rad, radmem[iv]);
255       pn[var[iv + 1]] = 1;
256       hIndSolve(pn, Npure + 1, rn, rad0, var, iv);
257       pn[var[iv + 1]] = 0;
258       b = rad0;
259       c = Nrad;
260       hElimR(rn, &rad0, b, c, var, iv);
261       hPure(rn, b, &c, var, iv, pn, &x);
262       hLex2R(rn, rad0, b, c, var, iv, hwork);
263       rad0 += (c - b);
264       hIndSolve(pn, Npure + x, rn, rad0, var, iv);
265     }
266     else
267     {
268       hIndSolve(pure, Npure, rad, Nrad, var, iv);
269     }
270   }
271   else
272   {
273     hCo = Npure + 1;
274     for (x=(currRing->N); x; x--)
275     {
276       if (pure[x])
277         hInd[x] = 0;
278       else
279         hInd[x] = 1;
280     }
281     hInd[var[iv]] = 0;
282   }
283 }
284 
scIndIntvec(ideal S,ideal Q)285 intvec * scIndIntvec(ideal S, ideal Q)
286 {
287   id_Test(S, currRing);
288   if( Q!=NULL ) id_Test(Q, currRing);
289 
290   intvec *Set=new intvec((currRing->N));
291   int  mc,i;
292   hexist = hInit(S, Q, &hNexist, currRing);
293   if (hNexist==0)
294   {
295     for(i=0; i<(currRing->N); i++)
296       (*Set)[i]=1;
297     return Set;
298   }
299   hwork = (scfmon)omAlloc(hNexist * sizeof(scmon));
300   hvar = (varset)omAlloc(((currRing->N) + 1) * sizeof(int));
301   hpure = (scmon)omAlloc((1 + ((currRing->N) * (currRing->N))) * sizeof(int));
302   hInd = (scmon)omAlloc0((1 + (currRing->N)) * sizeof(int));
303   mc = hisModule;
304   if (mc==0)
305   {
306     hrad = hexist;
307     hNrad = hNexist;
308   }
309   else
310     hrad = (scfmon)omAlloc(hNexist * sizeof(scmon));
311   radmem = hCreate((currRing->N) - 1);
312   hCo = (currRing->N) + 1;
313   loop
314   {
315     if (mc!=0)
316       hComp(hexist, hNexist, mc, hrad, &hNrad);
317     if (hNrad!=0)
318     {
319       hNvar = (currRing->N);
320       hRadical(hrad, &hNrad, hNvar);
321       hSupp(hrad, hNrad, hvar, &hNvar);
322       if (hNvar!=0)
323       {
324         memset(hpure, 0, ((currRing->N) + 1) * sizeof(int));
325         hPure(hrad, 0, &hNrad, hvar, hNvar, hpure, &hNpure);
326         hLexR(hrad, hNrad, hvar, hNvar);
327         hIndSolve(hpure, hNpure, hrad, hNrad, hvar, hNvar);
328       }
329     }
330     else
331     {
332       hCo = 0;
333       break;
334     }
335     mc--;
336     if (mc <= 0)
337       break;
338   }
339   for(i=0; i<(currRing->N); i++)
340     (*Set)[i] = hInd[i+1];
341   hKill(radmem, (currRing->N) - 1);
342   omFreeSize((ADDRESS)hpure, (1 + ((currRing->N) * (currRing->N))) * sizeof(int));
343   omFreeSize((ADDRESS)hInd, (1 + (currRing->N)) * sizeof(int));
344   omFreeSize((ADDRESS)hvar, ((currRing->N) + 1) * sizeof(int));
345   omFreeSize((ADDRESS)hwork, hNexist * sizeof(scmon));
346   hDelete(hexist, hNexist);
347   if (hisModule)
348     omFreeSize((ADDRESS)hrad, hNexist * sizeof(scmon));
349   return Set;
350 }
351 
352 VAR indset ISet, JSet;
353 
hNotZero(scfmon rad,int Nrad,varset var,int Nvar)354 static BOOLEAN hNotZero(scfmon rad, int Nrad, varset var, int Nvar)
355 {
356   int  k1, i;
357   k1 = var[Nvar];
358   i = 0;
359   loop
360   {
361     if (rad[i][k1]==0)
362       return FALSE;
363     i++;
364     if (i == Nrad)
365       return TRUE;
366   }
367 }
368 
hIndep(scmon pure)369 static void hIndep(scmon pure)
370 {
371   int iv;
372   intvec *Set;
373 
374   Set = ISet->set = new intvec((currRing->N));
375   for (iv=(currRing->N); iv!=0 ; iv--)
376   {
377     if (pure[iv])
378       (*Set)[iv-1] = 0;
379     else
380       (*Set)[iv-1] = 1;
381   }
382   ISet = ISet->nx = (indset)omAlloc0Bin(indlist_bin);
383   hMu++;
384 }
385 
hIndMult(scmon pure,int Npure,scfmon rad,int Nrad,varset var,int Nvar)386 void hIndMult(scmon pure, int Npure, scfmon rad, int Nrad,
387  varset var, int Nvar)
388 {
389   int  dn, iv, rad0, b, c, x;
390   scmon pn;
391   scfmon rn;
392   if (Nrad < 2)
393   {
394     dn = Npure + Nrad;
395     if (dn == hCo)
396     {
397       if (Nrad==0)
398         hIndep(pure);
399       else
400       {
401         pn = *rad;
402         for (iv = Nvar; iv!=0; iv--)
403         {
404           x = var[iv];
405           if (pn[x])
406           {
407             pure[x] = 1;
408             hIndep(pure);
409             pure[x] = 0;
410           }
411         }
412       }
413     }
414     return;
415   }
416   iv = Nvar;
417   dn = Npure+1;
418   if (dn >= hCo)
419   {
420     if (dn > hCo)
421       return;
422     loop
423     {
424       if(!pure[var[iv]])
425       {
426         if(hNotZero(rad, Nrad, var, iv))
427         {
428           pure[var[iv]] = 1;
429           hIndep(pure);
430           pure[var[iv]] = 0;
431         }
432       }
433       iv--;
434       if (!iv)
435         return;
436     }
437   }
438   while(pure[var[iv]]) iv--;
439   hStepR(rad, Nrad, var, iv, &rad0);
440   iv--;
441   if (rad0 < Nrad)
442   {
443     pn = hGetpure(pure);
444     rn = hGetmem(Nrad, rad, radmem[iv]);
445     pn[var[iv + 1]] = 1;
446     hIndMult(pn, Npure + 1, rn, rad0, var, iv);
447     pn[var[iv + 1]] = 0;
448     b = rad0;
449     c = Nrad;
450     hElimR(rn, &rad0, b, c, var, iv);
451     hPure(rn, b, &c, var, iv, pn, &x);
452     hLex2R(rn, rad0, b, c, var, iv, hwork);
453     rad0 += (c - b);
454     hIndMult(pn, Npure + x, rn, rad0, var, iv);
455   }
456   else
457   {
458     hIndMult(pure, Npure, rad, Nrad, var, iv);
459   }
460 }
461 
462 /*3
463 * consider indset x := !pure
464 * (for all i) (if(sm(i) > x) return FALSE)
465 * else return TRUE
466 */
hCheck1(indset sm,scmon pure)467 static BOOLEAN hCheck1(indset sm, scmon pure)
468 {
469   int iv;
470   intvec *Set;
471   while (sm->nx != NULL)
472   {
473     Set = sm->set;
474     iv=(currRing->N);
475     loop
476     {
477       if (((*Set)[iv-1] == 0) && (pure[iv] == 0))
478         break;
479       iv--;
480       if (iv == 0)
481         return FALSE;
482     }
483     sm = sm->nx;
484   }
485   return TRUE;
486 }
487 
488 /*3
489 * consider indset x := !pure
490 * (for all i) if(x > sm(i)) delete sm(i))
491 * return (place for x)
492 */
hCheck2(indset sm,scmon pure)493 static indset hCheck2(indset sm, scmon pure)
494 {
495   int iv;
496   intvec *Set;
497   indset be, a1 = NULL;
498   while (sm->nx != NULL)
499   {
500     Set = sm->set;
501     iv=(currRing->N);
502     loop
503     {
504       if ((pure[iv] == 1) && ((*Set)[iv-1] == 1))
505         break;
506       iv--;
507       if (iv == 0)
508       {
509         if (a1 == NULL)
510         {
511           a1 = sm;
512         }
513         else
514         {
515           hMu2--;
516           be->nx = sm->nx;
517           delete Set;
518           omFreeBin((ADDRESS)sm, indlist_bin);
519           sm = be;
520         }
521         break;
522       }
523     }
524     be = sm;
525     sm = sm->nx;
526   }
527   if (a1 != NULL)
528   {
529     return a1;
530   }
531   else
532   {
533     hMu2++;
534     sm->set = new intvec((currRing->N));
535     sm->nx = (indset)omAlloc0Bin(indlist_bin);
536     return sm;
537   }
538 }
539 
540 /*2
541 *  definition x >= y
542 *      x(i) == 0 => y(i) == 0
543 *      > ex. j with x(j) == 1 and y(j) == 0
544 */
hCheckIndep(scmon pure)545 static void hCheckIndep(scmon pure)
546 {
547   intvec *Set;
548   indset res;
549   int iv;
550   if (hCheck1(ISet, pure))
551   {
552     if (hCheck1(JSet, pure))
553     {
554       res = hCheck2(JSet,pure);
555       if (res == NULL)
556         return;
557       Set = res->set;
558       for (iv=(currRing->N); iv; iv--)
559       {
560         if (pure[iv])
561           (*Set)[iv-1] = 0;
562         else
563           (*Set)[iv-1] = 1;
564       }
565     }
566   }
567 }
568 
hIndAllMult(scmon pure,int Npure,scfmon rad,int Nrad,varset var,int Nvar)569 void hIndAllMult(scmon pure, int Npure, scfmon rad, int Nrad,
570  varset var, int Nvar)
571 {
572   int  dn, iv, rad0, b, c, x;
573   scmon pn;
574   scfmon rn;
575   if (Nrad < 2)
576   {
577     dn = Npure + Nrad;
578     if (dn > hCo)
579     {
580       if (!Nrad)
581         hCheckIndep(pure);
582       else
583       {
584         pn = *rad;
585         for (iv = Nvar; iv; iv--)
586         {
587           x = var[iv];
588           if (pn[x])
589           {
590             pure[x] = 1;
591             hCheckIndep(pure);
592             pure[x] = 0;
593           }
594         }
595       }
596     }
597     return;
598   }
599   iv = Nvar;
600   while(pure[var[iv]]) iv--;
601   hStepR(rad, Nrad, var, iv, &rad0);
602   iv--;
603   if (rad0 < Nrad)
604   {
605     pn = hGetpure(pure);
606     rn = hGetmem(Nrad, rad, radmem[iv]);
607     pn[var[iv + 1]] = 1;
608     hIndAllMult(pn, Npure + 1, rn, rad0, var, iv);
609     pn[var[iv + 1]] = 0;
610     b = rad0;
611     c = Nrad;
612     hElimR(rn, &rad0, b, c, var, iv);
613     hPure(rn, b, &c, var, iv, pn, &x);
614     hLex2R(rn, rad0, b, c, var, iv, hwork);
615     rad0 += (c - b);
616     hIndAllMult(pn, Npure + x, rn, rad0, var, iv);
617   }
618   else
619   {
620     hIndAllMult(pure, Npure, rad, Nrad, var, iv);
621   }
622 }
623 
624 // multiplicity
625 
hZeroMult(scmon pure,scfmon stc,int Nstc,varset var,int Nvar)626 static int hZeroMult(scmon pure, scfmon stc, int Nstc, varset var, int Nvar)
627 {
628   int  iv = Nvar -1, sum, a, a0, a1, b, i;
629   int  x, x0;
630   scmon pn;
631   scfmon sn;
632   if (!iv)
633     return pure[var[1]];
634   else if (!Nstc)
635   {
636     sum = 1;
637     for (i = Nvar; i; i--)
638       sum *= pure[var[i]];
639     return sum;
640   }
641   x = a = 0;
642   pn = hGetpure(pure);
643   sn = hGetmem(Nstc, stc, stcmem[iv]);
644   hStepS(sn, Nstc, var, Nvar, &a, &x);
645   if (a == Nstc)
646     return pure[var[Nvar]] * hZeroMult(pn, sn, a, var, iv);
647   else
648     sum = x * hZeroMult(pn, sn, a, var, iv);
649   b = a;
650   loop
651   {
652     a0 = a;
653     x0 = x;
654     hStepS(sn, Nstc, var, Nvar, &a, &x);
655     hElimS(sn, &b, a0, a, var, iv);
656     a1 = a;
657     hPure(sn, a0, &a1, var, iv, pn, &i);
658     hLex2S(sn, b, a0, a1, var, iv, hwork);
659     b += (a1 - a0);
660     if (a < Nstc)
661     {
662       sum += (x - x0) * hZeroMult(pn, sn, b, var, iv);
663     }
664     else
665     {
666       sum += (pure[var[Nvar]] - x0) * hZeroMult(pn, sn, b, var, iv);
667       return sum;
668     }
669   }
670 }
671 
hProject(scmon pure,varset sel)672 static void hProject(scmon pure, varset sel)
673 {
674   int  i, i0, k;
675   i0 = 0;
676   for (i = 1; i <= (currRing->N); i++)
677   {
678     if (pure[i])
679     {
680       i0++;
681       sel[i0] = i;
682     }
683   }
684   i = hNstc;
685   memcpy(hwork, hstc, i * sizeof(scmon));
686   hStaircase(hwork, &i, sel, i0);
687   if ((i0 > 2) && (i > 10))
688     hOrdSupp(hwork, i, sel, i0);
689   memset(hpur0, 0, ((currRing->N) + 1) * sizeof(int));
690   hPure(hwork, 0, &i, sel, i0, hpur0, &k);
691   hLexS(hwork, i, sel, i0);
692   hMu += hZeroMult(hpur0, hwork, i, sel, i0);
693 }
694 
hDimMult(scmon pure,int Npure,scfmon rad,int Nrad,varset var,int Nvar)695 static void hDimMult(scmon pure, int Npure, scfmon rad, int Nrad,
696  varset var, int Nvar)
697 {
698   int  dn, iv, rad0, b, c, x;
699   scmon pn;
700   scfmon rn;
701   if (Nrad < 2)
702   {
703     dn = Npure + Nrad;
704     if (dn == hCo)
705     {
706       if (!Nrad)
707         hProject(pure, hsel);
708       else
709       {
710         pn = *rad;
711         for (iv = Nvar; iv; iv--)
712         {
713           x = var[iv];
714           if (pn[x])
715           {
716             pure[x] = 1;
717             hProject(pure, hsel);
718             pure[x] = 0;
719           }
720         }
721       }
722     }
723     return;
724   }
725   iv = Nvar;
726   dn = Npure+1;
727   if (dn >= hCo)
728   {
729     if (dn > hCo)
730       return;
731     loop
732     {
733       if(!pure[var[iv]])
734       {
735         if(hNotZero(rad, Nrad, var, iv))
736         {
737           pure[var[iv]] = 1;
738           hProject(pure, hsel);
739           pure[var[iv]] = 0;
740         }
741       }
742       iv--;
743       if (!iv)
744         return;
745     }
746   }
747   while(pure[var[iv]]) iv--;
748   hStepR(rad, Nrad, var, iv, &rad0);
749   iv--;
750   if (rad0 < Nrad)
751   {
752     pn = hGetpure(pure);
753     rn = hGetmem(Nrad, rad, radmem[iv]);
754     pn[var[iv + 1]] = 1;
755     hDimMult(pn, Npure + 1, rn, rad0, var, iv);
756     pn[var[iv + 1]] = 0;
757     b = rad0;
758     c = Nrad;
759     hElimR(rn, &rad0, b, c, var, iv);
760     hPure(rn, b, &c, var, iv, pn, &x);
761     hLex2R(rn, rad0, b, c, var, iv, hwork);
762     rad0 += (c - b);
763     hDimMult(pn, Npure + x, rn, rad0, var, iv);
764   }
765   else
766   {
767     hDimMult(pure, Npure, rad, Nrad, var, iv);
768   }
769 }
770 
hDegree(ideal S,ideal Q)771 static void hDegree(ideal S, ideal Q)
772 {
773   id_Test(S, currRing);
774   if( Q!=NULL ) id_Test(Q, currRing);
775 
776   int  di;
777   int  mc;
778   hexist = hInit(S, Q, &hNexist, currRing);
779   if (!hNexist)
780   {
781     hCo = 0;
782     hMu = 1;
783     return;
784   }
785   //hWeight();
786   hwork = (scfmon)omAlloc(hNexist * sizeof(scmon));
787   hvar = (varset)omAlloc(((currRing->N) + 1) * sizeof(int));
788   hsel = (varset)omAlloc(((currRing->N) + 1) * sizeof(int));
789   hpure = (scmon)omAlloc((1 + ((currRing->N) * (currRing->N))) * sizeof(int));
790   hpur0 = (scmon)omAlloc((1 + ((currRing->N) * (currRing->N))) * sizeof(int));
791   mc = hisModule;
792   hrad = (scfmon)omAlloc(hNexist * sizeof(scmon));
793   if (!mc)
794   {
795     memcpy(hrad, hexist, hNexist * sizeof(scmon));
796     hstc = hexist;
797     hNrad = hNstc = hNexist;
798   }
799   else
800     hstc = (scfmon)omAlloc(hNexist * sizeof(scmon));
801   radmem = hCreate((currRing->N) - 1);
802   stcmem = hCreate((currRing->N) - 1);
803   hCo = (currRing->N) + 1;
804   di = hCo + 1;
805   loop
806   {
807     if (mc)
808     {
809       hComp(hexist, hNexist, mc, hrad, &hNrad);
810       hNstc = hNrad;
811       memcpy(hstc, hrad, hNrad * sizeof(scmon));
812     }
813     if (hNrad)
814     {
815       hNvar = (currRing->N);
816       hRadical(hrad, &hNrad, hNvar);
817       hSupp(hrad, hNrad, hvar, &hNvar);
818       if (hNvar)
819       {
820         hCo = hNvar;
821         memset(hpure, 0, ((currRing->N) + 1) * sizeof(int));
822         hPure(hrad, 0, &hNrad, hvar, hNvar, hpure, &hNpure);
823         hLexR(hrad, hNrad, hvar, hNvar);
824         hDimSolve(hpure, hNpure, hrad, hNrad, hvar, hNvar);
825       }
826     }
827     else
828     {
829       hNvar = 1;
830       hCo = 0;
831     }
832     if (hCo < di)
833     {
834       di = hCo;
835       hMu = 0;
836     }
837     if (hNvar && (hCo == di))
838     {
839       if (di && (di < (currRing->N)))
840         hDimMult(hpure, hNpure, hrad, hNrad, hvar, hNvar);
841       else if (!di)
842         hMu++;
843       else
844       {
845         hStaircase(hstc, &hNstc, hvar, hNvar);
846         if ((hNvar > 2) && (hNstc > 10))
847           hOrdSupp(hstc, hNstc, hvar, hNvar);
848         memset(hpur0, 0, ((currRing->N) + 1) * sizeof(int));
849         hPure(hstc, 0, &hNstc, hvar, hNvar, hpur0, &hNpure);
850         hLexS(hstc, hNstc, hvar, hNvar);
851         hMu += hZeroMult(hpur0, hstc, hNstc, hvar, hNvar);
852       }
853     }
854     mc--;
855     if (mc <= 0)
856       break;
857   }
858   hCo = di;
859   hKill(stcmem, (currRing->N) - 1);
860   hKill(radmem, (currRing->N) - 1);
861   omFreeSize((ADDRESS)hpur0, (1 + ((currRing->N) * (currRing->N))) * sizeof(int));
862   omFreeSize((ADDRESS)hpure, (1 + ((currRing->N) * (currRing->N))) * sizeof(int));
863   omFreeSize((ADDRESS)hsel, ((currRing->N) + 1) * sizeof(int));
864   omFreeSize((ADDRESS)hvar, ((currRing->N) + 1) * sizeof(int));
865   omFreeSize((ADDRESS)hwork, hNexist * sizeof(scmon));
866   omFreeSize((ADDRESS)hrad, hNexist * sizeof(scmon));
867   hDelete(hexist, hNexist);
868   if (hisModule)
869     omFreeSize((ADDRESS)hstc, hNexist * sizeof(scmon));
870 }
871 
scMultInt(ideal S,ideal Q)872 int  scMultInt(ideal S, ideal Q)
873 {
874   id_Test(S, currRing);
875   if( Q!=NULL ) id_Test(Q, currRing);
876 
877   hDegree(S, Q);
878   return hMu;
879 }
880 
scPrintDegree(int co,int mu)881 void scPrintDegree(int co, int mu)
882 {
883   int di = (currRing->N)-co;
884   if (currRing->OrdSgn == 1)
885   {
886     if (di>0)
887       Print("// dimension (proj.)  = %d\n// degree (proj.)   = %d\n", di-1, mu);
888     else
889       Print("// dimension (affine) = 0\n// degree (affine)  = %d\n",       mu);
890   }
891   else
892     Print("// dimension (local)   = %d\n// multiplicity = %d\n", di, mu);
893 }
894 
scDegree(ideal S,intvec * modulweight,ideal Q)895 void scDegree(ideal S, intvec *modulweight, ideal Q)
896 {
897   id_Test(S, currRing);
898   if( Q!=NULL ) id_Test(Q, currRing);
899 
900   int co, mu, l;
901   intvec *hseries2;
902   intvec *hseries1 = hFirstSeries(S, modulweight, Q);
903   l = hseries1->length()-1;
904   if (l > 1)
905     hseries2 = hSecondSeries(hseries1);
906   else
907     hseries2 = hseries1;
908   hDegreeSeries(hseries1, hseries2, &co, &mu);
909   if ((l == 1) &&(mu == 0))
910     scPrintDegree((currRing->N)+1, 0);
911   else
912     scPrintDegree(co, mu);
913   if (l>1)
914     delete hseries1;
915   delete hseries2;
916 }
917 
hDegree0(ideal S,ideal Q,const ring tailRing)918 static void hDegree0(ideal S, ideal Q, const ring tailRing)
919 {
920   id_TestTail(S, currRing, tailRing);
921   if (Q!=NULL) id_TestTail(Q, currRing, tailRing);
922 
923   int  mc;
924   hexist = hInit(S, Q, &hNexist, tailRing);
925   if (!hNexist)
926   {
927     hMu = -1;
928     return;
929   }
930   else
931     hMu = 0;
932 
933   const ring r = currRing;
934 
935   hwork = (scfmon)omAlloc(hNexist * sizeof(scmon));
936   hvar = (varset)omAlloc(((r->N) + 1) * sizeof(int));
937   hpur0 = (scmon)omAlloc((1 + ((r->N) * (r->N))) * sizeof(int));
938   mc = hisModule;
939   if (!mc)
940   {
941     hstc = hexist;
942     hNstc = hNexist;
943   }
944   else
945     hstc = (scfmon)omAlloc(hNexist * sizeof(scmon));
946   stcmem = hCreate((r->N) - 1);
947   loop
948   {
949     if (mc)
950     {
951       hComp(hexist, hNexist, mc, hstc, &hNstc);
952       if (!hNstc)
953       {
954         hMu = -1;
955         break;
956       }
957     }
958     hNvar = (r->N);
959     for (int i = hNvar; i; i--)
960       hvar[i] = i;
961     hStaircase(hstc, &hNstc, hvar, hNvar);
962     hSupp(hstc, hNstc, hvar, &hNvar);
963     if ((hNvar == (r->N)) && (hNstc >= (r->N)))
964     {
965       if ((hNvar > 2) && (hNstc > 10))
966         hOrdSupp(hstc, hNstc, hvar, hNvar);
967       memset(hpur0, 0, ((r->N) + 1) * sizeof(int));
968       hPure(hstc, 0, &hNstc, hvar, hNvar, hpur0, &hNpure);
969       if (hNpure == hNvar)
970       {
971         hLexS(hstc, hNstc, hvar, hNvar);
972         hMu += hZeroMult(hpur0, hstc, hNstc, hvar, hNvar);
973       }
974       else
975         hMu = -1;
976     }
977     else if (hNvar)
978       hMu = -1;
979     mc--;
980     if (mc <= 0 || hMu < 0)
981       break;
982   }
983   hKill(stcmem, (r->N) - 1);
984   omFreeSize((ADDRESS)hpur0, (1 + ((r->N) * (r->N))) * sizeof(int));
985   omFreeSize((ADDRESS)hvar, ((r->N) + 1) * sizeof(int));
986   omFreeSize((ADDRESS)hwork, hNexist * sizeof(scmon));
987   hDelete(hexist, hNexist);
988   if (hisModule)
989     omFreeSize((ADDRESS)hstc, hNexist * sizeof(scmon));
990 }
991 
scMult0Int(ideal S,ideal Q,const ring tailRing)992 int  scMult0Int(ideal S, ideal Q, const ring tailRing)
993 {
994   id_TestTail(S, currRing, tailRing);
995   if (Q!=NULL) id_TestTail(Q, currRing, tailRing);
996 
997   hDegree0(S, Q, tailRing);
998   return hMu;
999 }
1000 
1001 
1002 // HC
1003 
1004 STATIC_VAR poly pWork;
1005 
hHedge(poly hEdge)1006 static void hHedge(poly hEdge)
1007 {
1008   pSetm(pWork);
1009   if (pLmCmp(pWork, hEdge) == currRing->OrdSgn)
1010   {
1011     for (int i = hNvar; i>0; i--)
1012       pSetExp(hEdge,i, pGetExp(pWork,i));
1013     pSetm(hEdge);
1014   }
1015 }
1016 
1017 
hHedgeStep(scmon pure,scfmon stc,int Nstc,varset var,int Nvar,poly hEdge)1018 static void hHedgeStep(scmon pure, scfmon stc,
1019                        int Nstc, varset var, int Nvar,poly hEdge)
1020 {
1021   int  iv = Nvar -1, k = var[Nvar], a, a0, a1, b, i;
1022   int  x/*, x0*/;
1023   scmon pn;
1024   scfmon sn;
1025   if (iv==0)
1026   {
1027     pSetExp(pWork, k, pure[k]);
1028     hHedge(hEdge);
1029     return;
1030   }
1031   else if (Nstc==0)
1032   {
1033     for (i = Nvar; i>0; i--)
1034       pSetExp(pWork, var[i], pure[var[i]]);
1035     hHedge(hEdge);
1036     return;
1037   }
1038   x = a = 0;
1039   pn = hGetpure(pure);
1040   sn = hGetmem(Nstc, stc, stcmem[iv]);
1041   hStepS(sn, Nstc, var, Nvar, &a, &x);
1042   if (a == Nstc)
1043   {
1044     pSetExp(pWork, k, pure[k]);
1045     hHedgeStep(pn, sn, a, var, iv,hEdge);
1046     return;
1047   }
1048   else
1049   {
1050     pSetExp(pWork, k, x);
1051     hHedgeStep(pn, sn, a, var, iv,hEdge);
1052   }
1053   b = a;
1054   loop
1055   {
1056     a0 = a;
1057     // x0 = x;
1058     hStepS(sn, Nstc, var, Nvar, &a, &x);
1059     hElimS(sn, &b, a0, a, var, iv);
1060     a1 = a;
1061     hPure(sn, a0, &a1, var, iv, pn, &i);
1062     hLex2S(sn, b, a0, a1, var, iv, hwork);
1063     b += (a1 - a0);
1064     if (a < Nstc)
1065     {
1066       pSetExp(pWork, k, x);
1067       hHedgeStep(pn, sn, b, var, iv,hEdge);
1068     }
1069     else
1070     {
1071       pSetExp(pWork, k, pure[k]);
1072       hHedgeStep(pn, sn, b, var, iv,hEdge);
1073       return;
1074     }
1075   }
1076 }
1077 
scComputeHC(ideal S,ideal Q,int ak,poly & hEdge,ring tailRing)1078 void scComputeHC(ideal S, ideal Q, int ak, poly &hEdge, ring tailRing)
1079 {
1080   id_TestTail(S, currRing, tailRing);
1081   if (Q!=NULL) id_TestTail(Q, currRing, tailRing);
1082 
1083   int  i;
1084   int  k = ak;
1085   #ifdef HAVE_RINGS
1086   if (rField_is_Ring(currRing) && (currRing->OrdSgn == -1))
1087   {
1088     //consider just monic generators (over rings with zero-divisors)
1089     ideal SS=id_Copy(S,tailRing);
1090     for(i=0;i<=idElem(S);i++)
1091     {
1092       if((SS->m[i]!=NULL)
1093       && ((p_IsPurePower(SS->m[i],tailRing)==0)
1094         ||(!n_IsUnit(pGetCoeff(SS->m[i]), tailRing->cf))))
1095       {
1096         p_Delete(&SS->m[i],tailRing);
1097       }
1098     }
1099     S=id_Copy(SS,tailRing);
1100     idSkipZeroes(S);
1101   }
1102   #if 0
1103   printf("\nThis is HC:\n");
1104   for(int ii=0;ii<=idElem(S);ii++)
1105   {
1106     pWrite(S->m[ii]);
1107   }
1108   //getchar();
1109   #endif
1110   #endif
1111   if(idElem(S) == 0)
1112     return;
1113   hNvar = (currRing->N);
1114   hexist = hInit(S, Q, &hNexist, tailRing); // tailRing?
1115   if (k!=0)
1116     hComp(hexist, hNexist, k, hexist, &hNstc);
1117   else
1118     hNstc = hNexist;
1119   assume(hNexist > 0);
1120   hwork = (scfmon)omAlloc(hNexist * sizeof(scmon));
1121   hvar = (varset)omAlloc((hNvar + 1) * sizeof(int));
1122   hpure = (scmon)omAlloc((1 + (hNvar * hNvar)) * sizeof(int));
1123   stcmem = hCreate(hNvar - 1);
1124   for (i = hNvar; i>0; i--)
1125     hvar[i] = i;
1126   hStaircase(hexist, &hNstc, hvar, hNvar);
1127   if ((hNvar > 2) && (hNstc > 10))
1128     hOrdSupp(hexist, hNstc, hvar, hNvar);
1129   memset(hpure, 0, (hNvar + 1) * sizeof(int));
1130   hPure(hexist, 0, &hNstc, hvar, hNvar, hpure, &hNpure);
1131   hLexS(hexist, hNstc, hvar, hNvar);
1132   if (hEdge!=NULL)
1133     pLmFree(hEdge);
1134   hEdge = pInit();
1135   pWork = pInit();
1136   hHedgeStep(hpure, hexist, hNstc, hvar, hNvar,hEdge);
1137   pSetComp(hEdge,ak);
1138   hKill(stcmem, hNvar - 1);
1139   omFreeSize((ADDRESS)hwork, hNexist * sizeof(scmon));
1140   omFreeSize((ADDRESS)hvar, (hNvar + 1) * sizeof(int));
1141   omFreeSize((ADDRESS)hpure, (1 + (hNvar * hNvar)) * sizeof(int));
1142   hDelete(hexist, hNexist);
1143   pLmFree(pWork);
1144 }
1145 
1146 
1147 
1148 //  kbase
1149 
1150 STATIC_VAR poly last;
1151 STATIC_VAR scmon act;
1152 
scElKbase()1153 static void scElKbase()
1154 {
1155   poly q = pInit();
1156   pSetCoeff0(q,nInit(1));
1157   pSetExpV(q,act);
1158   pNext(q) = NULL;
1159   last = pNext(last) = q;
1160 }
1161 
scMax(int i,scfmon stc,int Nvar)1162 static int scMax( int i, scfmon stc, int Nvar)
1163 {
1164   int x, y=stc[0][Nvar];
1165   for (; i;)
1166   {
1167     i--;
1168     x = stc[i][Nvar];
1169     if (x > y) y = x;
1170   }
1171   return y;
1172 }
1173 
scMin(int i,scfmon stc,int Nvar)1174 static int scMin( int i, scfmon stc, int Nvar)
1175 {
1176   int x, y=stc[0][Nvar];
1177   for (; i;)
1178   {
1179     i--;
1180     x = stc[i][Nvar];
1181     if (x < y) y = x;
1182   }
1183   return y;
1184 }
1185 
scRestrict(int & Nstc,scfmon stc,int Nvar)1186 static int scRestrict( int &Nstc, scfmon stc, int Nvar)
1187 {
1188   int x, y;
1189   int i, j, Istc = Nstc;
1190 
1191   y = MAX_INT_VAL;
1192   for (i=Nstc-1; i>=0; i--)
1193   {
1194     j = Nvar-1;
1195     loop
1196     {
1197       if(stc[i][j] != 0) break;
1198       j--;
1199       if (j == 0)
1200       {
1201         Istc--;
1202         x = stc[i][Nvar];
1203         if (x < y) y = x;
1204         stc[i] = NULL;
1205         break;
1206       }
1207     }
1208   }
1209   if (Istc < Nstc)
1210   {
1211     for (i=Nstc-1; i>=0; i--)
1212     {
1213       if (stc[i] && (stc[i][Nvar] >= y))
1214       {
1215         Istc--;
1216         stc[i] = NULL;
1217       }
1218     }
1219     j = 0;
1220     while (stc[j]) j++;
1221     i = j+1;
1222     for(; i<Nstc; i++)
1223     {
1224       if (stc[i])
1225       {
1226         stc[j] = stc[i];
1227         j++;
1228       }
1229     }
1230     Nstc = Istc;
1231     return y;
1232   }
1233   else
1234     return -1;
1235 }
1236 
scAll(int Nvar,int deg)1237 static void scAll( int Nvar, int deg)
1238 {
1239   int i;
1240   int d = deg;
1241   if (d == 0)
1242   {
1243     for (i=Nvar; i; i--) act[i] = 0;
1244     scElKbase();
1245     return;
1246   }
1247   if (Nvar == 1)
1248   {
1249     act[1] = d;
1250     scElKbase();
1251     return;
1252   }
1253   do
1254   {
1255     act[Nvar] = d;
1256     scAll(Nvar-1, deg-d);
1257     d--;
1258   } while (d >= 0);
1259 }
1260 
scAllKbase(int Nvar,int ideg,int deg)1261 static void scAllKbase( int Nvar, int ideg, int deg)
1262 {
1263   do
1264   {
1265     act[Nvar] = ideg;
1266     scAll(Nvar-1, deg-ideg);
1267     ideg--;
1268   } while (ideg >= 0);
1269 }
1270 
scDegKbase(scfmon stc,int Nstc,int Nvar,int deg)1271 static void scDegKbase( scfmon stc, int Nstc, int Nvar, int deg)
1272 {
1273   int  Ivar, Istc, i, j;
1274   scfmon sn;
1275   int x, ideg;
1276 
1277   if (deg == 0)
1278   {
1279     for (i=Nstc-1; i>=0; i--)
1280     {
1281       for (j=Nvar;j;j--){ if(stc[i][j]) break; }
1282       if (j==0){return;}
1283     }
1284     for (i=Nvar; i; i--) act[i] = 0;
1285     scElKbase();
1286     return;
1287   }
1288   if (Nvar == 1)
1289   {
1290     for (i=Nstc-1; i>=0; i--) if(deg >= stc[i][1]) return;
1291     act[1] = deg;
1292     scElKbase();
1293     return;
1294   }
1295   Ivar = Nvar-1;
1296   sn = hGetmem(Nstc, stc, stcmem[Ivar]);
1297   x = scRestrict(Nstc, sn, Nvar);
1298   if (x <= 0)
1299   {
1300     if (x == 0) return;
1301     ideg = deg;
1302   }
1303   else
1304   {
1305     if (deg < x) ideg = deg;
1306     else ideg = x-1;
1307     if (Nstc == 0)
1308     {
1309       scAllKbase(Nvar, ideg, deg);
1310       return;
1311     }
1312   }
1313   loop
1314   {
1315     x = scMax(Nstc, sn, Nvar);
1316     while (ideg >= x)
1317     {
1318       act[Nvar] = ideg;
1319       scDegKbase(sn, Nstc, Ivar, deg-ideg);
1320       ideg--;
1321     }
1322     if (ideg < 0) return;
1323     Istc = Nstc;
1324     for (i=Nstc-1; i>=0; i--)
1325     {
1326       if (ideg < sn[i][Nvar])
1327       {
1328         Istc--;
1329         sn[i] = NULL;
1330       }
1331     }
1332     if (Istc == 0)
1333     {
1334       scAllKbase(Nvar, ideg, deg);
1335       return;
1336     }
1337     j = 0;
1338     while (sn[j]) j++;
1339     i = j+1;
1340     for (; i<Nstc; i++)
1341     {
1342       if (sn[i])
1343       {
1344         sn[j] = sn[i];
1345         j++;
1346       }
1347     }
1348     Nstc = Istc;
1349   }
1350 }
1351 
scInKbase(scfmon stc,int Nstc,int Nvar)1352 static void scInKbase( scfmon stc, int Nstc, int Nvar)
1353 {
1354   int  Ivar, Istc, i, j;
1355   scfmon sn;
1356   int x, ideg;
1357 
1358   if (Nvar == 1)
1359   {
1360     ideg = scMin(Nstc, stc, 1);
1361     while (ideg > 0)
1362     {
1363       ideg--;
1364       act[1] = ideg;
1365       scElKbase();
1366     }
1367     return;
1368   }
1369   Ivar = Nvar-1;
1370   sn = hGetmem(Nstc, stc, stcmem[Ivar]);
1371   x = scRestrict(Nstc, sn, Nvar);
1372   if (x == 0) return;
1373   ideg = x-1;
1374   loop
1375   {
1376     x = scMax(Nstc, sn, Nvar);
1377     while (ideg >= x)
1378     {
1379       act[Nvar] = ideg;
1380       scInKbase(sn, Nstc, Ivar);
1381       ideg--;
1382     }
1383     if (ideg < 0) return;
1384     Istc = Nstc;
1385     for (i=Nstc-1; i>=0; i--)
1386     {
1387       if (ideg < sn[i][Nvar])
1388       {
1389         Istc--;
1390         sn[i] = NULL;
1391       }
1392     }
1393     j = 0;
1394     while (sn[j]) j++;
1395     i = j+1;
1396     for (; i<Nstc; i++)
1397     {
1398       if (sn[i])
1399       {
1400         sn[j] = sn[i];
1401         j++;
1402       }
1403     }
1404     Nstc = Istc;
1405   }
1406 }
1407 
scIdKbase(poly q,const int rank)1408 static ideal scIdKbase(poly q, const int rank)
1409 {
1410   ideal res = idInit(pLength(q), rank);
1411   polyset mm = res->m;
1412   do
1413   {
1414     *mm = q; ++mm;
1415 
1416     const poly p = pNext(q);
1417     pNext(q) = NULL;
1418     q = p;
1419 
1420   } while (q!=NULL);
1421 
1422   id_Test(res, currRing);   // WRONG RANK!!!???
1423   return res;
1424 }
1425 
scKBase(int deg,ideal s,ideal Q,intvec * mv)1426 ideal scKBase(int deg, ideal s, ideal Q, intvec * mv)
1427 {
1428   if( Q!=NULL) id_Test(Q, currRing);
1429 
1430   int  i, di;
1431   poly p;
1432 
1433   if (deg < 0)
1434   {
1435     di = scDimInt(s, Q);
1436     if (di != 0)
1437     {
1438       //Werror("KBase not finite");
1439       return idInit(1,s->rank);
1440     }
1441   }
1442   stcmem = hCreate((currRing->N) - 1);
1443   hexist = hInit(s, Q, &hNexist, currRing);
1444   p = last = pInit();
1445   /*pNext(p) = NULL;*/
1446   act = (scmon)omAlloc(((currRing->N) + 1) * sizeof(int));
1447   *act = 0;
1448   if (!hNexist)
1449   {
1450     scAll((currRing->N), deg);
1451     goto ende;
1452   }
1453   if (!hisModule)
1454   {
1455     if (deg < 0) scInKbase(hexist, hNexist, (currRing->N));
1456     else scDegKbase(hexist, hNexist, (currRing->N), deg);
1457   }
1458   else
1459   {
1460     hstc = (scfmon)omAlloc(hNexist * sizeof(scmon));
1461     for (i = 1; i <= hisModule; i++)
1462     {
1463       *act = i;
1464       hComp(hexist, hNexist, i, hstc, &hNstc);
1465       int deg_ei=deg;
1466       if (mv!=NULL) deg_ei -= (*mv)[i-1];
1467       if ((deg < 0) || (deg_ei>=0))
1468       {
1469         if (hNstc)
1470         {
1471           if (deg < 0) scInKbase(hstc, hNstc, (currRing->N));
1472           else scDegKbase(hstc, hNstc, (currRing->N), deg_ei);
1473         }
1474         else
1475           scAll((currRing->N), deg_ei);
1476       }
1477     }
1478     omFreeSize((ADDRESS)hstc, hNexist * sizeof(scmon));
1479   }
1480 ende:
1481   hDelete(hexist, hNexist);
1482   omFreeSize((ADDRESS)act, ((currRing->N) + 1) * sizeof(int));
1483   hKill(stcmem, (currRing->N) - 1);
1484   pLmFree(&p);
1485   if (p == NULL)
1486     return idInit(1,s->rank);
1487 
1488   last = p;
1489   return scIdKbase(p, s->rank);
1490 }
1491 
1492 #if 0 //-- alternative implementation of scComputeHC
1493 /*
1494 void scComputeHCw(ideal ss, ideal Q, int ak, poly &hEdge, ring tailRing)
1495 {
1496   id_TestTail(ss, currRing, tailRing);
1497   if (Q!=NULL) id_TestTail(Q, currRing, tailRing);
1498 
1499   int  i, di;
1500   poly p;
1501 
1502   if (hEdge!=NULL)
1503     pLmFree(hEdge);
1504 
1505   ideal s=idInit(IDELEMS(ss),ak);
1506   for(i=IDELEMS(ss)-1;i>=0;i--)
1507   {
1508     if (ss->m[i]!=NULL) s->m[i]=pHead(ss->m[i]);
1509   }
1510   di = scDimInt(s, Q);
1511   stcmem = hCreate((currRing->N) - 1);
1512   hexist = hInit(s, Q, &hNexist, currRing);
1513   p = last = pInit();
1514   // pNext(p) = NULL;
1515   act = (scmon)omAlloc(((currRing->N) + 1) * sizeof(int));
1516   *act = 0;
1517   if (!hNexist)
1518   {
1519     scAll((currRing->N), -1);
1520     goto ende;
1521   }
1522   if (!hisModule)
1523   {
1524     scInKbase(hexist, hNexist, (currRing->N));
1525   }
1526   else
1527   {
1528     hstc = (scfmon)omAlloc(hNexist * sizeof(scmon));
1529     for (i = 1; i <= hisModule; i++)
1530     {
1531       *act = i;
1532       hComp(hexist, hNexist, i, hstc, &hNstc);
1533       if (hNstc)
1534       {
1535         scInKbase(hstc, hNstc, (currRing->N));
1536       }
1537       else
1538         scAll((currRing->N), -1);
1539     }
1540     omFreeSize((ADDRESS)hstc, hNexist * sizeof(scmon));
1541   }
1542 ende:
1543   hDelete(hexist, hNexist);
1544   omFreeSize((ADDRESS)act, ((currRing->N) + 1) * sizeof(int));
1545   hKill(stcmem, (currRing->N) - 1);
1546   pDeleteLm(&p);
1547   idDelete(&s);
1548   if (p == NULL)
1549   {
1550     return; // no HEdge
1551   }
1552   else
1553   {
1554     last = p;
1555     ideal res=scIdKbase(p, ss->rank);
1556     poly p_ind=res->m[0]; int ind=0;
1557     for(i=IDELEMS(res)-1;i>0;i--)
1558     {
1559       if (pCmp(res->m[i],p_ind)==-1) { p_ind=res->m[i]; ind=i; }
1560     }
1561     assume(p_ind!=NULL);
1562     assume(res->m[ind]==p_ind);
1563     hEdge=p_ind;
1564     res->m[ind]=NULL;
1565     nDelete(&pGetCoeff(hEdge));
1566     pGetCoeff(hEdge)=NULL;
1567     for(i=(currRing->N);i>0;i--)
1568       pIncrExp(hEdge,i);
1569     pSetm(hEdge);
1570 
1571     idDelete(&res);
1572     return;
1573   }
1574 }
1575  */
1576 #endif
1577 
1578 #ifdef HAVE_SHIFTBBA
1579 
1580 /*
1581  * Computation of the Gel'fand-Kirillov Dimension
1582  */
1583 
1584 #include "polys/shiftop.h"
1585 #include <vector>
1586 
countCycles(const intvec * _G,int v,std::vector<int> path,std::vector<BOOLEAN> visited,std::vector<BOOLEAN> cyclic,std::vector<int> cache)1587 static std::vector<int> countCycles(const intvec* _G, int v, std::vector<int> path, std::vector<BOOLEAN> visited, std::vector<BOOLEAN> cyclic, std::vector<int> cache)
1588 {
1589   intvec* G = ivCopy(_G); // modifications must be local
1590 
1591   if (cache[v] != -2) return cache; // value is already cached
1592 
1593   visited[v] = TRUE;
1594   path.push_back(v);
1595 
1596   int cycles = 0;
1597   for (int w = 0; w < G->cols(); w++)
1598   {
1599     if (IMATELEM(*G, v + 1, w + 1)) // edge v -> w exists in G
1600     {
1601       if (!visited[w])
1602       { // continue with w
1603         cache = countCycles(G, w, path, visited, cyclic, cache);
1604         if (cache[w] == -1)
1605         {
1606           cache[v] = -1;
1607           return cache;
1608         }
1609         cycles = si_max(cycles, cache[w]);
1610       }
1611       else
1612       { // found new cycle
1613         int pathIndexOfW = -1;
1614         for (int i = path.size() - 1; i >= 0; i--) {
1615           if (cyclic[path[i]] == 1) { // found an already cyclic vertex
1616             cache[v] = -1;
1617             return cache;
1618           }
1619           cyclic[path[i]] = TRUE;
1620 
1621           if (path[i] == w) { // end of the cycle
1622             assume(IMATELEM(*G, v + 1, w + 1) != 0);
1623             IMATELEM(*G, v + 1, w + 1) = 0; // remove edge v -> w
1624             pathIndexOfW = i;
1625             break;
1626           } else {
1627             assume(IMATELEM(*G, path[i - 1] + 1, path[i] + 1) != 0);
1628             IMATELEM(*G, path[i - 1] + 1, path[i] + 1) = 0; // remove edge vi-1 -> vi
1629           }
1630         }
1631         assume(pathIndexOfW != -1); // should never happen
1632         for (int i = path.size() - 1; i >= pathIndexOfW; i--) {
1633           cache = countCycles(G, path[i], path, visited, cyclic, cache);
1634           if (cache[path[i]] == -1)
1635           {
1636             cache[v] = -1;
1637             return cache;
1638           }
1639           cycles = si_max(cycles, cache[path[i]] + 1);
1640         }
1641       }
1642     }
1643   }
1644   cache[v] = cycles;
1645 
1646   delete G;
1647   return cache;
1648 }
1649 
1650 // -1 is infinity
graphGrowth(const intvec * G)1651 static int graphGrowth(const intvec* G)
1652 {
1653   // init
1654   int n = G->cols();
1655   std::vector<int> path;
1656   std::vector<BOOLEAN> visited;
1657   std::vector<BOOLEAN> cyclic;
1658   std::vector<int> cache;
1659   visited.resize(n, FALSE);
1660   cyclic.resize(n, FALSE);
1661   cache.resize(n, -2);
1662 
1663   // get max number of cycles
1664   int cycles = 0;
1665   for (int v = 0; v < n; v++)
1666   {
1667     cache = countCycles(G, v, path, visited, cyclic, cache);
1668     if (cache[v] == -1)
1669       return -1;
1670     cycles = si_max(cycles, cache[v]);
1671   }
1672   return cycles;
1673 }
1674 
1675 // ATTENTION:
1676 //  - `words` contains the words normal modulo M of length n
1677 //  - `numberOfNormalWords` contains the number of words normal modulo M of length 0 ... n
_lp_computeNormalWords(ideal words,int & numberOfNormalWords,int length,ideal M,int minDeg,int & last)1678 static void _lp_computeNormalWords(ideal words, int& numberOfNormalWords, int length, ideal M, int minDeg, int& last)
1679 {
1680   if (length <= 0){
1681     poly one = pOne();
1682     if (p_LPDivisibleBy(M, one, currRing)) // 1 \in M => no normal words at all
1683     {
1684       pDelete(&one);
1685       last = -1;
1686       numberOfNormalWords = 0;
1687     }
1688     else
1689     {
1690       words->m[0] = one;
1691       last = 0;
1692       numberOfNormalWords = 1;
1693     }
1694     return;
1695   }
1696 
1697   _lp_computeNormalWords(words, numberOfNormalWords, length - 1, M, minDeg, last);
1698 
1699   int nVars = currRing->isLPring - currRing->LPncGenCount;
1700   int numberOfNewNormalWords = 0;
1701 
1702   for (int j = nVars - 1; j >= 0; j--)
1703   {
1704     for (int i = last; i >= 0; i--)
1705     {
1706       int index = (j * (last + 1)) + i;
1707 
1708       if (words->m[i] != NULL)
1709       {
1710         if (j > 0) {
1711           words->m[index] = pCopy(words->m[i]);
1712         }
1713 
1714         int varOffset = ((length - 1) * currRing->isLPring) + 1;
1715         pSetExp(words->m[index], varOffset + j, 1);
1716         pSetm(words->m[index]);
1717         pTest(words->m[index]);
1718 
1719         if (length >= minDeg && p_LPDivisibleBy(M, words->m[index], currRing))
1720         {
1721           pDelete(&words->m[index]);
1722           words->m[index] = NULL;
1723         }
1724         else
1725         {
1726           numberOfNewNormalWords++;
1727         }
1728       }
1729     }
1730   }
1731 
1732   last = nVars * last + nVars - 1;
1733 
1734   numberOfNormalWords += numberOfNewNormalWords;
1735 }
1736 
lp_computeNormalWords(int length,ideal M)1737 static ideal lp_computeNormalWords(int length, ideal M)
1738 {
1739   long minDeg = IDELEMS(M) > 0 ? pTotaldegree(M->m[0]) : 0;
1740   for (int i = 1; i < IDELEMS(M); i++)
1741   {
1742     minDeg = si_min(minDeg, pTotaldegree(M->m[i]));
1743   }
1744 
1745   int nVars = currRing->isLPring - currRing->LPncGenCount;
1746 
1747   int maxElems = 1;
1748   for (int i = 0; i < length; i++) // maxElems = nVars^n
1749     maxElems *= nVars;
1750   ideal words = idInit(maxElems);
1751   int last, numberOfNormalWords;
1752   _lp_computeNormalWords(words, numberOfNormalWords, length, M, minDeg, last);
1753   idSkipZeroes(words);
1754   return words;
1755 }
1756 
lp_countNormalWords(int upToLength,ideal M)1757 static int lp_countNormalWords(int upToLength, ideal M)
1758 {
1759   long minDeg = IDELEMS(M) > 0 ? pTotaldegree(M->m[0]) : 0;
1760   for (int i = 1; i < IDELEMS(M); i++)
1761   {
1762     minDeg = si_min(minDeg, pTotaldegree(M->m[i]));
1763   }
1764 
1765   int nVars = currRing->isLPring - currRing->LPncGenCount;
1766 
1767   int maxElems = 1;
1768   for (int i = 0; i < upToLength; i++) // maxElems = nVars^n
1769     maxElems *= nVars;
1770   ideal words = idInit(maxElems);
1771   int last, numberOfNormalWords;
1772   _lp_computeNormalWords(words, numberOfNormalWords, upToLength, M, minDeg, last);
1773   idDelete(&words);
1774   return numberOfNormalWords;
1775 }
1776 
1777 // NULL if graph is undefined
lp_ufnarovskiGraph(ideal G,ideal & standardWords)1778 intvec* lp_ufnarovskiGraph(ideal G, ideal &standardWords)
1779 {
1780   long l = 0;
1781   for (int i = 0; i < IDELEMS(G); i++)
1782     l = si_max(pTotaldegree(G->m[i]), l);
1783   l--;
1784   if (l <= 0)
1785   {
1786     WerrorS("Ufnarovski graph not implemented for l <= 0");
1787     return NULL;
1788   }
1789   int lV = currRing->isLPring;
1790 
1791   standardWords = lp_computeNormalWords(l, G);
1792 
1793   int n = IDELEMS(standardWords);
1794   intvec* UG = new intvec(n, n, 0);
1795   for (int i = 0; i < n; i++)
1796   {
1797     for (int j = 0; j < n; j++)
1798     {
1799       poly v = standardWords->m[i];
1800       poly w = standardWords->m[j];
1801 
1802       // check whether v*x1 = x2*w (overlap)
1803       bool overlap = true;
1804       for (int k = 1; k <= (l - 1) * lV; k++)
1805       {
1806         if (pGetExp(v, k + lV) != pGetExp(w, k)) {
1807           overlap = false;
1808           break;
1809         }
1810       }
1811 
1812       if (overlap)
1813       {
1814         // create the overlap
1815         poly p = pMult(pCopy(v), p_LPVarAt(w, l, currRing));
1816 
1817         // check whether the overlap is normal
1818         bool normal = true;
1819         for (int k = 0; k < IDELEMS(G); k++)
1820         {
1821           if (p_LPDivisibleBy(G->m[k], p, currRing))
1822           {
1823             normal = false;
1824             break;
1825           }
1826         }
1827 
1828         if (normal)
1829         {
1830           IMATELEM(*UG, i + 1, j + 1) = 1;
1831         }
1832       }
1833     }
1834   }
1835   return UG;
1836 }
1837 
1838 // -1 is infinity, -2 is error
lp_gkDim(const ideal _G)1839 int lp_gkDim(const ideal _G)
1840 {
1841   id_Test(_G, currRing);
1842 
1843   if (rField_is_Ring(currRing)) {
1844       WerrorS("GK-Dim not implemented for rings");
1845       return -2;
1846   }
1847 
1848   for (int i=IDELEMS(_G)-1;i>=0; i--)
1849   {
1850     if (_G->m[i] != NULL)
1851     {
1852       if (pGetComp(_G->m[i]) != 0)
1853       {
1854         WerrorS("GK-Dim not implemented for modules");
1855         return -2;
1856       }
1857       if (pGetNCGen(_G->m[i]) != 0)
1858       {
1859         WerrorS("GK-Dim not implemented for bi-modules");
1860         return -2;
1861       }
1862     }
1863   }
1864 
1865   ideal G = id_Head(_G, currRing); // G = LM(G) (and copy)
1866   idSkipZeroes(G); // remove zeros
1867   id_DelLmEquals(G, currRing); // remove duplicates
1868 
1869   // check if G is the zero ideal
1870   if (IDELEMS(G) == 1 && G->m[0] == NULL)
1871   {
1872     // NOTE: this is needed because if the ideal is <0>, then idSkipZeroes keeps this element, and IDELEMS is still 1!
1873     int lV = currRing->isLPring;
1874     int ncGenCount = currRing->LPncGenCount;
1875     if (lV - ncGenCount == 0)
1876     {
1877       idDelete(&G);
1878       return 0;
1879     }
1880     if (lV - ncGenCount == 1)
1881     {
1882       idDelete(&G);
1883       return 1;
1884     }
1885     if (lV - ncGenCount >= 2)
1886     {
1887       idDelete(&G);
1888       return -1;
1889     }
1890   }
1891 
1892   // get the max deg
1893   long maxDeg = 0;
1894   for (int i = 0; i < IDELEMS(G); i++)
1895   {
1896     maxDeg = si_max(maxDeg, pTotaldegree(G->m[i]));
1897 
1898     // also check whether G = <1>
1899     if (pIsConstantComp(G->m[i]))
1900     {
1901       WerrorS("GK-Dim not defined for 0-ring");
1902       idDelete(&G);
1903       return -2;
1904     }
1905   }
1906 
1907   // early termination if G \subset X
1908   if (maxDeg <= 1)
1909   {
1910     int lV = currRing->isLPring;
1911     int ncGenCount = currRing->LPncGenCount;
1912     if (IDELEMS(G) == lV - ncGenCount) // V = {1} no edges
1913     {
1914       idDelete(&G);
1915       return 0;
1916     }
1917     if (IDELEMS(G) == lV - ncGenCount - 1) // V = {1} with loop
1918     {
1919       idDelete(&G);
1920       return 1;
1921     }
1922     if (IDELEMS(G) <= lV - ncGenCount - 2) // V = {1} with more than one loop
1923     {
1924       idDelete(&G);
1925       return -1;
1926     }
1927   }
1928 
1929   ideal standardWords;
1930   intvec* UG = lp_ufnarovskiGraph(G, standardWords);
1931   if (UG == NULL)
1932   {
1933     idDelete(&G);
1934     return -2;
1935   }
1936   if (errorreported)
1937   {
1938     delete UG;
1939     idDelete(&G);
1940     return -2;
1941   }
1942   int gkDim = graphGrowth(UG);
1943   delete UG;
1944   idDelete(&G);
1945   return gkDim;
1946 }
1947 
1948 // converts an intvec matrix to a vector<vector<int> >
iv2vv(intvec * M)1949 static std::vector<std::vector<int> > iv2vv(intvec* M)
1950 {
1951   int rows = M->rows();
1952   int cols = M->cols();
1953 
1954   std::vector<std::vector<int> > mat(rows, std::vector<int>(cols));
1955 
1956   for (int i = 0; i < rows; i++)
1957   {
1958     for (int j = 0; j < cols; j++)
1959     {
1960       mat[i][j] = IMATELEM(*M, i + 1, j + 1);
1961     }
1962   }
1963 
1964   return mat;
1965 }
1966 
vvPrint(const std::vector<std::vector<int>> & mat)1967 static void vvPrint(const std::vector<std::vector<int> >& mat)
1968 {
1969   for (int i = 0; i < mat.size(); i++)
1970   {
1971     for (int j = 0; j < mat[i].size(); j++)
1972     {
1973       Print("%d ", mat[i][j]);
1974     }
1975     PrintLn();
1976   }
1977 }
1978 
vvTest(const std::vector<std::vector<int>> & mat)1979 static void vvTest(const std::vector<std::vector<int> >& mat)
1980 {
1981   if (mat.size() > 0)
1982   {
1983     int cols = mat[0].size();
1984     for (int i = 1; i < mat.size(); i++)
1985     {
1986       if (cols != mat[i].size())
1987         WerrorS("number of cols in matrix inconsistent");
1988     }
1989   }
1990 }
1991 
vvDeleteRow(std::vector<std::vector<int>> & mat,int row)1992 static void vvDeleteRow(std::vector<std::vector<int> >& mat, int row)
1993 {
1994   mat.erase(mat.begin() + row);
1995 }
1996 
vvDeleteColumn(std::vector<std::vector<int>> & mat,int col)1997 static void vvDeleteColumn(std::vector<std::vector<int> >& mat, int col)
1998 {
1999   for (int i = 0; i < mat.size(); i++)
2000   {
2001     mat[i].erase(mat[i].begin() + col);
2002   }
2003 }
2004 
vvIsRowZero(const std::vector<std::vector<int>> & mat,int row)2005 static BOOLEAN vvIsRowZero(const std::vector<std::vector<int> >& mat, int row)
2006 {
2007   for (int i = 0; i < mat[row].size(); i++)
2008   {
2009     if (mat[row][i] != 0)
2010       return FALSE;
2011   }
2012   return TRUE;
2013 }
2014 
vvIsColumnZero(const std::vector<std::vector<int>> & mat,int col)2015 static BOOLEAN vvIsColumnZero(const std::vector<std::vector<int> >& mat, int col)
2016 {
2017   for (int i = 0; i < mat.size(); i++)
2018   {
2019     if (mat[i][col] != 0)
2020       return FALSE;
2021   }
2022   return TRUE;
2023 }
2024 
vvIsZero(const std::vector<std::vector<int>> & mat)2025 static BOOLEAN vvIsZero(const std::vector<std::vector<int> >& mat)
2026 {
2027   for (int i = 0; i < mat.size(); i++)
2028   {
2029     if (!vvIsRowZero(mat, i))
2030       return FALSE;
2031   }
2032   return TRUE;
2033 }
2034 
vvMult(const std::vector<std::vector<int>> & a,const std::vector<std::vector<int>> & b)2035 static std::vector<std::vector<int> > vvMult(const std::vector<std::vector<int> >& a, const std::vector<std::vector<int> >& b)
2036 {
2037   int ra = a.size();
2038   int rb = b.size();
2039   int ca = a.size() > 0 ? a[0].size() : 0;
2040   int cb = b.size() > 0 ? b[0].size() : 0;
2041 
2042   if (ca != rb)
2043   {
2044     WerrorS("matrix dimensions do not match");
2045     return std::vector<std::vector<int> >();
2046   }
2047 
2048   std::vector<std::vector<int> > res(ra, std::vector<int>(cb));
2049   for (int i = 0; i < ra; i++)
2050   {
2051     for (int j = 0; j < cb; j++)
2052     {
2053       int sum = 0;
2054       for (int k = 0; k < ca; k++)
2055         sum += a[i][k] * b[k][j];
2056       res[i][j] = sum;
2057     }
2058   }
2059   return res;
2060 }
2061 
isAcyclic(const intvec * G)2062 static BOOLEAN isAcyclic(const intvec* G)
2063 {
2064   // init
2065   int n = G->cols();
2066   std::vector<int> path;
2067   std::vector<BOOLEAN> visited;
2068   std::vector<BOOLEAN> cyclic;
2069   std::vector<int> cache;
2070   visited.resize(n, FALSE);
2071   cyclic.resize(n, FALSE);
2072   cache.resize(n, -2);
2073 
2074   for (int v = 0; v < n; v++)
2075   {
2076     cache = countCycles(G, v, path, visited, cyclic, cache);
2077     // check that there are 0 cycles from v
2078     if (cache[v] != 0)
2079       return FALSE;
2080   }
2081   return TRUE;
2082 }
2083 
2084 /*
2085  * Computation of the K-Dimension
2086  */
2087 
2088 // -1 is infinity, -2 is error
lp_kDim(const ideal _G)2089 int lp_kDim(const ideal _G)
2090 {
2091   if (rField_is_Ring(currRing)) {
2092       WerrorS("K-Dim not implemented for rings");
2093       return -2;
2094   }
2095 
2096   for (int i=IDELEMS(_G)-1;i>=0; i--)
2097   {
2098     if (_G->m[i] != NULL)
2099     {
2100       if (pGetComp(_G->m[i]) != 0)
2101       {
2102         WerrorS("K-Dim not implemented for modules");
2103         return -2;
2104       }
2105       if (pGetNCGen(_G->m[i]) != 0)
2106       {
2107         WerrorS("K-Dim not implemented for bi-modules");
2108         return -2;
2109       }
2110     }
2111   }
2112 
2113   ideal G = id_Head(_G, currRing); // G = LM(G) (and copy)
2114   if (TEST_OPT_PROT)
2115     Print("%d original generators\n", IDELEMS(G));
2116   idSkipZeroes(G); // remove zeros
2117   id_DelLmEquals(G, currRing); // remove duplicates
2118   if (TEST_OPT_PROT)
2119     Print("%d non-zero unique generators\n", IDELEMS(G));
2120 
2121   // check if G is the zero ideal
2122   if (IDELEMS(G) == 1 && G->m[0] == NULL)
2123   {
2124     // NOTE: this is needed because if the ideal is <0>, then idSkipZeroes keeps this element, and IDELEMS is still 1!
2125     int lV = currRing->isLPring;
2126     int ncGenCount = currRing->LPncGenCount;
2127     if (lV - ncGenCount == 0)
2128     {
2129       idDelete(&G);
2130       return 1;
2131     }
2132     if (lV - ncGenCount == 1)
2133     {
2134       idDelete(&G);
2135       return -1;
2136     }
2137     if (lV - ncGenCount >= 2)
2138     {
2139       idDelete(&G);
2140       return -1;
2141     }
2142   }
2143 
2144   // get the max deg
2145   long maxDeg = 0;
2146   for (int i = 0; i < IDELEMS(G); i++)
2147   {
2148     maxDeg = si_max(maxDeg, pTotaldegree(G->m[i]));
2149 
2150     // also check whether G = <1>
2151     if (pIsConstantComp(G->m[i]))
2152     {
2153       WerrorS("K-Dim not defined for 0-ring"); // TODO is it minus infinity ?
2154       idDelete(&G);
2155       return -2;
2156     }
2157   }
2158   if (TEST_OPT_PROT)
2159     Print("max deg: %ld\n", maxDeg);
2160 
2161 
2162   // for normal words of length minDeg ... maxDeg-1
2163   // brute-force the normal words
2164   if (TEST_OPT_PROT)
2165     PrintS("Computing normal words normally...\n");
2166   long numberOfNormalWords = lp_countNormalWords(maxDeg - 1, G);
2167 
2168   if (TEST_OPT_PROT)
2169     Print("%ld normal words up to length %ld\n", numberOfNormalWords, maxDeg - 1);
2170 
2171   // early termination if G \subset X
2172   if (maxDeg <= 1)
2173   {
2174     int lV = currRing->isLPring;
2175     int ncGenCount = currRing->LPncGenCount;
2176     if (IDELEMS(G) == lV - ncGenCount) // V = {1} no edges
2177     {
2178       idDelete(&G);
2179       return numberOfNormalWords;
2180     }
2181     if (IDELEMS(G) == lV - ncGenCount - 1) // V = {1} with loop
2182     {
2183       idDelete(&G);
2184       return -1;
2185     }
2186     if (IDELEMS(G) <= lV - ncGenCount - 2) // V = {1} with more than one loop
2187     {
2188       idDelete(&G);
2189       return -1;
2190     }
2191   }
2192 
2193   if (TEST_OPT_PROT)
2194     PrintS("Computing Ufnarovski graph...\n");
2195 
2196   ideal standardWords;
2197   intvec* UG = lp_ufnarovskiGraph(G, standardWords);
2198   if (UG == NULL)
2199   {
2200     idDelete(&G);
2201     return -2;
2202   }
2203   if (errorreported)
2204   {
2205     delete UG;
2206     idDelete(&G);
2207     return -2;
2208   }
2209 
2210   if (TEST_OPT_PROT)
2211     Print("Ufnarovski graph is %dx%d.\n", UG->rows(), UG->cols());
2212 
2213   if (TEST_OPT_PROT)
2214     PrintS("Checking whether Ufnarovski graph is acyclic...\n");
2215 
2216   if (!isAcyclic(UG))
2217   {
2218     // in this case we have infinitely many normal words
2219     return -1;
2220   }
2221 
2222   std::vector<std::vector<int> > vvUG = iv2vv(UG);
2223   for (int i = 0; i < vvUG.size(); i++)
2224   {
2225     if (vvIsRowZero(vvUG, i) && vvIsColumnZero(vvUG, i)) // i is isolated vertex
2226     {
2227       vvDeleteRow(vvUG, i);
2228       vvDeleteColumn(vvUG, i);
2229       i--;
2230     }
2231   }
2232   if (TEST_OPT_PROT)
2233     Print("Simplified Ufnarovski graph to %dx%d.\n", (int)vvUG.size(), (int)vvUG.size());
2234 
2235   // for normal words of length >= maxDeg
2236   // use Ufnarovski graph
2237   if (TEST_OPT_PROT)
2238     PrintS("Computing normal words via Ufnarovski graph...\n");
2239   std::vector<std::vector<int> > UGpower = vvUG;
2240   long nUGpower = 1;
2241   while (!vvIsZero(UGpower))
2242   {
2243     if (TEST_OPT_PROT)
2244       PrintS("Start count graph entries.\n");
2245     for (int i = 0; i < UGpower.size(); i++)
2246     {
2247       for (int j = 0; j < UGpower[i].size(); j++)
2248       {
2249         numberOfNormalWords += UGpower[i][j];
2250       }
2251     }
2252 
2253     if (TEST_OPT_PROT)
2254     {
2255       PrintS("Done count graph entries.\n");
2256       Print("%ld normal words up to length %ld\n", numberOfNormalWords, maxDeg - 1 + nUGpower);
2257     }
2258 
2259     if (TEST_OPT_PROT)
2260       PrintS("Start mat mult.\n");
2261     UGpower = vvMult(UGpower, vvUG); // TODO: avoid creation of new intvec
2262     if (TEST_OPT_PROT)
2263       PrintS("Done mat mult.\n");
2264     nUGpower++;
2265   }
2266 
2267   delete UG;
2268   idDelete(&G);
2269   return numberOfNormalWords;
2270 }
2271 #endif
2272