source: git/Singular/syz0.cc @ 63be42

spielwiese
Last change on this file since 63be42 was e0d91c, checked in by Hans Schönemann <hannes@…>, 26 years ago
* hannes: cosmetic changes git-svn-id: file:///usr/local/Singular/svn/trunk@2513 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 22.0 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/* $Id: syz0.cc,v 1.17 1998-09-22 14:09:02 Singular 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 "spolys.h"
18#include "stairc.h"
19#include "ipid.h"
20#include "cntrlc.h"
21#include "ipid.h"
22#include "intvec.h"
23#include "ipshell.h"
24#include "numbers.h"
25#include "ideals.h"
26#include "intvec.h"
27#include "ring.h"
28#include "syz.h"
29
30
31static polyset syInitSort(polyset oldF,int rkF,int Fmax,
32         int syComponentOrder,intvec **modcomp)
33{
34  int i,j,k,kk,kkk,jj;
35  polyset F;
36  int Fl=Fmax;
37
38  while ((Fl!=0) && (oldF[Fl-1]==NULL)) Fl--;
39  if (*modcomp!=NULL) delete modcomp;
40  *modcomp = new intvec(rkF+2);
41  F=(polyset)Alloc0(Fmax*sizeof(poly));
42  j=0;
43  for(i=0;i<=rkF;i++)
44  {
45    k=0;
46    jj = j;
47    (**modcomp)[i] = j;
48    while (k<Fl)
49    {
50      while ((k<Fl) && (pGetComp(oldF[k]) != i)) k++;
51      if (k<Fl)
52      {
53        kk=jj;
54        while ((kk<Fl) && (F[kk]) && (pComp0(oldF[k],F[kk])!=syComponentOrder))
55        {
56            kk++;
57        }
58        for (kkk=j;kkk>kk;kkk--)
59        {
60          F[kkk] = F[kkk-1];
61        }
62        F[kk] = oldF[k];
63//Print("Element %d: ",kk);pWrite(F[kk]);
64        j++;
65        k++;
66      }
67    }
68  }
69  (**modcomp)[rkF+1] = Fl;
70  return F;
71}
72
73static void syCreatePairs(polyset F,int lini,int wend,int k,int j,int i,
74           polyset pairs,int regularPairs=0,ideal mW=NULL)
75{
76  int l,ii=0,jj;
77  poly p,q;
78
79  while (((k<wend) && (pGetComp(F[k]) == i)) ||
80         ((currQuotient!=NULL) && (k<regularPairs+IDELEMS(currQuotient))))
81  {
82    p = pOne();
83    if ((k<wend) && (pGetComp(F[k]) == i) && (k!=j))
84      pLcm(F[j],F[k],p);
85    else if (ii<IDELEMS(currQuotient))
86    {
87      q = pHead(F[j]);
88      if (mW!=NULL)
89      {
90        for(jj=1;jj<=pVariables;jj++)
91          pSetExp(q,jj,pGetExp(q,jj) -pGetExp(mW->m[pGetComp(q)-1],jj));
92        pSetm(q);
93      }
94      pLcm(q,currQuotient->m[ii],p);
95      if (mW!=NULL)
96      {
97        for(jj=1;jj<=pVariables;jj++)
98          pSetExp(p,jj,pGetExp(p,jj) +pGetExp(mW->m[pGetComp(p)-1],jj));
99        pSetm(p);
100      }
101      pDelete(&q);
102      k = regularPairs+ii;
103      ii++;
104    }
105    l=lini;
106    while ((l<k) && ((pairs[l]==NULL) || (!pDivisibleBy(pairs[l],p))))
107    {
108      if ((pairs[l]!=NULL) && (pDivisibleBy(p,pairs[l])))
109        pDelete(&(pairs[l]));
110      l++;
111    }
112    if (l==k)
113    {
114      pSetm(p);
115      pairs[l] = p;
116    }
117    else
118      pDelete(&p);
119    k++;
120  }
121}
122
123static poly syRedtail2(poly p, polyset redWith, intvec *modcomp, spSpolyLoopProc SpolyLoop = NULL)
124{
125  poly h, hn;
126  int hncomp,nxt;
127  int j;
128
129  h = p;
130  hn = pNext(h);
131  while(hn != NULL)
132  {
133    hncomp = pGetComp(hn);
134    j = (*modcomp)[hncomp];
135    nxt = (*modcomp)[hncomp+1];
136    while (j < nxt)
137    {
138      if (pDivisibleBy2(redWith[j], hn))
139      {
140        //if (TEST_OPT_PROT) Print("r");
141        hn = spSpolyRed(redWith[j],hn,NULL, SpolyLoop);
142        if (hn == NULL)
143        {
144          pNext(h) = NULL;
145          return p;
146        }
147        hncomp = pGetComp(hn);
148        j = (*modcomp)[hncomp];
149        nxt = (*modcomp)[hncomp+1];
150      }
151      else
152      {
153        j++;
154      }
155    }
156    h = pNext(h) = hn;
157    hn = pNext(h);
158  }
159  return p;
160}
161
162/*2
163* computes the Schreyer syzygies in the local case
164* input: F, Fmax,  noSort: F is already ordered by: Schreyer-order
165*              (only allocated: Shdl, Smax)
166* output: Shdl, Smax
167*/
168void sySchreyersSyzygiesFM(polyset F,int Fmax,polyset* Shdl,int* Smax,
169  BOOLEAN noSort)
170{
171  int Fl=Fmax;
172  while ((Fl!=0) && (F[Fl-1]==NULL)) Fl--;
173  if (Fl==0) return;
174
175  int i,j,l,k,totalToRed,ecartToRed,kk,bestEcart,totalmax,rkF,
176    Sl=0,smax,tmax,tl;
177  int *ecartS, *ecartT, *totalS,
178    *totalT=NULL, *temp=NULL;
179  intvec * modcomp=NULL;
180  polyset pairs,S,T,ST,oldF;
181  poly p,q,toRed;
182  BOOLEAN notFound = FALSE;
183
184/*-------------initializing the sets--------------------*/
185  ideal idF=(ideal)Alloc(sizeof(ip_sideal));
186  ST=(polyset)Alloc0(Fl*sizeof(poly));
187  S=(polyset)Alloc0(Fl*sizeof(poly));
188  ecartS=(int*)Alloc(Fl*sizeof(int));
189  totalS=(int*)Alloc(Fl*sizeof(int));
190  T=(polyset)Alloc0(2*Fl*sizeof(poly));
191  ecartT=(int*)Alloc(2*Fl*sizeof(int));
192  totalT=(int*)Alloc(2*Fl*sizeof(int));
193  pairs=(polyset)Alloc0(Fl*sizeof(poly));
194
195  smax = Fl;
196  tmax = 2*Fl;
197  idF->m=F;IDELEMS(idF)=Fmax;
198  rkF=idRankFreeModule(idF);
199  Free((ADDRESS)idF,sizeof(ip_sideal));
200  spSet(currRing);
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 = spSpolyCreate(S[j],S[k],NULL, NULL);
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 = spSpolyRed(p,toRed,NULL, NULL);
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, spSpolyLoopProc SpolyLoop=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 = spSpolyRed(currQuotient->m[i],toNorm,NULL, SpolyLoop);
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;
480  intvec *newmodcomp;
481  polyset pairs,oldF,F=*FF;
482  poly p,q,toRed,syz,lastmonom,multWith;
483  ideal idF=(ideal)Alloc(sizeof(*idF)),null;
484  BOOLEAN isNotReduced=TRUE;
485
486  while ((Fl!=0) && (F[Fl-1]==NULL)) Fl--;
487  newmodcomp = new intvec(Fl+2);
488//for (j=0;j<Fl;j++) pWrite(F[j]);
489//PrintLn();
490  if (currQuotient==NULL)
491    pairs=(polyset)Alloc0(Fl*sizeof(poly));
492  else
493    pairs=(polyset)Alloc0((Fl+IDELEMS(currQuotient))*sizeof(poly));
494  idF->m=F;IDELEMS(idF)=Fmax;
495  rkF=idRankFreeModule(idF);
496  null = idInit(1,rkF);
497  Free((ADDRESS)idF,sizeof(*idF));
498  if (noSort)
499  {
500    oldF = *FF;
501    F=syInitSort(*FF,rkF,Fmax,syComponentOrder,modcomp);
502  }
503  else
504  {
505    F = *FF;
506  }
507  for(j=0;j<Fl;j++)
508  {
509    (*newmodcomp)[j+1] = Sl;
510    if (TEST_OPT_PROT)
511    {
512      Print("(%d)",Fl-j);
513      mflush();
514    }
515    i = pGetComp(F[j]);
516    if (syComponentOrder==1)
517    {
518      lini=k=j+1;
519      wend=Fl;
520    }
521    else
522    {
523      lini=k=0;
524      while ((k<j) && (pGetComp(F[k]) != i)) k++;
525      wend=j;
526    }
527    syCreatePairs(F,lini,wend,k,j,i,pairs,Fl,mW);
528    if (currQuotient!=NULL) wend = Fl+IDELEMS(currQuotient);
529    for (k=lini;k<wend;k++)
530    {
531      if (pairs[k]!=NULL)
532      {
533        if (TEST_OPT_PROT)
534        {
535           PrintS(".");
536           mflush();
537        }
538        //begins to construct the syzygy
539        if (k<Fl)
540        {
541          syz = pCopy(pairs[k]);
542          syz->coef = nCopy(F[k]->coef);
543          syz->coef = nNeg(syz->coef);
544          pNext(syz) = pairs[k];
545          lastmonom = pNext(syz);
546          lastmonom->coef = nCopy(F[j]->coef);
547          pSetComp(lastmonom,k+1);
548        }
549        else
550        {
551          syz = pairs[k];
552          syz->coef = nCopy(currQuotient->m[k-Fl]->coef);
553          lastmonom = syz;
554          multWith = pDivide(syz,F[j]);
555          multWith->coef = nCopy(currQuotient->m[k-Fl]->coef);
556        }
557        pSetComp(syz,j+1);
558        pairs[k] = NULL;
559        //the next term of the syzygy
560        //constructs the spoly
561        if (BTEST1(6))
562        {
563          if (k<Fl)
564          {
565            PrintS("pair: ");pWrite0(F[j]);PrintS("  ");pWrite(F[k]);
566          }
567          else
568          {
569            PrintS("pair: ");pWrite0(F[j]);PrintS("  ");pWrite(currQuotient->m[k-Fl]);
570          }
571        }
572        if (k<Fl)
573          toRed = spSpolyCreate(F[j],F[k],NULL, NULL);
574        else
575        {
576          q = pMultT(pCopy(F[j]),multWith);
577          toRed = sySpecNormalize(q,mW, NULL);
578          pDelete(&multWith);
579        }
580        isNotReduced = TRUE;
581        while (toRed!=NULL)
582        {
583          if (BTEST1(6))
584          {
585            PrintS("toRed: ");pWrite(toRed);
586          }
587//          l=0;
588//          while ((l<Fl) && (!pDivisibleBy(F[l],toRed))) l++;
589//          if (l>=Fl)
590          l = (**modcomp)[pGetComp(toRed)+1]-1;
591          kkk = (**modcomp)[pGetComp(toRed)];
592          while ((l>=kkk) && (!pDivisibleBy(F[l],toRed))) l--;
593          if (l<kkk)
594          {
595            if ((currQuotient!=NULL) && (isNotReduced))
596            {
597              toRed = sySpecNormalize(toRed,mW, NULL);
598              isNotReduced = FALSE;
599            }
600            else
601            {
602              //no polynom for reduction
603              WerrorS("ideal not a standardbasis");
604              pDelete(&toRed);
605              pDelete(&syz);
606              for(k=j;k<Fl;k++) pDelete(&(pairs[k]));
607              Free((ADDRESS)pairs,Fl*sizeof(poly));
608              if (noSort)
609              {
610                Free((ADDRESS)F,Fl*sizeof(poly));
611                F = oldF;
612              }
613              for(k=0;k<*Smax;k++) pDelete(&((*Shdl)[k]));
614              return;
615            }
616          }
617          else
618          {
619            //the next monom of the syzygy
620            isNotReduced = TRUE;
621            if (BTEST1(6))
622            {
623              PrintS("reduced with: ");pWrite(F[l]);
624            }
625            multWith = pDivide(toRed,F[l]);
626            multWith->coef = nDiv(toRed->coef,F[l]->coef);
627            multWith->coef = nNeg(multWith->coef);
628            pNext(lastmonom) = toRed;
629            pIter(lastmonom);
630            pIter(toRed);
631            pNext(lastmonom) = NULL;
632            lastmonom->coef = nDiv(lastmonom->coef,F[l]->coef);
633            lastmonom->coef = nNeg(lastmonom->coef);
634            pSetComp(lastmonom,l+1);
635            //computes the new toRed
636            p = pCopy(pNext(F[l]));
637            p = pMultT(p,multWith);
638            pDelete(&multWith);
639            toRed = pAdd(toRed,p);
640            //the module component of the new monom
641//pWrite(toRed);
642          }
643        }
644//PrintLn();
645        if (syz!=NULL)
646        {
647          if (Sl>=*Smax)
648          {
649            pEnlargeSet(Shdl,*Smax,16);
650            *Smax += 16;
651          }
652          pNorm(syz);
653          if (BTEST1(OPT_REDTAIL))
654          {
655            (*newmodcomp)[j+2] = Sl;
656            (*Shdl)[Sl] = syRedtail2(syz,*Shdl,newmodcomp, NULL);
657            (*newmodcomp)[j+2] = 0;
658          }
659          else
660            (*Shdl)[Sl] = syz;
661          Sl++;
662        }
663      }
664    }
665//    for(k=j;k<Fl;k++) pDelete(&(pairs[k]));
666  }
667  (*newmodcomp)[Fl+1] = Sl;
668  if (currQuotient==NULL)
669    Free((ADDRESS)pairs,Fl*sizeof(poly));
670  else
671    Free((ADDRESS)pairs,(Fl+IDELEMS(currQuotient))*sizeof(poly));
672  if (noSort)
673  {
674    Free((ADDRESS)oldF,Fmax*sizeof(poly));
675    *FF = F;
676  }
677  delete *modcomp;
678  *length = Fl+2;
679  *modcomp = newmodcomp;
680}
681
682void syReOrderResolventFB(resolvente res,int length, int initial)
683{
684  int syzIndex=length-1,i,j;
685  poly p;
686
687  while ((syzIndex!=0) && (res[syzIndex]==NULL)) syzIndex--;
688  while (syzIndex>=initial)
689  {
690    for(i=0;i<IDELEMS(res[syzIndex]);i++)
691    {
692      p = res[syzIndex]->m[i];
693      while (p!=NULL)
694      {
695        if (res[syzIndex-1]->m[pGetComp(p)-1]!=NULL)
696        {
697          for(j=1;j<=pVariables;j++)
698          {
699            pSetExp(p,j,pGetExp(p,j)
700                        -pGetExp(res[syzIndex-1]->m[pGetComp(p)-1],j));
701          }
702        }
703        else
704          PrintS("error in the resolvent\n");
705        pSetm(p);
706        pIter(p);
707      }
708    }
709    syzIndex--;
710  }
711}
712
713BOOLEAN syTestOrder(ideal M)
714{
715  int i=idRankFreeModule(M);
716  int j=0;
717
718  while ((currRing->order[j]!=ringorder_c) && (currRing->order[j]!=ringorder_C))
719    j++;
720  if ((i>0) && (currRing->order[j+1]!=0))
721  {
722    return TRUE;
723  }
724  return FALSE;
725}
726
727resolvente sySchreyerResolvente(ideal arg, int maxlength, int * length,
728  BOOLEAN isMonomial,BOOLEAN notReplace)
729{
730  ideal mW=NULL;
731  int i,syzIndex = 0,j=0,lgth,*ord=NULL,*bl0=NULL,*bl1=NULL;
732  intvec * modcomp=NULL,*w=NULL;
733  short ** wv=NULL;
734  BOOLEAN sort = TRUE;
735  tHomog hom=(tHomog)idHomModule(arg,NULL,&w);
736  ring origR = currRing;
737  sip_sring tmpR;
738
739  if ((!isMonomial) && syTestOrder(arg))
740  {
741    WerrorS("sres only implemented for modules with ordering  ..,c or ..,C");
742    return NULL;
743  }
744  *length = 4;
745  resolvente res = (resolvente)Alloc0(4*sizeof(ideal)),newres;
746  res[0] = idCopy(arg);
747  while ((!idIs0(res[syzIndex])) && ((maxlength==-1) || (syzIndex<maxlength)))
748  {
749    i = IDELEMS(res[syzIndex]);
750    //while ((i!=0) && (!res[syzIndex]->m[i-1])) i--;
751    if (syzIndex+1==*length)
752    {
753      newres = (resolvente)Alloc((*length+4)*sizeof(ideal));
754      for (j=0;j<*length+4;j++) newres[j] = NULL;
755      for (j=0;j<*length;j++) newres[j] = res[j];
756      Free((ADDRESS)res,*length*sizeof(ideal));
757      *length += 4;
758      res=newres;
759    }
760    res[syzIndex+1] = idInit(16,1);
761    if ((currRing->OrdSgn == 1) || (hom==isHomog))
762    {
763      sySchreyersSyzygiesFB(&(res[syzIndex]->m),i,&(res[syzIndex+1]->m),
764        &(IDELEMS(res[syzIndex+1])),sort,&modcomp,&lgth,mW);
765      mW = res[syzIndex];
766    }
767    else
768      sySchreyersSyzygiesFM(res[syzIndex]->m,i,&(res[syzIndex+1]->m),
769        &(IDELEMS(res[syzIndex+1])),sort);
770//idPrint(res[syzIndex+1]);
771    if ((syzIndex==0) && (currRing->OrdSgn==1))
772    {
773      j = 0;
774      while ((currRing->order[j]!=ringorder_c)
775              && (currRing->order[j]!=ringorder_C))
776        j++;
777      if ((!notReplace) && (currRing->order[j]!=0))
778      {
779        while (currRing->order[j]!=0) j++;
780        ord = (int*)Alloc0((j+2)*sizeof(int));
781        wv = (short**)Alloc0((j+2)*sizeof(short*));
782        bl0 = (int*)Alloc0((j+2)*sizeof(int));
783        bl1 = (int*)Alloc0((j+2)*sizeof(int));
784        j = 0;
785        while ((currRing->order[j]!=ringorder_c)
786                && (currRing->order[j]!=ringorder_C))
787        {
788          ord[j] = currRing->order[j];
789          bl0[j] = currRing->block0[j];
790          bl1[j] = currRing->block1[j];
791          wv[j] = currRing->wvhdl[j];
792          j++;
793        }
794        int m_order=j;
795        while (currRing->order[j+1]!=0)
796        {
797          ord[j] = currRing->order[j+1];
798          bl0[j] = currRing->block0[j+1];
799          bl1[j] = currRing->block1[j+1];
800          wv[j] = currRing->wvhdl[j+1];
801          j++;
802        }
803        ord[j] = currRing->order[m_order];
804        bl0[j] = currRing->block0[m_order];
805        bl1[j] = currRing->block1[m_order];
806        wv[j] = currRing->wvhdl[m_order];
807        tmpR = *currRing;
808        tmpR.order = ord;
809        tmpR.block0 = bl0;
810        tmpR.block1 = bl1;
811        tmpR.wvhdl = wv;
812        rComplete(&tmpR);
813        rChangeCurrRing(&tmpR, TRUE);
814      }
815    }
816    if (sort) sort=FALSE;
817    syzIndex++;
818    if (TEST_OPT_PROT) Print("[%d]\n",syzIndex);
819  }
820  if (currRing->OrdSgn == -1)
821    pSetSchreyerOrdM(NULL,0,0);
822  syReOrderResolventFB(res,*length);
823  syzIndex = 1;
824  if (/*ringOrderChanged:*/ ord!=NULL)
825  {
826    j = 0;
827    while (currRing->order[j]!=0) j++;
828    Free((ADDRESS)ord,(j+2)*sizeof(int));
829    Free((ADDRESS)bl0,(j+2)*sizeof(int));
830    Free((ADDRESS)bl1,(j+2)*sizeof(int));
831    Free((ADDRESS)wv,(j+2)*sizeof(short*));
832    rChangeCurrRing(origR, TRUE);
833  }
834  while ((syzIndex < *length) && (res[syzIndex]))
835  {
836    for (i=0;i<IDELEMS(res[syzIndex]);i++)
837    {
838      if (res[syzIndex]->m[i])
839        res[syzIndex]->m[i] = pOrdPolyMerge(res[syzIndex]->m[i]);
840    }
841    syzIndex++;
842  }
843  if (modcomp!=NULL) delete modcomp;
844  if (w!=NULL) delete w;
845  return res;
846}
847
848syStrategy sySchreyer(ideal arg, int maxlength)
849{
850  int typ0;
851  syStrategy result=(syStrategy)Alloc0(sizeof(ssyStrategy));
852
853  resolvente fr = sySchreyerResolvente(arg,maxlength,&(result->length));
854  result->fullres = (resolvente)Alloc0((result->length+1)*sizeof(ideal));
855  for (int i=result->length-1;i>=0;i--)
856  {
857    if (fr[i]!=NULL)
858      result->fullres[i] = fr[i];
859      fr[i] = NULL;
860  }
861  Free((ADDRESS)fr,(result->length)*sizeof(ideal));
862  return result;
863}
864
Note: See TracBrowser for help on using the repository browser.