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