source: git/kernel/syz0.cc @ 2a38d90

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