1 /****************************************
2 * Computer Algebra System SINGULAR *
3 ****************************************/
4 /*
5 * ABSTRACT: resolutions
6 */
7
8 #include "kernel/mod2.h"
9 #include "misc/options.h"
10 #include "kernel/polys.h"
11 #include "kernel/GBEngine/kstd1.h"
12 #include "kernel/GBEngine/kutil.h"
13 #include "kernel/combinatorics/stairc.h"
14 #include "misc/intvec.h"
15 #include "coeffs/numbers.h"
16 #include "kernel/ideals.h"
17 #include "misc/intvec.h"
18 #include "polys/monomials/ring.h"
19 #include "kernel/GBEngine/syz.h"
20 #include "polys/kbuckets.h"
21 #include "polys/prCopy.h"
22
syInitSort(ideal arg,intvec ** modcomp)23 static void syInitSort(ideal arg,intvec **modcomp)
24 {
25 int i,j,k,kk,kkk,jj;
26 idSkipZeroes(arg);
27 polyset F,oldF=arg->m;
28 int Fl=IDELEMS(arg);
29 int rkF=id_RankFreeModule(arg,currRing);
30 int syComponentOrder=currRing->ComponentOrder;
31
32 while ((Fl!=0) && (oldF[Fl-1]==NULL)) Fl--;
33 if (*modcomp!=NULL) delete modcomp;
34 *modcomp = new intvec(rkF+2);
35 F=(polyset)omAlloc0(IDELEMS(arg)*sizeof(poly));
36 j=0;
37 for(i=0;i<=rkF;i++)
38 {
39 k=0;
40 jj = j;
41 (**modcomp)[i] = j;
42 while (k<Fl)
43 {
44 while ((k<Fl) && (pGetComp(oldF[k]) != i)) k++;
45 if (k<Fl)
46 {
47 kk=jj;
48 while ((kk<Fl) && (F[kk]) && (pLmCmp(oldF[k],F[kk])!=syComponentOrder))
49 {
50 kk++;
51 }
52 for (kkk=j;kkk>kk;kkk--)
53 {
54 F[kkk] = F[kkk-1];
55 }
56 F[kk] = oldF[k];
57 //Print("Element %d: ",kk);pWrite(F[kk]);
58 j++;
59 k++;
60 }
61 }
62 }
63 (**modcomp)[rkF+1] = Fl;
64 arg->m = F;
65 omFreeSize((ADDRESS)oldF,IDELEMS(arg)*sizeof(poly));
66 }
67
syCreatePairs(polyset F,int lini,int wend,int k,int j,int i,polyset pairs,int regularPairs=0,ideal mW=NULL)68 static void syCreatePairs(polyset F,int lini,int wend,int k,int j,int i,
69 polyset pairs,int regularPairs=0,ideal mW=NULL)
70 {
71 int l,ii=0,jj;
72 poly p,q;
73
74 while (((k<wend) && (pGetComp(F[k]) == i)) ||
75 ((currRing->qideal!=NULL) && (k<regularPairs+IDELEMS(currRing->qideal))))
76 {
77 p = pOne();
78 if ((k<wend) && (pGetComp(F[k]) == i) && (k!=j))
79 pLcm(F[j],F[k],p);
80 else if (ii<IDELEMS(currRing->qideal))
81 {
82 q = pHead(F[j]);
83 if (mW!=NULL)
84 {
85 for(jj=1;jj<=(currRing->N);jj++)
86 pSetExp(q,jj,pGetExp(q,jj) -pGetExp(mW->m[pGetComp(q)-1],jj));
87 pSetm(q);
88 }
89 pLcm(q,currRing->qideal->m[ii],p);
90 if (mW!=NULL)
91 {
92 for(jj=1;jj<=(currRing->N);jj++)
93 pSetExp(p,jj,pGetExp(p,jj) +pGetExp(mW->m[pGetComp(p)-1],jj));
94 pSetm(p);
95 }
96 pDelete(&q);
97 k = regularPairs+ii;
98 ii++;
99 }
100 l=lini;
101 while ((l<k) && ((pairs[l]==NULL) || (!pDivisibleBy(pairs[l],p))))
102 {
103 if ((pairs[l]!=NULL) && (pDivisibleBy(p,pairs[l])))
104 pDelete(&(pairs[l]));
105 l++;
106 }
107 if (l==k)
108 {
109 pSetm(p);
110 pairs[l] = p;
111 }
112 else
113 pDelete(&p);
114 k++;
115 }
116 }
117
syRedtail2(poly p,polyset redWith,intvec * modcomp)118 static poly syRedtail2(poly p, polyset redWith, intvec *modcomp)
119 {
120 poly h, hn;
121 int hncomp,nxt;
122 int j;
123
124 h = p;
125 hn = pNext(h);
126 while(hn != NULL)
127 {
128 hncomp = pGetComp(hn);
129 j = (*modcomp)[hncomp];
130 nxt = (*modcomp)[hncomp+1];
131 while (j < nxt)
132 {
133 if (pLmDivisibleByNoComp(redWith[j], hn))
134 {
135 //if (TEST_OPT_PROT) PrintS("r");
136 hn = ksOldSpolyRed(redWith[j],hn);
137 if (hn == NULL)
138 {
139 pNext(h) = NULL;
140 return p;
141 }
142 hncomp = pGetComp(hn);
143 j = (*modcomp)[hncomp];
144 nxt = (*modcomp)[hncomp+1];
145 }
146 else
147 {
148 j++;
149 }
150 }
151 h = pNext(h) = hn;
152 hn = pNext(h);
153 }
154 return p;
155 }
156
157 /*2
158 * computes the Schreyer syzygies in the local case
159 * input: arg (only allocated: Shdl, Smax)
160 * output: Shdl, Smax
161 */
sySchreyersSyzygiesFM(ideal arg,intvec ** modcomp)162 static ideal sySchreyersSyzygiesFM(ideal arg,intvec ** modcomp)
163 {
164 int Fl=IDELEMS(arg);
165 while ((Fl!=0) && (arg->m[Fl-1]==NULL)) Fl--;
166 ideal result=idInit(16,arg->rank+Fl);
167 polyset F=arg->m,*Shdl=&(result->m);
168 if (Fl==0) return result;
169
170 int i,j,l,k,totalToRed,ecartToRed,kk;
171 int bestEcart,totalmax,rkF,Sl=0,smax,tmax,tl;
172 int *ecartS, *ecartT, *totalS,
173 *totalT=NULL, *temp=NULL;
174 polyset pairs,S,T,ST/*,oldF*/;
175 poly p,q,toRed;
176 BOOLEAN notFound = FALSE;
177 intvec * newmodcomp = new intvec(Fl+2);
178 intvec * tempcomp;
179
180 //Print("Naechster Modul\n");
181 //idPrint(arg);
182 /*-------------initializing the sets--------------------*/
183 ST=(polyset)omAlloc0(Fl*sizeof(poly));
184 S=(polyset)omAlloc0(Fl*sizeof(poly));
185 ecartS=(int*)omAlloc(Fl*sizeof(int));
186 totalS=(int*)omAlloc(Fl*sizeof(int));
187 T=(polyset)omAlloc0(2*Fl*sizeof(poly));
188 ecartT=(int*)omAlloc(2*Fl*sizeof(int));
189 totalT=(int*)omAlloc(2*Fl*sizeof(int));
190 pairs=(polyset)omAlloc0(Fl*sizeof(poly));
191
192 smax = Fl;
193 tmax = 2*Fl;
194 for (j=1;j<IDELEMS(arg);j++)
195 {
196 if (arg->m[j] != NULL)
197 {
198 assume (arg->m[j-1] != NULL);
199 assume (pGetComp(arg->m[j-1])-pGetComp(arg->m[j])<=0);
200 }
201 }
202 rkF=id_RankFreeModule(arg,currRing);
203 /*----------------construction of the new ordering----------*/
204 if (rkF>0)
205 rSetSyzComp(rkF, currRing);
206 else
207 rSetSyzComp(1, currRing);
208 /*----------------creating S--------------------------------*/
209 for(j=0;j<Fl;j++)
210 {
211 S[j] = pCopy(F[j]);
212 totalS[j] = p_LDeg(S[j],&k,currRing);
213 ecartS[j] = totalS[j]-p_FDeg(S[j],currRing);
214 //Print("%d", pGetComp(S[j]));PrintS(" ");
215 p = S[j];
216 if (rkF==0) pSetCompP(p,1);
217 while (pNext(p)!=NULL) pIter(p);
218 pNext(p) = pHead(F[j]);
219 pIter(p);
220 if (rkF==0)
221 pSetComp(p,j+2);
222 else
223 pSetComp(p,rkF+j+1);
224 pSetmComp(p);
225 }
226 //PrintLn();
227 if (rkF==0) rkF = 1;
228 /*---------------creating the initial for T----------------*/
229 j=0;
230 l=-1;
231 totalmax=-1;
232 for (k=0;k<smax;k++)
233 if (totalS[k]>totalmax) totalmax=totalS[k];
234 for (kk=1;kk<=rkF;kk++)
235 {
236 for (k=0;k<=totalmax;k++)
237 {
238 for (l=0;l<smax;l++)
239 {
240 if ((pGetComp(S[l])==kk) && (totalS[l]==k))
241 {
242 ST[j] = S[l];
243 totalT[j] = totalS[l];
244 ecartT[j] = ecartS[l];
245 //Print("%d", totalS[l]);PrintS(" ");
246 j++;
247 }
248 }
249 }
250 }
251 //PrintLn();
252 for (j=0;j<smax;j++)
253 {
254 totalS[j] = totalT[j];
255 ecartS[j] = ecartT[j];
256 }
257
258 /*---------------computing---------------------------------*/
259 for(j=0;j<smax;j++)
260 {
261 (*newmodcomp)[j+1] = Sl;
262 i = pGetComp(S[j]);
263 int syComponentOrder= currRing->ComponentOrder;
264 int lini,wend;
265 if (syComponentOrder==1)
266 {
267 lini=k=j+1;
268 wend=Fl;
269 }
270 else
271 {
272 lini=k=0;
273 while ((k<j) && (pGetComp(S[k]) != i)) k++;
274 wend=j;
275 }
276 if (TEST_OPT_PROT)
277 {
278 Print("(%d)",Fl-j);
279 mflush();
280 }
281 syCreatePairs(S,lini,wend,k,j,i,pairs);
282 for (k=lini;k<wend;k++)
283 {
284 if (pairs[k]!=NULL)
285 {
286 /*--------------creating T----------------------------------*/
287 for (l=0;l<smax;l++)
288 {
289 ecartT[l] = ecartS[l];
290 totalT[l] = totalS[l];
291 T[l] = ST[l];
292 }
293 tl = smax;
294 tempcomp = ivCopy(*modcomp);
295 /*--------------begin to reduce-----------------------------*/
296 toRed = ksOldCreateSpoly(S[j],S[k]);
297 ecartToRed = 1;
298 bestEcart = 1;
299 if (TEST_OPT_DEBUG)
300 {
301 PrintS("pair: ");pWrite0(S[j]);PrintS(" ");pWrite(S[k]);
302 }
303 if (TEST_OPT_PROT)
304 {
305 PrintS(".");
306 mflush();
307 }
308 //Print("Reduziere Paar %d,%d (ecart %d): \n",j,k,ecartToRed);
309 //Print("Poly %d: ",j);pWrite(S[j]);
310 //Print("Poly %d: ",k);pWrite(S[k]);
311 //Print("Spoly: ");pWrite(toRed);
312 while (pGetComp(toRed)<=rkF)
313 {
314 if (TEST_OPT_DEBUG)
315 {
316 PrintS("toRed: ");pWrite(toRed);
317 }
318 /*
319 * if ((bestEcart) || (ecartToRed!=0))
320 * {
321 */
322 totalToRed = p_LDeg(toRed,&kk,currRing);
323 ecartToRed = totalToRed-p_FDeg(toRed,currRing);
324 /*
325 * }
326 */
327 //Print("toRed now (neuer ecart %d): ",ecartToRed);pWrite(toRed);
328 notFound = TRUE;
329 bestEcart = 32000; //a very large integer
330 p = NULL;
331 int l=0;
332 #define OLD_SEARCH
333 #ifdef OLD_SEARCH
334 while ((l<tl) && (pGetComp(T[l])<pGetComp(toRed))) l++;
335 //assume (l==(**modcomp)[pGetComp(toRed)]);
336 while ((l<tl) && (notFound))
337 #else
338 l = (**modcomp)[pGetComp(toRed)];
339 int kkk = (**modcomp)[pGetComp(toRed)+1];
340 while ((l<kkk) && (notFound))
341 #endif
342 {
343 if ((ecartT[l]<bestEcart) && (pDivisibleBy(T[l],toRed)))
344 {
345 if (ecartT[l]<=ecartToRed) notFound = FALSE;
346 p = T[l];
347 bestEcart = ecartT[l];
348 }
349 l++;
350 }
351 if (p==NULL)
352 {
353 pDelete(&toRed);
354 for(k=j;k<Fl;k++) pDelete(&(pairs[k]));
355 omFreeSize((ADDRESS)pairs,Fl*sizeof(poly));
356 omFreeSize((ADDRESS)ST,Fl*sizeof(poly));
357 omFreeSize((ADDRESS)S,Fl*sizeof(poly));
358 omFreeSize((ADDRESS)T,tmax*sizeof(poly));
359 omFreeSize((ADDRESS)ecartT,tmax*sizeof(int));
360 omFreeSize((ADDRESS)totalT,tmax*sizeof(int));
361 omFreeSize((ADDRESS)ecartS,Fl*sizeof(int));
362 omFreeSize((ADDRESS)totalS,Fl*sizeof(int));
363 for(k=0;k<IDELEMS(result);k++) pDelete(&((*Shdl)[k]));
364 WerrorS("ideal not a standard basis");//no polynom for reduction
365 return result;
366 }
367 else
368 {
369 //Print("reduced with (ecart %d): ",bestEcart);wrp(p);PrintLn();
370 if (notFound)
371 {
372 if (tl>=tmax)
373 {
374 pEnlargeSet(&T,tmax,16);
375 tmax += 16;
376 temp = (int*)omAlloc((tmax+16)*sizeof(int));
377 for(l=0;l<tmax;l++) temp[l]=totalT[l];
378 totalT = temp;
379 temp = (int*)omAlloc((tmax+16)*sizeof(int));
380 for(l=0;l<tmax;l++) temp[l]=ecartT[l];
381 ecartT = temp;
382 }
383 //PrintS("t");
384 int comptR=pGetComp(toRed);
385 for (l=tempcomp->length()-1;l>comptR;l--)
386 {
387 if ((*tempcomp)[l]>0)
388 (*tempcomp)[l]++;
389 }
390 l=0;
391 while ((l<tl) && (comptR>pGetComp(T[l]))) l++;
392 while ((l<tl) && (totalT[l]<=totalToRed)) l++;
393 for (kk=tl;kk>l;kk--)
394 {
395 T[kk]=T[kk-1];
396 totalT[kk]=totalT[kk-1];
397 ecartT[kk]=ecartT[kk-1];
398 }
399 q = pCopy(toRed);
400 pNorm(q);
401 T[l] = q;
402 totalT[l] = totalToRed;
403 ecartT[l] = ecartToRed;
404 tl++;
405 }
406 toRed = ksOldSpolyRed(p,toRed);
407 }
408 }
409 //Print("toRed finally (neuer ecart %d): ",ecartToRed);pWrite(toRed);
410 //PrintS("s");
411 if (pGetComp(toRed)>rkF)
412 {
413 if (Sl>=IDELEMS(result))
414 {
415 pEnlargeSet(Shdl,IDELEMS(result),16);
416 IDELEMS(result) += 16;
417 }
418 //p_Shift(&toRed,-rkF,currRing);
419 pNorm(toRed);
420 (*Shdl)[Sl] = toRed;
421 Sl++;
422 }
423 /*----------------deleting all polys not from ST--------------*/
424 for(l=0;l<tl;l++)
425 {
426 kk=0;
427 while ((kk<smax) && (T[l] != S[kk])) kk++;
428 if (kk>=smax)
429 {
430 pDelete(&T[l]);
431 //Print ("#");
432 }
433 }
434 delete tempcomp;
435 }
436 }
437 for(k=lini;k<wend;k++) pDelete(&(pairs[k]));
438 }
439 (*newmodcomp)[Fl+1] = Sl;
440 omFreeSize((ADDRESS)pairs,Fl*sizeof(poly));
441 omFreeSize((ADDRESS)ST,Fl*sizeof(poly));
442 omFreeSize((ADDRESS)S,Fl*sizeof(poly));
443 omFreeSize((ADDRESS)T,tmax*sizeof(poly));
444 omFreeSize((ADDRESS)ecartT,tmax*sizeof(int));
445 omFreeSize((ADDRESS)totalT,tmax*sizeof(int));
446 omFreeSize((ADDRESS)ecartS,Fl*sizeof(int));
447 omFreeSize((ADDRESS)totalS,Fl*sizeof(int));
448 delete *modcomp;
449 *modcomp = newmodcomp;
450 return result;
451 }
452
453 /*3
454 *special Normalform for Schreyer in factor rings
455 */
sySpecNormalize(poly toNorm,ideal mW=NULL)456 poly sySpecNormalize(poly toNorm,ideal mW=NULL)
457 {
458 int j,i=0;
459 poly p;
460
461 if (toNorm==NULL) return NULL;
462 p = pHead(toNorm);
463 if (mW!=NULL)
464 {
465 for(j=1;j<=(currRing->N);j++)
466 pSetExp(p,j,pGetExp(p,j) -pGetExp(mW->m[pGetComp(p)-1],j));
467 }
468 while ((p!=NULL) && (i<IDELEMS(currRing->qideal)))
469 {
470 if (pDivisibleBy(currRing->qideal->m[i],p))
471 {
472 //pNorm(toNorm);
473 toNorm = ksOldSpolyRed(currRing->qideal->m[i],toNorm);
474 pDelete(&p);
475 if (toNorm==NULL) return NULL;
476 p = pHead(toNorm);
477 if (mW!=NULL)
478 {
479 for(j=1;j<=(currRing->N);j++)
480 pSetExp(p,j,pGetExp(p,j) -pGetExp(mW->m[pGetComp(p)-1],j));
481 }
482 i = 0;
483 }
484 else
485 {
486 i++;
487 }
488 }
489 pDelete(&p);
490 return toNorm;
491 }
492
493 /*2
494 * computes the Schreyer syzygies in the global case
495 * input: F
496 * output: Shdl, Smax
497 * modcomp, length stores the start position of the module comp. in arg
498 */
sySchreyersSyzygiesFB(ideal arg,intvec ** modcomp,ideal mW,BOOLEAN redTail=TRUE)499 static ideal sySchreyersSyzygiesFB(ideal arg,intvec ** modcomp,ideal mW,BOOLEAN redTail=TRUE)
500 {
501 kBucket_pt sy0buck = kBucketCreate(currRing);
502
503 int Fl=IDELEMS(arg);
504 while ((Fl!=0) && (arg->m[Fl-1]==NULL)) Fl--;
505 ideal result=idInit(16,Fl);
506 int i,j,l,k,kkk,/*rkF,*/Sl=0,syComponentOrder=currRing->ComponentOrder;
507 int /*fstart,*/wend,lini,ltR,gencQ=0;
508 intvec *newmodcomp;
509 int *Flength;
510 polyset pairs,F=arg->m,*Shdl=&(result->m);
511 poly /*p,*/q,toRed,syz,lastmonom,multWith;
512 BOOLEAN isNotReduced=TRUE;
513
514 //#define WRITE_BUCKETS
515 #ifdef WRITE_BUCKETS
516 PrintS("Input: \n");
517 ideal twr=idHead(arg);
518 idPrint(arg);
519 idDelete(&twr);
520 if (modcomp!=NULL) (*modcomp)->show(0,0);
521 #endif
522
523 newmodcomp = new intvec(Fl+2);
524 //for (j=0;j<Fl;j++) pWrite(F[j]);
525 //PrintLn();
526 if (currRing->qideal==NULL)
527 pairs=(polyset)omAlloc0(Fl*sizeof(poly));
528 else
529 {
530 gencQ = IDELEMS(currRing->qideal);
531 pairs=(polyset)omAlloc0((Fl+gencQ)*sizeof(poly));
532 }
533 // rkF=id_RankFreeModule(arg,currRing);
534 Flength = (int*)omAlloc0(Fl*sizeof(int));
535 for(j=0;j<Fl;j++)
536 {
537 Flength[j] = pLength(F[j]);
538 }
539 for(j=0;j<Fl;j++)
540 {
541 (*newmodcomp)[j+1] = Sl;
542 if (TEST_OPT_PROT)
543 {
544 Print("(%d)",Fl-j);
545 mflush();
546 }
547 i = pGetComp(F[j]);
548 if (syComponentOrder==1)
549 {
550 lini=k=j+1;
551 wend=Fl;
552 }
553 else
554 {
555 lini=k=0;
556 while ((k<j) && (pGetComp(F[k]) != i)) k++;
557 wend=j;
558 }
559 syCreatePairs(F,lini,wend,k,j,i,pairs,Fl,mW);
560 if (currRing->qideal!=NULL) wend = Fl+gencQ;
561 for (k=lini;k<wend;k++)
562 {
563 if (pairs[k]!=NULL)
564 {
565 if (TEST_OPT_PROT)
566 {
567 PrintS(".");
568 mflush();
569 }
570 //begins to construct the syzygy
571 if (k<Fl)
572 {
573 number an=nCopy(pGetCoeff(F[k])),bn=nCopy(pGetCoeff(F[j]));
574 /*int ct =*/ (void) ksCheckCoeff(&an, &bn, currRing->cf);
575 syz = pCopy(pairs[k]);
576 //syz->coef = nCopy(F[k]->coef);
577 syz->coef = an;
578 //syz->coef = nInpNeg(syz->coef);
579 pNext(syz) = pairs[k];
580 lastmonom = pNext(syz);
581 //lastmonom->coef = nCopy(F[j]->coef);
582 lastmonom->coef = bn;
583 lastmonom->coef = nInpNeg(lastmonom->coef);
584 pSetComp(lastmonom,k+1);
585 }
586 else
587 {
588 syz = pairs[k];
589 syz->coef = nCopy(currRing->qideal->m[k-Fl]->coef);
590 syz->coef = nInpNeg(syz->coef);
591 lastmonom = syz;
592 multWith = pMDivide(syz,F[j]);
593 multWith->coef = nCopy(currRing->qideal->m[k-Fl]->coef);
594 }
595 pSetComp(syz,j+1);
596 pairs[k] = NULL;
597 //the next term of the syzygy
598 //constructs the spoly
599 if (TEST_OPT_DEBUG)
600 {
601 if (k<Fl)
602 {
603 PrintS("pair: ");pWrite0(F[j]);PrintS(" ");pWrite(F[k]);
604 }
605 else
606 {
607 PrintS("pair: ");pWrite0(F[j]);PrintS(" ");pWrite(currRing->qideal->m[k-Fl]);
608 }
609 }
610 if (k<Fl)
611 toRed = ksOldCreateSpoly(F[j],F[k]);
612 else
613 {
614 q = pMult_mm(pCopy(F[j]),multWith);
615 toRed = sySpecNormalize(q,mW);
616 pDelete(&multWith);
617 }
618 kBucketInit(sy0buck,toRed,-1);
619 toRed = kBucketGetLm(sy0buck);
620 isNotReduced = TRUE;
621 while (toRed!=NULL)
622 {
623 if (TEST_OPT_DEBUG)
624 {
625 PrintS("toRed: ");pWrite(toRed);
626 }
627 // l=0;
628 // while ((l<Fl) && (!pDivisibleBy(F[l],toRed))) l++;
629 // if (l>=Fl)
630 l = (**modcomp)[pGetComp(toRed)+1]-1;
631 kkk = (**modcomp)[pGetComp(toRed)];
632 while ((l>=kkk) && (!pDivisibleBy(F[l],toRed))) l--;
633 #ifdef WRITE_BUCKETS
634 kBucketClear(sy0buck,&toRed,<R);
635 printf("toRed in Pair[%d, %d]:", j, k);
636 pWrite(toRed);
637 kBucketInit(sy0buck,toRed,-1);
638 #endif
639
640 if (l<kkk)
641 {
642 if ((currRing->qideal!=NULL) && (isNotReduced))
643 {
644 kBucketClear(sy0buck,&toRed,<R);
645 toRed = sySpecNormalize(toRed,mW);
646 #ifdef WRITE_BUCKETS
647 printf("toRed in Pair[%d, %d]:", j, k);
648 pWrite(toRed);
649 #endif
650 kBucketInit(sy0buck,toRed,-1);
651 toRed = kBucketGetLm(sy0buck);
652 isNotReduced = FALSE;
653 }
654 else
655 {
656 pDelete(&toRed);
657
658 pDelete(&syz);
659 for(k=j;k<Fl;k++) pDelete(&(pairs[k]));
660 omFreeSize((ADDRESS)pairs,(Fl + gencQ)*sizeof(poly));
661
662
663 for(k=0;k<IDELEMS(result);k++) pDelete(&((*Shdl)[k]));
664
665 kBucketDestroy(&(sy0buck));
666
667 //no polynom for reduction
668 WerrorS("ideal not a standard basis");
669
670 return result;
671 }
672 }
673 else
674 {
675 //the next monom of the syzygy
676 isNotReduced = TRUE;
677 if (TEST_OPT_DEBUG)
678 {
679 PrintS("reduced with: ");pWrite(F[l]);
680 }
681 pNext(lastmonom) = pHead(toRed);
682 pIter(lastmonom);
683 lastmonom->coef = nDiv(lastmonom->coef,F[l]->coef);
684 //lastmonom->coef = nInpNeg(lastmonom->coef);
685 pSetComp(lastmonom,l+1);
686 //computes the new toRed
687 number up = kBucketPolyRed(sy0buck,F[l],Flength[l],NULL);
688 if (! nIsOne(up))
689 {
690 // Thomas: Now do whatever you need to do
691 #ifdef WRITE_BUCKETS
692 PrintS("multiplied with: ");nWrite(up);PrintLn();
693 #endif
694 syz=__p_Mult_nn(syz,up,currRing);
695 }
696 nDelete(&up);
697
698 toRed = kBucketGetLm(sy0buck);
699 //the module component of the new monom
700 //pWrite(toRed);
701 }
702 }
703 kBucketClear(sy0buck,&toRed,<R); //Zur Sichereheit
704 //PrintLn();
705 if (syz!=NULL)
706 {
707 if (Sl>=IDELEMS(result))
708 {
709 pEnlargeSet(Shdl,IDELEMS(result),16);
710 IDELEMS(result) += 16;
711 }
712 pNorm(syz);
713 if (BTEST1(OPT_REDTAIL) && redTail)
714 {
715 (*newmodcomp)[j+2] = Sl;
716 (*Shdl)[Sl] = syRedtail2(syz,*Shdl,newmodcomp);
717 (*newmodcomp)[j+2] = 0;
718 }
719 else
720 (*Shdl)[Sl] = syz;
721 Sl++;
722 }
723 }
724 }
725 // for(k=j;k<Fl;k++) pDelete(&(pairs[k]));
726 }
727 (*newmodcomp)[Fl+1] = Sl;
728 if (currRing->qideal==NULL)
729 omFreeSize((ADDRESS)pairs,Fl*sizeof(poly));
730 else
731 omFreeSize((ADDRESS)pairs,(Fl+IDELEMS(currRing->qideal))*sizeof(poly));
732 omFreeSize((ADDRESS)Flength,Fl*sizeof(int));
733 delete *modcomp;
734 *modcomp = newmodcomp;
735
736 kBucketDestroy(&(sy0buck));
737 return result;
738 }
739
syReOrderResolventFB(resolvente res,int length,int initial)740 void syReOrderResolventFB(resolvente res,int length, int initial)
741 {
742 int syzIndex=length-1,i,j;
743 poly p;
744
745 while ((syzIndex!=0) && (res[syzIndex]==NULL)) syzIndex--;
746 while (syzIndex>=initial)
747 {
748 for(i=0;i<IDELEMS(res[syzIndex]);i++)
749 {
750 p = res[syzIndex]->m[i];
751
752 while (p!=NULL)
753 {
754 if (res[syzIndex-1]->m[pGetComp(p)-1]!=NULL)
755 {
756 for(j=1;j<=(currRing->N);j++)
757 {
758 pSetExp(p,j,pGetExp(p,j)
759 -pGetExp(res[syzIndex-1]->m[pGetComp(p)-1],j));
760 }
761 }
762 else
763 PrintS("error in the resolvent\n");
764 pSetm(p);
765 pIter(p);
766 }
767 }
768 syzIndex--;
769 }
770 }
771
772 #if 0
773 static void syMergeSortResolventFB(resolvente res,int length, int initial=1)
774 {
775 int syzIndex=length-1,i,j;
776 poly qq,pp,result=NULL;
777 poly p;
778
779 while ((syzIndex!=0) && (res[syzIndex]==NULL)) syzIndex--;
780 while (syzIndex>=initial)
781 {
782 for(i=0;i<IDELEMS(res[syzIndex]);i++)
783 {
784 p = res[syzIndex]->m[i];
785 if (p != NULL)
786 {
787 for (;;)
788 {
789 qq = p;
790 for(j=1;j<=(currRing->N);j++)
791 {
792 pSetExp(p,j,pGetExp(p,j)
793 -pGetExp(res[syzIndex-1]->m[pGetComp(p)-1],j));
794 }
795 pSetm(p);
796 for (;;)
797 {
798 if (pNext(p) == NULL)
799 {
800 pAdd(result, qq);
801 break;
802 }
803 pp = pNext(p);
804 for(j=1;j<=(currRing->N);j++)
805 {
806 pSetExp(pp,j,pGetExp(pp,j)
807 -pGetExp(res[syzIndex-1]->m[pGetComp(pp)-1],j));
808 }
809 pSetm(pp);
810 if (pCmp(p,pNext(p)) != 1)
811 {
812 pp = p;
813 pIter(p);
814 pNext(pp) = NULL;
815 result = pAdd(result, qq);
816 break;
817 }
818 pIter(p);
819 }
820 }
821 }
822 res[syzIndex]->m[i] = p;
823 }
824 syzIndex--;
825 }
826 }
827 #endif
828
syTestOrder(ideal M)829 BOOLEAN syTestOrder(ideal M)
830 {
831 int i=id_RankFreeModule(M,currRing);
832 if (i == 0) return FALSE;
833 int j=0;
834
835 while ((currRing->order[j]!=ringorder_c) && (currRing->order[j]!=ringorder_C))
836 j++;
837 if (currRing->order[j+1]!=0)
838 return TRUE;
839 return FALSE;
840 }
841
842 #if 0 /*debug only */
843 static void syPrintResolution(resolvente res,int start,int length)
844 {
845 while ((start < length) && (res[start]))
846 {
847 Print("Syz(%d): \n",start);
848 idTest(res[start]);
849 //idPrint(res[start]);
850 start++;
851 }
852 }
853 #endif
854
sySchreyerResolvente(ideal arg,int maxlength,int * length,BOOLEAN isMonomial,BOOLEAN)855 resolvente sySchreyerResolvente(ideal arg, int maxlength, int * length,
856 BOOLEAN isMonomial, BOOLEAN /*notReplace*/)
857 {
858 ideal mW=NULL;
859 int i,syzIndex = 0,j=0;
860 intvec * modcomp=NULL,*w=NULL;
861 // int ** wv=NULL;
862 tHomog hom=(tHomog)idHomModule(arg,NULL,&w);
863 ring origR = currRing;
864 ring syRing = NULL;
865
866 if ((!isMonomial) && syTestOrder(arg))
867 {
868 WerrorS("sres only implemented for modules with ordering ..,c or ..,C");
869 return NULL;
870 }
871 *length = 4;
872 resolvente res = (resolvente)omAlloc0(4*sizeof(ideal)),newres;
873 res[0] = idCopy(arg);
874
875 while ((!idIs0(res[syzIndex])) && ((maxlength==-1) || (syzIndex<maxlength)))
876 {
877 i = IDELEMS(res[syzIndex]);
878 //while ((i!=0) && (!res[syzIndex]->m[i-1])) i--;
879 if (syzIndex+1==*length)
880 {
881 newres = (resolvente)omAlloc0((*length+4)*sizeof(ideal));
882 // for (j=0;j<*length+4;j++) newres[j] = NULL;
883 for (j=0;j<*length;j++) newres[j] = res[j];
884 omFreeSize((ADDRESS)res,*length*sizeof(ideal));
885 *length += 4;
886 res=newres;
887 }
888
889 if ((hom==isHomog)|| (rHasGlobalOrdering(origR)))
890 {
891 if (syzIndex==0) syInitSort(res[0],&modcomp);
892
893 if ((syzIndex==0) && !rRing_has_CompLastBlock(currRing))
894 res[syzIndex+1] = sySchreyersSyzygiesFB(res[syzIndex],&modcomp,mW,FALSE);
895 else
896 res[syzIndex+1] = sySchreyersSyzygiesFB(res[syzIndex],&modcomp,mW);
897
898 if (errorreported)
899 {
900 for (j=0;j<*length;j++) idDelete( &res[j] );
901 omFreeSize((ADDRESS)res,*length*sizeof(ideal));
902 return NULL;
903 }
904
905 mW = res[syzIndex];
906 }
907 //idPrint(res[syzIndex+1]);
908
909 if ( /*(*/ syzIndex==0 /*)*/ )
910 {
911 if ((hom==isHomog)|| (rHasGlobalOrdering(origR)))
912 {
913 syRing = rAssure_CompLastBlock(origR, TRUE);
914 if (syRing != origR)
915 {
916 rChangeCurrRing(syRing);
917 for (i=0; i<IDELEMS(res[1]); i++)
918 {
919 res[1]->m[i] = prMoveR( res[1]->m[i], origR, syRing);
920 }
921 }
922 idTest(res[1]);
923 }
924 else
925 {
926 syRing = rAssure_SyzComp_CompLastBlock(origR);
927 if (syRing != origR)
928 {
929 rChangeCurrRing(syRing);
930 for (i=0; i<IDELEMS(res[0]); i++)
931 {
932 res[0]->m[i] = prMoveR( res[0]->m[i], origR, syRing);
933 }
934 }
935 idTest(res[0]);
936 }
937 }
938 if ((hom!=isHomog) && (rHasLocalOrMixedOrdering(origR)))
939 {
940 if (syzIndex==0) syInitSort(res[0],&modcomp);
941 res[syzIndex+1] = sySchreyersSyzygiesFM(res[syzIndex],&modcomp);
942 if (errorreported)
943 {
944 for (j=0;j<*length;j++) idDelete( &res[j] );
945 omFreeSize((ADDRESS)res,*length*sizeof(ideal));
946 return NULL;
947 }
948 }
949 syzIndex++;
950 if (TEST_OPT_PROT) Print("[%d]\n",syzIndex);
951 }
952 //syPrintResolution(res,1,*length);
953 if ((hom!=isHomog) && (rHasLocalOrMixedOrdering(origR)))
954 {
955 syzIndex = 1;
956 while ((syzIndex < *length) && (!idIs0(res[syzIndex])))
957 {
958 id_Shift(res[syzIndex],-rGetMaxSyzComp(syzIndex, currRing),currRing);
959 syzIndex++;
960 }
961 }
962 if ((hom==isHomog) || (rHasGlobalOrdering(origR)))
963 syzIndex = 1;
964 else
965 syzIndex = 0;
966 syReOrderResolventFB(res,*length,syzIndex+1);
967 if (/*ringOrderChanged:*/ origR!=syRing && syRing != NULL)
968 {
969 rChangeCurrRing(origR);
970 // Thomas: Here I assume that all (!) polys of res live in tmpR
971 while ((syzIndex < *length) && (res[syzIndex]))
972 {
973 for (i=0;i<IDELEMS(res[syzIndex]);i++)
974 {
975 if (res[syzIndex]->m[i])
976 {
977 res[syzIndex]->m[i] = prMoveR( res[syzIndex]->m[i], syRing, origR);
978 }
979 }
980 syzIndex++;
981 }
982 // j = 0; while (currRing->order[j]!=0) j++; // What was this for???!
983 rDelete(syRing);
984 }
985 else
986 {
987 // Thomas -- are you sure that you have to "reorder" here?
988 while ((syzIndex < *length) && (res[syzIndex]))
989 {
990 for (i=0;i<IDELEMS(res[syzIndex]);i++)
991 {
992 if (res[syzIndex]->m[i])
993 res[syzIndex]->m[i] = pSortCompCorrect(res[syzIndex]->m[i]);
994 }
995 syzIndex++;
996 }
997 }
998 if ((hom==isHomog) || (rHasGlobalOrdering(origR)))
999 {
1000 if (res[1]!=NULL)
1001 {
1002 syReOrderResolventFB(res,2,1);
1003 for (i=0;i<IDELEMS(res[1]);i++)
1004 {
1005 if (res[1]->m[i])
1006 res[1]->m[i] = pSort(res[1]->m[i]);
1007 }
1008 }
1009 }
1010 //syPrintResolution(res,0,*length);
1011
1012 //syMergeSortResolventFB(res,*length);
1013 if (modcomp!=NULL) delete modcomp;
1014 if (w!=NULL) delete w;
1015 return res;
1016 }
1017
sySchreyer(ideal arg,int maxlength)1018 syStrategy sySchreyer(ideal arg, int maxlength)
1019 {
1020 int rl;
1021 resolvente fr = sySchreyerResolvente(arg,maxlength,&(rl));
1022 if (fr==NULL) return NULL;
1023
1024 // int typ0;
1025 syStrategy result=(syStrategy)omAlloc0(sizeof(ssyStrategy));
1026 result->length=rl;
1027 result->fullres = (resolvente)omAlloc0((rl /*result->length*/+1)*sizeof(ideal));
1028 for (int i=rl /*result->length*/-1;i>=0;i--)
1029 {
1030 if (fr[i]!=NULL)
1031 {
1032 idSkipZeroes(fr[i]);
1033 result->fullres[i] = fr[i];
1034 fr[i] = NULL;
1035 }
1036 }
1037 if (currRing->qideal!=NULL)
1038 {
1039 for (int i=0; i<rl; i++)
1040 {
1041 if (result->fullres[i]!=NULL)
1042 {
1043 ideal t=kNF(currRing->qideal,NULL,result->fullres[i]);
1044 idDelete(&result->fullres[i]);
1045 result->fullres[i]=t;
1046 if (i<rl-1)
1047 {
1048 for(int j=IDELEMS(t)-1;j>=0; j--)
1049 {
1050 if ((t->m[j]==NULL) && (result->fullres[i+1]!=NULL))
1051 {
1052 for(int k=IDELEMS(result->fullres[i+1])-1;k>=0; k--)
1053 {
1054 if (result->fullres[i+1]->m[k]!=NULL)
1055 {
1056 pDeleteComp(&(result->fullres[i+1]->m[k]),j+1);
1057 }
1058 }
1059 }
1060 }
1061 }
1062 idSkipZeroes(result->fullres[i]);
1063 }
1064 }
1065 if ((rl>maxlength) && (result->fullres[rl-1]!=NULL))
1066 {
1067 idDelete(&result->fullres[rl-1]);
1068 }
1069 }
1070 omFreeSize((ADDRESS)fr,(rl /*result->length*/)*sizeof(ideal));
1071 return result;
1072 }
1073
1074