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