source: git/kernel/syz0.cc @ 0f7420d

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