source: git/kernel/syz0.cc @ bc669c

fieker-DuValspielwiese
Last change on this file since bc669c was 690e21e, checked in by Hans Schönemann <hannes@…>, 14 years ago
moved option marcos to options.h git-svn-id: file:///usr/local/Singular/svn/trunk@12466 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 26.7 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/* $Id$ */
5/*
6* ABSTRACT: resolutions
7*/
8
9
10#include "mod2.h"
11#include "options.h"
12#include "omalloc.h"
13#include "polys.h"
14#include "febase.h"
15#include "kstd1.h"
16#include "kutil.h"
17#include "stairc.h"
18//#include "cntrlc.h"
19#include "intvec.h"
20#include "numbers.h"
21#include "ideals.h"
22#include "intvec.h"
23#include "ring.h"
24#include "syz.h"
25#include "kbuckets.h"
26#include "prCopy.h"
27
28static kBucket_pt sy0buck;
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=idRankFreeModule(arg);
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<=pVariables;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<=pVariables;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,kkk;
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=idRankFreeModule(arg);
210/*----------------construction of the new ordering----------*/
211  if (rkF>0)
212    rSetSyzComp(rkF);
213  else
214    rSetSyzComp(1);
215/*----------------creating S--------------------------------*/
216  for(j=0;j<Fl;j++)
217  {
218    S[j] = pCopy(F[j]);
219    totalS[j] = pLDeg(S[j],&k,currRing);
220    ecartS[j] = totalS[j]-pFDeg(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 = pLDeg(toRed,&kk,currRing);
330            ecartToRed = totalToRed-pFDeg(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          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          //pShift(&toRed,-rkF);
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<=pVariables;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<=pVariables;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  int Fl=IDELEMS(arg);
509  while ((Fl!=0) && (arg->m[Fl-1]==NULL)) Fl--;
510  ideal result=idInit(16,Fl);
511  int i,j,l,k,kkk,rkF,Sl=0,syComponentOrder=currRing->ComponentOrder;
512  int fstart,wend,lini,ltR,gencQ=0;
513  intvec *newmodcomp;
514  int *Flength;
515  polyset pairs,F=arg->m,*Shdl=&(result->m);
516  poly p,q,toRed,syz,lastmonom,multWith;
517  BOOLEAN isNotReduced=TRUE;
518
519//#define WRITE_BUCKETS
520#ifdef WRITE_BUCKETS
521PrintS("Input: \n");
522ideal twr=idHead(arg);
523idPrint(arg);
524idDelete(&twr);
525if (modcomp!=NULL) (*modcomp)->show(0,0);
526#endif
527  newmodcomp = new intvec(Fl+2);
528//for (j=0;j<Fl;j++) pWrite(F[j]);
529//PrintLn();
530  if (currQuotient==NULL)
531    pairs=(polyset)omAlloc0(Fl*sizeof(poly));
532  else
533  {
534    gencQ = IDELEMS(currQuotient);
535    pairs=(polyset)omAlloc0((Fl+gencQ)*sizeof(poly));
536  }
537  rkF=idRankFreeModule(arg);
538  Flength = (int*)omAlloc0(Fl*sizeof(int));
539  for(j=0;j<Fl;j++)
540  {
541    Flength[j] = pLength(F[j]);
542  }
543  for(j=0;j<Fl;j++)
544  {
545    (*newmodcomp)[j+1] = Sl;
546    if (TEST_OPT_PROT)
547    {
548      Print("(%d)",Fl-j);
549      mflush();
550    }
551    i = pGetComp(F[j]);
552    if (syComponentOrder==1)
553    {
554      lini=k=j+1;
555      wend=Fl;
556    }
557    else
558    {
559      lini=k=0;
560      while ((k<j) && (pGetComp(F[k]) != i)) k++;
561      wend=j;
562    }
563    syCreatePairs(F,lini,wend,k,j,i,pairs,Fl,mW);
564    if (currQuotient!=NULL) wend = Fl+gencQ;
565    for (k=lini;k<wend;k++)
566    {
567      if (pairs[k]!=NULL)
568      {
569        if (TEST_OPT_PROT)
570        {
571           PrintS(".");
572           mflush();
573        }
574        //begins to construct the syzygy
575        if (k<Fl)
576        {
577          number an=nCopy(pGetCoeff(F[k])),bn=nCopy(pGetCoeff(F[j]));
578          int ct = ksCheckCoeff(&an, &bn);
579          syz = pCopy(pairs[k]);
580          //syz->coef = nCopy(F[k]->coef);
581          syz->coef = an;
582          //syz->coef = nNeg(syz->coef);
583          pNext(syz) = pairs[k];
584          lastmonom = pNext(syz);
585          //lastmonom->coef = nCopy(F[j]->coef);
586          lastmonom->coef = bn;
587          lastmonom->coef = nNeg(lastmonom->coef);
588          pSetComp(lastmonom,k+1);
589        }
590        else
591        {
592          syz = pairs[k];
593          syz->coef = nCopy(currQuotient->m[k-Fl]->coef);
594          syz->coef = nNeg(syz->coef);
595          lastmonom = syz;
596          multWith = pDivide(syz,F[j]);
597          multWith->coef = nCopy(currQuotient->m[k-Fl]->coef);
598        }
599        pSetComp(syz,j+1);
600        pairs[k] = NULL;
601        //the next term of the syzygy
602        //constructs the spoly
603        if (TEST_OPT_DEBUG)
604        {
605          if (k<Fl)
606          {
607            PrintS("pair: ");pWrite0(F[j]);PrintS("  ");pWrite(F[k]);
608          }
609          else
610          {
611            PrintS("pair: ");pWrite0(F[j]);PrintS("  ");pWrite(currQuotient->m[k-Fl]);
612          }
613        }
614        if (k<Fl)
615          toRed = ksOldCreateSpoly(F[j],F[k]);
616        else
617        {
618          q = pMult_mm(pCopy(F[j]),multWith);
619          toRed = sySpecNormalize(q,mW);
620          pDelete(&multWith);
621        }
622        kBucketInit(sy0buck,toRed,-1);
623        toRed = kBucketGetLm(sy0buck);
624        isNotReduced = TRUE;
625        while (toRed!=NULL)
626        {
627          if (TEST_OPT_DEBUG)
628          {
629            PrintS("toRed: ");pWrite(toRed);
630          }
631//          l=0;
632//          while ((l<Fl) && (!pDivisibleBy(F[l],toRed))) l++;
633//          if (l>=Fl)
634          l = (**modcomp)[pGetComp(toRed)+1]-1;
635          kkk = (**modcomp)[pGetComp(toRed)];
636          while ((l>=kkk) && (!pDivisibleBy(F[l],toRed))) l--;
637#ifdef WRITE_BUCKETS
638          kBucketClear(sy0buck,&toRed,&ltR);
639          printf("toRed in Pair[%d, %d]:", j, k);
640          pWrite(toRed);
641          kBucketInit(sy0buck,toRed,-1);
642#endif
643
644          if (l<kkk)
645          {
646            if ((currQuotient!=NULL) && (isNotReduced))
647            {
648              kBucketClear(sy0buck,&toRed,&ltR);
649              toRed = sySpecNormalize(toRed,mW);
650#ifdef WRITE_BUCKETS
651              printf("toRed in Pair[%d, %d]:", j, k);
652              pWrite(toRed);
653#endif
654              kBucketInit(sy0buck,toRed,-1);
655              toRed = kBucketGetLm(sy0buck);
656              isNotReduced = FALSE;
657            }
658            else
659            {
660              //no polynom for reduction
661              WerrorS("ideal not a standard basis");
662              pDelete(&toRed);
663              pDelete(&syz);
664              for(k=j;k<Fl;k++) pDelete(&(pairs[k]));
665              omFreeSize((ADDRESS)pairs,(Fl + gencQ)*sizeof(poly));
666              for(k=0;k<IDELEMS(result);k++) pDelete(&((*Shdl)[k]));
667              return result;
668            }
669          }
670          else
671          {
672            //the next monom of the syzygy
673            isNotReduced = TRUE;
674            if (TEST_OPT_DEBUG)
675            {
676              PrintS("reduced with: ");pWrite(F[l]);
677            }
678            pNext(lastmonom) = pHead(toRed);
679            pIter(lastmonom);
680            lastmonom->coef = nDiv(lastmonom->coef,F[l]->coef);
681            //lastmonom->coef = nNeg(lastmonom->coef);
682            pSetComp(lastmonom,l+1);
683            //computes the new toRed
684            number up = kBucketPolyRed(sy0buck,F[l],Flength[l],NULL);
685            if (! nIsOne(up))
686            {
687              // Thomas: Now do whatever you need to do
688#ifdef WRITE_BUCKETS
689              PrintS("multiplied with: ");nWrite(up);PrintLn();
690#endif
691              pMult_nn(syz,up);
692            }
693            nDelete(&up);
694
695            toRed = kBucketGetLm(sy0buck);
696            //the module component of the new monom
697//pWrite(toRed);
698          }
699        }
700        kBucketClear(sy0buck,&toRed,&ltR); //Zur Sichereheit
701//PrintLn();
702        if (syz!=NULL)
703        {
704          if (Sl>=IDELEMS(result))
705          {
706            pEnlargeSet(Shdl,IDELEMS(result),16);
707            IDELEMS(result) += 16;
708          }
709          pNorm(syz);
710          if (BTEST1(OPT_REDTAIL) && redTail)
711          {
712            (*newmodcomp)[j+2] = Sl;
713            (*Shdl)[Sl] = syRedtail2(syz,*Shdl,newmodcomp);
714            (*newmodcomp)[j+2] = 0;
715          }
716          else
717            (*Shdl)[Sl] = syz;
718          Sl++;
719        }
720      }
721    }
722//    for(k=j;k<Fl;k++) pDelete(&(pairs[k]));
723  }
724  (*newmodcomp)[Fl+1] = Sl;
725  if (currQuotient==NULL)
726    omFreeSize((ADDRESS)pairs,Fl*sizeof(poly));
727  else
728    omFreeSize((ADDRESS)pairs,(Fl+IDELEMS(currQuotient))*sizeof(poly));
729  omFreeSize((ADDRESS)Flength,Fl*sizeof(int));
730  delete *modcomp;
731  *modcomp = newmodcomp;
732  return result;
733}
734
735void syReOrderResolventFB(resolvente res,int length, int initial)
736{
737  int syzIndex=length-1,i,j;
738  poly p;
739
740  while ((syzIndex!=0) && (res[syzIndex]==NULL)) syzIndex--;
741  while (syzIndex>=initial)
742  {
743    for(i=0;i<IDELEMS(res[syzIndex]);i++)
744    {
745      p = res[syzIndex]->m[i];
746
747      while (p!=NULL)
748      {
749        if (res[syzIndex-1]->m[pGetComp(p)-1]!=NULL)
750        {
751          for(j=1;j<=pVariables;j++)
752          {
753            pSetExp(p,j,pGetExp(p,j)
754                        -pGetExp(res[syzIndex-1]->m[pGetComp(p)-1],j));
755          }
756        }
757        else
758          PrintS("error in the resolvent\n");
759        pSetm(p);
760        pIter(p);
761      }
762    }
763    syzIndex--;
764  }
765}
766
767static void syMergeSortResolventFB(resolvente res,int length, int initial=1)
768{
769  int syzIndex=length-1,i,j;
770  poly qq,pp,result=NULL;
771  poly p;
772
773  while ((syzIndex!=0) && (res[syzIndex]==NULL)) syzIndex--;
774  while (syzIndex>=initial)
775  {
776    for(i=0;i<IDELEMS(res[syzIndex]);i++)
777    {
778      p = res[syzIndex]->m[i];
779      if (p != NULL)
780      {
781        for (;;)
782        {
783          qq = p;
784          for(j=1;j<=pVariables;j++)
785          {
786            pSetExp(p,j,pGetExp(p,j)
787                        -pGetExp(res[syzIndex-1]->m[pGetComp(p)-1],j));
788          }
789          pSetm(p);
790          for (;;)
791          {
792            if (pNext(p) == NULL)
793            {
794              pAdd(result, qq);
795              break;
796            }
797            pp = pNext(p);
798            for(j=1;j<=pVariables;j++)
799            {
800              pSetExp(pp,j,pGetExp(pp,j)
801                          -pGetExp(res[syzIndex-1]->m[pGetComp(pp)-1],j));
802            }
803            pSetm(pp);
804            if (pCmp(p,pNext(p)) != 1)
805            {
806              pp = p;
807              pIter(p);
808              pNext(pp) = NULL;
809              result = pAdd(result, qq);
810              break;
811            }
812            pIter(p);
813          }
814        }
815      }
816      res[syzIndex]->m[i] = p;
817    }
818    syzIndex--;
819  }
820}
821
822BOOLEAN syTestOrder(ideal M)
823{
824  int i=idRankFreeModule(M);
825  if (i == 0) return FALSE;
826  int j=0;
827
828  while ((currRing->order[j]!=ringorder_c) && (currRing->order[j]!=ringorder_C))
829    j++;
830  if (currRing->order[j+1]!=0)
831    return TRUE;
832  return FALSE;
833}
834
835static void idShift(ideal arg,int index)
836{
837  int i,j=rGetMaxSyzComp(index);
838  for (i=0;i<IDELEMS(arg);i++)
839  {
840    if (arg->m[i]!=NULL)
841      pShift(&arg->m[i],-j);
842  }
843}
844
845static void syPrintResolution(resolvente res,int start,int length)
846{
847  while ((start < length) && (res[start]))
848  {
849    Print("Syz(%d): \n",start);
850    idTest(res[start]);
851    //idPrint(res[start]);
852    start++;
853  }
854}
855
856resolvente sySchreyerResolvente(ideal arg, int maxlength, int * length,
857                                BOOLEAN isMonomial, BOOLEAN notReplace)
858{
859  ideal mW=NULL;
860  int i,syzIndex = 0,j=0;
861  intvec * modcomp=NULL,*w=NULL;
862  int ** wv=NULL;
863  tHomog hom=(tHomog)idHomModule(arg,NULL,&w);
864  ring origR = currRing;
865  ring syRing=NULL;
866
867  if ((!isMonomial) && syTestOrder(arg))
868  {
869    WerrorS("sres only implemented for modules with ordering  ..,c or ..,C");
870    return NULL;
871  }
872  *length = 4;
873  resolvente res = (resolvente)omAlloc0(4*sizeof(ideal)),newres;
874  res[0] = idCopy(arg);
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    sy0buck = kBucketCreate();
880    if (syzIndex+1==*length)
881    {
882      newres = (resolvente)omAlloc0((*length+4)*sizeof(ideal));
883      // for (j=0;j<*length+4;j++) newres[j] = NULL;
884      for (j=0;j<*length;j++) newres[j] = res[j];
885      omFreeSize((ADDRESS)res,*length*sizeof(ideal));
886      *length += 4;
887      res=newres;
888    }
889    if ((hom==isHomog)|| (rHasGlobalOrdering(origR)))
890    {
891      if (syzIndex==0) syInitSort(res[0],&modcomp);
892      if ((syzIndex==0) && !rRing_has_CompLastBlock(currRing))
893        res[syzIndex+1] = sySchreyersSyzygiesFB(res[syzIndex],&modcomp,mW,FALSE);
894      else
895        res[syzIndex+1] = sySchreyersSyzygiesFB(res[syzIndex],&modcomp,mW);
896      mW = res[syzIndex];
897    }
898//idPrint(res[syzIndex+1]);
899
900    if ((syzIndex==0))
901    {
902      if ((hom==isHomog)|| (rHasGlobalOrdering(origR)))
903      {
904        syRing = rCurrRingAssure_CompLastBlock();
905        if (syRing != origR)
906        {
907          for (i=0; i<IDELEMS(res[1]); i++)
908          {
909            res[1]->m[i] = prMoveR( res[1]->m[i], origR);
910          }
911        }
912        idTest(res[1]);
913      }
914      else
915      {
916        syRing = rCurrRingAssure_SyzComp_CompLastBlock();
917        if (syRing != origR)
918        {
919          for (i=0; i<IDELEMS(res[0]); i++)
920          {
921            res[0]->m[i] = prMoveR( res[0]->m[i], origR);
922          }
923        }
924        idTest(res[0]);
925      }
926    }
927    if ((hom!=isHomog) && (rHasLocalOrMixedOrdering(origR)))
928    {
929      if (syzIndex==0) syInitSort(res[0],&modcomp);
930      res[syzIndex+1] = sySchreyersSyzygiesFM(res[syzIndex],&modcomp);
931    }
932    syzIndex++;
933    if (TEST_OPT_PROT) Print("[%d]\n",syzIndex);
934    kBucketDestroy(&(sy0buck));
935  }
936  //syPrintResolution(res,1,*length);
937  if ((hom!=isHomog) && (rHasLocalOrMixedOrdering(origR)))
938  {
939    syzIndex = 1;
940    while ((syzIndex < *length) && (!idIs0(res[syzIndex])))
941    {
942      idShift(res[syzIndex],syzIndex);
943      syzIndex++;
944    }
945  }
946  if ((hom==isHomog) || (rHasGlobalOrdering(origR)))
947    syzIndex = 1;
948  else
949    syzIndex = 0;
950  syReOrderResolventFB(res,*length,syzIndex+1);
951  if (/*ringOrderChanged:*/ origR!=syRing && syRing != NULL)
952  {
953    rChangeCurrRing(origR);
954    // Thomas: Here I assume that all (!) polys of res live in tmpR
955    while ((syzIndex < *length) && (res[syzIndex]))
956    {
957      for (i=0;i<IDELEMS(res[syzIndex]);i++)
958      {
959        if (res[syzIndex]->m[i])
960        {
961          res[syzIndex]->m[i] = prMoveR( res[syzIndex]->m[i], syRing);
962        }
963      }
964      syzIndex++;
965    }
966    j = 0;
967    while (currRing->order[j]!=0) j++;
968    rKill(syRing);
969  }
970  else
971  {
972    // Thomas -- are you sure that you have to "reorder" here?
973    while ((syzIndex < *length) && (res[syzIndex]))
974    {
975      for (i=0;i<IDELEMS(res[syzIndex]);i++)
976      {
977        if (res[syzIndex]->m[i])
978          res[syzIndex]->m[i] = pSortCompCorrect(res[syzIndex]->m[i]);
979      }
980      syzIndex++;
981    }
982  }
983  if ((hom==isHomog) || (rHasGlobalOrdering(origR)))
984  {
985    if (res[1]!=NULL)
986    {
987      syReOrderResolventFB(res,2,1);
988      for (i=0;i<IDELEMS(res[1]);i++)
989      {
990        if (res[1]->m[i])
991          res[1]->m[i] = pSort(res[1]->m[i]);
992      }
993    }
994  }
995  //syPrintResolution(res,0,*length);
996
997  //syMergeSortResolventFB(res,*length);
998  if (modcomp!=NULL) delete modcomp;
999  if (w!=NULL) delete w;
1000  return res;
1001}
1002
1003syStrategy sySchreyer(ideal arg, int maxlength)
1004{
1005  int rl;
1006  resolvente fr = sySchreyerResolvente(arg,maxlength,&(rl));
1007  if (fr==NULL) return NULL;
1008
1009  int typ0;
1010  syStrategy result=(syStrategy)omAlloc0(sizeof(ssyStrategy));
1011  result->length=rl;
1012  result->fullres = (resolvente)omAlloc0((rl /*result->length*/+1)*sizeof(ideal));
1013  for (int i=rl /*result->length*/-1;i>=0;i--)
1014  {
1015    if (fr[i]!=NULL)
1016      result->fullres[i] = fr[i];
1017      fr[i] = NULL;
1018  }
1019  if (currQuotient!=NULL)
1020  {
1021    for (int i=0; i<rl; i++)
1022    {
1023      if (result->fullres[i]!=NULL)
1024      {
1025        ideal t=kNF(currQuotient,NULL,result->fullres[i]);
1026        idDelete(&result->fullres[i]);
1027        result->fullres[i]=t;
1028        if (i<rl-1)
1029        {
1030          for(int j=IDELEMS(t)-1;j>=0; j--)
1031          {
1032            if ((t->m[j]==NULL) && (result->fullres[i+1]!=NULL))
1033            {
1034              for(int k=IDELEMS(result->fullres[i+1])-1;k>=0; k--)
1035              {
1036                if (result->fullres[i+1]->m[k]!=NULL)
1037                {
1038                  pDeleteComp(&(result->fullres[i+1]->m[k]),j+1);
1039                }
1040              }
1041            }
1042          }
1043        }
1044        idSkipZeroes(result->fullres[i]);
1045      }
1046    }
1047    if ((rl>maxlength) && (result->fullres[rl-1]!=NULL))
1048    {
1049      idDelete(&result->fullres[rl-1]);
1050    }
1051  }
1052  omFreeSize((ADDRESS)fr,(rl /*result->length*/)*sizeof(ideal));
1053  return result;
1054}
1055
Note: See TracBrowser for help on using the repository browser.