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