source: git/kernel/syz0.cc @ 601105

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