My Project  debian-1:4.1.2-p1+ds-2
hilb.cc
Go to the documentation of this file.
1 /****************************************
2 * Computer Algebra System SINGULAR *
3 ****************************************/
4 /*
5 * ABSTRACT - Hilbert series
6 */
7 
8 #include "kernel/mod2.h"
9 
10 #include "misc/mylimits.h"
11 #include "misc/intvec.h"
12 
16 
17 #include "polys/monomials/ring.h"
19 #include "polys/simpleideals.h"
20 
21 
22 // #include "kernel/structs.h"
23 // #include "kernel/polys.h"
24 //ADICHANGES:
25 #include "kernel/ideals.h"
26 // #include "kernel/GBEngine/kstd1.h"
27 
28 # define omsai 1
29 #if omsai
30 
32 #include "coeffs/coeffs.h"
34 #include "coeffs/numbers.h"
35 #include <vector>
36 #include "Singular/ipshell.h"
37 
38 
39 #include <Singular/ipshell.h>
40 #include <ctime>
41 #include <iostream>
42 #endif
43 
45 STATIC_VAR int *Q0, *Ql;
47 
48 
49 static int hMinModulweight(intvec *modulweight)
50 {
51  int i,j,k;
52 
53  if(modulweight==NULL) return 0;
54  j=(*modulweight)[0];
55  for(i=modulweight->rows()-1;i!=0;i--)
56  {
57  k=(*modulweight)[i];
58  if(k<j) j=k;
59  }
60  return j;
61 }
62 
63 static void hHilbEst(scfmon stc, int Nstc, varset var, int Nvar)
64 {
65  int i, j;
66  int x, y, z = 1;
67  int *p;
68  for (i = Nvar; i>0; i--)
69  {
70  x = 0;
71  for (j = 0; j < Nstc; j++)
72  {
73  y = stc[j][var[i]];
74  if (y > x)
75  x = y;
76  }
77  z += x;
78  j = i - 1;
79  if (z > Ql[j])
80  {
81  if (z>(MAX_INT_VAL)/2)
82  {
83  WerrorS("internal arrays too big");
84  return;
85  }
86  p = (int *)omAlloc((unsigned long)z * sizeof(int));
87  if (Ql[j]!=0)
88  {
89  if (j==0)
90  memcpy(p, Qpol[j], Ql[j] * sizeof(int));
91  omFreeSize((ADDRESS)Qpol[j], Ql[j] * sizeof(int));
92  }
93  if (j==0)
94  {
95  for (x = Ql[j]; x < z; x++)
96  p[x] = 0;
97  }
98  Ql[j] = z;
99  Qpol[j] = p;
100  }
101  }
102 }
103 
104 static int *hAddHilb(int Nv, int x, int *pol, int *lp)
105 {
106  int l = *lp, ln, i;
107  int *pon;
108  *lp = ln = l + x;
109  pon = Qpol[Nv];
110  memcpy(pon, pol, l * sizeof(int));
111  if (l > x)
112  {/*pon[i] -= pol[i - x];*/
113  for (i = x; i < l; i++)
114  { int64 t=pon[i];
115  int64 t2=pol[i - x];
116  t-=t2;
117  if ((t>=INT_MIN)&&(t<=INT_MAX)) pon[i]=t;
118  else if (!errorreported) WerrorS("int overflow in hilb 1");
119  }
120  for (i = l; i < ln; i++)
121  { /*pon[i] = -pol[i - x];*/
122  int64 t= -pol[i - x];
123  if ((t>=INT_MIN)&&(t<=INT_MAX)) pon[i]=t;
124  else if (!errorreported) WerrorS("int overflow in hilb 2");
125  }
126  }
127  else
128  {
129  for (i = l; i < x; i++)
130  pon[i] = 0;
131  for (i = x; i < ln; i++)
132  pon[i] = -pol[i - x];
133  }
134  return pon;
135 }
136 
137 static void hLastHilb(scmon pure, int Nv, varset var, int *pol, int lp)
138 {
139  int l = lp, x, i, j;
140  int *pl;
141  int *p;
142  p = pol;
143  for (i = Nv; i>0; i--)
144  {
145  x = pure[var[i + 1]];
146  if (x!=0)
147  p = hAddHilb(i, x, p, &l);
148  }
149  pl = *Qpol;
150  j = Q0[Nv + 1];
151  for (i = 0; i < l; i++)
152  { /* pl[i + j] += p[i];*/
153  int64 t=pl[i+j];
154  int64 t2=p[i];
155  t+=t2;
156  if ((t>=INT_MIN)&&(t<=INT_MAX)) pl[i+j]=t;
157  else if (!errorreported) WerrorS("int overflow in hilb 3");
158  }
159  x = pure[var[1]];
160  if (x!=0)
161  {
162  j += x;
163  for (i = 0; i < l; i++)
164  { /* pl[i + j] -= p[i];*/
165  int64 t=pl[i+j];
166  int64 t2=p[i];
167  t-=t2;
168  if ((t>=INT_MIN)&&(t<=INT_MAX)) pl[i+j]=t;
169  else if (!errorreported) WerrorS("int overflow in hilb 4");
170  }
171  }
172  j += l;
173  if (j > hLength)
174  hLength = j;
175 }
176 
177 static void hHilbStep(scmon pure, scfmon stc, int Nstc, varset var,
178  int Nvar, int *pol, int Lpol)
179 {
180  int iv = Nvar -1, ln, a, a0, a1, b, i;
181  int x, x0;
182  scmon pn;
183  scfmon sn;
184  int *pon;
185  if (Nstc==0)
186  {
187  hLastHilb(pure, iv, var, pol, Lpol);
188  return;
189  }
190  x = a = 0;
191  pn = hGetpure(pure);
192  sn = hGetmem(Nstc, stc, stcmem[iv]);
193  hStepS(sn, Nstc, var, Nvar, &a, &x);
194  Q0[iv] = Q0[Nvar];
195  ln = Lpol;
196  pon = pol;
197  if (a == Nstc)
198  {
199  x = pure[var[Nvar]];
200  if (x!=0)
201  pon = hAddHilb(iv, x, pon, &ln);
202  hHilbStep(pn, sn, a, var, iv, pon, ln);
203  return;
204  }
205  else
206  {
207  pon = hAddHilb(iv, x, pon, &ln);
208  hHilbStep(pn, sn, a, var, iv, pon, ln);
209  }
210  b = a;
211  x0 = 0;
212  loop
213  {
214  Q0[iv] += (x - x0);
215  a0 = a;
216  x0 = x;
217  hStepS(sn, Nstc, var, Nvar, &a, &x);
218  hElimS(sn, &b, a0, a, var, iv);
219  a1 = a;
220  hPure(sn, a0, &a1, var, iv, pn, &i);
221  hLex2S(sn, b, a0, a1, var, iv, hwork);
222  b += (a1 - a0);
223  ln = Lpol;
224  if (a < Nstc)
225  {
226  pon = hAddHilb(iv, x - x0, pol, &ln);
227  hHilbStep(pn, sn, b, var, iv, pon, ln);
228  }
229  else
230  {
231  x = pure[var[Nvar]];
232  if (x!=0)
233  pon = hAddHilb(iv, x - x0, pol, &ln);
234  else
235  pon = pol;
236  hHilbStep(pn, sn, b, var, iv, pon, ln);
237  return;
238  }
239  }
240 }
241 
242 /*
243 *basic routines
244 */
245 static void hWDegree(intvec *wdegree)
246 {
247  int i, k;
248  int x;
249 
250  for (i=(currRing->N); i; i--)
251  {
252  x = (*wdegree)[i-1];
253  if (x != 1)
254  {
255  for (k=hNexist-1; k>=0; k--)
256  {
257  hexist[k][i] *= x;
258  }
259  }
260  }
261 }
262 // ---------------------------------- ADICHANGES ---------------------------------------------
263 //!!!!!!!!!!!!!!!!!!!!! Just for Monomial Ideals !!!!!!!!!!!!!!!!!!!!!!!!!!!!
264 
265 #if 0 // unused
266 //Tests if the ideal is sorted by degree
267 static bool idDegSortTest(ideal I)
268 {
269  if((I == NULL)||(idIs0(I)))
270  {
271  return(TRUE);
272  }
273  for(int i = 0; i<IDELEMS(I)-1; i++)
274  {
275  if(p_Totaldegree(I->m[i],currRing)>p_Totaldegree(I->m[i+1],currRing))
276  {
277  idPrint(I);
278  WerrorS("Ideal is not deg sorted!!");
279  return(FALSE);
280  }
281  }
282  return(TRUE);
283 }
284 #endif
285 
286 //adds the new polynomial at the coresponding position
287 //and simplifies the ideal, destroys p
288 static void SortByDeg_p(ideal I, poly p)
289 {
290  int i,j;
291  if(idIs0(I))
292  {
293  I->m[0] = p;
294  return;
295  }
296  idSkipZeroes(I);
297  #if 1
298  for(i = 0; (i<IDELEMS(I)) && (p_Totaldegree(I->m[i],currRing)<=p_Totaldegree(p,currRing)); i++)
299  {
300  if(p_DivisibleBy( I->m[i],p, currRing))
301  {
302  p_Delete(&p,currRing);
303  return;
304  }
305  }
306  for(i = IDELEMS(I)-1; (i>=0) && (p_Totaldegree(I->m[i],currRing)>=p_Totaldegree(p,currRing)); i--)
307  {
308  if(p_DivisibleBy(p,I->m[i], currRing))
309  {
310  p_Delete(&I->m[i],currRing);
311  }
312  }
313  if(idIs0(I))
314  {
315  idSkipZeroes(I);
316  I->m[0] = p;
317  return;
318  }
319  #endif
320  idSkipZeroes(I);
321  //First I take the case when all generators have the same degree
322  if(p_Totaldegree(I->m[0],currRing) == p_Totaldegree(I->m[IDELEMS(I)-1],currRing))
323  {
325  {
326  idInsertPoly(I,p);
327  idSkipZeroes(I);
328  for(i=IDELEMS(I)-1;i>=1; i--)
329  {
330  I->m[i] = I->m[i-1];
331  }
332  I->m[0] = p;
333  return;
334  }
336  {
337  idInsertPoly(I,p);
338  idSkipZeroes(I);
339  return;
340  }
341  }
343  {
344  idInsertPoly(I,p);
345  idSkipZeroes(I);
346  for(i=IDELEMS(I)-1;i>=1; i--)
347  {
348  I->m[i] = I->m[i-1];
349  }
350  I->m[0] = p;
351  return;
352  }
354  {
355  idInsertPoly(I,p);
356  idSkipZeroes(I);
357  return;
358  }
359  for(i = IDELEMS(I)-2; ;)
360  {
362  {
363  idInsertPoly(I,p);
364  idSkipZeroes(I);
365  for(j = IDELEMS(I)-1; j>=i+1;j--)
366  {
367  I->m[j] = I->m[j-1];
368  }
369  I->m[i] = p;
370  return;
371  }
373  {
374  idInsertPoly(I,p);
375  idSkipZeroes(I);
376  for(j = IDELEMS(I)-1; j>=i+2;j--)
377  {
378  I->m[j] = I->m[j-1];
379  }
380  I->m[i+1] = p;
381  return;
382  }
383  i--;
384  }
385 }
386 
387 //it sorts the ideal by the degrees
388 static ideal SortByDeg(ideal I)
389 {
390  if(idIs0(I))
391  {
392  return id_Copy(I,currRing);
393  }
394  int i;
395  ideal res;
396  idSkipZeroes(I);
397  res = idInit(1,1);
398  for(i = 0; i<=IDELEMS(I)-1;i++)
399  {
400  SortByDeg_p(res, I->m[i]);
401  I->m[i]=NULL; // I->m[i] is now in res
402  }
403  idSkipZeroes(res);
404  //idDegSortTest(res);
405  return(res);
406 }
407 
408 //idQuot(I,p) for I monomial ideal, p a ideal with a single monomial.
409 ideal idQuotMon(ideal Iorig, ideal p)
410 {
411  if(idIs0(Iorig))
412  {
413  ideal res = idInit(1,1);
414  res->m[0] = poly(0);
415  return(res);
416  }
417  if(idIs0(p))
418  {
419  ideal res = idInit(1,1);
420  res->m[0] = pOne();
421  return(res);
422  }
423  ideal I = id_Head(Iorig,currRing);
424  ideal res = idInit(IDELEMS(I),1);
425  int i,j;
426  int dummy;
427  for(i = 0; i<IDELEMS(I); i++)
428  {
429  res->m[i] = p_Head(I->m[i], currRing);
430  for(j = 1; (j<=currRing->N) ; j++)
431  {
432  dummy = p_GetExp(p->m[0], j, currRing);
433  if(dummy > 0)
434  {
435  if(p_GetExp(I->m[i], j, currRing) < dummy)
436  {
437  p_SetExp(res->m[i], j, 0, currRing);
438  }
439  else
440  {
441  p_SetExp(res->m[i], j, p_GetExp(I->m[i], j, currRing) - dummy, currRing);
442  }
443  }
444  }
445  p_Setm(res->m[i], currRing);
446  if(p_Totaldegree(res->m[i],currRing) == p_Totaldegree(I->m[i],currRing))
447  {
448  p_Delete(&res->m[i],currRing);
449  }
450  else
451  {
452  p_Delete(&I->m[i],currRing);
453  }
454  }
455  idSkipZeroes(res);
456  idSkipZeroes(I);
457  if(!idIs0(res))
458  {
459  for(i = 0; i<=IDELEMS(res)-1; i++)
460  {
461  SortByDeg_p(I,res->m[i]);
462  res->m[i]=NULL; // is now in I
463  }
464  }
466  //idDegSortTest(I);
467  return(I);
468 }
469 
470 //id_Add for monomials
471 static void idAddMon(ideal I, ideal p)
472 {
473  SortByDeg_p(I,p->m[0]);
474  p->m[0]=NULL; // is now in I
475  //idSkipZeroes(I);
476 }
477 
478 //searches for a variable that is not yet used (assumes that the ideal is sqrfree)
479 static poly ChoosePVar (ideal I)
480 {
481  bool flag=TRUE;
482  int i,j;
483  poly res;
484  for(i=1;i<=currRing->N;i++)
485  {
486  flag=TRUE;
487  for(j=IDELEMS(I)-1;(j>=0)&&(flag);j--)
488  {
489  if(p_GetExp(I->m[j], i, currRing)>0)
490  {
491  flag=FALSE;
492  }
493  }
494 
495  if(flag == TRUE)
496  {
497  res = p_ISet(1, currRing);
498  p_SetExp(res, i, 1, currRing);
500  return(res);
501  }
502  else
503  {
504  p_Delete(&res, currRing);
505  }
506  }
507  return(NULL); //i.e. it is the maximal ideal
508 }
509 
510 #if 0 // unused
511 //choice XL: last entry divided by x (xy10z15 -> y9z14)
512 static poly ChoosePXL(ideal I)
513 {
514  int i,j,dummy=0;
515  poly m;
516  for(i = IDELEMS(I)-1; (i>=0) && (dummy == 0); i--)
517  {
518  for(j = 1; (j<=currRing->N) && (dummy == 0); j++)
519  {
520  if(p_GetExp(I->m[i],j, currRing)>1)
521  {
522  dummy = 1;
523  }
524  }
525  }
526  m = p_Copy(I->m[i+1],currRing);
527  for(j = 1; j<=currRing->N; j++)
528  {
529  dummy = p_GetExp(m,j,currRing);
530  if(dummy >= 1)
531  {
532  p_SetExp(m, j, dummy-1, currRing);
533  }
534  }
535  if(!p_IsOne(m, currRing))
536  {
537  p_Setm(m, currRing);
538  return(m);
539  }
540  m = ChoosePVar(I);
541  return(m);
542 }
543 #endif
544 
545 #if 0 // unused
546 //choice XF: first entry divided by x (xy10z15 -> y9z14)
547 static poly ChoosePXF(ideal I)
548 {
549  int i,j,dummy=0;
550  poly m;
551  for(i =0 ; (i<=IDELEMS(I)-1) && (dummy == 0); i++)
552  {
553  for(j = 1; (j<=currRing->N) && (dummy == 0); j++)
554  {
555  if(p_GetExp(I->m[i],j, currRing)>1)
556  {
557  dummy = 1;
558  }
559  }
560  }
561  m = p_Copy(I->m[i-1],currRing);
562  for(j = 1; j<=currRing->N; j++)
563  {
564  dummy = p_GetExp(m,j,currRing);
565  if(dummy >= 1)
566  {
567  p_SetExp(m, j, dummy-1, currRing);
568  }
569  }
570  if(!p_IsOne(m, currRing))
571  {
572  p_Setm(m, currRing);
573  return(m);
574  }
575  m = ChoosePVar(I);
576  return(m);
577 }
578 #endif
579 
580 #if 0 // unused
581 //choice OL: last entry the first power (xy10z15 -> xy9z15)
582 static poly ChoosePOL(ideal I)
583 {
584  int i,j,dummy;
585  poly m;
586  for(i = IDELEMS(I)-1;i>=0;i--)
587  {
588  m = p_Copy(I->m[i],currRing);
589  for(j=1;j<=currRing->N;j++)
590  {
591  dummy = p_GetExp(m,j,currRing);
592  if(dummy > 0)
593  {
594  p_SetExp(m,j,dummy-1,currRing);
595  p_Setm(m,currRing);
596  }
597  }
598  if(!p_IsOne(m, currRing))
599  {
600  return(m);
601  }
602  else
603  {
604  p_Delete(&m,currRing);
605  }
606  }
607  m = ChoosePVar(I);
608  return(m);
609 }
610 #endif
611 
612 #if 0 // unused
613 //choice OF: first entry the first power (xy10z15 -> xy9z15)
614 static poly ChoosePOF(ideal I)
615 {
616  int i,j,dummy;
617  poly m;
618  for(i = 0 ;i<=IDELEMS(I)-1;i++)
619  {
620  m = p_Copy(I->m[i],currRing);
621  for(j=1;j<=currRing->N;j++)
622  {
623  dummy = p_GetExp(m,j,currRing);
624  if(dummy > 0)
625  {
626  p_SetExp(m,j,dummy-1,currRing);
627  p_Setm(m,currRing);
628  }
629  }
630  if(!p_IsOne(m, currRing))
631  {
632  return(m);
633  }
634  else
635  {
636  p_Delete(&m,currRing);
637  }
638  }
639  m = ChoosePVar(I);
640  return(m);
641 }
642 #endif
643 
644 #if 0 // unused
645 //choice VL: last entry the first variable with power (xy10z15 -> y)
646 static poly ChoosePVL(ideal I)
647 {
648  int i,j,dummy;
649  bool flag = TRUE;
650  poly m = p_ISet(1,currRing);
651  for(i = IDELEMS(I)-1;(i>=0) && (flag);i--)
652  {
653  flag = TRUE;
654  for(j=1;(j<=currRing->N) && (flag);j++)
655  {
656  dummy = p_GetExp(I->m[i],j,currRing);
657  if(dummy >= 2)
658  {
659  p_SetExp(m,j,1,currRing);
660  p_Setm(m,currRing);
661  flag = FALSE;
662  }
663  }
664  if(!p_IsOne(m, currRing))
665  {
666  return(m);
667  }
668  }
669  m = ChoosePVar(I);
670  return(m);
671 }
672 #endif
673 
674 #if 0 // unused
675 //choice VF: first entry the first variable with power (xy10z15 -> y)
676 static poly ChoosePVF(ideal I)
677 {
678  int i,j,dummy;
679  bool flag = TRUE;
680  poly m = p_ISet(1,currRing);
681  for(i = 0;(i<=IDELEMS(I)-1) && (flag);i++)
682  {
683  flag = TRUE;
684  for(j=1;(j<=currRing->N) && (flag);j++)
685  {
686  dummy = p_GetExp(I->m[i],j,currRing);
687  if(dummy >= 2)
688  {
689  p_SetExp(m,j,1,currRing);
690  p_Setm(m,currRing);
691  flag = FALSE;
692  }
693  }
694  if(!p_IsOne(m, currRing))
695  {
696  return(m);
697  }
698  }
699  m = ChoosePVar(I);
700  return(m);
701 }
702 #endif
703 
704 //choice JL: last entry just variable with power (xy10z15 -> y10)
705 static poly ChoosePJL(ideal I)
706 {
707  int i,j,dummy;
708  bool flag = TRUE;
709  poly m = p_ISet(1,currRing);
710  for(i = IDELEMS(I)-1;(i>=0) && (flag);i--)
711  {
712  flag = TRUE;
713  for(j=1;(j<=currRing->N) && (flag);j++)
714  {
715  dummy = p_GetExp(I->m[i],j,currRing);
716  if(dummy >= 2)
717  {
718  p_SetExp(m,j,dummy-1,currRing);
719  p_Setm(m,currRing);
720  flag = FALSE;
721  }
722  }
723  if(!p_IsOne(m, currRing))
724  {
725  return(m);
726  }
727  }
728  p_Delete(&m,currRing);
729  m = ChoosePVar(I);
730  return(m);
731 }
732 
733 #if 0
734 //choice JF: last entry just variable with power -1 (xy10z15 -> y9)
735 static poly ChoosePJF(ideal I)
736 {
737  int i,j,dummy;
738  bool flag = TRUE;
739  poly m = p_ISet(1,currRing);
740  for(i = 0;(i<=IDELEMS(I)-1) && (flag);i++)
741  {
742  flag = TRUE;
743  for(j=1;(j<=currRing->N) && (flag);j++)
744  {
745  dummy = p_GetExp(I->m[i],j,currRing);
746  if(dummy >= 2)
747  {
748  p_SetExp(m,j,dummy-1,currRing);
749  p_Setm(m,currRing);
750  flag = FALSE;
751  }
752  }
753  if(!p_IsOne(m, currRing))
754  {
755  return(m);
756  }
757  }
758  m = ChoosePVar(I);
759  return(m);
760 }
761 #endif
762 
763 //chooses 1 \neq p \not\in S. This choice should be made optimal
764 static poly ChooseP(ideal I)
765 {
766  poly m;
767  // TEST TO SEE WHICH ONE IS BETTER
768  //m = ChoosePXL(I);
769  //m = ChoosePXF(I);
770  //m = ChoosePOL(I);
771  //m = ChoosePOF(I);
772  //m = ChoosePVL(I);
773  //m = ChoosePVF(I);
774  m = ChoosePJL(I);
775  //m = ChoosePJF(I);
776  return(m);
777 }
778 
779 ///searches for a monomial of degree d>=2 and divides it by a variable (result monomial of deg d-1)
780 static poly SearchP(ideal I)
781 {
782  int i,j,exp;
783  poly res;
784  if(p_Totaldegree(I->m[IDELEMS(I)-1],currRing)<=1)
785  {
786  res = ChoosePVar(I);
787  return(res);
788  }
789  i = IDELEMS(I)-1;
790  res = p_Copy(I->m[i], currRing);
791  for(j=1;j<=currRing->N;j++)
792  {
793  exp = p_GetExp(I->m[i], j, currRing);
794  if(exp > 0)
795  {
796  p_SetExp(res, j, exp - 1, currRing);
798  break;
799  }
800  }
801  assume( j <= currRing->N );
802  return(res);
803 }
804 
805 //test if the ideal is of the form (x1, ..., xr)
806 static bool JustVar(ideal I)
807 {
808  #if 0
809  int i,j;
810  bool foundone;
811  for(i=0;i<=IDELEMS(I)-1;i++)
812  {
813  foundone = FALSE;
814  for(j = 1;j<=currRing->N;j++)
815  {
816  if(p_GetExp(I->m[i], j, currRing)>0)
817  {
818  if(foundone == TRUE)
819  {
820  return(FALSE);
821  }
822  foundone = TRUE;
823  }
824  }
825  }
826  return(TRUE);
827  #else
828  if(p_Totaldegree(I->m[IDELEMS(I)-1],currRing)>1)
829  {
830  return(FALSE);
831  }
832  return(TRUE);
833  #endif
834 }
835 
836 //computes the Euler Characteristic of the ideal
837 static void eulerchar (ideal I, int variables, mpz_ptr ec)
838 {
839  loop
840  {
841  mpz_t dummy;
842  if(JustVar(I) == TRUE)
843  {
844  if(IDELEMS(I) == variables)
845  {
846  mpz_init(dummy);
847  if((variables % 2) == 0)
848  mpz_set_ui(dummy, 1);
849  else
850  mpz_set_si(dummy, -1);
851  mpz_add(ec, ec, dummy);
852  mpz_clear(dummy);
853  }
854  return;
855  }
856  ideal p = idInit(1,1);
857  p->m[0] = SearchP(I);
858  //idPrint(I);
859  //idPrint(p);
860  //printf("\nNow get in idQuotMon\n");
861  ideal Ip = idQuotMon(I,p);
862  //idPrint(Ip);
863  //Ip = SortByDeg(Ip);
864  int i,howmanyvarinp = 0;
865  for(i = 1;i<=currRing->N;i++)
866  {
867  if(p_GetExp(p->m[0],i,currRing)>0)
868  {
869  howmanyvarinp++;
870  }
871  }
872  eulerchar(Ip, variables-howmanyvarinp, ec);
873  id_Delete(&Ip, currRing);
874  idAddMon(I,p);
875  id_Delete(&p, currRing);
876  }
877 }
878 
879 //tests if an ideal is Square Free, if no, returns the variable which appears at powers >1
880 static poly SqFree (ideal I)
881 {
882  int i,j;
883  bool flag=TRUE;
884  poly notsqrfree = NULL;
885  if(p_Totaldegree(I->m[IDELEMS(I)-1],currRing)<=1)
886  {
887  return(notsqrfree);
888  }
889  for(i=IDELEMS(I)-1;(i>=0)&&(flag);i--)
890  {
891  for(j=1;(j<=currRing->N)&&(flag);j++)
892  {
893  if(p_GetExp(I->m[i],j,currRing)>1)
894  {
895  flag=FALSE;
896  notsqrfree = p_ISet(1,currRing);
897  p_SetExp(notsqrfree,j,1,currRing);
898  }
899  }
900  }
901  if(notsqrfree != NULL)
902  {
903  p_Setm(notsqrfree,currRing);
904  }
905  return(notsqrfree);
906 }
907 
908 //checks if a polynomial is in I
909 static bool IsIn(poly p, ideal I)
910 {
911  //assumes that I is ordered by degree
912  if(idIs0(I))
913  {
914  if(p==poly(0))
915  {
916  return(TRUE);
917  }
918  else
919  {
920  return(FALSE);
921  }
922  }
923  if(p==poly(0))
924  {
925  return(FALSE);
926  }
927  int i,j;
928  bool flag;
929  for(i = 0;i<IDELEMS(I);i++)
930  {
931  flag = TRUE;
932  for(j = 1;(j<=currRing->N) &&(flag);j++)
933  {
934  if(p_GetExp(p, j, currRing)<p_GetExp(I->m[i], j, currRing))
935  {
936  flag = FALSE;
937  }
938  }
939  if(flag)
940  {
941  return(TRUE);
942  }
943  }
944  return(FALSE);
945 }
946 
947 //computes the lcm of min I, I monomial ideal
948 static poly LCMmon(ideal I)
949 {
950  if(idIs0(I))
951  {
952  return(NULL);
953  }
954  poly m;
955  int dummy,i,j;
956  m = p_ISet(1,currRing);
957  for(i=1;i<=currRing->N;i++)
958  {
959  dummy=0;
960  for(j=IDELEMS(I)-1;j>=0;j--)
961  {
962  if(p_GetExp(I->m[j],i,currRing) > dummy)
963  {
964  dummy = p_GetExp(I->m[j],i,currRing);
965  }
966  }
967  p_SetExp(m,i,dummy,currRing);
968  }
969  p_Setm(m,currRing);
970  return(m);
971 }
972 
973 //the Roune Slice Algorithm
974 void rouneslice(ideal I, ideal S, poly q, poly x, int &prune, int &moreprune, int &steps, int &NNN, mpz_ptr &hilbertcoef, int* &hilbpower)
975 {
976  loop
977  {
978  (steps)++;
979  int i,j;
980  int dummy;
981  poly m;
982  ideal p;
983  //----------- PRUNING OF S ---------------
984  //S SHOULD IN THIS POINT BE ORDERED BY DEGREE
985  for(i=IDELEMS(S)-1;i>=0;i--)
986  {
987  if(IsIn(S->m[i],I))
988  {
989  p_Delete(&S->m[i],currRing);
990  prune++;
991  }
992  }
993  idSkipZeroes(S);
994  //----------------------------------------
995  for(i=IDELEMS(I)-1;i>=0;i--)
996  {
997  m = p_Head(I->m[i],currRing);
998  for(j=1;j<=currRing->N;j++)
999  {
1000  dummy = p_GetExp(m,j,currRing);
1001  if(dummy > 0)
1002  {
1003  p_SetExp(m,j,dummy-1,currRing);
1004  }
1005  }
1006  p_Setm(m, currRing);
1007  if(IsIn(m,S))
1008  {
1009  p_Delete(&I->m[i],currRing);
1010  //printf("\n Deleted, since pi(m) is in S\n");pWrite(m);
1011  }
1012  p_Delete(&m,currRing);
1013  }
1014  idSkipZeroes(I);
1015  //----------- MORE PRUNING OF S ------------
1016  m = LCMmon(I);
1017  if(m != NULL)
1018  {
1019  for(i=0;i<IDELEMS(S);i++)
1020  {
1021  if(!(p_DivisibleBy(S->m[i], m, currRing)))
1022  {
1023  S->m[i] = NULL;
1024  j++;
1025  moreprune++;
1026  }
1027  else
1028  {
1029  if(pLmEqual(S->m[i],m))
1030  {
1031  S->m[i] = NULL;
1032  moreprune++;
1033  }
1034  }
1035  }
1036  idSkipZeroes(S);
1037  }
1038  p_Delete(&m,currRing);
1039  /*printf("\n---------------------------\n");
1040  printf("\n I\n");idPrint(I);
1041  printf("\n S\n");idPrint(S);
1042  printf("\n q\n");pWrite(q);
1043  getchar();*/
1044 
1045  if(idIs0(I))
1046  {
1047  id_Delete(&I, currRing);
1048  id_Delete(&S, currRing);
1049  break;
1050  }
1051  m = LCMmon(I);
1052  if(!p_DivisibleBy(x,m, currRing))
1053  {
1054  //printf("\nx does not divide lcm(I)");
1055  //printf("\nEmpty set");pWrite(q);
1056  id_Delete(&I, currRing);
1057  id_Delete(&S, currRing);
1058  p_Delete(&m, currRing);
1059  break;
1060  }
1061  p_Delete(&m, currRing);
1062  m = SqFree(I);
1063  if(m==NULL)
1064  {
1065  //printf("\n Corner: ");
1066  //pWrite(q);
1067  //printf("\n With the facets of the dual simplex:\n");
1068  //idPrint(I);
1069  mpz_t ec;
1070  mpz_init(ec);
1071  mpz_ptr ec_ptr = ec;
1072  eulerchar(I, currRing->N, ec_ptr);
1073  bool flag = FALSE;
1074  if(NNN==0)
1075  {
1076  hilbertcoef = (mpz_ptr)omAlloc((NNN+1)*sizeof(mpz_t));
1077  hilbpower = (int*)omAlloc((NNN+1)*sizeof(int));
1078  mpz_init_set( &hilbertcoef[NNN], ec);
1079  hilbpower[NNN] = p_Totaldegree(q,currRing);
1080  NNN++;
1081  }
1082  else
1083  {
1084  //I look if the power appears already
1085  for(i = 0;(i<NNN)&&(flag == FALSE)&&(p_Totaldegree(q,currRing)>=hilbpower[i]);i++)
1086  {
1087  if((hilbpower[i]) == (p_Totaldegree(q,currRing)))
1088  {
1089  flag = TRUE;
1090  mpz_add(&hilbertcoef[i],&hilbertcoef[i],ec_ptr);
1091  }
1092  }
1093  if(flag == FALSE)
1094  {
1095  hilbertcoef = (mpz_ptr)omRealloc(hilbertcoef, (NNN+1)*sizeof(mpz_t));
1096  hilbpower = (int*)omRealloc(hilbpower, (NNN+1)*sizeof(int));
1097  mpz_init(&hilbertcoef[NNN]);
1098  for(j = NNN; j>i; j--)
1099  {
1100  mpz_set(&hilbertcoef[j],&hilbertcoef[j-1]);
1101  hilbpower[j] = hilbpower[j-1];
1102  }
1103  mpz_set( &hilbertcoef[i], ec);
1104  hilbpower[i] = p_Totaldegree(q,currRing);
1105  NNN++;
1106  }
1107  }
1108  mpz_clear(ec);
1109  id_Delete(&I, currRing);
1110  id_Delete(&S, currRing);
1111  break;
1112  }
1113  else
1114  p_Delete(&m, currRing);
1115  m = ChooseP(I);
1116  p = idInit(1,1);
1117  p->m[0] = m;
1118  ideal Ip = idQuotMon(I,p);
1119  ideal Sp = idQuotMon(S,p);
1120  poly pq = pp_Mult_mm(q,m,currRing);
1121  rouneslice(Ip, Sp, pq, x, prune, moreprune, steps, NNN, hilbertcoef,hilbpower);
1122  idAddMon(S,p);
1123  p->m[0]=NULL;
1124  id_Delete(&p, currRing); // p->m[0] was also in S
1125  p_Delete(&pq,currRing);
1126  }
1127 }
1128 
1129 //it computes the first hilbert series by means of Roune Slice Algorithm
1130 void slicehilb(ideal I)
1131 {
1132  //printf("Adi changes are here: \n");
1133  int i, NNN = 0;
1134  int steps = 0, prune = 0, moreprune = 0;
1135  mpz_ptr hilbertcoef;
1136  int *hilbpower;
1137  ideal S = idInit(1,1);
1138  poly q = p_One(currRing);
1139  ideal X = idInit(1,1);
1140  X->m[0]=p_One(currRing);
1141  for(i=1;i<=currRing->N;i++)
1142  {
1143  p_SetExp(X->m[0],i,1,currRing);
1144  }
1145  p_Setm(X->m[0],currRing);
1146  I = id_Mult(I,X,currRing);
1147  ideal Itmp = SortByDeg(I);
1148  id_Delete(&I,currRing);
1149  I = Itmp;
1150  //printf("\n-------------RouneSlice--------------\n");
1151  rouneslice(I,S,q,X->m[0],prune, moreprune, steps, NNN, hilbertcoef, hilbpower);
1152  id_Delete(&X,currRing);
1153  p_Delete(&q,currRing);
1154  //printf("\nIn total Prune got rid of %i elements\n",prune);
1155  //printf("\nIn total More Prune got rid of %i elements\n",moreprune);
1156  //printf("\nSteps of rouneslice: %i\n\n", steps);
1157  printf("\n// %8d t^0",1);
1158  for(i = 0; i<NNN; i++)
1159  {
1160  if(mpz_sgn(&hilbertcoef[i])!=0)
1161  {
1162  gmp_printf("\n// %8Zd t^%d",&hilbertcoef[i],hilbpower[i]);
1163  }
1164  }
1165  PrintLn();
1166  omFreeSize(hilbertcoef, (NNN)*sizeof(mpz_t));
1167  omFreeSize(hilbpower, (NNN)*sizeof(int));
1168  //printf("\n-------------------------------------\n");
1169 }
1170 
1171 // -------------------------------- END OF CHANGES -------------------------------------------
1172 static intvec * hSeries(ideal S, intvec *modulweight,
1173  int /*notstc*/, intvec *wdegree, ideal Q, ring tailRing)
1174 {
1175 // id_TestTail(S, currRing, tailRing);
1176 
1177  intvec *work, *hseries1=NULL;
1178  int mc;
1179  int p0;
1180  int i, j, k, l, ii, mw;
1181  hexist = hInit(S, Q, &hNexist, tailRing);
1182  if (hNexist==0)
1183  {
1184  hseries1=new intvec(2);
1185  (*hseries1)[0]=1;
1186  (*hseries1)[1]=0;
1187  return hseries1;
1188  }
1189 
1190  #if 0
1191  if (wdegree == NULL)
1192  hWeight();
1193  else
1194  hWDegree(wdegree);
1195  #else
1196  if (wdegree != NULL) hWDegree(wdegree);
1197  #endif
1198 
1199  p0 = 1;
1200  hwork = (scfmon)omAlloc(hNexist * sizeof(scmon));
1201  hvar = (varset)omAlloc(((currRing->N) + 1) * sizeof(int));
1202  hpure = (scmon)omAlloc((1 + ((currRing->N) * (currRing->N))) * sizeof(int));
1203  stcmem = hCreate((currRing->N) - 1);
1204  Qpol = (int **)omAlloc(((currRing->N) + 1) * sizeof(int *));
1205  Ql = (int *)omAlloc0(((currRing->N) + 1) * sizeof(int));
1206  Q0 = (int *)omAlloc(((currRing->N) + 1) * sizeof(int));
1207  *Qpol = NULL;
1208  hLength = k = j = 0;
1209  mc = hisModule;
1210  if (mc!=0)
1211  {
1212  mw = hMinModulweight(modulweight);
1213  hstc = (scfmon)omAlloc(hNexist * sizeof(scmon));
1214  }
1215  else
1216  {
1217  mw = 0;
1218  hstc = hexist;
1219  hNstc = hNexist;
1220  }
1221  loop
1222  {
1223  if (mc!=0)
1224  {
1225  hComp(hexist, hNexist, mc, hstc, &hNstc);
1226  if (modulweight != NULL)
1227  j = (*modulweight)[mc-1]-mw;
1228  }
1229  if (hNstc!=0)
1230  {
1231  hNvar = (currRing->N);
1232  for (i = hNvar; i>=0; i--)
1233  hvar[i] = i;
1234  //if (notstc) // TODO: no mon divides another
1236  hSupp(hstc, hNstc, hvar, &hNvar);
1237  if (hNvar!=0)
1238  {
1239  if ((hNvar > 2) && (hNstc > 10))
1242  memset(hpure, 0, ((currRing->N) + 1) * sizeof(int));
1243  hPure(hstc, 0, &hNstc, hvar, hNvar, hpure, &hNpure);
1244  hLexS(hstc, hNstc, hvar, hNvar);
1245  Q0[hNvar] = 0;
1246  hHilbStep(hpure, hstc, hNstc, hvar, hNvar, &p0, 1);
1247  }
1248  }
1249  else
1250  {
1251  if(*Qpol!=NULL)
1252  (**Qpol)++;
1253  else
1254  {
1255  *Qpol = (int *)omAlloc(sizeof(int));
1256  hLength = *Ql = **Qpol = 1;
1257  }
1258  }
1259  if (*Qpol!=NULL)
1260  {
1261  i = hLength;
1262  while ((i > 0) && ((*Qpol)[i - 1] == 0))
1263  i--;
1264  if (i > 0)
1265  {
1266  l = i + j;
1267  if (l > k)
1268  {
1269  work = new intvec(l);
1270  for (ii=0; ii<k; ii++)
1271  (*work)[ii] = (*hseries1)[ii];
1272  if (hseries1 != NULL)
1273  delete hseries1;
1274  hseries1 = work;
1275  k = l;
1276  }
1277  while (i > 0)
1278  {
1279  (*hseries1)[i + j - 1] += (*Qpol)[i - 1];
1280  (*Qpol)[i - 1] = 0;
1281  i--;
1282  }
1283  }
1284  }
1285  mc--;
1286  if (mc <= 0)
1287  break;
1288  }
1289  if (k==0)
1290  {
1291  hseries1=new intvec(2);
1292  (*hseries1)[0]=0;
1293  (*hseries1)[1]=0;
1294  }
1295  else
1296  {
1297  l = k+1;
1298  while ((*hseries1)[l-2]==0) l--;
1299  if (l!=k)
1300  {
1301  work = new intvec(l);
1302  for (ii=l-2; ii>=0; ii--)
1303  (*work)[ii] = (*hseries1)[ii];
1304  delete hseries1;
1305  hseries1 = work;
1306  }
1307  (*hseries1)[l-1] = mw;
1308  }
1309  for (i = 0; i <= (currRing->N); i++)
1310  {
1311  if (Ql[i]!=0)
1312  omFreeSize((ADDRESS)Qpol[i], Ql[i] * sizeof(int));
1313  }
1314  omFreeSize((ADDRESS)Q0, ((currRing->N) + 1) * sizeof(int));
1315  omFreeSize((ADDRESS)Ql, ((currRing->N) + 1) * sizeof(int));
1316  omFreeSize((ADDRESS)Qpol, ((currRing->N) + 1) * sizeof(int *));
1317  hKill(stcmem, (currRing->N) - 1);
1318  omFreeSize((ADDRESS)hpure, (1 + ((currRing->N) * (currRing->N))) * sizeof(int));
1319  omFreeSize((ADDRESS)hvar, ((currRing->N) + 1) * sizeof(int));
1320  omFreeSize((ADDRESS)hwork, hNexist * sizeof(scmon));
1322  if (hisModule!=0)
1323  omFreeSize((ADDRESS)hstc, hNexist * sizeof(scmon));
1324  return hseries1;
1325 }
1326 
1327 
1328 intvec * hHstdSeries(ideal S, intvec *modulweight, intvec *wdegree, ideal Q, ring tailRing)
1329 {
1330  id_TestTail(S, currRing, tailRing);
1331  if (Q!=NULL) id_TestTail(Q, currRing, tailRing);
1332  return hSeries(S, modulweight, 0, wdegree, Q, tailRing);
1333 }
1334 
1335 intvec * hFirstSeries(ideal S, intvec *modulweight, ideal Q, intvec *wdegree, ring tailRing)
1336 {
1337  id_TestTail(S, currRing, tailRing);
1338  if (Q!= NULL) id_TestTail(Q, currRing, tailRing);
1339 
1340  intvec *hseries1= hSeries(S, modulweight, 1, wdegree, Q, tailRing);
1341  if (errorreported) { delete hseries1; hseries1=NULL; }
1342  return hseries1;
1343 }
1344 
1346 {
1347  intvec *work, *hseries2;
1348  int i, j, k, t, l;
1349  int s;
1350  if (hseries1 == NULL)
1351  return NULL;
1352  work = new intvec(hseries1);
1353  k = l = work->length()-1;
1354  s = 0;
1355  for (i = k-1; i >= 0; i--)
1356  s += (*work)[i];
1357  loop
1358  {
1359  if ((s != 0) || (k == 1))
1360  break;
1361  s = 0;
1362  t = (*work)[k-1];
1363  k--;
1364  for (i = k-1; i >= 0; i--)
1365  {
1366  j = (*work)[i];
1367  (*work)[i] = -t;
1368  s += t;
1369  t += j;
1370  }
1371  }
1372  hseries2 = new intvec(k+1);
1373  for (i = k-1; i >= 0; i--)
1374  (*hseries2)[i] = (*work)[i];
1375  (*hseries2)[k] = (*work)[l];
1376  delete work;
1377  return hseries2;
1378 }
1379 
1380 void hDegreeSeries(intvec *s1, intvec *s2, int *co, int *mu)
1381 {
1382  int i, j, k;
1383  int m;
1384  *co = *mu = 0;
1385  if ((s1 == NULL) || (s2 == NULL))
1386  return;
1387  i = s1->length();
1388  j = s2->length();
1389  if (j > i)
1390  return;
1391  m = 0;
1392  for(k=j-2; k>=0; k--)
1393  m += (*s2)[k];
1394  *mu = m;
1395  *co = i - j;
1396 }
1397 
1398 static void hPrintHilb(intvec *hseries)
1399 {
1400  int i, j, l, k;
1401  if (hseries == NULL)
1402  return;
1403  l = hseries->length()-1;
1404  k = (*hseries)[l];
1405  for (i = 0; i < l; i++)
1406  {
1407  j = (*hseries)[i];
1408  if (j != 0)
1409  {
1410  Print("// %8d t^%d\n", j, i+k);
1411  }
1412  }
1413 }
1414 
1415 /*
1416 *caller
1417 */
1418 void hLookSeries(ideal S, intvec *modulweight, ideal Q, intvec *wdegree, ring tailRing)
1419 {
1420  id_TestTail(S, currRing, tailRing);
1421 
1422  intvec *hseries1 = hFirstSeries(S, modulweight, Q, wdegree, tailRing);
1423  if (errorreported) return;
1424 
1425  hPrintHilb(hseries1);
1426 
1427  const int l = hseries1->length()-1;
1428 
1429  intvec *hseries2 = (l > 1) ? hSecondSeries(hseries1) : hseries1;
1430 
1431  int co, mu;
1432  hDegreeSeries(hseries1, hseries2, &co, &mu);
1433 
1434  PrintLn();
1435  hPrintHilb(hseries2);
1436  if ((l == 1) &&(mu == 0))
1437  scPrintDegree(rVar(currRing)+1, 0);
1438  else
1439  scPrintDegree(co, mu);
1440  if (l>1)
1441  delete hseries1;
1442  delete hseries2;
1443 }
1444 
1445 /***********************************************************************
1446  Computation of Hilbert series of non-commutative monomial algebras
1447 ************************************************************************/
1448 
1450 static void idInsertMonomial(ideal I, poly p)
1451 {
1452  /*
1453  * It adds monomial in I and if required,
1454  * enlarge the size of poly-set by 16.
1455  * It does not make copy of p.
1456  */
1457 
1458  if(I == NULL)
1459  {
1460  return;
1461  }
1462 
1463  int j = IDELEMS(I) - 1;
1464  while ((j >= 0) && (I->m[j] == NULL))
1465  {
1466  j--;
1467  }
1468  j++;
1469  if (j == IDELEMS(I))
1470  {
1471  pEnlargeSet(&(I->m), IDELEMS(I), 16);
1472  IDELEMS(I) +=16;
1473  }
1474  I->m[j] = p;
1475 }
1477 static int comapreMonoIdBases(ideal J, ideal Ob)
1478 {
1479  /*
1480  * Monomials of J and Ob are assumed to
1481  * be already sorted. J and Ob are
1482  * represented by the minimal generating set.
1483  */
1484  int i, s;
1485  s = 1;
1486  int JCount = IDELEMS(J);
1487  int ObCount = IDELEMS(Ob);
1488 
1489  if(idIs0(J))
1490  {
1491  return(1);
1492  }
1493  if(JCount != ObCount)
1494  {
1495  return(0);
1496  }
1497 
1498  for(i = 0; i < JCount; i++)
1499  {
1500  if(!(p_LmEqual(J->m[i], Ob->m[i], currRing)))
1501  {
1502  return(0);
1503  }
1504  }
1505  return(s);
1506 }
1508 static int CountOnIdUptoTruncationIndex(ideal I, int tr)
1509 {
1510  /*
1511  * The ideal I must be sorted in increasing total degree.
1512  * It counts the number of monomials in I up to
1513  * degree less than or equal to tr.
1514  */
1515 
1516  //case when I=1;
1517  if(p_Totaldegree(I->m[0], currRing) == 0)
1518  {
1519  return(1);
1520  }
1521 
1522  int count = 0;
1523  for(int i = 0; i < IDELEMS(I); i++)
1524  {
1525  if(p_Totaldegree(I->m[i], currRing) > tr)
1526  {
1527  return (count);
1528  }
1529  count = count + 1;
1530  }
1531 
1532  return(count);
1533 }
1535 static int comapreMonoIdBases_IG_Case(ideal J, int JCount, ideal Ob, int ObCount)
1536 {
1537  /*
1538  * Monomials of J and Ob are assumed to
1539  * be already sorted in increasing degrees.
1540  * J and Ob are represented by the minimal
1541  * generating set. It checks if J and Ob have
1542  * same monomials up to deg <=tr.
1543  */
1544 
1545  int i, s;
1546  s = 1;
1547  //when J is null
1548  //
1549  if(JCount != ObCount)
1550  {
1551  return(0);
1552  }
1553 
1554  if(JCount == 0)
1555  {
1556  return(1);
1557  }
1558 
1559  for(i = 0; i< JCount; i++)
1560  {
1561  if(!(p_LmEqual(J->m[i], Ob->m[i], currRing)))
1562  {
1563  return(0);
1564  }
1565  }
1566 
1567  return(s);
1568 }
1570 static int positionInOrbit_IG_Case(ideal I, poly w, std::vector<ideal> idorb, std::vector<poly> polist, int trInd, int trunDegHs)
1571 {
1572  /*
1573  * It compares the ideal I with ideals in the set 'idorb'
1574  * up to total degree =
1575  * trInd - max(deg of w, deg of word in polist) polynomials.
1576  *
1577  * It returns 0 if I is not equal to any ideal in the
1578  * 'idorb' else returns position of the matched ideal.
1579  */
1580 
1581  int ps = 0;
1582  int i, s = 0;
1583  int orbCount = idorb.size();
1584 
1585  if(idIs0(I))
1586  {
1587  return(1);
1588  }
1589 
1590  int degw = p_Totaldegree(w, currRing);
1591  int degp;
1592  int dtr;
1593  int dtrp;
1594 
1595  dtr = trInd - degw;
1596  int IwCount;
1597 
1598  IwCount = CountOnIdUptoTruncationIndex(I, dtr);
1599 
1600  if(IwCount == 0)
1601  {
1602  return(1);
1603  }
1604 
1605  int ObCount;
1606 
1607  bool flag2 = FALSE;
1608 
1609  for(i = 1;i < orbCount; i++)
1610  {
1611  degp = p_Totaldegree(polist[i], currRing);
1612  if(degw > degp)
1613  {
1614  dtr = trInd - degw;
1615 
1616  ObCount = 0;
1617  ObCount = CountOnIdUptoTruncationIndex(idorb[i], dtr);
1618  if(ObCount == 0)
1619  {continue;}
1620  if(flag2)
1621  {
1622  IwCount = 0;
1623  IwCount = CountOnIdUptoTruncationIndex(I, dtr);
1624  flag2 = FALSE;
1625  }
1626  }
1627  else
1628  {
1629  flag2 = TRUE;
1630  dtrp = trInd - degp;
1631  ObCount = 0;
1632  ObCount = CountOnIdUptoTruncationIndex(idorb[i], dtrp);
1633  IwCount = 0;
1634  IwCount = CountOnIdUptoTruncationIndex(I, dtrp);
1635  }
1636 
1637  s = comapreMonoIdBases_IG_Case(I, IwCount, idorb[i], ObCount);
1638 
1639  if(s)
1640  {
1641  ps = i + 1;
1642  break;
1643  }
1644  }
1645  return(ps);
1646 }
1648 static int positionInOrbit_FG_Case(ideal I, poly, std::vector<ideal> idorb, std::vector<poly>, int, int)
1649 {
1650  /*
1651  * It compares the ideal I with ideals in the set 'idorb'.
1652  * I and ideals of 'idorb' are sorted.
1653  *
1654  * It returns 0 if I is not equal to any ideal of 'idorb'
1655  * else returns position of the matched ideal.
1656  */
1657  int ps = 0;
1658  int i, s = 0;
1659  int OrbCount = idorb.size();
1660 
1661  if(idIs0(I))
1662  {
1663  return(1);
1664  }
1665 
1666  for(i = 1; i < OrbCount; i++)
1667  {
1668  s = comapreMonoIdBases(I, idorb[i]);
1669  if(s)
1670  {
1671  ps = i + 1;
1672  break;
1673  }
1674  }
1675 
1676  return(ps);
1677 }
1679 static int positionInOrbitTruncationCase(ideal I, poly w, std::vector<ideal> idorb, std::vector<poly> polist, int , int trunDegHs)
1680 {
1681  /*
1682  * It compares the ideal I with ideals in the set 'idorb'.
1683  * I and ideals in 'idorb' are sorted.
1684 
1685  * returns 0 if I is not equal to any ideal of 'idorb'
1686  * else returns position of the matched ideal.
1687  */
1688 
1689  int ps = 0;
1690  int i, s = 0;
1691  int OrbCount = idorb.size();
1692  int dtr=0; int IwCount, ObCount;
1693  dtr = trunDegHs - 1 - p_Totaldegree(w, currRing);
1694 
1695  if(idIs0(I))
1696  {
1697  for(i = 1; i < OrbCount; i++)
1698  {
1699  if(p_Totaldegree(w, currRing) == p_Totaldegree(polist[i], currRing))
1700  {
1701  if(idIs0(idorb[i]))
1702  return(i+1);
1703  ObCount=0;
1704  ObCount = CountOnIdUptoTruncationIndex(idorb[i], dtr);
1705  if(ObCount==0)
1706  {
1707  ps = i + 1;
1708  break;
1709  }
1710  }
1711  }
1712 
1713  return(ps);
1714  }
1715 
1716  IwCount = CountOnIdUptoTruncationIndex(I, dtr);
1717 
1718  if(p_Totaldegree(I->m[0], currRing)==0)
1719  {
1720  for(i = 1; i < OrbCount; i++)
1721  {
1722  if(idIs0(idorb[i]))
1723  continue;
1724  if(p_Totaldegree(idorb[i]->m[0], currRing)==0)
1725  {
1726  ps = i + 1;
1727  break;
1728  }
1729  }
1730  return(ps);
1731  }
1732 
1733  for(i = 1; i < OrbCount; i++)
1734  {
1735  if(p_Totaldegree(w, currRing) == p_Totaldegree(polist[i], currRing))
1736  {
1737  if(idIs0(idorb[i]))
1738  continue;
1739  ObCount=0;
1740  ObCount = CountOnIdUptoTruncationIndex(idorb[i], dtr);
1741  s = comapreMonoIdBases_IG_Case(I, IwCount, idorb[i], ObCount);
1742  if(s)
1743  {
1744  ps = i + 1;
1745  break;
1746  }
1747  }
1748  }
1749 
1750  return(ps);
1751 }
1753 static int monCompare( const void *m, const void *n)
1754 {
1755  /* compares monomials */
1756 
1757  return(p_Compare(*(poly*) m, *(poly*)n, currRing));
1758 }
1760 void sortMonoIdeal_pCompare(ideal I)
1761 {
1762  /*
1763  * sorts monomial ideal in ascending order
1764  * order must be a total degree
1765  */
1766 
1767  qsort(I->m, IDELEMS(I), sizeof(poly), monCompare);
1768 
1769 }
1770 
1771 
1773 static ideal minimalMonomialGenSet(ideal I)
1774 {
1775  /*
1776  * eliminates monomials which
1777  * can be generated by others in I
1778  */
1779  //first sort monomials of the ideal
1780 
1781  idSkipZeroes(I);
1782 
1784 
1785  int i, k;
1786  int ICount = IDELEMS(I);
1787 
1788  for(k = ICount - 1; k >=1; k--)
1789  {
1790  for(i = 0; i < k; i++)
1791  {
1792 
1793  if(p_LmDivisibleBy(I->m[i], I->m[k], currRing))
1794  {
1795  pDelete(&(I->m[k]));
1796  break;
1797  }
1798  }
1799  }
1800 
1801  idSkipZeroes(I);
1802  return(I);
1803 }
1805 static poly shiftInMon(poly p, int i, int lV, const ring r)
1806 {
1807  /*
1808  * shifts the varibles of monomial p in the i^th layer,
1809  * p remains unchanged,
1810  * creates new poly and returns it for the colon ideal
1811  */
1812  poly smon = p_One(r);
1813  int j, sh, cnt;
1814  cnt = r->N;
1815  sh = i*lV;
1816  int *e=(int *)omAlloc((r->N+1)*sizeof(int));
1817  int *s=(int *)omAlloc0((r->N+1)*sizeof(int));
1818  p_GetExpV(p, e, r);
1819 
1820  for(j = 1; j <= cnt; j++)
1821  {
1822  if(e[j] == 1)
1823  {
1824  s[j+sh] = e[j];
1825  }
1826  }
1827 
1828  p_SetExpV(smon, s, currRing);
1829  omFree(e);
1830  omFree(s);
1831 
1833  p_Setm(smon, currRing);
1834 
1835  return(smon);
1836 }
1838 static poly deleteInMon(poly w, int i, int lV, const ring r)
1839 {
1840  /*
1841  * deletes the variables upto i^th layer of monomial w
1842  * w remains unchanged
1843  * creates new poly and returns it for the colon ideal
1844  */
1845 
1846  poly dw = p_One(currRing);
1847  int *e = (int *)omAlloc((r->N+1)*sizeof(int));
1848  int *s=(int *)omAlloc0((r->N+1)*sizeof(int));
1849  p_GetExpV(w, e, r);
1850  int j, cnt;
1851  cnt = i*lV;
1852  /*
1853  for(j=1;j<=cnt;j++)
1854  {
1855  e[j]=0;
1856  }*/
1857  for(j = (cnt+1); j < (r->N+1); j++)
1858  {
1859  s[j] = e[j];
1860  }
1861 
1862  p_SetExpV(dw, s, currRing);//new exponents
1863  omFree(e);
1864  omFree(s);
1865 
1867  p_Setm(dw, currRing);
1868 
1869  return(dw);
1870 }
1872 static void TwordMap(poly p, poly w, int lV, int d, ideal Jwi, bool &flag)
1873 {
1874  /*
1875  * computes T_w(p) in a new poly object and places it
1876  * in Jwi which stores elements of colon ideal of I,
1877  * p and w remain unchanged,
1878  * the new polys for Jwi are constructed by sub-routines
1879  * deleteInMon, shiftInMon, p_MDivide,
1880  * places the result in Jwi and deletes the new polys
1881  * coming in dw, smon, qmon
1882  */
1883  int i;
1884  poly smon, dw;
1885  poly qmonp = NULL;
1886  bool del;
1887 
1888  for(i = 0;i <= d - 1; i++)
1889  {
1890  dw = deleteInMon(w, i, lV, currRing);
1891  smon = shiftInMon(p, i, lV, currRing);
1892  del = TRUE;
1893 
1894  if(pLmDivisibleBy(smon, w))
1895  {
1896  flag = TRUE;
1897  del = FALSE;
1898 
1899  pDelete(&dw);
1900  pDelete(&smon);
1901 
1902  //delete all monomials of Jwi
1903  //and make Jwi =1
1904 
1905  for(int j = 0;j < IDELEMS(Jwi); j++)
1906  {
1907  pDelete(&Jwi->m[j]);
1908  }
1909 
1911  break;
1912  }
1913 
1914  if(pLmDivisibleBy(dw, smon))
1915  {
1916  del = FALSE;
1917  qmonp = p_MDivide(smon, dw, currRing);
1918  idInsertMonomial(Jwi, shiftInMon(qmonp, -d, lV, currRing));
1919  pLmFree(&qmonp);
1920  pDelete(&dw);
1921  pDelete(&smon);
1922  }
1923  //in case both if are false, delete dw and smon
1924  if(del)
1925  {
1926  pDelete(&dw);
1927  pDelete(&smon);
1928  }
1929  }
1930 
1931 }
1933 static ideal colonIdeal(ideal S, poly w, int lV, ideal Jwi, int trunDegHs)
1934 {
1935  /*
1936  * It computes the right colon ideal of a two-sided ideal S
1937  * w.r.t. word w and save it in a new object Jwi.
1938  * It keeps S and w unchanged.
1939  */
1940 
1941  if(idIs0(S))
1942  {
1943  return(S);
1944  }
1945 
1946  int i, d;
1947  d = p_Totaldegree(w, currRing);
1948  if(trunDegHs !=0 && d >= trunDegHs)
1949  {
1951  return(Jwi);
1952  }
1953  bool flag = FALSE;
1954  int SCount = IDELEMS(S);
1955  for(i = 0; i < SCount; i++)
1956  {
1957  TwordMap(S->m[i], w, lV, d, Jwi, flag);
1958  if(flag)
1959  {
1960  break;
1961  }
1962  }
1963 
1964  Jwi = minimalMonomialGenSet(Jwi);
1965  return(Jwi);
1966 }
1968 void HilbertSeries_OrbitData(ideal S, int lV, bool IG_CASE, bool mgrad, bool odp, int trunDegHs)
1969 {
1970  /*
1971  * This is based on iterative right colon operations on a
1972  * two-sided monomial ideal of the free associative algebra.
1973  * The algorithm terminates for those monomial ideals
1974  * whose monomials define "regular formal languages",
1975  * that is, all monomials of the input ideal can be obtained
1976  * from finite languages by applying finite number of
1977  * rational operations.
1978  */
1979 
1980  int trInd;
1981  S = minimalMonomialGenSet(S);
1982  if( !idIs0(S) && p_Totaldegree(S->m[0], currRing)==0)
1983  {
1984  PrintS("Hilbert Series:\n 0\n");
1985  return;
1986  }
1987  int (*POS)(ideal, poly, std::vector<ideal>, std::vector<poly>, int, int);
1988  if(trunDegHs != 0)
1989  {
1990  Print("\nTruncation degree = %d\n",trunDegHs);
1992  }
1993  else
1994  {
1995  if(IG_CASE)
1996  {
1997  if(idIs0(S))
1998  {
1999  WerrorS("wrong input: it is not an infinitely gen. case");
2000  return;
2001  }
2002  trInd = p_Totaldegree(S->m[IDELEMS(S)-1], currRing);
2003  POS = &positionInOrbit_IG_Case;
2004  }
2005  else
2006  POS = &positionInOrbit_FG_Case;
2007  }
2008  std::vector<ideal > idorb;
2009  std::vector< poly > polist;
2010 
2011  ideal orb_init = idInit(1, 1);
2012  idorb.push_back(orb_init);
2013 
2014  polist.push_back( p_One(currRing));
2015 
2016  std::vector< std::vector<int> > posMat;
2017  std::vector<int> posRow(lV,0);
2018  std::vector<int> C;
2019 
2020  int ds, is, ps;
2021  int lpcnt = 0;
2022 
2023  poly w, wi;
2024  ideal Jwi;
2025 
2026  while(lpcnt < idorb.size())
2027  {
2028  w = NULL;
2029  w = polist[lpcnt];
2030  if(lpcnt >= 1 && idIs0(idorb[lpcnt]) == FALSE)
2031  {
2032  if(p_Totaldegree(idorb[lpcnt]->m[0], currRing) != 0)
2033  {
2034  C.push_back(1);
2035  }
2036  else
2037  C.push_back(0);
2038  }
2039  else
2040  {
2041  C.push_back(1);
2042  }
2043 
2044  ds = p_Totaldegree(w, currRing);
2045  lpcnt++;
2046 
2047  for(is = 1; is <= lV; is++)
2048  {
2049  wi = NULL;
2050  //make new copy 'wi' of word w=polist[lpcnt]
2051  //and update it (for the colon operation).
2052  //if corresponding to wi, right colon operation gives
2053  //a new (right colon) ideal of S,
2054  //keep 'wi' in the polist else delete it
2055 
2056  wi = pCopy(w);
2057  p_SetExp(wi, (ds*lV)+is, 1, currRing);
2058  p_Setm(wi, currRing);
2059  Jwi = NULL;
2060  //Jwi stores (right) colon ideal of S w.r.t. word
2061  //wi if colon operation gives a new ideal place it
2062  //in the vector of ideals 'idorb'
2063  //otherwise delete it
2064 
2065  Jwi = idInit(1,1);
2066 
2067  Jwi = colonIdeal(S, wi, lV, Jwi, trunDegHs);
2068  ps = (*POS)(Jwi, wi, idorb, polist, trInd, trunDegHs);
2069 
2070  if(ps == 0) // finds a new ideal
2071  {
2072  posRow[is-1] = idorb.size();
2073 
2074  idorb.push_back(Jwi);
2075  polist.push_back(wi);
2076  }
2077  else // ideal is already there in the set
2078  {
2079  posRow[is-1]=ps-1;
2080  idDelete(&Jwi);
2081  pDelete(&wi);
2082  }
2083  }
2084  posMat.push_back(posRow);
2085  posRow.resize(lV,0);
2086  }
2087  int lO = C.size();//size of the orbit
2088  PrintLn();
2089  Print("maximal length of words = %ld\n", p_Totaldegree(polist[lO-1], currRing));
2090  Print("\nlength of the Orbit = %d", lO);
2091  PrintLn();
2092 
2093  if(odp)
2094  {
2095  Print("words description of the Orbit: \n");
2096  for(is = 0; is < lO; is++)
2097  {
2098  pWrite0(polist[is]);
2099  PrintS(" ");
2100  }
2101  PrintLn();
2102  PrintS("\nmaximal degree, #(sum_j R(w,w_j))");
2103  PrintLn();
2104  for(is = 0; is < lO; is++)
2105  {
2106  if(idIs0(idorb[is]))
2107  {
2108  PrintS("NULL\n");
2109  }
2110  else
2111  {
2112  Print("%ld, %d \n",p_Totaldegree(idorb[is]->m[IDELEMS(idorb[is])-1], currRing),IDELEMS(idorb[is]));
2113  }
2114  }
2115  }
2116 
2117  for(is = idorb.size()-1; is >= 0; is--)
2118  {
2119  idDelete(&idorb[is]);
2120  }
2121  for(is = polist.size()-1; is >= 0; is--)
2122  {
2123  pDelete(&polist[is]);
2124  }
2125 
2126  idorb.resize(0);
2127  polist.resize(0);
2128 
2129  int adjMatrix[lO][lO];
2130  memset(adjMatrix, 0, lO*lO*sizeof(int));
2131  int rowCount, colCount;
2132  int tm = 0;
2133  if(!mgrad)
2134  {
2135  for(rowCount = 0; rowCount < lO; rowCount++)
2136  {
2137  for(colCount = 0; colCount < lV; colCount++)
2138  {
2139  tm = posMat[rowCount][colCount];
2140  adjMatrix[rowCount][tm] = adjMatrix[rowCount][tm] + 1;
2141  }
2142  }
2143  }
2144 
2145  ring r = currRing;
2146  int npar;
2147  char** tt;
2148  TransExtInfo p;
2149  if(!mgrad)
2150  {
2151  tt=(char**)omAlloc(sizeof(char*));
2152  tt[0] = omStrDup("t");
2153  npar = 1;
2154  }
2155  else
2156  {
2157  tt=(char**)omalloc(lV*sizeof(char*));
2158  for(is = 0; is < lV; is++)
2159  {
2160  tt[is] = (char*)omAlloc(7*sizeof(char)); //if required enlarge it later
2161  sprintf (tt[is], "t%d", is+1);
2162  }
2163  npar = lV;
2164  }
2165 
2166  p.r = rDefault(0, npar, tt);
2168  char** xx = (char**)omAlloc(sizeof(char*));
2169  xx[0] = omStrDup("x");
2170  ring R = rDefault(cf, 1, xx);
2171  rChangeCurrRing(R);//rWrite(R);
2172  /*
2173  * matrix corresponding to the orbit of the ideal
2174  */
2175  matrix mR = mpNew(lO, lO);
2176  matrix cMat = mpNew(lO,1);
2177  poly rc;
2178 
2179  if(!mgrad)
2180  {
2181  for(rowCount = 0; rowCount < lO; rowCount++)
2182  {
2183  for(colCount = 0; colCount < lO; colCount++)
2184  {
2185  if(adjMatrix[rowCount][colCount] != 0)
2186  {
2187  MATELEM(mR, rowCount + 1, colCount + 1) = p_ISet(adjMatrix[rowCount][colCount], R);
2188  p_SetCoeff(MATELEM(mR, rowCount + 1, colCount + 1), n_Mult(pGetCoeff(mR->m[lO*rowCount+colCount]),n_Param(1, R->cf), R->cf), R);
2189  }
2190  }
2191  }
2192  }
2193  else
2194  {
2195  for(rowCount = 0; rowCount < lO; rowCount++)
2196  {
2197  for(colCount = 0; colCount < lV; colCount++)
2198  {
2199  rc=NULL;
2200  rc=p_One(R);
2201  p_SetCoeff(rc, n_Mult(pGetCoeff(rc), n_Param(colCount+1, R->cf),R->cf), R);
2202  MATELEM(mR, rowCount +1, posMat[rowCount][colCount]+1)=p_Add_q(rc,MATELEM(mR, rowCount +1, posMat[rowCount][colCount]+1), R);
2203  }
2204  }
2205  }
2206 
2207  for(rowCount = 0; rowCount < lO; rowCount++)
2208  {
2209  if(C[rowCount] != 0)
2210  {
2211  MATELEM(cMat, rowCount + 1, 1) = p_ISet(C[rowCount], R);
2212  }
2213  }
2214 
2215  matrix u;
2216  unitMatrix(lO, u); //unit matrix
2217  matrix gMat = mp_Sub(u, mR, R);
2218 
2219  char* s;
2220 
2221  if(odp)
2222  {
2223  PrintS("\nlinear system:\n");
2224  if(!mgrad)
2225  {
2226  for(rowCount = 0; rowCount < lO; rowCount++)
2227  {
2228  Print("H(%d) = ", rowCount+1);
2229  for(colCount = 0; colCount < lV; colCount++)
2230  {
2231  StringSetS(""); nWrite(n_Param(1, R->cf));
2232  s = StringEndS(); PrintS(s);
2233  Print("*"); omFree(s);
2234  Print("H(%d) + ", posMat[rowCount][colCount] + 1);
2235  }
2236  Print(" %d\n", C[rowCount] );
2237  }
2238  PrintS("where H(1) represents the series corresp. to input ideal\n");
2239  PrintS("and i^th summand in the rhs of an eqn. is according\n");
2240  PrintS("to the right colon map corresp. to the i^th variable\n");
2241  }
2242  else
2243  {
2244  for(rowCount = 0; rowCount < lO; rowCount++)
2245  {
2246  Print("H(%d) = ", rowCount+1);
2247  for(colCount = 0; colCount < lV; colCount++)
2248  {
2249  StringSetS(""); nWrite(n_Param(colCount+1, R->cf));
2250  s = StringEndS(); PrintS(s);
2251  Print("*");omFree(s);
2252  Print("H(%d) + ", posMat[rowCount][colCount] + 1);
2253  }
2254  Print(" %d\n", C[rowCount] );
2255  }
2256  PrintS("where H(1) represents the series corresp. to input ideal\n");
2257  }
2258  }
2259  PrintLn();
2260  posMat.resize(0);
2261  C.resize(0);
2262  matrix pMat;
2263  matrix lMat;
2264  matrix uMat;
2265  matrix H_serVec = mpNew(lO, 1);
2266  matrix Hnot;
2267 
2268  //std::clock_t start;
2269  //start = std::clock();
2270 
2271  luDecomp(gMat, pMat, lMat, uMat, R);
2272  luSolveViaLUDecomp(pMat, lMat, uMat, cMat, H_serVec, Hnot);
2273 
2274  //to print system solving time
2275  //if(odp){
2276  //std::cout<<"solving time of the system = "<<(std::clock()-start)/(double)(CLOCKS_PER_SEC / 1000)<<" ms"<<std::endl;}
2277 
2278  mp_Delete(&mR, R);
2279  mp_Delete(&u, R);
2280  mp_Delete(&pMat, R);
2281  mp_Delete(&lMat, R);
2282  mp_Delete(&uMat, R);
2283  mp_Delete(&cMat, R);
2284  mp_Delete(&gMat, R);
2285  mp_Delete(&Hnot, R);
2286  //print the Hilbert series and length of the Orbit
2287  PrintLn();
2288  Print("Hilbert series:");
2289  PrintLn();
2290  pWrite(H_serVec->m[0]);
2291  if(!mgrad)
2292  {
2293  omFree(tt[0]);
2294  }
2295  else
2296  {
2297  for(is = lV-1; is >= 0; is--)
2298 
2299  omFree( tt[is]);
2300  }
2301  omFree(tt);
2302  omFree(xx[0]);
2303  omFree(xx);
2304  rChangeCurrRing(r);
2305  rKill(R);
2306 }
2308 ideal RightColonOperation(ideal S, poly w, int lV)
2309 {
2310  /*
2311  * This returns right colon ideal of a monomial two-sided ideal of
2312  * the free associative algebra with respect to a monomial 'w'
2313  * (S:_R w).
2314  */
2315  S = minimalMonomialGenSet(S);
2316  ideal Iw = idInit(1,1);
2317  Iw = colonIdeal(S, w, lV, Iw, 0);
2318  return (Iw);
2319 }
2320 
hStaircase
void hStaircase(scfmon stc, int *Nstc, varset var, int Nvar)
Definition: hutil.cc:315
FALSE
#define FALSE
Definition: auxiliary.h:96
hutil.h
hMinModulweight
static int hMinModulweight(intvec *modulweight)
Definition: hilb.cc:49
hElimS
void hElimS(scfmon stc, int *e1, int a2, int e2, varset var, int Nvar)
Definition: hutil.cc:674
p_GetComp
#define p_GetComp(p, r)
Definition: monomials.h:61
hPrintHilb
static void hPrintHilb(intvec *hseries)
Definition: hilb.cc:1398
ip_smatrix
Definition: matpol.h:13
pLmEqual
#define pLmEqual(p1, p2)
Definition: polys.h:107
TransExtInfo
struct for passing initialization parameters to naInitChar
Definition: transext.h:87
hvar
VAR varset hvar
Definition: hutil.cc:17
hexist
VAR scfmon hexist
Definition: hutil.cc:15
monCompare
static int monCompare(const void *m, const void *n)
Definition: hilb.cc:1752
p_GetExp
static long p_GetExp(const poly p, const unsigned long iBitmask, const int VarOffset)
get a single variable exponent @Note: the integer VarOffset encodes:
Definition: p_polys.h:456
j
int j
Definition: facHensel.cc:105
omFree
#define omFree(addr)
Definition: omAllocDecl.h:259
Ql
STATIC_VAR int * Ql
Definition: hilb.cc:45
pWrite0
void pWrite0(poly p)
Definition: polys.h:293
k
int k
Definition: cfEzgcd.cc:92
idDelete
#define idDelete(H)
delete an ideal
Definition: ideals.h:28
pLmFree
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:68
hpure
VAR scmon hpure
Definition: hutil.cc:16
x
Variable x
Definition: cfModGcd.cc:4023
MATELEM
#define MATELEM(mat, i, j)
1-based access to matrix
Definition: matpol.h:29
y
const CanonicalForm int const CFList const Variable & y
Definition: facAbsFact.cc:57
minimalMonomialGenSet
static ideal minimalMonomialGenSet(ideal I)
Definition: hilb.cc:1772
rChangeCurrRing
void rChangeCurrRing(ring r)
Definition: polys.cc:15
nWrite
#define nWrite(n)
Definition: numbers.h:28
SearchP
static poly SearchP(ideal I)
searches for a monomial of degree d>=2 and divides it by a variable (result monomial of deg d-1)
Definition: hilb.cc:780
unitMatrix
bool unitMatrix(const int n, matrix &unitMat, const ring R)
Creates a new matrix which is the (nxn) unit matrix, and returns true in case of success.
Definition: linearAlgebra.cc:95
idInsertMonomial
static void idInsertMonomial(ideal I, poly p)
Definition: hilb.cc:1449
pEnlargeSet
void pEnlargeSet(poly **p, int l, int increment)
Definition: p_polys.cc:3653
hWDegree
static void hWDegree(intvec *wdegree)
Definition: hilb.cc:245
positionInOrbit_FG_Case
static int positionInOrbit_FG_Case(ideal I, poly, std::vector< ideal > idorb, std::vector< poly >, int, int)
Definition: hilb.cc:1647
ADDRESS
void * ADDRESS
Definition: auxiliary.h:135
p_Head
static poly p_Head(poly p, const ring r)
copy the i(leading) term of p
Definition: p_polys.h:810
cf
CanonicalForm cf
Definition: cfModGcd.cc:4024
simpleideals.h
shiftInMon
static poly shiftInMon(poly p, int i, int lV, const ring r)
Definition: hilb.cc:1804
STATIC_VAR
#define STATIC_VAR
Definition: globaldefs.h:7
idInsertPoly
BOOLEAN idInsertPoly(ideal h1, poly h2)
insert h2 into h1 (if h2 is not the zero polynomial) return TRUE iff h2 was indeed inserted
Definition: simpleideals.cc:648
rKill
void rKill(ring r)
Definition: ipshell.cc:6123
luSolveViaLUDecomp
bool luSolveViaLUDecomp(const matrix pMat, const matrix lMat, const matrix uMat, const matrix bVec, matrix &xVec, matrix &H)
Solves the linear system A * x = b, where A is an (m x n)-matrix which is given by its LU-decompositi...
Definition: linearAlgebra.cc:377
hNexist
VAR int hNexist
Definition: hutil.cc:18
omStrDup
#define omStrDup(s)
Definition: omAllocDecl.h:261
ChooseP
static poly ChooseP(ideal I)
Definition: hilb.cc:764
p_IsOne
static BOOLEAN p_IsOne(const poly p, const ring R)
either poly(1) or gen(k)?!
Definition: p_polys.h:1917
nInitChar
coeffs nInitChar(n_coeffType t, void *parameter)
one-time initialisations for new coeffs in case of an error return NULL
Definition: numbers.cc:349
pp_Mult_mm
static poly pp_Mult_mm(poly p, poly m, const ring r)
Definition: p_polys.h:975
comapreMonoIdBases
static int comapreMonoIdBases(ideal J, ideal Ob)
Definition: hilb.cc:1476
hGetmem
scfmon hGetmem(int lm, scfmon old, monp monmem)
Definition: hutil.cc:1025
pDelete
#define pDelete(p_ptr)
Definition: polys.h:175
N
const CanonicalForm CFMap CFMap & N
Definition: cfEzgcd.cc:48
InternalPoly::var
Variable var
Definition: int_poly.h:74
IsIn
static bool IsIn(poly p, ideal I)
Definition: hilb.cc:909
scfmon
scmon * scfmon
Definition: hutil.h:14
n_Param
static FORCE_INLINE number n_Param(const int iParameter, const coeffs r)
return the (iParameter^th) parameter as a NEW number NOTE: parameter numbering: 1....
Definition: coeffs.h:804
StringEndS
char * StringEndS()
Definition: reporter.cc:150
idIs0
BOOLEAN idIs0(ideal h)
returns true if h is the zero ideal
Definition: simpleideals.cc:776
loop
#define loop
Definition: structs.h:79
w
const CanonicalForm & w
Definition: facAbsFact.cc:55
b
CanonicalForm b
Definition: cfModGcd.cc:4044
p_LmDivisibleBy
static BOOLEAN p_LmDivisibleBy(poly a, poly b, const ring r)
Definition: p_polys.h:1802
hHilbStep
static void hHilbStep(scmon pure, scfmon stc, int Nstc, varset var, int Nvar, int *pol, int Lpol)
Definition: hilb.cc:177
p_LmEqual
#define p_LmEqual(p1, p2, r)
Definition: p_polys.h:1631
idAddMon
static void idAddMon(ideal I, ideal p)
Definition: hilb.cc:471
hNvar
VAR int hNvar
Definition: hutil.cc:18
p_SetExp
static unsigned long p_SetExp(poly p, const unsigned long e, const unsigned long iBitmask, const int VarOffset)
set a single variable exponent @Note: VarOffset encodes the position in p->exp
Definition: p_polys.h:475
mu
void mu(int **points, int sizePoints)
Definition: cfNewtonPolygon.cc:467
HilbertSeries_OrbitData
void HilbertSeries_OrbitData(ideal S, int lV, bool IG_CASE, bool mgrad, bool odp, int trunDegHs)
Definition: hilb.cc:1967
rDefault
ring rDefault(const coeffs cf, int N, char **n, int ord_size, rRingOrder_t *ord, int *block0, int *block1, int **wvhdl, unsigned long bitmask)
Definition: ring.cc:101
deleteInMon
static poly deleteInMon(poly w, int i, int lV, const ring r)
Definition: hilb.cc:1837
omRealloc
#define omRealloc(addr, size)
Definition: omAllocDecl.h:223
p_SetExpV
static void p_SetExpV(poly p, int *ev, const ring r)
Definition: p_polys.h:1456
p_Copy
static poly p_Copy(poly p, const ring r)
returns a copy of p
Definition: p_polys.h:796
stairc.h
rVar
static short rVar(const ring r)
#define rVar(r) (r->N)
Definition: ring.h:586
hNstc
VAR int hNstc
Definition: hutil.cc:18
TRUE
#define TRUE
Definition: auxiliary.h:100
i
int i
Definition: cfEzgcd.cc:125
id_Delete
void id_Delete(ideal *h, ring r)
deletes an ideal/module/matrix
Definition: simpleideals.cc:122
res
CanonicalForm res
Definition: facAbsFact.cc:64
hSeries
static intvec * hSeries(ideal S, intvec *modulweight, int, intvec *wdegree, ideal Q, ring tailRing)
Definition: hilb.cc:1172
CountOnIdUptoTruncationIndex
static int CountOnIdUptoTruncationIndex(ideal I, int tr)
Definition: hilb.cc:1507
colonIdeal
static ideal colonIdeal(ideal S, poly w, int lV, ideal Jwi, int trunDegHs)
Definition: hilb.cc:1932
PrintS
void PrintS(const char *s)
Definition: reporter.cc:283
variables
STATIC_VAR int variables
Definition: interpolation.cc:87
omFreeSize
#define omFreeSize(addr, size)
Definition: omAllocDecl.h:258
hilb.h
p_MDivide
poly p_MDivide(poly a, poly b, const ring r)
Definition: p_polys.cc:1471
idSkipZeroes
void idSkipZeroes(ideal ide)
gives an ideal/module the minimal possible size
Definition: simpleideals.cc:180
TwordMap
static void TwordMap(poly p, poly w, int lV, int d, ideal Jwi, bool &flag)
Definition: hilb.cc:1871
RightColonOperation
ideal RightColonOperation(ideal S, poly w, int lV)
Definition: hilb.cc:2307
hComp
void hComp(scfmon exist, int Nexist, int ak, scfmon stc, int *Nstc)
Definition: hutil.cc:156
mod2.h
ip_smatrix::m
poly * m
Definition: matpol.h:18
hKill
void hKill(monf xmem, int Nvar)
Definition: hutil.cc:1012
coeffs
pOne
#define pOne()
Definition: polys.h:299
Qpol
STATIC_VAR int ** Qpol
Definition: hilb.cc:44
prune
void prune(Variable &alpha)
Definition: variable.cc:261
intvec
Definition: intvec.h:18
positionInOrbitTruncationCase
static int positionInOrbitTruncationCase(ideal I, poly w, std::vector< ideal > idorb, std::vector< poly > polist, int, int trunDegHs)
Definition: hilb.cc:1678
Q
STATIC_VAR jList * Q
Definition: janet.cc:30
SortByDeg
static ideal SortByDeg(ideal I)
Definition: hilb.cc:388
n_Mult
static FORCE_INLINE number n_Mult(number a, number b, const coeffs r)
return the product of 'a' and 'b', i.e., a*b
Definition: coeffs.h:635
omAlloc
#define omAlloc(size)
Definition: omAllocDecl.h:208
p_polys.h
p_Compare
int p_Compare(const poly a, const poly b, const ring R)
Definition: p_polys.cc:4812
hHstdSeries
intvec * hHstdSeries(ideal S, intvec *modulweight, intvec *wdegree, ideal Q, ring tailRing)
Definition: hilb.cc:1328
hDelete
void hDelete(scfmon ev, int ev_length)
Definition: hutil.cc:142
p_DivisibleBy
static BOOLEAN p_DivisibleBy(poly a, poly b, const ring r)
Definition: p_polys.h:1811
ChoosePJL
static poly ChoosePJL(ideal I)
Definition: hilb.cc:705
p_GetExpV
static void p_GetExpV(poly p, int *ev, const ring r)
Definition: p_polys.h:1441
rouneslice
void rouneslice(ideal I, ideal S, poly q, poly x, int &prune, int &moreprune, int &steps, int &NNN, mpz_ptr &hilbertcoef, int *&hilbpower)
Definition: hilb.cc:974
mp_Delete
void mp_Delete(matrix *a, const ring r)
Definition: matpol.cc:879
SqFree
static poly SqFree(ideal I)
Definition: hilb.cc:880
hStepS
void hStepS(scfmon stc, int Nstc, varset var, int Nvar, int *a, int *x)
Definition: hutil.cc:951
eulerchar
static void eulerchar(ideal I, int variables, mpz_ptr ec)
Definition: hilb.cc:837
intvec.h
hLastHilb
static void hLastHilb(scmon pure, int Nv, varset var, int *pol, int lp)
Definition: hilb.cc:137
JustVar
static bool JustVar(ideal I)
Definition: hilb.cc:806
scmon
int * scmon
Definition: hutil.h:13
mpNew
matrix mpNew(int r, int c)
create a r x c zero-matrix
Definition: matpol.cc:36
n_transExt
used for all transcendental extensions, i.e., the top-most extension in an extension tower is transce...
Definition: coeffs.h:38
hstc
VAR scfmon hstc
Definition: hutil.cc:15
hSecondSeries
intvec * hSecondSeries(intvec *hseries1)
Definition: hilb.cc:1345
hisModule
VAR int hisModule
Definition: hutil.cc:19
exp
gmp_float exp(const gmp_float &a)
Definition: mpr_complex.cc:356
hPure
void hPure(scfmon stc, int a, int *Nstc, varset var, int Nvar, scmon pure, int *Npure)
Definition: hutil.cc:623
hLength
STATIC_VAR int hLength
Definition: hilb.cc:46
hLex2S
void hLex2S(scfmon rad, int e1, int a2, int e2, varset var, int Nvar, scfmon w)
Definition: hutil.cc:814
ring.h
varset
int * varset
Definition: hutil.h:15
transext.h
hLookSeries
void hLookSeries(ideal S, intvec *modulweight, ideal Q, intvec *wdegree, ring tailRing)
Definition: hilb.cc:1418
p_Delete
static void p_Delete(poly *p, const ring r)
Definition: p_polys.h:845
p_Add_q
static poly p_Add_q(poly p, poly q, const ring r)
Definition: p_polys.h:880
hInit
scfmon hInit(ideal S, ideal Q, int *Nexist, ring tailRing)
Definition: hutil.cc:30
stcmem
VAR monf stcmem
Definition: hutil.cc:20
hFirstSeries
intvec * hFirstSeries(ideal S, intvec *modulweight, ideal Q, intvec *wdegree, ring tailRing)
Definition: hilb.cc:1335
p_One
poly p_One(const ring r)
Definition: p_polys.cc:1300
idPrint
#define idPrint(id)
Definition: ideals.h:45
LCMmon
static poly LCMmon(ideal I)
Definition: hilb.cc:948
StringSetS
void StringSetS(const char *st)
Definition: reporter.cc:127
hOrdSupp
void hOrdSupp(scfmon stc, int Nstc, varset var, int Nvar)
Definition: hutil.cc:204
sortMonoIdeal_pCompare
void sortMonoIdeal_pCompare(ideal I)
Definition: hilb.cc:1759
Print
#define Print
Definition: emacs.cc:79
omalloc
#define omalloc(size)
Definition: omAllocDecl.h:226
mylimits.h
hSupp
void hSupp(scfmon stc, int Nstc, varset var, int *Nvar)
Definition: hutil.cc:176
comapreMonoIdBases_IG_Case
static int comapreMonoIdBases_IG_Case(ideal J, int JCount, ideal Ob, int ObCount)
Definition: hilb.cc:1534
int64
long int64
Definition: auxiliary.h:68
idInit
ideal idInit(int idsize, int rank)
initialise an ideal / module
Definition: simpleideals.cc:34
hHilbEst
static void hHilbEst(scfmon stc, int Nstc, varset var, int Nvar)
Definition: hilb.cc:63
p_SetCoeff
static number p_SetCoeff(poly p, number n, ring r)
Definition: p_polys.h:399
id_TestTail
#define id_TestTail(A, lR, tR)
Definition: simpleideals.h:77
ChoosePVar
static poly ChoosePVar(ideal I)
Definition: hilb.cc:479
SortByDeg_p
static void SortByDeg_p(ideal I, poly p)
!!!!!!!!!!!!!!!!!!!! Just for Monomial Ideals !!!!!!!!!!!!!!!!!!!!!!!!!!!!
Definition: hilb.cc:288
WerrorS
void WerrorS(const char *s)
Definition: feFopen.cc:24
m
int m
Definition: cfEzgcd.cc:121
assume
#define assume(x)
Definition: mod2.h:384
NULL
#define NULL
Definition: omList.c:11
hGetpure
scmon hGetpure(scmon p)
Definition: hutil.cc:1054
ideals.h
p_SetComp
static unsigned long p_SetComp(poly p, unsigned long c, ring r)
Definition: p_polys.h:237
l
int l
Definition: cfEzgcd.cc:93
errorreported
VAR short errorreported
Definition: feFopen.cc:23
R
#define R
Definition: sirandom.c:27
id_Mult
ideal id_Mult(ideal h1, ideal h2, const ring R)
h1 * h2 one h_i must be an ideal (with at least one column) the other h_i may be a module (with no co...
Definition: simpleideals.cc:735
intvec::rows
int rows() const
Definition: intvec.h:96
p_Setm
static void p_Setm(poly p, const ring r)
Definition: p_polys.h:223
positionInOrbit_IG_Case
static int positionInOrbit_IG_Case(ideal I, poly w, std::vector< ideal > idorb, std::vector< poly > polist, int trInd, int trunDegHs)
Definition: hilb.cc:1569
hwork
VAR scfmon hwork
Definition: hutil.cc:15
pLmDivisibleBy
#define pLmDivisibleBy(a, b)
like pDivisibleBy, except that it is assumed that a!=NULL, b!=NULL
Definition: polys.h:134
p_Totaldegree
static long p_Totaldegree(poly p, const ring r)
Definition: p_polys.h:1428
p
int p
Definition: cfModGcd.cc:4019
currRing
VAR ring currRing
Widely used global variable which specifies the current polynomial ring for Singular interpreter and ...
Definition: polys.cc:13
s
const CanonicalForm int s
Definition: facAbsFact.cc:55
count
int status int void size_t count
Definition: si_signals.h:58
pCopy
#define pCopy(p)
return a copy of the poly
Definition: polys.h:174
hDegreeSeries
void hDegreeSeries(intvec *s1, intvec *s2, int *co, int *mu)
Definition: hilb.cc:1380
ipshell.h
IDELEMS
#define IDELEMS(i)
Definition: simpleideals.h:23
p_ISet
poly p_ISet(long i, const ring r)
returns the poly representing the integer i
Definition: p_polys.cc:1284
id_Copy
ideal id_Copy(ideal h1, const ring r)
copy an ideal
Definition: simpleideals.cc:412
hNpure
VAR int hNpure
Definition: hutil.cc:18
pGetCoeff
static number & pGetCoeff(poly p)
return an alias to the leading coefficient of p assumes that p != NULL NOTE: not copy
Definition: monomials.h:41
PrintLn
void PrintLn()
Definition: reporter.cc:309
Q0
STATIC_VAR int * Q0
Definition: hilb.cc:45
intvec::length
int length() const
Definition: intvec.h:94
scPrintDegree
void scPrintDegree(int co, int mu)
Definition: hdegree.cc:807
id_Head
ideal id_Head(ideal h, const ring r)
returns the ideals of initial terms
Definition: simpleideals.cc:1183
idQuotMon
ideal idQuotMon(ideal Iorig, ideal p)
Definition: hilb.cc:409
mp_Sub
matrix mp_Sub(matrix a, matrix b, const ring R)
Definition: matpol.cc:195
numbers.h
hAddHilb
static int * hAddHilb(int Nv, int x, int *pol, int *lp)
Definition: hilb.cc:104
omAlloc0
#define omAlloc0(size)
Definition: omAllocDecl.h:209
luDecomp
void luDecomp(const matrix aMat, matrix &pMat, matrix &lMat, matrix &uMat, const ring R)
LU-decomposition of a given (m x n)-matrix.
Definition: linearAlgebra.cc:103
slicehilb
void slicehilb(ideal I)
Definition: hilb.cc:1130
MAX_INT_VAL
const int MAX_INT_VAL
Definition: mylimits.h:11
hLexS
void hLexS(scfmon stc, int Nstc, varset var, int Nvar)
Definition: hutil.cc:508
linearAlgebra.h
coeffs.h
hCreate
monf hCreate(int Nvar)
Definition: hutil.cc:998
pWrite
void pWrite(poly p)
Definition: polys.h:292