source: git/Singular/syz0.cc @ 0fcb82

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