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
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
67static 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
167static 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 }
183 //idDegSortTest(res);
184 return(res);
185}
186
187//idQuot(I,p) for I monomial ideal, p a ideal with a single monomial.
188ideal 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);
226 {
227 p_Delete(&res->m[i],currRing);
228 }
229 else
230 {
231 p_Delete(&I->m[i],currRing);
232 }
233 }
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
250static 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)
258static 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 {
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)
290static 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);
305 flag = FALSE;
306 }
307 }
308 if(!p_IsOne(m, currRing))
309 {
310 return(m);
311 }
312 }
314 m = ChoosePVar(I);
315 return(m);
316}
317
318//chooses 1 \neq p \not\in S. This choice should be made optimal
319static 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)
327static 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)
353static 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
363static 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);
402 }
403}
404
405//tests if an ideal is Square Free, if no, returns the variable which appears at powers >1
406static 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
435static 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
474static 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 }
496 return(m);
497}
498
499//the Roune Slice Algorithm
500static 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 }
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 }
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);
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);
585 break;
586 }
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
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
656void 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);
675 I = Itmp;
676 //printf("\n-------------RouneSlice--------------\n");
677 rouneslice(I,S,q,X->m[0],prune, moreprune, steps, NNN, hilbertcoef, hilbpower);
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
732void 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
750static 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*/
776void hLookSeries(ideal S, intvec *modulweight, ideal Q, intvec *wdegree)
777{
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))
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
808static 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
835static 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
866static 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
893static 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
928static 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
1006static 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
1037static 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
1111static int monCompare( const void *m, const void *n)
1112{
1113 /* compares monomials */
1114
1115 return(p_Compare(*(poly*) m, *(poly*)n, currRing));
1116}
1117
1118static 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
1131static 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
1163static 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
1196static 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
1230static 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
1291static 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
1326void 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);
1369 }
1370 else
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;
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
1673ideal 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"
1702#include "polys/simpleideals.h"
1703#include "coeffs/coeffs.h"
1704
1705#include "kernel/ideals.h"
1706
1707static 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__)
1721static int compare_rp(void *arg,const void *pp1, const void *pp2)
1722#else
1723static 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
1739static 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
1753poly 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
1840static 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
1863intvec* 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
1920static 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
1931static 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
1947intvec* 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/* ------------------------------------------------------------------ */
1982static int hMinModulweight(intvec *modulweight)
1983{
1984 if(modulweight==NULL) return 0;
1985 return modulweight->min_in();
1986}
1987
1988static 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}
2005static 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
2053static 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
2109static 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
2150static 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
2215static 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);
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
2361intvec * 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
intvec * hSecondSeries(intvec *hseries1)
Definition: hilb.cc:697
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
static intvec * hSeries(ideal S, intvec *modulweight, intvec *wdegree, ideal Q)
Definition: hilb.cc:2215
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
intvec * hFirstSeries(ideal A, intvec *module_w, ideal Q, intvec *wdegree)
Definition: hilb.cc:1947
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
intvec * hFirstSeries0(ideal A, ideal Q, intvec *wdegree, const ring src, const ring Qt)
Definition: hilb.cc:1863
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
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 int64 * hAddHilb(int Nv, int x, int64 *pol, int *lp)
Definition: hilb.cc:2005
static ideal minimalMonomialGenSet(ideal I)
Definition: hilb.cc:1131
intvec * hFirstSeries1(ideal S, intvec *modulweight, ideal Q, intvec *wdegree)
Definition: hilb.cc:2361
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
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
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