source: git/kernel/syz0.cc @ 68349d

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