source: git/Singular/syz0.cc @ c54075

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