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