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