My Project
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,intvec *modul_weight)
1399 {
1400  int i, j, l, k;
1401  if (hseries == NULL)
1402  return;
1403  l = hseries->length()-1;
1404  k = (*hseries)[l];
1405  if ((modul_weight!=NULL)&&(modul_weight->compare(0)!=0))
1406  {
1407  char *s=modul_weight->ivString(1,0,1);
1408  Print("module weights:%s\n",s);
1409  omFree(s);
1410  }
1411  for (i = 0; i < l; i++)
1412  {
1413  j = (*hseries)[i];
1414  if (j != 0)
1415  {
1416  Print("// %8d t^%d\n", j, i+k);
1417  }
1418  }
1419 }
1420 
1421 /*
1422 *caller
1423 */
1424 void hLookSeries(ideal S, intvec *modulweight, ideal Q, intvec *wdegree, ring tailRing)
1425 {
1426  id_TestTail(S, currRing, tailRing);
1427 
1428  intvec *hseries1 = hFirstSeries(S, modulweight, Q, wdegree, tailRing);
1429  if (errorreported) return;
1430 
1431  hPrintHilb(hseries1,modulweight);
1432 
1433  const int l = hseries1->length()-1;
1434 
1435  intvec *hseries2 = (l > 1) ? hSecondSeries(hseries1) : hseries1;
1436 
1437  int co, mu;
1438  hDegreeSeries(hseries1, hseries2, &co, &mu);
1439 
1440  PrintLn();
1441  hPrintHilb(hseries2,modulweight);
1442  if ((l == 1) &&(mu == 0))
1443  scPrintDegree(rVar(currRing)+1, 0);
1444  else
1445  scPrintDegree(co, mu);
1446  if (l>1)
1447  delete hseries1;
1448  delete hseries2;
1449 }
1450 
1451 /***********************************************************************
1452  Computation of Hilbert series of non-commutative monomial algebras
1453 ************************************************************************/
1454 
1455 
1456 static void idInsertMonomial(ideal I, poly p)
1457 {
1458  /*
1459  * It adds monomial in I and if required,
1460  * enlarge the size of poly-set by 16.
1461  * It does not make copy of p.
1462  */
1463 
1464  if(I == NULL)
1465  {
1466  return;
1467  }
1468 
1469  int j = IDELEMS(I) - 1;
1470  while ((j >= 0) && (I->m[j] == NULL))
1471  {
1472  j--;
1473  }
1474  j++;
1475  if (j == IDELEMS(I))
1476  {
1477  pEnlargeSet(&(I->m), IDELEMS(I), 16);
1478  IDELEMS(I) +=16;
1479  }
1480  I->m[j] = p;
1481 }
1482 
1483 static int comapreMonoIdBases(ideal J, ideal Ob)
1484 {
1485  /*
1486  * Monomials of J and Ob are assumed to
1487  * be already sorted. J and Ob are
1488  * represented by the minimal generating set.
1489  */
1490  int i, s;
1491  s = 1;
1492  int JCount = IDELEMS(J);
1493  int ObCount = IDELEMS(Ob);
1494 
1495  if(idIs0(J))
1496  {
1497  return(1);
1498  }
1499  if(JCount != ObCount)
1500  {
1501  return(0);
1502  }
1503 
1504  for(i = 0; i < JCount; i++)
1505  {
1506  if(!(p_LmEqual(J->m[i], Ob->m[i], currRing)))
1507  {
1508  return(0);
1509  }
1510  }
1511  return(s);
1512 }
1513 
1514 static int CountOnIdUptoTruncationIndex(ideal I, int tr)
1515 {
1516  /*
1517  * The ideal I must be sorted in increasing total degree.
1518  * It counts the number of monomials in I up to
1519  * degree less than or equal to tr.
1520  */
1521 
1522  //case when I=1;
1523  if(p_Totaldegree(I->m[0], currRing) == 0)
1524  {
1525  return(1);
1526  }
1527 
1528  int count = 0;
1529  for(int i = 0; i < IDELEMS(I); i++)
1530  {
1531  if(p_Totaldegree(I->m[i], currRing) > tr)
1532  {
1533  return (count);
1534  }
1535  count = count + 1;
1536  }
1537 
1538  return(count);
1539 }
1540 
1541 static int comapreMonoIdBases_IG_Case(ideal J, int JCount, ideal Ob, int ObCount)
1542 {
1543  /*
1544  * Monomials of J and Ob are assumed to
1545  * be already sorted in increasing degrees.
1546  * J and Ob are represented by the minimal
1547  * generating set. It checks if J and Ob have
1548  * same monomials up to deg <=tr.
1549  */
1550 
1551  int i, s;
1552  s = 1;
1553  //when J is null
1554  //
1555  if(JCount != ObCount)
1556  {
1557  return(0);
1558  }
1559 
1560  if(JCount == 0)
1561  {
1562  return(1);
1563  }
1564 
1565  for(i = 0; i< JCount; i++)
1566  {
1567  if(!(p_LmEqual(J->m[i], Ob->m[i], currRing)))
1568  {
1569  return(0);
1570  }
1571  }
1572 
1573  return(s);
1574 }
1575 
1576 static int positionInOrbit_IG_Case(ideal I, poly w, std::vector<ideal> idorb, std::vector<poly> polist, int trInd, int trunDegHs)
1577 {
1578  /*
1579  * It compares the ideal I with ideals in the set 'idorb'
1580  * up to total degree =
1581  * trInd - max(deg of w, deg of word in polist) polynomials.
1582  *
1583  * It returns 0 if I is not equal to any ideal in the
1584  * 'idorb' else returns position of the matched ideal.
1585  */
1586 
1587  int ps = 0;
1588  int i, s = 0;
1589  int orbCount = idorb.size();
1590 
1591  if(idIs0(I))
1592  {
1593  return(1);
1594  }
1595 
1596  int degw = p_Totaldegree(w, currRing);
1597  int degp;
1598  int dtr;
1599  int dtrp;
1600 
1601  dtr = trInd - degw;
1602  int IwCount;
1603 
1604  IwCount = CountOnIdUptoTruncationIndex(I, dtr);
1605 
1606  if(IwCount == 0)
1607  {
1608  return(1);
1609  }
1610 
1611  int ObCount;
1612 
1613  bool flag2 = FALSE;
1614 
1615  for(i = 1;i < orbCount; i++)
1616  {
1617  degp = p_Totaldegree(polist[i], currRing);
1618  if(degw > degp)
1619  {
1620  dtr = trInd - degw;
1621 
1622  ObCount = 0;
1623  ObCount = CountOnIdUptoTruncationIndex(idorb[i], dtr);
1624  if(ObCount == 0)
1625  {continue;}
1626  if(flag2)
1627  {
1628  IwCount = 0;
1629  IwCount = CountOnIdUptoTruncationIndex(I, dtr);
1630  flag2 = FALSE;
1631  }
1632  }
1633  else
1634  {
1635  flag2 = TRUE;
1636  dtrp = trInd - degp;
1637  ObCount = 0;
1638  ObCount = CountOnIdUptoTruncationIndex(idorb[i], dtrp);
1639  IwCount = 0;
1640  IwCount = CountOnIdUptoTruncationIndex(I, dtrp);
1641  }
1642 
1643  s = comapreMonoIdBases_IG_Case(I, IwCount, idorb[i], ObCount);
1644 
1645  if(s)
1646  {
1647  ps = i + 1;
1648  break;
1649  }
1650  }
1651  return(ps);
1652 }
1653 
1654 static int positionInOrbit_FG_Case(ideal I, poly, std::vector<ideal> idorb, std::vector<poly>, int, int)
1655 {
1656  /*
1657  * It compares the ideal I with ideals in the set 'idorb'.
1658  * I and ideals of 'idorb' are sorted.
1659  *
1660  * It returns 0 if I is not equal to any ideal of 'idorb'
1661  * else returns position of the matched ideal.
1662  */
1663  int ps = 0;
1664  int i, s = 0;
1665  int OrbCount = idorb.size();
1666 
1667  if(idIs0(I))
1668  {
1669  return(1);
1670  }
1671 
1672  for(i = 1; i < OrbCount; i++)
1673  {
1674  s = comapreMonoIdBases(I, idorb[i]);
1675  if(s)
1676  {
1677  ps = i + 1;
1678  break;
1679  }
1680  }
1681 
1682  return(ps);
1683 }
1684 
1685 static int positionInOrbitTruncationCase(ideal I, poly w, std::vector<ideal> idorb, std::vector<poly> polist, int , int trunDegHs)
1686 {
1687  /*
1688  * It compares the ideal I with ideals in the set 'idorb'.
1689  * I and ideals in 'idorb' are sorted.
1690 
1691  * returns 0 if I is not equal to any ideal of 'idorb'
1692  * else returns position of the matched ideal.
1693  */
1694 
1695  int ps = 0;
1696  int i, s = 0;
1697  int OrbCount = idorb.size();
1698  int dtr=0; int IwCount, ObCount;
1699  dtr = trunDegHs - 1 - p_Totaldegree(w, currRing);
1700 
1701  if(idIs0(I))
1702  {
1703  for(i = 1; i < OrbCount; i++)
1704  {
1705  if(p_Totaldegree(w, currRing) == p_Totaldegree(polist[i], currRing))
1706  {
1707  if(idIs0(idorb[i]))
1708  return(i+1);
1709  ObCount=0;
1710  ObCount = CountOnIdUptoTruncationIndex(idorb[i], dtr);
1711  if(ObCount==0)
1712  {
1713  ps = i + 1;
1714  break;
1715  }
1716  }
1717  }
1718 
1719  return(ps);
1720  }
1721 
1722  IwCount = CountOnIdUptoTruncationIndex(I, dtr);
1723 
1724  if(p_Totaldegree(I->m[0], currRing)==0)
1725  {
1726  for(i = 1; i < OrbCount; i++)
1727  {
1728  if(idIs0(idorb[i]))
1729  continue;
1730  if(p_Totaldegree(idorb[i]->m[0], currRing)==0)
1731  {
1732  ps = i + 1;
1733  break;
1734  }
1735  }
1736  return(ps);
1737  }
1738 
1739  for(i = 1; i < OrbCount; i++)
1740  {
1741  if(p_Totaldegree(w, currRing) == p_Totaldegree(polist[i], currRing))
1742  {
1743  if(idIs0(idorb[i]))
1744  continue;
1745  ObCount=0;
1746  ObCount = CountOnIdUptoTruncationIndex(idorb[i], dtr);
1747  s = comapreMonoIdBases_IG_Case(I, IwCount, idorb[i], ObCount);
1748  if(s)
1749  {
1750  ps = i + 1;
1751  break;
1752  }
1753  }
1754  }
1755 
1756  return(ps);
1757 }
1758 
1759 static int monCompare( const void *m, const void *n)
1760 {
1761  /* compares monomials */
1762 
1763  return(p_Compare(*(poly*) m, *(poly*)n, currRing));
1764 }
1765 
1767 {
1768  /*
1769  * sorts monomial ideal in ascending order
1770  * order must be a total degree
1771  */
1772 
1773  qsort(I->m, IDELEMS(I), sizeof(poly), monCompare);
1774 
1775 }
1776 
1777 
1778 
1779 static ideal minimalMonomialGenSet(ideal I)
1780 {
1781  /*
1782  * eliminates monomials which
1783  * can be generated by others in I
1784  */
1785  //first sort monomials of the ideal
1786 
1787  idSkipZeroes(I);
1788 
1790 
1791  int i, k;
1792  int ICount = IDELEMS(I);
1793 
1794  for(k = ICount - 1; k >=1; k--)
1795  {
1796  for(i = 0; i < k; i++)
1797  {
1798 
1799  if(p_LmDivisibleBy(I->m[i], I->m[k], currRing))
1800  {
1801  pDelete(&(I->m[k]));
1802  break;
1803  }
1804  }
1805  }
1806 
1807  idSkipZeroes(I);
1808  return(I);
1809 }
1810 
1811 static poly shiftInMon(poly p, int i, int lV, const ring r)
1812 {
1813  /*
1814  * shifts the varibles of monomial p in the i^th layer,
1815  * p remains unchanged,
1816  * creates new poly and returns it for the colon ideal
1817  */
1818  poly smon = p_One(r);
1819  int j, sh, cnt;
1820  cnt = r->N;
1821  sh = i*lV;
1822  int *e=(int *)omAlloc((r->N+1)*sizeof(int));
1823  int *s=(int *)omAlloc0((r->N+1)*sizeof(int));
1824  p_GetExpV(p, e, r);
1825 
1826  for(j = 1; j <= cnt; j++)
1827  {
1828  if(e[j] == 1)
1829  {
1830  s[j+sh] = e[j];
1831  }
1832  }
1833 
1834  p_SetExpV(smon, s, currRing);
1835  omFree(e);
1836  omFree(s);
1837 
1839  p_Setm(smon, currRing);
1840 
1841  return(smon);
1842 }
1843 
1844 static poly deleteInMon(poly w, int i, int lV, const ring r)
1845 {
1846  /*
1847  * deletes the variables upto i^th layer of monomial w
1848  * w remains unchanged
1849  * creates new poly and returns it for the colon ideal
1850  */
1851 
1852  poly dw = p_One(currRing);
1853  int *e = (int *)omAlloc((r->N+1)*sizeof(int));
1854  int *s=(int *)omAlloc0((r->N+1)*sizeof(int));
1855  p_GetExpV(w, e, r);
1856  int j, cnt;
1857  cnt = i*lV;
1858  /*
1859  for(j=1;j<=cnt;j++)
1860  {
1861  e[j]=0;
1862  }*/
1863  for(j = (cnt+1); j < (r->N+1); j++)
1864  {
1865  s[j] = e[j];
1866  }
1867 
1868  p_SetExpV(dw, s, currRing);//new exponents
1869  omFree(e);
1870  omFree(s);
1871 
1873  p_Setm(dw, currRing);
1874 
1875  return(dw);
1876 }
1877 
1878 static void TwordMap(poly p, poly w, int lV, int d, ideal Jwi, bool &flag)
1879 {
1880  /*
1881  * computes T_w(p) in a new poly object and places it
1882  * in Jwi which stores elements of colon ideal of I,
1883  * p and w remain unchanged,
1884  * the new polys for Jwi are constructed by sub-routines
1885  * deleteInMon, shiftInMon, p_MDivide,
1886  * places the result in Jwi and deletes the new polys
1887  * coming in dw, smon, qmon
1888  */
1889  int i;
1890  poly smon, dw;
1891  poly qmonp = NULL;
1892  bool del;
1893 
1894  for(i = 0;i <= d - 1; i++)
1895  {
1896  dw = deleteInMon(w, i, lV, currRing);
1897  smon = shiftInMon(p, i, lV, currRing);
1898  del = TRUE;
1899 
1900  if(pLmDivisibleBy(smon, w))
1901  {
1902  flag = TRUE;
1903  del = FALSE;
1904 
1905  pDelete(&dw);
1906  pDelete(&smon);
1907 
1908  //delete all monomials of Jwi
1909  //and make Jwi =1
1910 
1911  for(int j = 0;j < IDELEMS(Jwi); j++)
1912  {
1913  pDelete(&Jwi->m[j]);
1914  }
1915 
1917  break;
1918  }
1919 
1920  if(pLmDivisibleBy(dw, smon))
1921  {
1922  del = FALSE;
1923  qmonp = p_MDivide(smon, dw, currRing);
1924  idInsertMonomial(Jwi, shiftInMon(qmonp, -d, lV, currRing));
1925  pLmFree(&qmonp);
1926  pDelete(&dw);
1927  pDelete(&smon);
1928  }
1929  //in case both if are false, delete dw and smon
1930  if(del)
1931  {
1932  pDelete(&dw);
1933  pDelete(&smon);
1934  }
1935  }
1936 
1937 }
1938 
1939 static ideal colonIdeal(ideal S, poly w, int lV, ideal Jwi, int trunDegHs)
1940 {
1941  /*
1942  * It computes the right colon ideal of a two-sided ideal S
1943  * w.r.t. word w and save it in a new object Jwi.
1944  * It keeps S and w unchanged.
1945  */
1946 
1947  if(idIs0(S))
1948  {
1949  return(S);
1950  }
1951 
1952  int i, d;
1953  d = p_Totaldegree(w, currRing);
1954  if(trunDegHs !=0 && d >= trunDegHs)
1955  {
1957  return(Jwi);
1958  }
1959  bool flag = FALSE;
1960  int SCount = IDELEMS(S);
1961  for(i = 0; i < SCount; i++)
1962  {
1963  TwordMap(S->m[i], w, lV, d, Jwi, flag);
1964  if(flag)
1965  {
1966  break;
1967  }
1968  }
1969 
1970  Jwi = minimalMonomialGenSet(Jwi);
1971  return(Jwi);
1972 }
1973 
1974 void HilbertSeries_OrbitData(ideal S, int lV, bool IG_CASE, bool mgrad, bool odp, int trunDegHs)
1975 {
1976 
1977  /* new story:
1978  no lV is needed, i.e. it is to be determined
1979  the rest is extracted from the interface input list in extra.cc and makes the input of this proc
1980  called from extra.cc
1981  */
1982 
1983  /*
1984  * This is based on iterative right colon operations on a
1985  * two-sided monomial ideal of the free associative algebra.
1986  * The algorithm terminates for those monomial ideals
1987  * whose monomials define "regular formal languages",
1988  * that is, all monomials of the input ideal can be obtained
1989  * from finite languages by applying finite number of
1990  * rational operations.
1991  */
1992 
1993  int trInd;
1994  S = minimalMonomialGenSet(S);
1995  if( !idIs0(S) && p_Totaldegree(S->m[0], currRing)==0)
1996  {
1997  PrintS("Hilbert Series:\n 0\n");
1998  return;
1999  }
2000  int (*POS)(ideal, poly, std::vector<ideal>, std::vector<poly>, int, int);
2001  if(trunDegHs != 0)
2002  {
2003  Print("\nTruncation degree = %d\n",trunDegHs);
2005  }
2006  else
2007  {
2008  if(IG_CASE)
2009  {
2010  if(idIs0(S))
2011  {
2012  WerrorS("wrong input: it is not an infinitely gen. case");
2013  return;
2014  }
2015  trInd = p_Totaldegree(S->m[IDELEMS(S)-1], currRing);
2016  POS = &positionInOrbit_IG_Case;
2017  }
2018  else
2019  POS = &positionInOrbit_FG_Case;
2020  }
2021  std::vector<ideal > idorb;
2022  std::vector< poly > polist;
2023 
2024  ideal orb_init = idInit(1, 1);
2025  idorb.push_back(orb_init);
2026 
2027  polist.push_back( p_One(currRing));
2028 
2029  std::vector< std::vector<int> > posMat;
2030  std::vector<int> posRow(lV,0);
2031  std::vector<int> C;
2032 
2033  int ds, is, ps;
2034  unsigned long lpcnt = 0;
2035 
2036  poly w, wi;
2037  ideal Jwi;
2038 
2039  while(lpcnt < idorb.size())
2040  {
2041  w = NULL;
2042  w = polist[lpcnt];
2043  if(lpcnt >= 1 && idIs0(idorb[lpcnt]) == FALSE)
2044  {
2045  if(p_Totaldegree(idorb[lpcnt]->m[0], currRing) != 0)
2046  {
2047  C.push_back(1);
2048  }
2049  else
2050  C.push_back(0);
2051  }
2052  else
2053  {
2054  C.push_back(1);
2055  }
2056 
2057  ds = p_Totaldegree(w, currRing);
2058  lpcnt++;
2059 
2060  for(is = 1; is <= lV; is++)
2061  {
2062  wi = NULL;
2063  //make new copy 'wi' of word w=polist[lpcnt]
2064  //and update it (for the colon operation).
2065  //if corresponding to wi, right colon operation gives
2066  //a new (right colon) ideal of S,
2067  //keep 'wi' in the polist else delete it
2068 
2069  wi = pCopy(w);
2070  p_SetExp(wi, (ds*lV)+is, 1, currRing);
2071  p_Setm(wi, currRing);
2072  Jwi = NULL;
2073  //Jwi stores (right) colon ideal of S w.r.t. word
2074  //wi if colon operation gives a new ideal place it
2075  //in the vector of ideals 'idorb'
2076  //otherwise delete it
2077 
2078  Jwi = idInit(1,1);
2079 
2080  Jwi = colonIdeal(S, wi, lV, Jwi, trunDegHs);
2081  ps = (*POS)(Jwi, wi, idorb, polist, trInd, trunDegHs);
2082 
2083  if(ps == 0) // finds a new ideal
2084  {
2085  posRow[is-1] = idorb.size();
2086 
2087  idorb.push_back(Jwi);
2088  polist.push_back(wi);
2089  }
2090  else // ideal is already there in the set
2091  {
2092  posRow[is-1]=ps-1;
2093  idDelete(&Jwi);
2094  pDelete(&wi);
2095  }
2096  }
2097  posMat.push_back(posRow);
2098  posRow.resize(lV,0);
2099  }
2100  int lO = C.size();//size of the orbit
2101  PrintLn();
2102  Print("maximal length of words = %ld\n", p_Totaldegree(polist[lO-1], currRing));
2103  Print("\nlength of the Orbit = %d", lO);
2104  PrintLn();
2105 
2106  if(odp)
2107  {
2108  Print("words description of the Orbit: \n");
2109  for(is = 0; is < lO; is++)
2110  {
2111  pWrite0(polist[is]);
2112  PrintS(" ");
2113  }
2114  PrintLn();
2115  PrintS("\nmaximal degree, #(sum_j R(w,w_j))");
2116  PrintLn();
2117  for(is = 0; is < lO; is++)
2118  {
2119  if(idIs0(idorb[is]))
2120  {
2121  PrintS("NULL\n");
2122  }
2123  else
2124  {
2125  Print("%ld, %d \n",p_Totaldegree(idorb[is]->m[IDELEMS(idorb[is])-1], currRing),IDELEMS(idorb[is]));
2126  }
2127  }
2128  }
2129 
2130  for(is = idorb.size()-1; is >= 0; is--)
2131  {
2132  idDelete(&idorb[is]);
2133  }
2134  for(is = polist.size()-1; is >= 0; is--)
2135  {
2136  pDelete(&polist[is]);
2137  }
2138 
2139  idorb.resize(0);
2140  polist.resize(0);
2141 
2142  int adjMatrix[lO][lO];
2143  memset(adjMatrix, 0, lO*lO*sizeof(int));
2144  int rowCount, colCount;
2145  int tm = 0;
2146  if(!mgrad)
2147  {
2148  for(rowCount = 0; rowCount < lO; rowCount++)
2149  {
2150  for(colCount = 0; colCount < lV; colCount++)
2151  {
2152  tm = posMat[rowCount][colCount];
2153  adjMatrix[rowCount][tm] = adjMatrix[rowCount][tm] + 1;
2154  }
2155  }
2156  }
2157 
2158  ring r = currRing;
2159  int npar;
2160  char** tt;
2161  TransExtInfo p;
2162  if(!mgrad)
2163  {
2164  tt=(char**)omAlloc(sizeof(char*));
2165  tt[0] = omStrDup("t");
2166  npar = 1;
2167  }
2168  else
2169  {
2170  tt=(char**)omalloc(lV*sizeof(char*));
2171  for(is = 0; is < lV; is++)
2172  {
2173  tt[is] = (char*)omAlloc(7*sizeof(char)); //if required enlarge it later
2174  sprintf (tt[is], "t%d", is+1);
2175  }
2176  npar = lV;
2177  }
2178 
2179  p.r = rDefault(0, npar, tt);
2181  char** xx = (char**)omAlloc(sizeof(char*));
2182  xx[0] = omStrDup("x");
2183  ring R = rDefault(cf, 1, xx);
2184  rChangeCurrRing(R);//rWrite(R);
2185  /*
2186  * matrix corresponding to the orbit of the ideal
2187  */
2188  matrix mR = mpNew(lO, lO);
2189  matrix cMat = mpNew(lO,1);
2190  poly rc;
2191 
2192  if(!mgrad)
2193  {
2194  for(rowCount = 0; rowCount < lO; rowCount++)
2195  {
2196  for(colCount = 0; colCount < lO; colCount++)
2197  {
2198  if(adjMatrix[rowCount][colCount] != 0)
2199  {
2200  MATELEM(mR, rowCount + 1, colCount + 1) = p_ISet(adjMatrix[rowCount][colCount], R);
2201  p_SetCoeff(MATELEM(mR, rowCount + 1, colCount + 1), n_Mult(pGetCoeff(mR->m[lO*rowCount+colCount]),n_Param(1, R->cf), R->cf), R);
2202  }
2203  }
2204  }
2205  }
2206  else
2207  {
2208  for(rowCount = 0; rowCount < lO; rowCount++)
2209  {
2210  for(colCount = 0; colCount < lV; colCount++)
2211  {
2212  rc=NULL;
2213  rc=p_One(R);
2214  p_SetCoeff(rc, n_Mult(pGetCoeff(rc), n_Param(colCount+1, R->cf),R->cf), R);
2215  MATELEM(mR, rowCount +1, posMat[rowCount][colCount]+1)=p_Add_q(rc,MATELEM(mR, rowCount +1, posMat[rowCount][colCount]+1), R);
2216  }
2217  }
2218  }
2219 
2220  for(rowCount = 0; rowCount < lO; rowCount++)
2221  {
2222  if(C[rowCount] != 0)
2223  {
2224  MATELEM(cMat, rowCount + 1, 1) = p_ISet(C[rowCount], R);
2225  }
2226  }
2227 
2228  matrix u;
2229  unitMatrix(lO, u); //unit matrix
2230  matrix gMat = mp_Sub(u, mR, R);
2231 
2232  char* s;
2233 
2234  if(odp)
2235  {
2236  PrintS("\nlinear system:\n");
2237  if(!mgrad)
2238  {
2239  for(rowCount = 0; rowCount < lO; rowCount++)
2240  {
2241  Print("H(%d) = ", rowCount+1);
2242  for(colCount = 0; colCount < lV; colCount++)
2243  {
2244  StringSetS(""); nWrite(n_Param(1, R->cf));
2245  s = StringEndS(); PrintS(s);
2246  Print("*"); omFree(s);
2247  Print("H(%d) + ", posMat[rowCount][colCount] + 1);
2248  }
2249  Print(" %d\n", C[rowCount] );
2250  }
2251  PrintS("where H(1) represents the series corresp. to input ideal\n");
2252  PrintS("and i^th summand in the rhs of an eqn. is according\n");
2253  PrintS("to the right colon map corresp. to the i^th variable\n");
2254  }
2255  else
2256  {
2257  for(rowCount = 0; rowCount < lO; rowCount++)
2258  {
2259  Print("H(%d) = ", rowCount+1);
2260  for(colCount = 0; colCount < lV; colCount++)
2261  {
2262  StringSetS(""); nWrite(n_Param(colCount+1, R->cf));
2263  s = StringEndS(); PrintS(s);
2264  Print("*");omFree(s);
2265  Print("H(%d) + ", posMat[rowCount][colCount] + 1);
2266  }
2267  Print(" %d\n", C[rowCount] );
2268  }
2269  PrintS("where H(1) represents the series corresp. to input ideal\n");
2270  }
2271  }
2272  PrintLn();
2273  posMat.resize(0);
2274  C.resize(0);
2275  matrix pMat;
2276  matrix lMat;
2277  matrix uMat;
2278  matrix H_serVec = mpNew(lO, 1);
2279  matrix Hnot;
2280 
2281  //std::clock_t start;
2282  //start = std::clock();
2283 
2284  luDecomp(gMat, pMat, lMat, uMat, R);
2285  luSolveViaLUDecomp(pMat, lMat, uMat, cMat, H_serVec, Hnot);
2286 
2287  //to print system solving time
2288  //if(odp){
2289  //std::cout<<"solving time of the system = "<<(std::clock()-start)/(double)(CLOCKS_PER_SEC / 1000)<<" ms"<<std::endl;}
2290 
2291  mp_Delete(&mR, R);
2292  mp_Delete(&u, R);
2293  mp_Delete(&pMat, R);
2294  mp_Delete(&lMat, R);
2295  mp_Delete(&uMat, R);
2296  mp_Delete(&cMat, R);
2297  mp_Delete(&gMat, R);
2298  mp_Delete(&Hnot, R);
2299  //print the Hilbert series and length of the Orbit
2300  PrintLn();
2301  Print("Hilbert series:");
2302  PrintLn();
2303  pWrite(H_serVec->m[0]);
2304  if(!mgrad)
2305  {
2306  omFree(tt[0]);
2307  }
2308  else
2309  {
2310  for(is = lV-1; is >= 0; is--)
2311 
2312  omFree( tt[is]);
2313  }
2314  omFree(tt);
2315  omFree(xx[0]);
2316  omFree(xx);
2317  rChangeCurrRing(r);
2318  rKill(R);
2319 }
2320 
2321 ideal RightColonOperation(ideal S, poly w, int lV)
2322 {
2323  /*
2324  * This returns right colon ideal of a monomial two-sided ideal of
2325  * the free associative algebra with respect to a monomial 'w'
2326  * (S:_R w).
2327  */
2328  S = minimalMonomialGenSet(S);
2329  ideal Iw = idInit(1,1);
2330  Iw = colonIdeal(S, w, lV, Iw, 0);
2331  return (Iw);
2332 }
2333 
long int64
Definition: auxiliary.h:68
#define TRUE
Definition: auxiliary.h:100
#define FALSE
Definition: auxiliary.h:96
void * ADDRESS
Definition: auxiliary.h:119
const CanonicalForm CFMap CFMap & N
Definition: cfEzgcd.cc:56
int l
Definition: cfEzgcd.cc:100
int m
Definition: cfEzgcd.cc:128
int i
Definition: cfEzgcd.cc:132
int k
Definition: cfEzgcd.cc:99
Variable x
Definition: cfModGcd.cc:4084
int p
Definition: cfModGcd.cc:4080
CanonicalForm cf
Definition: cfModGcd.cc:4085
CanonicalForm b
Definition: cfModGcd.cc:4105
void mu(int **points, int sizePoints)
Definition: intvec.h:23
int length() const
Definition: intvec.h:94
int compare(const intvec *o) const
Definition: intvec.cc:206
char * ivString(int not_mat=1, int spaces=0, int dim=2) const
Definition: intvec.cc:58
int rows() const
Definition: intvec.h:96
poly * m
Definition: matpol.h:18
Coefficient rings, fields and other domains suitable for Singular polynomials.
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:637
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:807
@ n_transExt
used for all transcendental extensions, i.e., the top-most extension in an extension tower is transce...
Definition: coeffs.h:39
coeffs nInitChar(n_coeffType t, void *parameter)
one-time initialisations for new coeffs in case of an error return NULL
Definition: numbers.cc:353
#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
int j
Definition: facHensel.cc:110
void FACTORY_PUBLIC prune(Variable &alpha)
Definition: variable.cc:261
VAR short errorreported
Definition: feFopen.cc:23
void WerrorS(const char *s)
Definition: feFopen.cc:24
#define STATIC_VAR
Definition: globaldefs.h:7
void scPrintDegree(int co, int mu)
Definition: hdegree.cc:881
static void idInsertMonomial(ideal I, poly p)
Definition: hilb.cc:1456
static int comapreMonoIdBases_IG_Case(ideal J, int JCount, ideal Ob, int ObCount)
Definition: hilb.cc:1541
STATIC_VAR int * Q0
Definition: hilb.cc:45
static poly SqFree(ideal I)
Definition: hilb.cc:880
static void idAddMon(ideal I, ideal p)
Definition: hilb.cc:471
static int comapreMonoIdBases(ideal J, ideal Ob)
Definition: hilb.cc:1483
static void TwordMap(poly p, poly w, int lV, int d, ideal Jwi, bool &flag)
Definition: hilb.cc:1878
static poly ChooseP(ideal I)
Definition: hilb.cc:764
static poly deleteInMon(poly w, int i, int lV, const ring r)
Definition: hilb.cc:1844
STATIC_VAR int hLength
Definition: hilb.cc:46
static int CountOnIdUptoTruncationIndex(ideal I, int tr)
Definition: hilb.cc:1514
static poly ChoosePJL(ideal I)
Definition: hilb.cc:705
static int monCompare(const void *m, const void *n)
Definition: hilb.cc:1759
static void hHilbEst(scfmon stc, int Nstc, varset var, int Nvar)
Definition: hilb.cc:63
static void hLastHilb(scmon pure, int Nv, varset var, int *pol, int lp)
Definition: hilb.cc:137
static void hPrintHilb(intvec *hseries, intvec *modul_weight)
Definition: hilb.cc:1398
static int positionInOrbitTruncationCase(ideal I, poly w, std::vector< ideal > idorb, std::vector< poly > polist, int, int trunDegHs)
Definition: hilb.cc:1685
static poly LCMmon(ideal I)
Definition: hilb.cc:948
void HilbertSeries_OrbitData(ideal S, int lV, bool IG_CASE, bool mgrad, bool odp, int trunDegHs)
Definition: hilb.cc:1974
void hLookSeries(ideal S, intvec *modulweight, ideal Q, intvec *wdegree, ring tailRing)
Definition: hilb.cc:1424
static ideal colonIdeal(ideal S, poly w, int lV, ideal Jwi, int trunDegHs)
Definition: hilb.cc:1939
static int hMinModulweight(intvec *modulweight)
Definition: hilb.cc:49
static poly shiftInMon(poly p, int i, int lV, const ring r)
Definition: hilb.cc:1811
static poly ChoosePVar(ideal I)
Definition: hilb.cc:479
static int positionInOrbit_FG_Case(ideal I, poly, std::vector< ideal > idorb, std::vector< poly >, int, int)
Definition: hilb.cc:1654
static ideal SortByDeg(ideal I)
Definition: hilb.cc:388
static bool IsIn(poly p, ideal I)
Definition: hilb.cc:909
static void eulerchar(ideal I, int variables, mpz_ptr ec)
Definition: hilb.cc:837
ideal RightColonOperation(ideal S, poly w, int lV)
Definition: hilb.cc:2321
static void hWDegree(intvec *wdegree)
Definition: hilb.cc:245
void hDegreeSeries(intvec *s1, intvec *s2, int *co, int *mu)
Definition: hilb.cc:1380
static intvec * hSeries(ideal S, intvec *modulweight, int, intvec *wdegree, ideal Q, ring tailRing)
Definition: hilb.cc:1172
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
static ideal minimalMonomialGenSet(ideal I)
Definition: hilb.cc:1779
intvec * hSecondSeries(intvec *hseries1)
Definition: hilb.cc:1345
static int * hAddHilb(int Nv, int x, int *pol, int *lp)
Definition: hilb.cc:104
intvec * hHstdSeries(ideal S, intvec *modulweight, intvec *wdegree, ideal Q, ring tailRing)
Definition: hilb.cc:1328
STATIC_VAR int * Ql
Definition: hilb.cc:45
void sortMonoIdeal_pCompare(ideal I)
Definition: hilb.cc:1766
ideal idQuotMon(ideal Iorig, ideal p)
Definition: hilb.cc:409
static int positionInOrbit_IG_Case(ideal I, poly w, std::vector< ideal > idorb, std::vector< poly > polist, int trInd, int trunDegHs)
Definition: hilb.cc:1576
static void SortByDeg_p(ideal I, poly p)
!!!!!!!!!!!!!!!!!!!! Just for Monomial Ideals !!!!!!!!!!!!!!!!!!!!!!!!!!!!
Definition: hilb.cc:288
void slicehilb(ideal I)
Definition: hilb.cc:1130
static bool JustVar(ideal I)
Definition: hilb.cc:806
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
STATIC_VAR int ** Qpol
Definition: hilb.cc:44
static void hHilbStep(scmon pure, scfmon stc, int Nstc, varset var, int Nvar, int *pol, int Lpol)
Definition: hilb.cc:177
intvec * hFirstSeries(ideal S, intvec *modulweight, ideal Q, intvec *wdegree, ring tailRing)
Definition: hilb.cc:1335
monf hCreate(int Nvar)
Definition: hutil.cc:999
void hComp(scfmon exist, int Nexist, int ak, scfmon stc, int *Nstc)
Definition: hutil.cc:157
scfmon hInit(ideal S, ideal Q, int *Nexist, ring tailRing)
Definition: hutil.cc:31
void hLex2S(scfmon rad, int e1, int a2, int e2, varset var, int Nvar, scfmon w)
Definition: hutil.cc:815
VAR scfmon hstc
Definition: hutil.cc:16
VAR varset hvar
Definition: hutil.cc:18
void hKill(monf xmem, int Nvar)
Definition: hutil.cc:1013
VAR int hNexist
Definition: hutil.cc:19
void hElimS(scfmon stc, int *e1, int a2, int e2, varset var, int Nvar)
Definition: hutil.cc:675
void hLexS(scfmon stc, int Nstc, varset var, int Nvar)
Definition: hutil.cc:509
void hDelete(scfmon ev, int ev_length)
Definition: hutil.cc:143
VAR 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
VAR scmon hpure
Definition: hutil.cc:17
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 hOrdSupp(scfmon stc, int Nstc, varset var, int Nvar)
Definition: hutil.cc:205
VAR int hNpure
Definition: hutil.cc:19
VAR scfmon hexist
Definition: hutil.cc:16
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
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
BOOLEAN idIs0(ideal h)
returns true if h is the zero ideal
#define idPrint(id)
Definition: ideals.h:46
STATIC_VAR int variables
void rKill(ring r)
Definition: ipshell.cc:6179
STATIC_VAR jList * Q
Definition: janet.cc:30
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.
void luDecomp(const matrix aMat, matrix &pMat, matrix &lMat, matrix &uMat, const ring R)
LU-decomposition of a given (m x n)-matrix.
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...
void mp_Delete(matrix *a, const ring r)
Definition: matpol.cc:880
matrix mp_Sub(matrix a, matrix b, const ring R)
Definition: matpol.cc:196
matrix mpNew(int r, int c)
create a r x c zero-matrix
Definition: matpol.cc:37
#define MATELEM(mat, i, j)
1-based access to matrix
Definition: matpol.h:29
#define assume(x)
Definition: mod2.h:387
#define p_GetComp(p, r)
Definition: monomials.h:64
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
gmp_float exp(const gmp_float &a)
Definition: mpr_complex.cc:357
const int MAX_INT_VAL
Definition: mylimits.h:12
The main handler for Singular numbers which are suitable for Singular polynomials.
#define nWrite(n)
Definition: numbers.h:29
#define omStrDup(s)
Definition: omAllocDecl.h:263
#define omFreeSize(addr, size)
Definition: omAllocDecl.h:260
#define omAlloc(size)
Definition: omAllocDecl.h:210
#define omalloc(size)
Definition: omAllocDecl.h:228
#define omRealloc(addr, size)
Definition: omAllocDecl.h:225
#define omFree(addr)
Definition: omAllocDecl.h:261
#define omAlloc0(size)
Definition: omAllocDecl.h:211
#define NULL
Definition: omList.c:12
poly p_ISet(long i, const ring r)
returns the poly representing the integer i
Definition: p_polys.cc:1292
poly p_MDivide(poly a, poly b, const ring r)
Definition: p_polys.cc:1479
int p_Compare(const poly a, const poly b, const ring R)
Definition: p_polys.cc:4932
poly p_One(const ring r)
Definition: p_polys.cc:1308
void pEnlargeSet(poly **p, int l, int increment)
Definition: p_polys.cc:3766
static poly p_Add_q(poly p, poly q, const ring r)
Definition: p_polys.h:896
static poly p_Head(poly p, const ring r)
copy the i(leading) term of p
Definition: p_polys.h:826
#define p_LmEqual(p1, p2, r)
Definition: p_polys.h:1691
static void p_SetExpV(poly p, int *ev, const ring r)
Definition: p_polys.h:1504
static poly pp_Mult_mm(poly p, poly m, const ring r)
Definition: p_polys.h:991
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:488
static unsigned long p_SetComp(poly p, unsigned long c, ring r)
Definition: p_polys.h:247
static void p_Setm(poly p, const ring r)
Definition: p_polys.h:233
static number p_SetCoeff(poly p, number n, ring r)
Definition: p_polys.h:412
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:469
static BOOLEAN p_IsOne(const poly p, const ring R)
either poly(1) or gen(k)?!
Definition: p_polys.h:1979
static BOOLEAN p_LmDivisibleBy(poly a, poly b, const ring r)
Definition: p_polys.h:1863
static BOOLEAN p_DivisibleBy(poly a, poly b, const ring r)
Definition: p_polys.h:1872
static void p_Delete(poly *p, const ring r)
Definition: p_polys.h:861
static void p_GetExpV(poly p, int *ev, const ring r)
Definition: p_polys.h:1480
static poly p_Copy(poly p, const ring r)
returns a copy of p
Definition: p_polys.h:812
static long p_Totaldegree(poly p, const ring r)
Definition: p_polys.h:1467
void rChangeCurrRing(ring r)
Definition: polys.cc:15
VAR ring currRing
Widely used global variable which specifies the current polynomial ring for Singular interpreter and ...
Definition: polys.cc:13
#define pDelete(p_ptr)
Definition: polys.h:186
#define pLmEqual(p1, p2)
Definition: polys.h:111
void pWrite0(poly p)
Definition: polys.h:309
#define pLmDivisibleBy(a, b)
like pDivisibleBy, except that it is assumed that a!=NULL, b!=NULL
Definition: polys.h:140
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 pCopy(p)
return a copy of the poly
Definition: polys.h:185
#define pOne()
Definition: polys.h:315
void StringSetS(const char *st)
Definition: reporter.cc:128
void PrintS(const char *s)
Definition: reporter.cc:284
char * StringEndS()
Definition: reporter.cc:151
void PrintLn()
Definition: reporter.cc:310
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:102
static short rVar(const ring r)
#define rVar(r) (r->N)
Definition: ring.h:597
int status int void size_t count
Definition: si_signals.h:59
ideal idInit(int idsize, int rank)
initialise an ideal / module
Definition: simpleideals.cc:35
void id_Delete(ideal *h, ring r)
deletes an ideal/module/matrix
ideal id_Head(ideal h, const ring r)
returns the ideals of initial terms
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...
void idSkipZeroes(ideal ide)
gives an ideal/module the minimal possible size
#define IDELEMS(i)
Definition: simpleideals.h:23
#define id_TestTail(A, lR, tR)
Definition: simpleideals.h:77
#define R
Definition: sirandom.c:27
#define loop
Definition: structs.h:80
struct for passing initialization parameters to naInitChar
Definition: transext.h:88