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 "misc/options.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  if (errorreported) return;
904  l = hseries1->length()-1;
905  if (l > 1)
906  hseries2 = hSecondSeries(hseries1);
907  else
908  hseries2 = hseries1;
909  hDegreeSeries(hseries1, hseries2, &co, &mu);
910  if ((l == 1) &&(mu == 0))
911  scPrintDegree((currRing->N)+1, 0);
912  else
913  scPrintDegree(co, mu);
914  if (l>1)
915  delete hseries1;
916  delete hseries2;
917 }
918 
919 static void hDegree0(ideal S, ideal Q, const ring tailRing)
920 {
921  id_TestTail(S, currRing, tailRing);
922  if (Q!=NULL) id_TestTail(Q, currRing, tailRing);
923 
924  int mc;
925  hexist = hInit(S, Q, &hNexist, tailRing);
926  if (!hNexist)
927  {
928  hMu = -1;
929  return;
930  }
931  else
932  hMu = 0;
933 
934  const ring r = currRing;
935 
936  hwork = (scfmon)omAlloc(hNexist * sizeof(scmon));
937  hvar = (varset)omAlloc(((r->N) + 1) * sizeof(int));
938  hpur0 = (scmon)omAlloc((1 + ((r->N) * (r->N))) * sizeof(int));
939  mc = hisModule;
940  if (!mc)
941  {
942  hstc = hexist;
943  hNstc = hNexist;
944  }
945  else
946  hstc = (scfmon)omAlloc(hNexist * sizeof(scmon));
947  stcmem = hCreate((r->N) - 1);
948  loop
949  {
950  if (mc)
951  {
952  hComp(hexist, hNexist, mc, hstc, &hNstc);
953  if (!hNstc)
954  {
955  hMu = -1;
956  break;
957  }
958  }
959  hNvar = (r->N);
960  for (int i = hNvar; i; i--)
961  hvar[i] = i;
963  hSupp(hstc, hNstc, hvar, &hNvar);
964  if ((hNvar == (r->N)) && (hNstc >= (r->N)))
965  {
966  if ((hNvar > 2) && (hNstc > 10))
968  memset(hpur0, 0, ((r->N) + 1) * sizeof(int));
969  hPure(hstc, 0, &hNstc, hvar, hNvar, hpur0, &hNpure);
970  if (hNpure == hNvar)
971  {
972  hLexS(hstc, hNstc, hvar, hNvar);
974  }
975  else
976  hMu = -1;
977  }
978  else if (hNvar)
979  hMu = -1;
980  mc--;
981  if (mc <= 0 || hMu < 0)
982  break;
983  }
984  hKill(stcmem, (r->N) - 1);
985  omFreeSize((ADDRESS)hpur0, (1 + ((r->N) * (r->N))) * sizeof(int));
986  omFreeSize((ADDRESS)hvar, ((r->N) + 1) * sizeof(int));
987  omFreeSize((ADDRESS)hwork, hNexist * sizeof(scmon));
989  if (hisModule)
990  omFreeSize((ADDRESS)hstc, hNexist * sizeof(scmon));
991 }
992 
993 int scMult0Int(ideal S, ideal Q, const ring tailRing)
994 {
995  id_TestTail(S, currRing, tailRing);
996  if (Q!=NULL) id_TestTail(Q, currRing, tailRing);
997 
998  hDegree0(S, Q, tailRing);
999  return hMu;
1000 }
1001 
1002 
1003 // HC
1004 
1006 
1007 static void hHedge(poly hEdge)
1008 {
1009  pSetm(pWork);
1010  if (pLmCmp(pWork, hEdge) == currRing->OrdSgn)
1011  {
1012  for (int i = hNvar; i>0; i--)
1013  pSetExp(hEdge,i, pGetExp(pWork,i));
1014  pSetm(hEdge);
1015  }
1016 }
1017 
1018 
1019 static void hHedgeStep(scmon pure, scfmon stc,
1020  int Nstc, varset var, int Nvar,poly hEdge)
1021 {
1022  int iv = Nvar -1, k = var[Nvar], a, a0, a1, b, i;
1023  int x/*, x0*/;
1024  scmon pn;
1025  scfmon sn;
1026  if (iv==0)
1027  {
1028  pSetExp(pWork, k, pure[k]);
1029  hHedge(hEdge);
1030  return;
1031  }
1032  else if (Nstc==0)
1033  {
1034  for (i = Nvar; i>0; i--)
1035  pSetExp(pWork, var[i], pure[var[i]]);
1036  hHedge(hEdge);
1037  return;
1038  }
1039  x = a = 0;
1040  pn = hGetpure(pure);
1041  sn = hGetmem(Nstc, stc, stcmem[iv]);
1042  hStepS(sn, Nstc, var, Nvar, &a, &x);
1043  if (a == Nstc)
1044  {
1045  pSetExp(pWork, k, pure[k]);
1046  hHedgeStep(pn, sn, a, var, iv,hEdge);
1047  return;
1048  }
1049  else
1050  {
1051  pSetExp(pWork, k, x);
1052  hHedgeStep(pn, sn, a, var, iv,hEdge);
1053  }
1054  b = a;
1055  loop
1056  {
1057  a0 = a;
1058  // x0 = x;
1059  hStepS(sn, Nstc, var, Nvar, &a, &x);
1060  hElimS(sn, &b, a0, a, var, iv);
1061  a1 = a;
1062  hPure(sn, a0, &a1, var, iv, pn, &i);
1063  hLex2S(sn, b, a0, a1, var, iv, hwork);
1064  b += (a1 - a0);
1065  if (a < Nstc)
1066  {
1067  pSetExp(pWork, k, x);
1068  hHedgeStep(pn, sn, b, var, iv,hEdge);
1069  }
1070  else
1071  {
1072  pSetExp(pWork, k, pure[k]);
1073  hHedgeStep(pn, sn, b, var, iv,hEdge);
1074  return;
1075  }
1076  }
1077 }
1078 
1079 void scComputeHC(ideal S, ideal Q, int ak, poly &hEdge, ring tailRing)
1080 {
1081  id_TestTail(S, currRing, tailRing);
1082  if (Q!=NULL) id_TestTail(Q, currRing, tailRing);
1083 
1084  int i;
1085  int k = ak;
1086  #ifdef HAVE_RINGS
1087  if (rField_is_Ring(currRing) && (currRing->OrdSgn == -1))
1088  {
1089  //consider just monic generators (over rings with zero-divisors)
1090  ideal SS=id_Copy(S,tailRing);
1091  for(i=0;i<=idElem(S);i++)
1092  {
1093  if((SS->m[i]!=NULL)
1094  && ((p_IsPurePower(SS->m[i],tailRing)==0)
1095  ||(!n_IsUnit(pGetCoeff(SS->m[i]), tailRing->cf))))
1096  {
1097  p_Delete(&SS->m[i],tailRing);
1098  }
1099  }
1100  S=id_Copy(SS,tailRing);
1101  idSkipZeroes(S);
1102  }
1103  #if 0
1104  printf("\nThis is HC:\n");
1105  for(int ii=0;ii<=idElem(S);ii++)
1106  {
1107  pWrite(S->m[ii]);
1108  }
1109  //getchar();
1110  #endif
1111  #endif
1112  if(idElem(S) == 0)
1113  return;
1114  hNvar = (currRing->N);
1115  hexist = hInit(S, Q, &hNexist, tailRing); // tailRing?
1116  if (k!=0)
1117  hComp(hexist, hNexist, k, hexist, &hNstc);
1118  else
1119  hNstc = hNexist;
1120  assume(hNexist > 0);
1121  hwork = (scfmon)omAlloc(hNexist * sizeof(scmon));
1122  hvar = (varset)omAlloc((hNvar + 1) * sizeof(int));
1123  hpure = (scmon)omAlloc((1 + (hNvar * hNvar)) * sizeof(int));
1124  stcmem = hCreate(hNvar - 1);
1125  for (i = hNvar; i>0; i--)
1126  hvar[i] = i;
1128  if ((hNvar > 2) && (hNstc > 10))
1130  memset(hpure, 0, (hNvar + 1) * sizeof(int));
1131  hPure(hexist, 0, &hNstc, hvar, hNvar, hpure, &hNpure);
1132  hLexS(hexist, hNstc, hvar, hNvar);
1133  if (hEdge!=NULL)
1134  pLmFree(hEdge);
1135  hEdge = pInit();
1136  pWork = pInit();
1137  hHedgeStep(hpure, hexist, hNstc, hvar, hNvar,hEdge);
1138  pSetComp(hEdge,ak);
1139  hKill(stcmem, hNvar - 1);
1140  omFreeSize((ADDRESS)hwork, hNexist * sizeof(scmon));
1141  omFreeSize((ADDRESS)hvar, (hNvar + 1) * sizeof(int));
1142  omFreeSize((ADDRESS)hpure, (1 + (hNvar * hNvar)) * sizeof(int));
1144  pLmFree(pWork);
1145 }
1146 
1147 
1148 
1149 // kbase
1150 
1153 
1154 static void scElKbase()
1155 {
1156  poly q = pInit();
1157  pSetCoeff0(q,nInit(1));
1158  pSetExpV(q,act);
1159  pNext(q) = NULL;
1160  last = pNext(last) = q;
1161 }
1162 
1163 static int scMax( int i, scfmon stc, int Nvar)
1164 {
1165  int x, y=stc[0][Nvar];
1166  for (; i;)
1167  {
1168  i--;
1169  x = stc[i][Nvar];
1170  if (x > y) y = x;
1171  }
1172  return y;
1173 }
1174 
1175 static int scMin( int i, scfmon stc, int Nvar)
1176 {
1177  int x, y=stc[0][Nvar];
1178  for (; i;)
1179  {
1180  i--;
1181  x = stc[i][Nvar];
1182  if (x < y) y = x;
1183  }
1184  return y;
1185 }
1186 
1187 static int scRestrict( int &Nstc, scfmon stc, int Nvar)
1188 {
1189  int x, y;
1190  int i, j, Istc = Nstc;
1191 
1192  y = MAX_INT_VAL;
1193  for (i=Nstc-1; i>=0; i--)
1194  {
1195  j = Nvar-1;
1196  loop
1197  {
1198  if(stc[i][j] != 0) break;
1199  j--;
1200  if (j == 0)
1201  {
1202  Istc--;
1203  x = stc[i][Nvar];
1204  if (x < y) y = x;
1205  stc[i] = NULL;
1206  break;
1207  }
1208  }
1209  }
1210  if (Istc < Nstc)
1211  {
1212  for (i=Nstc-1; i>=0; i--)
1213  {
1214  if (stc[i] && (stc[i][Nvar] >= y))
1215  {
1216  Istc--;
1217  stc[i] = NULL;
1218  }
1219  }
1220  j = 0;
1221  while (stc[j]) j++;
1222  i = j+1;
1223  for(; i<Nstc; i++)
1224  {
1225  if (stc[i])
1226  {
1227  stc[j] = stc[i];
1228  j++;
1229  }
1230  }
1231  Nstc = Istc;
1232  return y;
1233  }
1234  else
1235  return -1;
1236 }
1237 
1238 static void scAll( int Nvar, int deg)
1239 {
1240  int i;
1241  int d = deg;
1242  if (d == 0)
1243  {
1244  for (i=Nvar; i; i--) act[i] = 0;
1245  scElKbase();
1246  return;
1247  }
1248  if (Nvar == 1)
1249  {
1250  act[1] = d;
1251  scElKbase();
1252  return;
1253  }
1254  do
1255  {
1256  act[Nvar] = d;
1257  scAll(Nvar-1, deg-d);
1258  d--;
1259  } while (d >= 0);
1260 }
1261 
1262 static void scAllKbase( int Nvar, int ideg, int deg)
1263 {
1264  do
1265  {
1266  act[Nvar] = ideg;
1267  scAll(Nvar-1, deg-ideg);
1268  ideg--;
1269  } while (ideg >= 0);
1270 }
1271 
1272 static void scDegKbase( scfmon stc, int Nstc, int Nvar, int deg)
1273 {
1274  int Ivar, Istc, i, j;
1275  scfmon sn;
1276  int x, ideg;
1277 
1278  if (deg == 0)
1279  {
1280  for (i=Nstc-1; i>=0; i--)
1281  {
1282  for (j=Nvar;j;j--){ if(stc[i][j]) break; }
1283  if (j==0){return;}
1284  }
1285  for (i=Nvar; i; i--) act[i] = 0;
1286  scElKbase();
1287  return;
1288  }
1289  if (Nvar == 1)
1290  {
1291  for (i=Nstc-1; i>=0; i--) if(deg >= stc[i][1]) return;
1292  act[1] = deg;
1293  scElKbase();
1294  return;
1295  }
1296  Ivar = Nvar-1;
1297  sn = hGetmem(Nstc, stc, stcmem[Ivar]);
1298  x = scRestrict(Nstc, sn, Nvar);
1299  if (x <= 0)
1300  {
1301  if (x == 0) return;
1302  ideg = deg;
1303  }
1304  else
1305  {
1306  if (deg < x) ideg = deg;
1307  else ideg = x-1;
1308  if (Nstc == 0)
1309  {
1310  scAllKbase(Nvar, ideg, deg);
1311  return;
1312  }
1313  }
1314  loop
1315  {
1316  x = scMax(Nstc, sn, Nvar);
1317  while (ideg >= x)
1318  {
1319  act[Nvar] = ideg;
1320  scDegKbase(sn, Nstc, Ivar, deg-ideg);
1321  ideg--;
1322  }
1323  if (ideg < 0) return;
1324  Istc = Nstc;
1325  for (i=Nstc-1; i>=0; i--)
1326  {
1327  if (ideg < sn[i][Nvar])
1328  {
1329  Istc--;
1330  sn[i] = NULL;
1331  }
1332  }
1333  if (Istc == 0)
1334  {
1335  scAllKbase(Nvar, ideg, deg);
1336  return;
1337  }
1338  j = 0;
1339  while (sn[j]) j++;
1340  i = j+1;
1341  for (; i<Nstc; i++)
1342  {
1343  if (sn[i])
1344  {
1345  sn[j] = sn[i];
1346  j++;
1347  }
1348  }
1349  Nstc = Istc;
1350  }
1351 }
1352 
1353 static void scInKbase( scfmon stc, int Nstc, int Nvar)
1354 {
1355  int Ivar, Istc, i, j;
1356  scfmon sn;
1357  int x, ideg;
1358 
1359  if (Nvar == 1)
1360  {
1361  ideg = scMin(Nstc, stc, 1);
1362  while (ideg > 0)
1363  {
1364  ideg--;
1365  act[1] = ideg;
1366  scElKbase();
1367  }
1368  return;
1369  }
1370  Ivar = Nvar-1;
1371  sn = hGetmem(Nstc, stc, stcmem[Ivar]);
1372  x = scRestrict(Nstc, sn, Nvar);
1373  if (x == 0) return;
1374  ideg = x-1;
1375  loop
1376  {
1377  x = scMax(Nstc, sn, Nvar);
1378  while (ideg >= x)
1379  {
1380  act[Nvar] = ideg;
1381  scInKbase(sn, Nstc, Ivar);
1382  ideg--;
1383  }
1384  if (ideg < 0) return;
1385  Istc = Nstc;
1386  for (i=Nstc-1; i>=0; i--)
1387  {
1388  if (ideg < sn[i][Nvar])
1389  {
1390  Istc--;
1391  sn[i] = NULL;
1392  }
1393  }
1394  j = 0;
1395  while (sn[j]) j++;
1396  i = j+1;
1397  for (; i<Nstc; i++)
1398  {
1399  if (sn[i])
1400  {
1401  sn[j] = sn[i];
1402  j++;
1403  }
1404  }
1405  Nstc = Istc;
1406  }
1407 }
1408 
1409 static ideal scIdKbase(poly q, const int rank)
1410 {
1411  ideal res = idInit(pLength(q), rank);
1412  polyset mm = res->m;
1413  do
1414  {
1415  *mm = q; ++mm;
1416 
1417  const poly p = pNext(q);
1418  pNext(q) = NULL;
1419  q = p;
1420 
1421  } while (q!=NULL);
1422 
1423  id_Test(res, currRing); // WRONG RANK!!!???
1424  return res;
1425 }
1426 
1427 ideal scKBase(int deg, ideal s, ideal Q, intvec * mv)
1428 {
1429  if( Q!=NULL) id_Test(Q, currRing);
1430 
1431  int i, di;
1432  poly p;
1433 
1434  if (deg < 0)
1435  {
1436  di = scDimInt(s, Q);
1437  if (di != 0)
1438  {
1439  //Werror("KBase not finite");
1440  return idInit(1,s->rank);
1441  }
1442  }
1443  stcmem = hCreate((currRing->N) - 1);
1444  hexist = hInit(s, Q, &hNexist, currRing);
1445  p = last = pInit();
1446  /*pNext(p) = NULL;*/
1447  act = (scmon)omAlloc(((currRing->N) + 1) * sizeof(int));
1448  *act = 0;
1449  if (!hNexist)
1450  {
1451  scAll((currRing->N), deg);
1452  goto ende;
1453  }
1454  if (!hisModule)
1455  {
1456  if (deg < 0) scInKbase(hexist, hNexist, (currRing->N));
1457  else scDegKbase(hexist, hNexist, (currRing->N), deg);
1458  }
1459  else
1460  {
1461  hstc = (scfmon)omAlloc(hNexist * sizeof(scmon));
1462  for (i = 1; i <= hisModule; i++)
1463  {
1464  *act = i;
1465  hComp(hexist, hNexist, i, hstc, &hNstc);
1466  int deg_ei=deg;
1467  if (mv!=NULL) deg_ei -= (*mv)[i-1];
1468  if ((deg < 0) || (deg_ei>=0))
1469  {
1470  if (hNstc)
1471  {
1472  if (deg < 0) scInKbase(hstc, hNstc, (currRing->N));
1473  else scDegKbase(hstc, hNstc, (currRing->N), deg_ei);
1474  }
1475  else
1476  scAll((currRing->N), deg_ei);
1477  }
1478  }
1479  omFreeSize((ADDRESS)hstc, hNexist * sizeof(scmon));
1480  }
1481 ende:
1483  omFreeSize((ADDRESS)act, ((currRing->N) + 1) * sizeof(int));
1484  hKill(stcmem, (currRing->N) - 1);
1485  pLmFree(&p);
1486  if (p == NULL)
1487  return idInit(1,s->rank);
1488 
1489  last = p;
1490  return scIdKbase(p, s->rank);
1491 }
1492 
1493 #if 0 //-- alternative implementation of scComputeHC
1494 /*
1495 void scComputeHCw(ideal ss, ideal Q, int ak, poly &hEdge, ring tailRing)
1496 {
1497  id_TestTail(ss, currRing, tailRing);
1498  if (Q!=NULL) id_TestTail(Q, currRing, tailRing);
1499 
1500  int i, di;
1501  poly p;
1502 
1503  if (hEdge!=NULL)
1504  pLmFree(hEdge);
1505 
1506  ideal s=idInit(IDELEMS(ss),ak);
1507  for(i=IDELEMS(ss)-1;i>=0;i--)
1508  {
1509  if (ss->m[i]!=NULL) s->m[i]=pHead(ss->m[i]);
1510  }
1511  di = scDimInt(s, Q);
1512  stcmem = hCreate((currRing->N) - 1);
1513  hexist = hInit(s, Q, &hNexist, currRing);
1514  p = last = pInit();
1515  // pNext(p) = NULL;
1516  act = (scmon)omAlloc(((currRing->N) + 1) * sizeof(int));
1517  *act = 0;
1518  if (!hNexist)
1519  {
1520  scAll((currRing->N), -1);
1521  goto ende;
1522  }
1523  if (!hisModule)
1524  {
1525  scInKbase(hexist, hNexist, (currRing->N));
1526  }
1527  else
1528  {
1529  hstc = (scfmon)omAlloc(hNexist * sizeof(scmon));
1530  for (i = 1; i <= hisModule; i++)
1531  {
1532  *act = i;
1533  hComp(hexist, hNexist, i, hstc, &hNstc);
1534  if (hNstc)
1535  {
1536  scInKbase(hstc, hNstc, (currRing->N));
1537  }
1538  else
1539  scAll((currRing->N), -1);
1540  }
1541  omFreeSize((ADDRESS)hstc, hNexist * sizeof(scmon));
1542  }
1543 ende:
1544  hDelete(hexist, hNexist);
1545  omFreeSize((ADDRESS)act, ((currRing->N) + 1) * sizeof(int));
1546  hKill(stcmem, (currRing->N) - 1);
1547  pDeleteLm(&p);
1548  idDelete(&s);
1549  if (p == NULL)
1550  {
1551  return; // no HEdge
1552  }
1553  else
1554  {
1555  last = p;
1556  ideal res=scIdKbase(p, ss->rank);
1557  poly p_ind=res->m[0]; int ind=0;
1558  for(i=IDELEMS(res)-1;i>0;i--)
1559  {
1560  if (pCmp(res->m[i],p_ind)==-1) { p_ind=res->m[i]; ind=i; }
1561  }
1562  assume(p_ind!=NULL);
1563  assume(res->m[ind]==p_ind);
1564  hEdge=p_ind;
1565  res->m[ind]=NULL;
1566  nDelete(&pGetCoeff(hEdge));
1567  pGetCoeff(hEdge)=NULL;
1568  for(i=(currRing->N);i>0;i--)
1569  pIncrExp(hEdge,i);
1570  pSetm(hEdge);
1571 
1572  idDelete(&res);
1573  return;
1574  }
1575 }
1576  */
1577 #endif
1578 
1579 #ifdef HAVE_SHIFTBBA
1580 
1581 /*
1582  * Computation of the Gel'fand-Kirillov Dimension
1583  */
1584 
1585 #include "polys/shiftop.h"
1586 #include <vector>
1587 
1588 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)
1589 {
1590  intvec* G = ivCopy(_G); // modifications must be local
1591 
1592  if (cache[v] != -2) return cache; // value is already cached
1593 
1594  visited[v] = TRUE;
1595  path.push_back(v);
1596 
1597  int cycles = 0;
1598  for (int w = 0; w < G->cols(); w++)
1599  {
1600  if (IMATELEM(*G, v + 1, w + 1)) // edge v -> w exists in G
1601  {
1602  if (!visited[w])
1603  { // continue with w
1604  cache = countCycles(G, w, path, visited, cyclic, cache);
1605  if (cache[w] == -1)
1606  {
1607  cache[v] = -1;
1608  return cache;
1609  }
1610  cycles = si_max(cycles, cache[w]);
1611  }
1612  else
1613  { // found new cycle
1614  int pathIndexOfW = -1;
1615  for (int i = path.size() - 1; i >= 0; i--) {
1616  if (cyclic[path[i]] == 1) { // found an already cyclic vertex
1617  cache[v] = -1;
1618  return cache;
1619  }
1620  cyclic[path[i]] = TRUE;
1621 
1622  if (path[i] == w) { // end of the cycle
1623  assume(IMATELEM(*G, v + 1, w + 1) != 0);
1624  IMATELEM(*G, v + 1, w + 1) = 0; // remove edge v -> w
1625  pathIndexOfW = i;
1626  break;
1627  } else {
1628  assume(IMATELEM(*G, path[i - 1] + 1, path[i] + 1) != 0);
1629  IMATELEM(*G, path[i - 1] + 1, path[i] + 1) = 0; // remove edge vi-1 -> vi
1630  }
1631  }
1632  assume(pathIndexOfW != -1); // should never happen
1633  for (int i = path.size() - 1; i >= pathIndexOfW; i--) {
1634  cache = countCycles(G, path[i], path, visited, cyclic, cache);
1635  if (cache[path[i]] == -1)
1636  {
1637  cache[v] = -1;
1638  return cache;
1639  }
1640  cycles = si_max(cycles, cache[path[i]] + 1);
1641  }
1642  }
1643  }
1644  }
1645  cache[v] = cycles;
1646 
1647  delete G;
1648  return cache;
1649 }
1650 
1651 // -1 is infinity
1652 static int graphGrowth(const intvec* G)
1653 {
1654  // init
1655  int n = G->cols();
1656  std::vector<int> path;
1657  std::vector<BOOLEAN> visited;
1658  std::vector<BOOLEAN> cyclic;
1659  std::vector<int> cache;
1660  visited.resize(n, FALSE);
1661  cyclic.resize(n, FALSE);
1662  cache.resize(n, -2);
1663 
1664  // get max number of cycles
1665  int cycles = 0;
1666  for (int v = 0; v < n; v++)
1667  {
1668  cache = countCycles(G, v, path, visited, cyclic, cache);
1669  if (cache[v] == -1)
1670  return -1;
1671  cycles = si_max(cycles, cache[v]);
1672  }
1673  return cycles;
1674 }
1675 
1676 // ATTENTION:
1677 // - `words` contains the words normal modulo M of length n
1678 // - `numberOfNormalWords` contains the number of words normal modulo M of length 0 ... n
1679 static void _lp_computeNormalWords(ideal words, int& numberOfNormalWords, int length, ideal M, int minDeg, int& last)
1680 {
1681  if (length <= 0){
1682  poly one = pOne();
1683  if (p_LPDivisibleBy(M, one, currRing)) // 1 \in M => no normal words at all
1684  {
1685  pDelete(&one);
1686  last = -1;
1687  numberOfNormalWords = 0;
1688  }
1689  else
1690  {
1691  words->m[0] = one;
1692  last = 0;
1693  numberOfNormalWords = 1;
1694  }
1695  return;
1696  }
1697 
1698  _lp_computeNormalWords(words, numberOfNormalWords, length - 1, M, minDeg, last);
1699 
1700  int nVars = currRing->isLPring - currRing->LPncGenCount;
1701  int numberOfNewNormalWords = 0;
1702 
1703  for (int j = nVars - 1; j >= 0; j--)
1704  {
1705  for (int i = last; i >= 0; i--)
1706  {
1707  int index = (j * (last + 1)) + i;
1708 
1709  if (words->m[i] != NULL)
1710  {
1711  if (j > 0) {
1712  words->m[index] = pCopy(words->m[i]);
1713  }
1714 
1715  int varOffset = ((length - 1) * currRing->isLPring) + 1;
1716  pSetExp(words->m[index], varOffset + j, 1);
1717  pSetm(words->m[index]);
1718  pTest(words->m[index]);
1719 
1720  if (length >= minDeg && p_LPDivisibleBy(M, words->m[index], currRing))
1721  {
1722  pDelete(&words->m[index]);
1723  words->m[index] = NULL;
1724  }
1725  else
1726  {
1727  numberOfNewNormalWords++;
1728  }
1729  }
1730  }
1731  }
1732 
1733  last = nVars * last + nVars - 1;
1734 
1735  numberOfNormalWords += numberOfNewNormalWords;
1736 }
1737 
1738 static ideal lp_computeNormalWords(int length, ideal M)
1739 {
1740  long minDeg = IDELEMS(M) > 0 ? pTotaldegree(M->m[0]) : 0;
1741  for (int i = 1; i < IDELEMS(M); i++)
1742  {
1743  minDeg = si_min(minDeg, pTotaldegree(M->m[i]));
1744  }
1745 
1746  int nVars = currRing->isLPring - currRing->LPncGenCount;
1747 
1748  int maxElems = 1;
1749  for (int i = 0; i < length; i++) // maxElems = nVars^n
1750  maxElems *= nVars;
1751  ideal words = idInit(maxElems);
1752  int last, numberOfNormalWords;
1753  _lp_computeNormalWords(words, numberOfNormalWords, length, M, minDeg, last);
1754  idSkipZeroes(words);
1755  return words;
1756 }
1757 
1758 static int lp_countNormalWords(int upToLength, ideal M)
1759 {
1760  long minDeg = IDELEMS(M) > 0 ? pTotaldegree(M->m[0]) : 0;
1761  for (int i = 1; i < IDELEMS(M); i++)
1762  {
1763  minDeg = si_min(minDeg, pTotaldegree(M->m[i]));
1764  }
1765 
1766  int nVars = currRing->isLPring - currRing->LPncGenCount;
1767 
1768  int maxElems = 1;
1769  for (int i = 0; i < upToLength; i++) // maxElems = nVars^n
1770  maxElems *= nVars;
1771  ideal words = idInit(maxElems);
1772  int last, numberOfNormalWords;
1773  _lp_computeNormalWords(words, numberOfNormalWords, upToLength, M, minDeg, last);
1774  idDelete(&words);
1775  return numberOfNormalWords;
1776 }
1777 
1778 // NULL if graph is undefined
1779 intvec* lp_ufnarovskiGraph(ideal G, ideal &standardWords)
1780 {
1781  long l = 0;
1782  for (int i = 0; i < IDELEMS(G); i++)
1783  l = si_max(pTotaldegree(G->m[i]), l);
1784  l--;
1785  if (l <= 0)
1786  {
1787  WerrorS("Ufnarovski graph not implemented for l <= 0");
1788  return NULL;
1789  }
1790  int lV = currRing->isLPring;
1791 
1792  standardWords = lp_computeNormalWords(l, G);
1793 
1794  int n = IDELEMS(standardWords);
1795  intvec* UG = new intvec(n, n, 0);
1796  for (int i = 0; i < n; i++)
1797  {
1798  for (int j = 0; j < n; j++)
1799  {
1800  poly v = standardWords->m[i];
1801  poly w = standardWords->m[j];
1802 
1803  // check whether v*x1 = x2*w (overlap)
1804  bool overlap = true;
1805  for (int k = 1; k <= (l - 1) * lV; k++)
1806  {
1807  if (pGetExp(v, k + lV) != pGetExp(w, k)) {
1808  overlap = false;
1809  break;
1810  }
1811  }
1812 
1813  if (overlap)
1814  {
1815  // create the overlap
1816  poly p = pMult(pCopy(v), p_LPVarAt(w, l, currRing));
1817 
1818  // check whether the overlap is normal
1819  bool normal = true;
1820  for (int k = 0; k < IDELEMS(G); k++)
1821  {
1822  if (p_LPDivisibleBy(G->m[k], p, currRing))
1823  {
1824  normal = false;
1825  break;
1826  }
1827  }
1828 
1829  if (normal)
1830  {
1831  IMATELEM(*UG, i + 1, j + 1) = 1;
1832  }
1833  }
1834  }
1835  }
1836  return UG;
1837 }
1838 
1839 // -1 is infinity, -2 is error
1840 int lp_gkDim(const ideal _G)
1841 {
1842  id_Test(_G, currRing);
1843 
1844  if (rField_is_Ring(currRing)) {
1845  WerrorS("GK-Dim not implemented for rings");
1846  return -2;
1847  }
1848 
1849  for (int i=IDELEMS(_G)-1;i>=0; i--)
1850  {
1851  if (_G->m[i] != NULL)
1852  {
1853  if (pGetComp(_G->m[i]) != 0)
1854  {
1855  WerrorS("GK-Dim not implemented for modules");
1856  return -2;
1857  }
1858  if (pGetNCGen(_G->m[i]) != 0)
1859  {
1860  WerrorS("GK-Dim not implemented for bi-modules");
1861  return -2;
1862  }
1863  }
1864  }
1865 
1866  ideal G = id_Head(_G, currRing); // G = LM(G) (and copy)
1867  idSkipZeroes(G); // remove zeros
1868  id_DelLmEquals(G, currRing); // remove duplicates
1869 
1870  // check if G is the zero ideal
1871  if (IDELEMS(G) == 1 && G->m[0] == NULL)
1872  {
1873  // NOTE: this is needed because if the ideal is <0>, then idSkipZeroes keeps this element, and IDELEMS is still 1!
1874  int lV = currRing->isLPring;
1875  int ncGenCount = currRing->LPncGenCount;
1876  if (lV - ncGenCount == 0)
1877  {
1878  idDelete(&G);
1879  return 0;
1880  }
1881  if (lV - ncGenCount == 1)
1882  {
1883  idDelete(&G);
1884  return 1;
1885  }
1886  if (lV - ncGenCount >= 2)
1887  {
1888  idDelete(&G);
1889  return -1;
1890  }
1891  }
1892 
1893  // get the max deg
1894  long maxDeg = 0;
1895  for (int i = 0; i < IDELEMS(G); i++)
1896  {
1897  maxDeg = si_max(maxDeg, pTotaldegree(G->m[i]));
1898 
1899  // also check whether G = <1>
1900  if (pIsConstantComp(G->m[i]))
1901  {
1902  WerrorS("GK-Dim not defined for 0-ring");
1903  idDelete(&G);
1904  return -2;
1905  }
1906  }
1907 
1908  // early termination if G \subset X
1909  if (maxDeg <= 1)
1910  {
1911  int lV = currRing->isLPring;
1912  int ncGenCount = currRing->LPncGenCount;
1913  if (IDELEMS(G) == lV - ncGenCount) // V = {1} no edges
1914  {
1915  idDelete(&G);
1916  return 0;
1917  }
1918  if (IDELEMS(G) == lV - ncGenCount - 1) // V = {1} with loop
1919  {
1920  idDelete(&G);
1921  return 1;
1922  }
1923  if (IDELEMS(G) <= lV - ncGenCount - 2) // V = {1} with more than one loop
1924  {
1925  idDelete(&G);
1926  return -1;
1927  }
1928  }
1929 
1930  ideal standardWords;
1931  intvec* UG = lp_ufnarovskiGraph(G, standardWords);
1932  if (UG == NULL)
1933  {
1934  idDelete(&G);
1935  return -2;
1936  }
1937  if (errorreported)
1938  {
1939  delete UG;
1940  idDelete(&G);
1941  return -2;
1942  }
1943  int gkDim = graphGrowth(UG);
1944  delete UG;
1945  idDelete(&G);
1946  return gkDim;
1947 }
1948 
1949 // converts an intvec matrix to a vector<vector<int> >
1950 static std::vector<std::vector<int> > iv2vv(intvec* M)
1951 {
1952  int rows = M->rows();
1953  int cols = M->cols();
1954 
1955  std::vector<std::vector<int> > mat(rows, std::vector<int>(cols));
1956 
1957  for (int i = 0; i < rows; i++)
1958  {
1959  for (int j = 0; j < cols; j++)
1960  {
1961  mat[i][j] = IMATELEM(*M, i + 1, j + 1);
1962  }
1963  }
1964 
1965  return mat;
1966 }
1967 
1968 static void vvPrint(const std::vector<std::vector<int> >& mat)
1969 {
1970  for (int i = 0; i < mat.size(); i++)
1971  {
1972  for (int j = 0; j < mat[i].size(); j++)
1973  {
1974  Print("%d ", mat[i][j]);
1975  }
1976  PrintLn();
1977  }
1978 }
1979 
1980 static void vvTest(const std::vector<std::vector<int> >& mat)
1981 {
1982  if (mat.size() > 0)
1983  {
1984  int cols = mat[0].size();
1985  for (int i = 1; i < mat.size(); i++)
1986  {
1987  if (cols != mat[i].size())
1988  WerrorS("number of cols in matrix inconsistent");
1989  }
1990  }
1991 }
1992 
1993 static void vvDeleteRow(std::vector<std::vector<int> >& mat, int row)
1994 {
1995  mat.erase(mat.begin() + row);
1996 }
1997 
1998 static void vvDeleteColumn(std::vector<std::vector<int> >& mat, int col)
1999 {
2000  for (int i = 0; i < mat.size(); i++)
2001  {
2002  mat[i].erase(mat[i].begin() + col);
2003  }
2004 }
2005 
2006 static BOOLEAN vvIsRowZero(const std::vector<std::vector<int> >& mat, int row)
2007 {
2008  for (int i = 0; i < mat[row].size(); i++)
2009  {
2010  if (mat[row][i] != 0)
2011  return FALSE;
2012  }
2013  return TRUE;
2014 }
2015 
2016 static BOOLEAN vvIsColumnZero(const std::vector<std::vector<int> >& mat, int col)
2017 {
2018  for (int i = 0; i < mat.size(); i++)
2019  {
2020  if (mat[i][col] != 0)
2021  return FALSE;
2022  }
2023  return TRUE;
2024 }
2025 
2026 static BOOLEAN vvIsZero(const std::vector<std::vector<int> >& mat)
2027 {
2028  for (int i = 0; i < mat.size(); i++)
2029  {
2030  if (!vvIsRowZero(mat, i))
2031  return FALSE;
2032  }
2033  return TRUE;
2034 }
2035 
2036 static std::vector<std::vector<int> > vvMult(const std::vector<std::vector<int> >& a, const std::vector<std::vector<int> >& b)
2037 {
2038  int ra = a.size();
2039  int rb = b.size();
2040  int ca = a.size() > 0 ? a[0].size() : 0;
2041  int cb = b.size() > 0 ? b[0].size() : 0;
2042 
2043  if (ca != rb)
2044  {
2045  WerrorS("matrix dimensions do not match");
2046  return std::vector<std::vector<int> >();
2047  }
2048 
2049  std::vector<std::vector<int> > res(ra, std::vector<int>(cb));
2050  for (int i = 0; i < ra; i++)
2051  {
2052  for (int j = 0; j < cb; j++)
2053  {
2054  int sum = 0;
2055  for (int k = 0; k < ca; k++)
2056  sum += a[i][k] * b[k][j];
2057  res[i][j] = sum;
2058  }
2059  }
2060  return res;
2061 }
2062 
2063 static BOOLEAN isAcyclic(const intvec* G)
2064 {
2065  // init
2066  int n = G->cols();
2067  std::vector<int> path;
2068  std::vector<BOOLEAN> visited;
2069  std::vector<BOOLEAN> cyclic;
2070  std::vector<int> cache;
2071  visited.resize(n, FALSE);
2072  cyclic.resize(n, FALSE);
2073  cache.resize(n, -2);
2074 
2075  for (int v = 0; v < n; v++)
2076  {
2077  cache = countCycles(G, v, path, visited, cyclic, cache);
2078  // check that there are 0 cycles from v
2079  if (cache[v] != 0)
2080  return FALSE;
2081  }
2082  return TRUE;
2083 }
2084 
2085 /*
2086  * Computation of the K-Dimension
2087  */
2088 
2089 // -1 is infinity, -2 is error
2090 int lp_kDim(const ideal _G)
2091 {
2092  if (rField_is_Ring(currRing)) {
2093  WerrorS("K-Dim not implemented for rings");
2094  return -2;
2095  }
2096 
2097  for (int i=IDELEMS(_G)-1;i>=0; i--)
2098  {
2099  if (_G->m[i] != NULL)
2100  {
2101  if (pGetComp(_G->m[i]) != 0)
2102  {
2103  WerrorS("K-Dim not implemented for modules");
2104  return -2;
2105  }
2106  if (pGetNCGen(_G->m[i]) != 0)
2107  {
2108  WerrorS("K-Dim not implemented for bi-modules");
2109  return -2;
2110  }
2111  }
2112  }
2113 
2114  ideal G = id_Head(_G, currRing); // G = LM(G) (and copy)
2115  if (TEST_OPT_PROT)
2116  Print("%d original generators\n", IDELEMS(G));
2117  idSkipZeroes(G); // remove zeros
2118  id_DelLmEquals(G, currRing); // remove duplicates
2119  if (TEST_OPT_PROT)
2120  Print("%d non-zero unique generators\n", IDELEMS(G));
2121 
2122  // check if G is the zero ideal
2123  if (IDELEMS(G) == 1 && G->m[0] == NULL)
2124  {
2125  // NOTE: this is needed because if the ideal is <0>, then idSkipZeroes keeps this element, and IDELEMS is still 1!
2126  int lV = currRing->isLPring;
2127  int ncGenCount = currRing->LPncGenCount;
2128  if (lV - ncGenCount == 0)
2129  {
2130  idDelete(&G);
2131  return 1;
2132  }
2133  if (lV - ncGenCount == 1)
2134  {
2135  idDelete(&G);
2136  return -1;
2137  }
2138  if (lV - ncGenCount >= 2)
2139  {
2140  idDelete(&G);
2141  return -1;
2142  }
2143  }
2144 
2145  // get the max deg
2146  long maxDeg = 0;
2147  for (int i = 0; i < IDELEMS(G); i++)
2148  {
2149  maxDeg = si_max(maxDeg, pTotaldegree(G->m[i]));
2150 
2151  // also check whether G = <1>
2152  if (pIsConstantComp(G->m[i]))
2153  {
2154  WerrorS("K-Dim not defined for 0-ring"); // TODO is it minus infinity ?
2155  idDelete(&G);
2156  return -2;
2157  }
2158  }
2159  if (TEST_OPT_PROT)
2160  Print("max deg: %ld\n", maxDeg);
2161 
2162 
2163  // for normal words of length minDeg ... maxDeg-1
2164  // brute-force the normal words
2165  if (TEST_OPT_PROT)
2166  PrintS("Computing normal words normally...\n");
2167  long numberOfNormalWords = lp_countNormalWords(maxDeg - 1, G);
2168 
2169  if (TEST_OPT_PROT)
2170  Print("%ld normal words up to length %ld\n", numberOfNormalWords, maxDeg - 1);
2171 
2172  // early termination if G \subset X
2173  if (maxDeg <= 1)
2174  {
2175  int lV = currRing->isLPring;
2176  int ncGenCount = currRing->LPncGenCount;
2177  if (IDELEMS(G) == lV - ncGenCount) // V = {1} no edges
2178  {
2179  idDelete(&G);
2180  return numberOfNormalWords;
2181  }
2182  if (IDELEMS(G) == lV - ncGenCount - 1) // V = {1} with loop
2183  {
2184  idDelete(&G);
2185  return -1;
2186  }
2187  if (IDELEMS(G) <= lV - ncGenCount - 2) // V = {1} with more than one loop
2188  {
2189  idDelete(&G);
2190  return -1;
2191  }
2192  }
2193 
2194  if (TEST_OPT_PROT)
2195  PrintS("Computing Ufnarovski graph...\n");
2196 
2197  ideal standardWords;
2198  intvec* UG = lp_ufnarovskiGraph(G, standardWords);
2199  if (UG == NULL)
2200  {
2201  idDelete(&G);
2202  return -2;
2203  }
2204  if (errorreported)
2205  {
2206  delete UG;
2207  idDelete(&G);
2208  return -2;
2209  }
2210 
2211  if (TEST_OPT_PROT)
2212  Print("Ufnarovski graph is %dx%d.\n", UG->rows(), UG->cols());
2213 
2214  if (TEST_OPT_PROT)
2215  PrintS("Checking whether Ufnarovski graph is acyclic...\n");
2216 
2217  if (!isAcyclic(UG))
2218  {
2219  // in this case we have infinitely many normal words
2220  return -1;
2221  }
2222 
2223  std::vector<std::vector<int> > vvUG = iv2vv(UG);
2224  for (int i = 0; i < vvUG.size(); i++)
2225  {
2226  if (vvIsRowZero(vvUG, i) && vvIsColumnZero(vvUG, i)) // i is isolated vertex
2227  {
2228  vvDeleteRow(vvUG, i);
2229  vvDeleteColumn(vvUG, i);
2230  i--;
2231  }
2232  }
2233  if (TEST_OPT_PROT)
2234  Print("Simplified Ufnarovski graph to %dx%d.\n", (int)vvUG.size(), (int)vvUG.size());
2235 
2236  // for normal words of length >= maxDeg
2237  // use Ufnarovski graph
2238  if (TEST_OPT_PROT)
2239  PrintS("Computing normal words via Ufnarovski graph...\n");
2240  std::vector<std::vector<int> > UGpower = vvUG;
2241  long nUGpower = 1;
2242  while (!vvIsZero(UGpower))
2243  {
2244  if (TEST_OPT_PROT)
2245  PrintS("Start count graph entries.\n");
2246  for (int i = 0; i < UGpower.size(); i++)
2247  {
2248  for (int j = 0; j < UGpower[i].size(); j++)
2249  {
2250  numberOfNormalWords += UGpower[i][j];
2251  }
2252  }
2253 
2254  if (TEST_OPT_PROT)
2255  {
2256  PrintS("Done count graph entries.\n");
2257  Print("%ld normal words up to length %ld\n", numberOfNormalWords, maxDeg - 1 + nUGpower);
2258  }
2259 
2260  if (TEST_OPT_PROT)
2261  PrintS("Start mat mult.\n");
2262  UGpower = vvMult(UGpower, vvUG); // TODO: avoid creation of new intvec
2263  if (TEST_OPT_PROT)
2264  PrintS("Done mat mult.\n");
2265  nUGpower++;
2266  }
2267 
2268  delete UG;
2269  idDelete(&G);
2270  return numberOfNormalWords;
2271 }
2272 #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:4082
int p
Definition: cfModGcd.cc:4078
CanonicalForm b
Definition: cfModGcd.cc:4103
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:515
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:753
#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:993
static ideal lp_computeNormalWords(int length, ideal M)
Definition: hdegree.cc:1738
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:1019
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:1950
ideal scKBase(int deg, ideal s, ideal Q, intvec *mv)
Definition: hdegree.cc:1427
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:1779
intvec * scIndIntvec(ideal S, ideal Q)
Definition: hdegree.cc:285
static int scMin(int i, scfmon stc, int Nvar)
Definition: hdegree.cc:1175
static void vvDeleteRow(std::vector< std::vector< int > > &mat, int row)
Definition: hdegree.cc:1993
static indset hCheck2(indset sm, scmon pure)
Definition: hdegree.cc:493
STATIC_VAR poly last
Definition: hdegree.cc:1151
void scComputeHC(ideal S, ideal Q, int ak, poly &hEdge, ring tailRing)
Definition: hdegree.cc:1079
static BOOLEAN hCheck1(indset sm, scmon pure)
Definition: hdegree.cc:467
static int graphGrowth(const intvec *G)
Definition: hdegree.cc:1652
VAR int hMu
Definition: hdegree.cc:27
static BOOLEAN vvIsColumnZero(const std::vector< std::vector< int > > &mat, int col)
Definition: hdegree.cc:2016
VAR omBin indlist_bin
Definition: hdegree.cc:28
STATIC_VAR poly pWork
Definition: hdegree.cc:1005
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:1998
static BOOLEAN hNotZero(scfmon rad, int Nrad, varset var, int Nvar)
Definition: hdegree.cc:354
int lp_kDim(const ideal _G)
Definition: hdegree.cc:2090
static void hDegree0(ideal S, ideal Q, const ring tailRing)
Definition: hdegree.cc:919
static void scElKbase()
Definition: hdegree.cc:1154
static void hHedge(poly hEdge)
Definition: hdegree.cc:1007
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:1187
int lp_gkDim(const ideal _G)
Definition: hdegree.cc:1840
VAR indset ISet
Definition: hdegree.cc:352
static void vvPrint(const std::vector< std::vector< int > > &mat)
Definition: hdegree.cc:1968
static void vvTest(const std::vector< std::vector< int > > &mat)
Definition: hdegree.cc:1980
static void scAllKbase(int Nvar, int ideg, int deg)
Definition: hdegree.cc:1262
static void scAll(int Nvar, int deg)
Definition: hdegree.cc:1238
int scMultInt(ideal S, ideal Q)
Definition: hdegree.cc:872
static void scDegKbase(scfmon stc, int Nstc, int Nvar, int deg)
Definition: hdegree.cc:1272
STATIC_VAR scmon act
Definition: hdegree.cc:1152
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:1758
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:2063
static int scMax(int i, scfmon stc, int Nvar)
Definition: hdegree.cc:1163
static ideal scIdKbase(poly q, const int rank)
Definition: hdegree.cc:1409
static void hIndep(scmon pure)
Definition: hdegree.cc:369
static void scInKbase(scfmon stc, int Nstc, int Nvar)
Definition: hdegree.cc:1353
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:2026
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:2006
static std::vector< std::vector< int > > vvMult(const std::vector< std::vector< int > > &a, const std::vector< std::vector< int > > &b)
Definition: hdegree.cc:2036
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:1588
static void _lp_computeNormalWords(ideal words, int &numberOfNormalWords, int length, ideal M, int minDeg, int &last)
Definition: hdegree.cc:1679
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:1418
intvec * hSecondSeries(intvec *hseries1)
Definition: hilb.cc:1383
intvec * hFirstSeries(ideal S, intvec *modulweight, ideal Q, intvec *wdegree, ring tailRing)
Definition: hilb.cc:1373
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:103
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:1222
static void p_Delete(poly *p, const ring r)
Definition: p_polys.h:873
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_Z(const ring r)
Definition: ring.h:510
#define rField_is_Ring(R)
Definition: ring.h:486
BOOLEAN p_LPDivisibleBy(poly a, poly b, const ring r)
Definition: shiftop.cc:774
poly p_LPVarAt(poly p, int pos, const ring r)
Definition: shiftop.cc:843
#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:75