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