source: git/Singular/syz0.cc @ a3bc95e

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