source: git/kernel/syz0.cc @ 854405

spielwiese
Last change on this file since 854405 was 762407, checked in by Oleksandr Motsak <motsak@…>, 12 years ago
config.h is for sources files only FIX: config.h should only be used by source (not from inside kernel/mod2.h!) NOTE: each source file should better include mod2.h right after config.h, while headers should better not include mod2.h.
  • Property mode set to 100644
File size: 27.2 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/* $Id$ */
5/*
6* ABSTRACT: resolutions
7*/
8
9
10#include "config.h"
11#include <kernel/mod2.h>
12#include <misc/options.h>
13#include <omalloc/omalloc.h>
14#include <kernel/polys.h>
15#include <kernel/febase.h>
16#include <kernel/kstd1.h>
17#include <kernel/kutil.h>
18#include <kernel/stairc.h>
19//#include "cntrlc.h"
20#include <misc/intvec.h>
21#include <coeffs/numbers.h>
22#include <kernel/ideals.h>
23#include <misc/intvec.h>
24#include <polys/monomials/ring.h>
25#include <kernel/syz.h>
26#include <polys/kbuckets.h>
27#include <polys/prCopy.h>
28
29static void syInitSort(ideal arg,intvec **modcomp)
30{
31  int i,j,k,kk,kkk,jj;
32  idSkipZeroes(arg);
33  polyset F,oldF=arg->m;
34  int Fl=IDELEMS(arg);
35  int rkF=id_RankFreeModule(arg,currRing);
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}
73
74static void syCreatePairs(polyset F,int lini,int wend,int k,int j,int i,
75           polyset pairs,int regularPairs=0,ideal mW=NULL)
76{
77  int l,ii=0,jj;
78  poly p,q;
79
80  while (((k<wend) && (pGetComp(F[k]) == i)) ||
81         ((currQuotient!=NULL) && (k<regularPairs+IDELEMS(currQuotient))))
82  {
83    p = pOne();
84    if ((k<wend) && (pGetComp(F[k]) == i) && (k!=j))
85      pLcm(F[j],F[k],p);
86    else if (ii<IDELEMS(currQuotient))
87    {
88      q = pHead(F[j]);
89      if (mW!=NULL)
90      {
91        for(jj=1;jj<=(currRing->N);jj++)
92          pSetExp(q,jj,pGetExp(q,jj) -pGetExp(mW->m[pGetComp(q)-1],jj));
93        pSetm(q);
94      }
95      pLcm(q,currQuotient->m[ii],p);
96      if (mW!=NULL)
97      {
98        for(jj=1;jj<=(currRing->N);jj++)
99          pSetExp(p,jj,pGetExp(p,jj) +pGetExp(mW->m[pGetComp(p)-1],jj));
100        pSetm(p);
101      }
102      pDelete(&q);
103      k = regularPairs+ii;
104      ii++;
105    }
106    l=lini;
107    while ((l<k) && ((pairs[l]==NULL) || (!pDivisibleBy(pairs[l],p))))
108    {
109      if ((pairs[l]!=NULL) && (pDivisibleBy(p,pairs[l])))
110        pDelete(&(pairs[l]));
111      l++;
112    }
113    if (l==k)
114    {
115      pSetm(p);
116      pairs[l] = p;
117    }
118    else
119      pDelete(&p);
120    k++;
121  }
122}
123
124static poly syRedtail2(poly p, polyset redWith, intvec *modcomp)
125{
126  poly h, hn;
127  int hncomp,nxt;
128  int j;
129
130  h = p;
131  hn = pNext(h);
132  while(hn != NULL)
133  {
134    hncomp = pGetComp(hn);
135    j = (*modcomp)[hncomp];
136    nxt = (*modcomp)[hncomp+1];
137    while (j < nxt)
138    {
139      if (pLmDivisibleByNoComp(redWith[j], hn))
140      {
141        //if (TEST_OPT_PROT) PrintS("r");
142        hn = ksOldSpolyRed(redWith[j],hn);
143        if (hn == NULL)
144        {
145          pNext(h) = NULL;
146          return p;
147        }
148        hncomp = pGetComp(hn);
149        j = (*modcomp)[hncomp];
150        nxt = (*modcomp)[hncomp+1];
151      }
152      else
153      {
154        j++;
155      }
156    }
157    h = pNext(h) = hn;
158    hn = pNext(h);
159  }
160  return p;
161}
162
163/*2
164* computes the Schreyer syzygies in the local case
165* input: arg   (only allocated: Shdl, Smax)
166* output: Shdl, Smax
167*/
168static ideal sySchreyersSyzygiesFM(ideal arg,intvec ** modcomp)
169{
170  int Fl=IDELEMS(arg);
171  while ((Fl!=0) && (arg->m[Fl-1]==NULL)) Fl--;
172  ideal result=idInit(16,arg->rank+Fl);
173  polyset F=arg->m,*Shdl=&(result->m);
174  if (Fl==0) return result;
175
176  int i,j,l,k,totalToRed,ecartToRed,kk,kkk;
177  int bestEcart,totalmax,rkF,Sl=0,smax,tmax,tl;
178  int *ecartS, *ecartT, *totalS,
179    *totalT=NULL, *temp=NULL;
180  polyset pairs,S,T,ST,oldF;
181  poly p,q,toRed;
182  BOOLEAN notFound = FALSE;
183  intvec * newmodcomp = new intvec(Fl+2);
184  intvec * tempcomp;
185
186//Print("Naechster Modul\n");
187//idPrint(arg);
188/*-------------initializing the sets--------------------*/
189  ST=(polyset)omAlloc0(Fl*sizeof(poly));
190  S=(polyset)omAlloc0(Fl*sizeof(poly));
191  ecartS=(int*)omAlloc(Fl*sizeof(int));
192  totalS=(int*)omAlloc(Fl*sizeof(int));
193  T=(polyset)omAlloc0(2*Fl*sizeof(poly));
194  ecartT=(int*)omAlloc(2*Fl*sizeof(int));
195  totalT=(int*)omAlloc(2*Fl*sizeof(int));
196  pairs=(polyset)omAlloc0(Fl*sizeof(poly));
197
198  smax = Fl;
199  tmax = 2*Fl;
200  for (j=1;j<IDELEMS(arg);j++)
201  {
202    if (arg->m[j] != NULL)
203    {
204      assume (arg->m[j-1] != NULL);
205      assume (pGetComp(arg->m[j-1])-pGetComp(arg->m[j])<=0);
206    }
207  }
208  rkF=id_RankFreeModule(arg,currRing);
209/*----------------construction of the new ordering----------*/
210  if (rkF>0)
211    rSetSyzComp(rkF, currRing);
212  else
213    rSetSyzComp(1, currRing);
214/*----------------creating S--------------------------------*/
215  for(j=0;j<Fl;j++)
216  {
217    S[j] = pCopy(F[j]);
218    totalS[j] = p_LDeg(S[j],&k,currRing);
219    ecartS[j] = totalS[j]-p_FDeg(S[j],currRing);
220//Print("%d", pGetComp(S[j]));PrintS("  ");
221    p = S[j];
222    if (rkF==0) pSetCompP(p,1);
223    while (pNext(p)!=NULL) pIter(p);
224    pNext(p) = pHead(F[j]);
225    pIter(p);
226    if (rkF==0)
227      pSetComp(p,j+2);
228    else
229      pSetComp(p,rkF+j+1);
230    pSetmComp(p);
231  }
232//PrintLn();
233  if (rkF==0) rkF = 1;
234/*---------------creating the initial for T----------------*/
235  j=0;
236  l=-1;
237  totalmax=-1;
238  for (k=0;k<smax;k++)
239    if (totalS[k]>totalmax) totalmax=totalS[k];
240  for (kk=1;kk<=rkF;kk++)
241  {
242    for (k=0;k<=totalmax;k++)
243    {
244      for (l=0;l<smax;l++)
245      {
246        if ((pGetComp(S[l])==kk) && (totalS[l]==k))
247        {
248          ST[j] = S[l];
249          totalT[j] = totalS[l];
250          ecartT[j] = ecartS[l];
251//Print("%d", totalS[l]);PrintS("  ");
252          j++;
253        }
254      }
255    }
256  }
257//PrintLn();
258  for (j=0;j<smax;j++)
259  {
260     totalS[j] = totalT[j];
261     ecartS[j] = ecartT[j];
262  }
263
264/*---------------computing---------------------------------*/
265  for(j=0;j<smax;j++)
266  {
267    (*newmodcomp)[j+1] = Sl;
268    i = pGetComp(S[j]);
269    int syComponentOrder= currRing->ComponentOrder;
270    int lini,wend;
271    if (syComponentOrder==1)
272    {
273      lini=k=j+1;
274      wend=Fl;
275    }
276    else
277    {
278      lini=k=0;
279      while ((k<j) && (pGetComp(S[k]) != i)) k++;
280      wend=j;
281    }
282    if (TEST_OPT_PROT)
283    {
284      Print("(%d)",Fl-j);
285      mflush();
286    }
287    syCreatePairs(S,lini,wend,k,j,i,pairs);
288    for (k=lini;k<wend;k++)
289    {
290      if (pairs[k]!=NULL)
291      {
292/*--------------creating T----------------------------------*/
293        for (l=0;l<smax;l++)
294        {
295          ecartT[l] = ecartS[l];
296          totalT[l] = totalS[l];
297          T[l] = ST[l];
298        }
299        tl = smax;
300        tempcomp = ivCopy(*modcomp);
301/*--------------begin to reduce-----------------------------*/
302        toRed = ksOldCreateSpoly(S[j],S[k]);
303        ecartToRed = 1;
304        bestEcart = 1;
305        if (TEST_OPT_DEBUG)
306        {
307          PrintS("pair: ");pWrite0(S[j]);PrintS(" ");pWrite(S[k]);
308        }
309        if (TEST_OPT_PROT)
310        {
311           PrintS(".");
312           mflush();
313        }
314//Print("Reduziere Paar %d,%d (ecart %d): \n",j,k,ecartToRed);
315//Print("Poly %d: ",j);pWrite(S[j]);
316//Print("Poly %d: ",k);pWrite(S[k]);
317//Print("Spoly: ");pWrite(toRed);
318        while (pGetComp(toRed)<=rkF)
319        {
320          if (TEST_OPT_DEBUG)
321          {
322            PrintS("toRed: ");pWrite(toRed);
323          }
324/*
325*         if ((bestEcart) || (ecartToRed!=0))
326*         {
327*/
328            totalToRed = p_LDeg(toRed,&kk,currRing);
329            ecartToRed = totalToRed-p_FDeg(toRed,currRing);
330/*
331*         }
332*/
333//Print("toRed now (neuer ecart %d): ",ecartToRed);pWrite(toRed);
334          notFound = TRUE;
335          bestEcart = 32000;  //a very large integer
336          p = NULL;
337          int l=0;
338#define OLD_SEARCH
339#ifdef OLD_SEARCH
340          while ((l<tl) && (pGetComp(T[l])<pGetComp(toRed))) l++;
341          //assume (l==(**modcomp)[pGetComp(toRed)]);
342          while ((l<tl) && (notFound))
343#else
344          l = (**modcomp)[pGetComp(toRed)];
345          kkk = (**modcomp)[pGetComp(toRed)+1];
346          while ((l<kkk) && (notFound))
347#endif
348          {
349            if ((ecartT[l]<bestEcart) && (pDivisibleBy(T[l],toRed)))
350            {
351              if (ecartT[l]<=ecartToRed) notFound = FALSE;
352              p = T[l];
353              bestEcart = ecartT[l];
354            }
355            l++;
356          }
357          if (p==NULL)
358          {
359            WerrorS("ideal not a standard basis");//no polynom for reduction
360            pDelete(&toRed);
361            for(k=j;k<Fl;k++) pDelete(&(pairs[k]));
362            omFreeSize((ADDRESS)pairs,Fl*sizeof(poly));
363            omFreeSize((ADDRESS)ST,Fl*sizeof(poly));
364            omFreeSize((ADDRESS)S,Fl*sizeof(poly));
365            omFreeSize((ADDRESS)T,tmax*sizeof(poly));
366            omFreeSize((ADDRESS)ecartT,tmax*sizeof(int));
367            omFreeSize((ADDRESS)totalT,tmax*sizeof(int));
368            omFreeSize((ADDRESS)ecartS,Fl*sizeof(int));
369            omFreeSize((ADDRESS)totalS,Fl*sizeof(int));
370            for(k=0;k<IDELEMS(result);k++) pDelete(&((*Shdl)[k]));
371            return result;
372          }
373          else
374          {
375//Print("reduced with (ecart %d): ",bestEcart);wrp(p);PrintLn();
376            if (notFound)
377            {
378              if (tl>=tmax)
379              {
380                pEnlargeSet(&T,tmax,16);
381                tmax += 16;
382                temp = (int*)omAlloc((tmax+16)*sizeof(int));
383                for(l=0;l<tmax;l++) temp[l]=totalT[l];
384                totalT = temp;
385                temp = (int*)omAlloc((tmax+16)*sizeof(int));
386                for(l=0;l<tmax;l++) temp[l]=ecartT[l];
387                ecartT = temp;
388              }
389//PrintS("t");
390              int comptR=pGetComp(toRed);
391              for (l=tempcomp->length()-1;l>comptR;l--)
392              {
393                if ((*tempcomp)[l]>0)
394                  (*tempcomp)[l]++;
395              }
396              l=0;
397              while ((l<tl) && (comptR>pGetComp(T[l]))) l++;
398              while ((l<tl) && (totalT[l]<=totalToRed)) l++;
399              for (kk=tl;kk>l;kk--)
400              {
401                T[kk]=T[kk-1];
402                totalT[kk]=totalT[kk-1];
403                ecartT[kk]=ecartT[kk-1];
404              }
405              q = pCopy(toRed);
406              pNorm(q);
407              T[l] = q;
408              totalT[l] = totalToRed;
409              ecartT[l] = ecartToRed;
410              tl++;
411            }
412            toRed = ksOldSpolyRed(p,toRed);
413          }
414        }
415//Print("toRed finally (neuer ecart %d): ",ecartToRed);pWrite(toRed);
416//PrintS("s");
417        if (pGetComp(toRed)>rkF)
418        {
419          if (Sl>=IDELEMS(result))
420          {
421            pEnlargeSet(Shdl,IDELEMS(result),16);
422            IDELEMS(result) += 16;
423          }
424          //p_Shift(&toRed,-rkF,currRing);
425          pNorm(toRed);
426          (*Shdl)[Sl] = toRed;
427          Sl++;
428        }
429/*----------------deleting all polys not from ST--------------*/
430        for(l=0;l<tl;l++)
431        {
432          kk=0;
433          while ((kk<smax) && (T[l] != S[kk])) kk++;
434          if (kk>=smax)
435          {
436            pDelete(&T[l]);
437//Print ("#");
438          }
439        }
440        delete tempcomp;
441      }
442    }
443    for(k=lini;k<wend;k++) pDelete(&(pairs[k]));
444  }
445  (*newmodcomp)[Fl+1] = Sl;
446  omFreeSize((ADDRESS)pairs,Fl*sizeof(poly));
447  omFreeSize((ADDRESS)ST,Fl*sizeof(poly));
448  omFreeSize((ADDRESS)S,Fl*sizeof(poly));
449  omFreeSize((ADDRESS)T,tmax*sizeof(poly));
450  omFreeSize((ADDRESS)ecartT,tmax*sizeof(int));
451  omFreeSize((ADDRESS)totalT,tmax*sizeof(int));
452  omFreeSize((ADDRESS)ecartS,Fl*sizeof(int));
453  omFreeSize((ADDRESS)totalS,Fl*sizeof(int));
454  delete *modcomp;
455  *modcomp = newmodcomp;
456  return result;
457}
458
459/*3
460*special Normalform for Schreyer in factor rings
461*/
462poly sySpecNormalize(poly toNorm,ideal mW=NULL)
463{
464  int j,i=0;
465  poly p;
466
467  if (toNorm==NULL) return NULL;
468  p = pHead(toNorm);
469  if (mW!=NULL)
470  {
471    for(j=1;j<=(currRing->N);j++)
472      pSetExp(p,j,pGetExp(p,j) -pGetExp(mW->m[pGetComp(p)-1],j));
473  }
474  while ((p!=NULL) && (i<IDELEMS(currQuotient)))
475  {
476    if (pDivisibleBy(currQuotient->m[i],p))
477    {
478      //pNorm(toNorm);
479      toNorm = ksOldSpolyRed(currQuotient->m[i],toNorm);
480      pDelete(&p);
481      if (toNorm==NULL) return NULL;
482      p = pHead(toNorm);
483      if (mW!=NULL)
484      {
485        for(j=1;j<=(currRing->N);j++)
486          pSetExp(p,j,pGetExp(p,j) -pGetExp(mW->m[pGetComp(p)-1],j));
487      }
488      i = 0;
489    }
490    else
491    {
492      i++;
493    }
494  }
495  pDelete(&p);
496  return toNorm;
497}
498
499/*2
500* computes the Schreyer syzygies in the global case
501* input: F
502* output: Shdl, Smax
503* modcomp, length stores the start position of the module comp. in arg
504*/
505static ideal sySchreyersSyzygiesFB(ideal arg,intvec ** modcomp,ideal mW,BOOLEAN redTail=TRUE)
506{
507  kBucket_pt sy0buck = kBucketCreate(currRing);
508
509  int Fl=IDELEMS(arg);
510  while ((Fl!=0) && (arg->m[Fl-1]==NULL)) Fl--;
511  ideal result=idInit(16,Fl);
512  int i,j,l,k,kkk,rkF,Sl=0,syComponentOrder=currRing->ComponentOrder;
513  int fstart,wend,lini,ltR,gencQ=0;
514  intvec *newmodcomp;
515  int *Flength;
516  polyset pairs,F=arg->m,*Shdl=&(result->m);
517  poly p,q,toRed,syz,lastmonom,multWith;
518  BOOLEAN isNotReduced=TRUE;
519
520//#define WRITE_BUCKETS
521#ifdef WRITE_BUCKETS
522  PrintS("Input: \n");
523  ideal twr=idHead(arg);
524  idPrint(arg);
525  idDelete(&twr);
526  if (modcomp!=NULL) (*modcomp)->show(0,0);
527#endif
528   
529  newmodcomp = new intvec(Fl+2);
530  //for (j=0;j<Fl;j++) pWrite(F[j]);
531  //PrintLn();
532  if (currQuotient==NULL)
533    pairs=(polyset)omAlloc0(Fl*sizeof(poly));
534  else
535  {
536    gencQ = IDELEMS(currQuotient);
537    pairs=(polyset)omAlloc0((Fl+gencQ)*sizeof(poly));
538  }
539  rkF=id_RankFreeModule(arg,currRing);
540  Flength = (int*)omAlloc0(Fl*sizeof(int));
541  for(j=0;j<Fl;j++)
542  {
543    Flength[j] = pLength(F[j]);
544  }
545  for(j=0;j<Fl;j++)
546  {
547    (*newmodcomp)[j+1] = Sl;
548    if (TEST_OPT_PROT)
549    {
550      Print("(%d)",Fl-j);
551      mflush();
552    }
553    i = pGetComp(F[j]);
554    if (syComponentOrder==1)
555    {
556      lini=k=j+1;
557      wend=Fl;
558    }
559    else
560    {
561      lini=k=0;
562      while ((k<j) && (pGetComp(F[k]) != i)) k++;
563      wend=j;
564    }
565    syCreatePairs(F,lini,wend,k,j,i,pairs,Fl,mW);
566    if (currQuotient!=NULL) wend = Fl+gencQ;
567    for (k=lini;k<wend;k++)
568    {
569      if (pairs[k]!=NULL)
570      {
571        if (TEST_OPT_PROT)
572        {
573           PrintS(".");
574           mflush();
575        }
576        //begins to construct the syzygy
577        if (k<Fl)
578        {
579          number an=nCopy(pGetCoeff(F[k])),bn=nCopy(pGetCoeff(F[j]));
580          int ct = ksCheckCoeff(&an, &bn, currRing->cf);
581          syz = pCopy(pairs[k]);
582          //syz->coef = nCopy(F[k]->coef);
583          syz->coef = an;
584          //syz->coef = nNeg(syz->coef);
585          pNext(syz) = pairs[k];
586          lastmonom = pNext(syz);
587          //lastmonom->coef = nCopy(F[j]->coef);
588          lastmonom->coef = bn;
589          lastmonom->coef = nNeg(lastmonom->coef);
590          pSetComp(lastmonom,k+1);
591        }
592        else
593        {
594          syz = pairs[k];
595          syz->coef = nCopy(currQuotient->m[k-Fl]->coef);
596          syz->coef = nNeg(syz->coef);
597          lastmonom = syz;
598          multWith = pDivide(syz,F[j]);
599          multWith->coef = nCopy(currQuotient->m[k-Fl]->coef);
600        }
601        pSetComp(syz,j+1);
602        pairs[k] = NULL;
603        //the next term of the syzygy
604        //constructs the spoly
605        if (TEST_OPT_DEBUG)
606        {
607          if (k<Fl)
608          {
609            PrintS("pair: ");pWrite0(F[j]);PrintS("  ");pWrite(F[k]);
610          }
611          else
612          {
613            PrintS("pair: ");pWrite0(F[j]);PrintS("  ");pWrite(currQuotient->m[k-Fl]);
614          }
615        }
616        if (k<Fl)
617          toRed = ksOldCreateSpoly(F[j],F[k]);
618        else
619        {
620          q = pMult_mm(pCopy(F[j]),multWith);
621          toRed = sySpecNormalize(q,mW);
622          pDelete(&multWith);
623        }
624        kBucketInit(sy0buck,toRed,-1);
625        toRed = kBucketGetLm(sy0buck);
626        isNotReduced = TRUE;
627        while (toRed!=NULL)
628        {
629          if (TEST_OPT_DEBUG)
630          {
631            PrintS("toRed: ");pWrite(toRed);
632          }
633//          l=0;
634//          while ((l<Fl) && (!pDivisibleBy(F[l],toRed))) l++;
635//          if (l>=Fl)
636          l = (**modcomp)[pGetComp(toRed)+1]-1;
637          kkk = (**modcomp)[pGetComp(toRed)];
638          while ((l>=kkk) && (!pDivisibleBy(F[l],toRed))) l--;
639#ifdef WRITE_BUCKETS
640          kBucketClear(sy0buck,&toRed,&ltR);
641          printf("toRed in Pair[%d, %d]:", j, k);
642          pWrite(toRed);
643          kBucketInit(sy0buck,toRed,-1);
644#endif
645
646          if (l<kkk)
647          {
648            if ((currQuotient!=NULL) && (isNotReduced))
649            {
650              kBucketClear(sy0buck,&toRed,&ltR);
651              toRed = sySpecNormalize(toRed,mW);
652#ifdef WRITE_BUCKETS
653              printf("toRed in Pair[%d, %d]:", j, k);
654              pWrite(toRed);
655#endif
656              kBucketInit(sy0buck,toRed,-1);
657              toRed = kBucketGetLm(sy0buck);
658              isNotReduced = FALSE;
659            }
660            else
661            {
662              //no polynom for reduction
663              WerrorS("ideal not a standard basis");
664              pDelete(&toRed);
665              pDelete(&syz);
666              for(k=j;k<Fl;k++) pDelete(&(pairs[k]));
667              omFreeSize((ADDRESS)pairs,(Fl + gencQ)*sizeof(poly));
668              for(k=0;k<IDELEMS(result);k++) pDelete(&((*Shdl)[k]));
669               
670              kBucketDestroy(&(sy0buck));
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   
737  kBucketDestroy(&(sy0buck));
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<=(currRing->N);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<=(currRing->N);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<=(currRing->N);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=id_RankFreeModule(M,currRing);
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, currRing);
844  for (i=0;i<IDELEMS(arg);i++)
845  {
846    if (arg->m[i]!=NULL)
847      p_Shift(&arg->m[i],-j,currRing);
848  }
849}
850
851#if 0 /*debug only */
852static void syPrintResolution(resolvente res,int start,int length)
853{
854  while ((start < length) && (res[start]))
855  {
856    Print("Syz(%d): \n",start);
857    idTest(res[start]);
858    //idPrint(res[start]);
859    start++;
860  }
861}
862#endif
863
864resolvente sySchreyerResolvente(ideal arg, int maxlength, int * length,
865                                BOOLEAN isMonomial, BOOLEAN notReplace)
866{
867  ideal mW=NULL;
868  int i,syzIndex = 0,j=0;
869  intvec * modcomp=NULL,*w=NULL;
870  int ** wv=NULL;
871  tHomog hom=(tHomog)idHomModule(arg,NULL,&w);
872  ring origR = currRing;
873  ring syRing = NULL;
874
875  if ((!isMonomial) && syTestOrder(arg))
876  {
877    WerrorS("sres only implemented for modules with ordering  ..,c or ..,C");
878    return NULL;
879  }
880  *length = 4;
881  resolvente res = (resolvente)omAlloc0(4*sizeof(ideal)),newres;
882  res[0] = idCopy(arg);
883
884  while ((!idIs0(res[syzIndex])) && ((maxlength==-1) || (syzIndex<maxlength)))
885  {
886    i = IDELEMS(res[syzIndex]);
887    //while ((i!=0) && (!res[syzIndex]->m[i-1])) i--;
888    if (syzIndex+1==*length)
889    {
890      newres = (resolvente)omAlloc0((*length+4)*sizeof(ideal));
891      // for (j=0;j<*length+4;j++) newres[j] = NULL;
892      for (j=0;j<*length;j++) newres[j] = res[j];
893      omFreeSize((ADDRESS)res,*length*sizeof(ideal));
894      *length += 4;
895      res=newres;
896    }
897     
898    if ((hom==isHomog)|| (rHasGlobalOrdering(origR)))
899    {
900      if (syzIndex==0) syInitSort(res[0],&modcomp);
901       
902      if ((syzIndex==0) && !rRing_has_CompLastBlock(currRing))
903        res[syzIndex+1] = sySchreyersSyzygiesFB(res[syzIndex],&modcomp,mW,FALSE);
904      else
905        res[syzIndex+1] = sySchreyersSyzygiesFB(res[syzIndex],&modcomp,mW);
906       
907      mW = res[syzIndex];
908    }
909//idPrint(res[syzIndex+1]);
910
911    if ((syzIndex==0))
912    {
913      if ((hom==isHomog)|| (rHasGlobalOrdering(origR)))
914      {
915        syRing = rAssure_CompLastBlock(origR, TRUE);
916        if (syRing != origR)
917        {
918          rChangeCurrRing(syRing);
919          for (i=0; i<IDELEMS(res[1]); i++)
920          {
921            res[1]->m[i] = prMoveR( res[1]->m[i], origR, syRing);
922          }
923        }
924        idTest(res[1]);
925      }
926      else
927      {
928        syRing = rAssure_SyzComp_CompLastBlock(origR, TRUE);
929        if (syRing != origR)
930        {
931          rChangeCurrRing(syRing);
932          for (i=0; i<IDELEMS(res[0]); i++)
933          {
934            res[0]->m[i] = prMoveR( res[0]->m[i], origR, syRing);
935          }
936        }
937        idTest(res[0]);
938      }
939    }
940    if ((hom!=isHomog) && (rHasLocalOrMixedOrdering(origR)))
941    {
942      if (syzIndex==0) syInitSort(res[0],&modcomp);
943      res[syzIndex+1] = sySchreyersSyzygiesFM(res[syzIndex],&modcomp);
944    }
945    syzIndex++;
946    if (TEST_OPT_PROT) Print("[%d]\n",syzIndex);
947  }
948  //syPrintResolution(res,1,*length);
949  if ((hom!=isHomog) && (rHasLocalOrMixedOrdering(origR)))
950  {
951    syzIndex = 1;
952    while ((syzIndex < *length) && (!idIs0(res[syzIndex])))
953    {
954      idShift(res[syzIndex],syzIndex);
955      syzIndex++;
956    }
957  }
958  if ((hom==isHomog) || (rHasGlobalOrdering(origR)))
959    syzIndex = 1;
960  else
961    syzIndex = 0;
962  syReOrderResolventFB(res,*length,syzIndex+1);
963  if (/*ringOrderChanged:*/ origR!=syRing && syRing != NULL)
964  {
965    rChangeCurrRing(origR);
966    // Thomas: Here I assume that all (!) polys of res live in tmpR
967    while ((syzIndex < *length) && (res[syzIndex]))
968    {
969      for (i=0;i<IDELEMS(res[syzIndex]);i++)
970      {
971        if (res[syzIndex]->m[i])
972        {
973          res[syzIndex]->m[i] = prMoveR( res[syzIndex]->m[i], syRing, origR);
974        }
975      }
976      syzIndex++;
977    }
978//    j = 0; while (currRing->order[j]!=0) j++; // What was this for???!
979    rDelete(syRing);
980  }
981  else
982  {
983    // Thomas -- are you sure that you have to "reorder" here?
984    while ((syzIndex < *length) && (res[syzIndex]))
985    {
986      for (i=0;i<IDELEMS(res[syzIndex]);i++)
987      {
988        if (res[syzIndex]->m[i])
989          res[syzIndex]->m[i] = pSortCompCorrect(res[syzIndex]->m[i]);
990      }
991      syzIndex++;
992    }
993  }
994  if ((hom==isHomog) || (rHasGlobalOrdering(origR)))
995  {
996    if (res[1]!=NULL)
997    {
998      syReOrderResolventFB(res,2,1);
999      for (i=0;i<IDELEMS(res[1]);i++)
1000      {
1001        if (res[1]->m[i])
1002          res[1]->m[i] = pSort(res[1]->m[i]);
1003      }
1004    }
1005  }
1006  //syPrintResolution(res,0,*length);
1007
1008  //syMergeSortResolventFB(res,*length);
1009  if (modcomp!=NULL) delete modcomp;
1010  if (w!=NULL) delete w;
1011  return res;
1012}
1013
1014syStrategy sySchreyer(ideal arg, int maxlength)
1015{
1016  int rl;
1017  resolvente fr = sySchreyerResolvente(arg,maxlength,&(rl));
1018  if (fr==NULL) return NULL;
1019
1020  int typ0;
1021  syStrategy result=(syStrategy)omAlloc0(sizeof(ssyStrategy));
1022  result->length=rl;
1023  result->fullres = (resolvente)omAlloc0((rl /*result->length*/+1)*sizeof(ideal));
1024  for (int i=rl /*result->length*/-1;i>=0;i--)
1025  {
1026    if (fr[i]!=NULL)
1027      result->fullres[i] = fr[i];
1028      fr[i] = NULL;
1029  }
1030  if (currQuotient!=NULL)
1031  {
1032    for (int i=0; i<rl; i++)
1033    {
1034      if (result->fullres[i]!=NULL)
1035      {
1036        ideal t=kNF(currQuotient,NULL,result->fullres[i]);
1037        idDelete(&result->fullres[i]);
1038        result->fullres[i]=t;
1039        if (i<rl-1)
1040        {
1041          for(int j=IDELEMS(t)-1;j>=0; j--)
1042          {
1043            if ((t->m[j]==NULL) && (result->fullres[i+1]!=NULL))
1044            {
1045              for(int k=IDELEMS(result->fullres[i+1])-1;k>=0; k--)
1046              {
1047                if (result->fullres[i+1]->m[k]!=NULL)
1048                {
1049                  pDeleteComp(&(result->fullres[i+1]->m[k]),j+1);
1050                }
1051              }
1052            }
1053          }
1054        }
1055        idSkipZeroes(result->fullres[i]);
1056      }
1057    }
1058    if ((rl>maxlength) && (result->fullres[rl-1]!=NULL))
1059    {
1060      idDelete(&result->fullres[rl-1]);
1061    }
1062  }
1063  omFreeSize((ADDRESS)fr,(rl /*result->length*/)*sizeof(ideal));
1064  return result;
1065}
1066
Note: See TracBrowser for help on using the repository browser.