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