My Project
syz0.cc
Go to the documentation of this file.
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"
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 
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 
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 
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 */
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 */
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 */
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,&ltR);
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,&ltR);
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,&ltR); //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 
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  {
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 
830 {
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 
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  {
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 
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;
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 
int BOOLEAN
Definition: auxiliary.h:87
#define TRUE
Definition: auxiliary.h:100
#define FALSE
Definition: auxiliary.h:96
void * ADDRESS
Definition: auxiliary.h:119
CanonicalForm FACTORY_PUBLIC pp(const CanonicalForm &)
CanonicalForm pp ( const CanonicalForm & f )
Definition: cf_gcd.cc:676
int l
Definition: cfEzgcd.cc:100
int m
Definition: cfEzgcd.cc:128
int i
Definition: cfEzgcd.cc:132
int k
Definition: cfEzgcd.cc:99
int p
Definition: cfModGcd.cc:4078
Definition: intvec.h:23
int length() const
Definition: intvec.h:94
#define Print
Definition: emacs.cc:80
return result
Definition: facAbsBiFact.cc:75
CanonicalForm res
Definition: facAbsFact.cc:60
const CanonicalForm & w
Definition: facAbsFact.cc:51
int j
Definition: facHensel.cc:110
VAR short errorreported
Definition: feFopen.cc:23
void WerrorS(const char *s)
Definition: feFopen.cc:24
template CanonicalForm tmax(const CanonicalForm &, const CanonicalForm &)
ideal idHead(ideal h)
#define idDelete(H)
delete an ideal
Definition: ideals.h:29
BOOLEAN idIs0(ideal h)
returns true if h is the zero ideal
#define idPrint(id)
Definition: ideals.h:46
static BOOLEAN idHomModule(ideal m, ideal Q, intvec **w)
Definition: ideals.h:96
#define idTest(id)
Definition: ideals.h:47
ideal idCopy(ideal A)
Definition: ideals.h:60
ideal * resolvente
Definition: ideals.h:18
poly initial(const poly p, const ring r, const gfan::ZVector &w)
Returns the initial form of p with respect to w.
Definition: initial.cc:30
static BOOLEAN length(leftv result, leftv arg)
Definition: interval.cc:257
intvec * ivCopy(const intvec *o)
Definition: intvec.h:145
STATIC_VAR jList * T
Definition: janet.cc:30
STATIC_VAR Poly * h
Definition: janet.cc:971
KINLINE poly ksOldCreateSpoly(poly p1, poly p2, poly spNoether, ring r)
Definition: kInline.h:1196
KINLINE poly ksOldSpolyRed(poly p1, poly p2, poly spNoether)
Definition: kInline.h:1176
void kBucketClear(kBucket_pt bucket, poly *p, int *length)
Definition: kbuckets.cc:521
void kBucketDestroy(kBucket_pt *bucket_pt)
Definition: kbuckets.cc:216
void kBucketInit(kBucket_pt bucket, poly lm, int length)
Definition: kbuckets.cc:493
int ksCheckCoeff(number *a, number *b, const coeffs r)
Definition: kbuckets.cc:1504
kBucket_pt kBucketCreate(const ring bucket_ring)
Creation/Destruction of buckets.
Definition: kbuckets.cc:209
number kBucketPolyRed(kBucket_pt bucket, poly p1, int l1, poly spNoether)
Definition: kbuckets.cc:1071
const poly kBucketGetLm(kBucket_pt bucket)
Definition: kbuckets.cc:506
poly kNF(ideal F, ideal Q, poly p, int syzComp, int lazyReduce)
Definition: kstd1.cc:3182
void pairs()
#define assume(x)
Definition: mod2.h:389
#define pIter(p)
Definition: monomials.h:37
#define pNext(p)
Definition: monomials.h:36
static number & pGetCoeff(poly p)
return an alias to the leading coefficient of p assumes that p != NULL NOTE: not copy
Definition: monomials.h:44
#define nDiv(a, b)
Definition: numbers.h:32
#define nWrite(n)
Definition: numbers.h:29
#define nDelete(n)
Definition: numbers.h:16
#define nInpNeg(n)
Definition: numbers.h:21
#define nCopy(n)
Definition: numbers.h:15
#define nIsOne(n)
Definition: numbers.h:25
#define omFreeSize(addr, size)
Definition: omAllocDecl.h:260
#define omAlloc(size)
Definition: omAllocDecl.h:210
#define omAlloc0(size)
Definition: omAllocDecl.h:211
#define NULL
Definition: omList.c:12
#define OPT_REDTAIL
Definition: options.h:92
#define TEST_OPT_PROT
Definition: options.h:104
#define BTEST1(a)
Definition: options.h:34
#define TEST_OPT_DEBUG
Definition: options.h:109
void pEnlargeSet(poly **p, int l, int increment)
Definition: p_polys.cc:3692
static int pLength(poly a)
Definition: p_polys.h:188
static long p_FDeg(const poly p, const ring r)
Definition: p_polys.h:378
static long p_LDeg(const poly p, int *l, const ring r)
Definition: p_polys.h:379
#define __p_Mult_nn(p, n, r)
Definition: p_polys.h:969
void rChangeCurrRing(ring r)
Definition: polys.cc:15
VAR ring currRing
Widely used global variable which specifies the current polynomial ring for Singular interpreter and ...
Definition: polys.cc:13
Compatibility layer for legacy polynomial operations (over currRing)
#define pAdd(p, q)
Definition: polys.h:203
#define pSort(p)
Definition: polys.h:218
#define pDelete(p_ptr)
Definition: polys.h:186
#define pHead(p)
returns newly allocated copy of Lm(p), coef is copied, next=NULL, p might be NULL
Definition: polys.h:67
#define pSetm(p)
Definition: polys.h:271
#define pSetCompP(a, i)
Definition: polys.h:303
#define pDeleteComp(p, k)
Definition: polys.h:360
#define pGetComp(p)
Component.
Definition: polys.h:37
void pNorm(poly p)
Definition: polys.h:362
#define pCmp(p1, p2)
pCmp: args may be NULL returns: (p2==NULL ? 1 : (p1 == NULL ? -1 : p_LmCmp(p1, p2)))
Definition: polys.h:115
#define pSetComp(p, v)
Definition: polys.h:38
void pWrite0(poly p)
Definition: polys.h:309
void pWrite(poly p)
Definition: polys.h:308
#define pGetExp(p, i)
Exponent.
Definition: polys.h:41
#define pSetmComp(p)
TODO:
Definition: polys.h:273
#define pMult_mm(p, m)
Definition: polys.h:202
#define pDivisibleBy(a, b)
returns TRUE, if leading monom of a divides leading monom of b i.e., if there exists a expvector c > ...
Definition: polys.h:138
#define pSetExp(p, i, v)
Definition: polys.h:42
#define pLmCmp(p, q)
returns 0|1|-1 if p=q|p>q|p<q w.r.t monomial ordering
Definition: polys.h:105
#define pMDivide(a, b)
Definition: polys.h:293
#define pCopy(p)
return a copy of the poly
Definition: polys.h:185
#define pOne()
Definition: polys.h:315
poly * polyset
Definition: polys.h:259
#define pSortCompCorrect(p)
Assume: If considered only as poly in any component of p (say, monomials of other components of p are...
Definition: polys.h:227
#define pLcm(a, b, m)
Definition: polys.h:295
#define pLmDivisibleByNoComp(a, b)
like pLmDivisibleBy, does not check components
Definition: polys.h:142
poly prMoveR(poly &p, ring src_r, ring dest_r)
Definition: prCopy.cc:90
void PrintS(const char *s)
Definition: reporter.cc:284
void PrintLn()
Definition: reporter.cc:310
#define mflush()
Definition: reporter.h:58
BOOLEAN rRing_has_CompLastBlock(const ring r)
Definition: ring.cc:5185
int rGetMaxSyzComp(int i, const ring r)
return the max-comonent wchich has syzIndex i Assume: i<= syzIndex_limit
Definition: ring.cc:5158
ring rAssure_SyzComp_CompLastBlock(const ring r)
makes sure that c/C ordering is last ordering and SyzIndex is first
Definition: ring.cc:4749
void rDelete(ring r)
unconditionally deletes fields in r
Definition: ring.cc:450
ring rAssure_CompLastBlock(ring r, BOOLEAN complete)
makes sure that c/C ordering is last ordering
Definition: ring.cc:4694
void rSetSyzComp(int k, const ring r)
Definition: ring.cc:5086
@ ringorder_C
Definition: ring.h:73
@ ringorder_c
Definition: ring.h:72
BOOLEAN rHasGlobalOrdering(const ring r)
Definition: ring.h:759
BOOLEAN rHasLocalOrMixedOrdering(const ring r)
Definition: ring.h:760
ideal idInit(int idsize, int rank)
initialise an ideal / module
Definition: simpleideals.cc:35
long id_RankFreeModule(ideal s, ring lmRing, ring tailRing)
return the maximal component number found in any polynomial in s
void idSkipZeroes(ideal ide)
gives an ideal/module the minimal possible size
void id_Shift(ideal M, int s, const ring r)
#define IDELEMS(i)
Definition: simpleideals.h:23
#define M
Definition: sirandom.c:25
tHomog
Definition: structs.h:35
@ isHomog
Definition: structs.h:37
resolvente sySchreyerResolvente(ideal arg, int maxlength, int *length, BOOLEAN isMonomial, BOOLEAN)
Definition: syz0.cc:855
static ideal sySchreyersSyzygiesFB(ideal arg, intvec **modcomp, ideal mW, BOOLEAN redTail=TRUE)
Definition: syz0.cc:499
syStrategy sySchreyer(ideal arg, int maxlength)
Definition: syz0.cc:1018
static void syCreatePairs(polyset F, int lini, int wend, int k, int j, int i, polyset pairs, int regularPairs=0, ideal mW=NULL)
Definition: syz0.cc:68
static void syInitSort(ideal arg, intvec **modcomp)
Definition: syz0.cc:23
poly sySpecNormalize(poly toNorm, ideal mW=NULL)
Definition: syz0.cc:456
static ideal sySchreyersSyzygiesFM(ideal arg, intvec **modcomp)
Definition: syz0.cc:162
static poly syRedtail2(poly p, polyset redWith, intvec *modcomp)
Definition: syz0.cc:118
void syReOrderResolventFB(resolvente res, int length, int initial)
Definition: syz0.cc:740
BOOLEAN syTestOrder(ideal M)
Definition: syz0.cc:829
ssyStrategy * syStrategy
Definition: syz.h:35