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