source: git/kernel/GBEngine/syz.cc @ 4f8fd1d

spielwiese
Last change on this file since 4f8fd1d was 4f8fd1d, checked in by Frédéric Chapoton <chapoton@…>, 17 months ago
trying to add a codespell linter for kernel/
  • Property mode set to 100644
File size: 29.6 KB
RevLine 
[35aab3]1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4
5/*
6* ABSTRACT: resolutions
7*/
8
[86ac3fc]9#include "kernel/mod2.h"
10#include "misc/options.h"
11#include "kernel/polys.h"
12#include "kernel/GBEngine/kstd1.h"
13#include "kernel/GBEngine/kutil.h"
14#include "kernel/combinatorics/stairc.h"
15#include "misc/intvec.h"
16#include "coeffs/numbers.h"
17#include "kernel/ideals.h"
18#include "misc/intvec.h"
19#include "polys/monomials/ring.h"
20#include "kernel/GBEngine/syz.h"
21#include "polys/prCopy.h"
[35aab3]22
[86ac3fc]23#include "polys/nc/sca.h"
[35aab3]24
25static intvec * syPrepareModComp(ideal arg,intvec ** w)
26{
27  intvec *w1 = NULL;
28  int i;
29  BOOLEAN isIdeal=FALSE;
30
31  if ((w==NULL) || (*w==NULL)) return w1;
32  int maxxx = (*w)->length();
33  if (maxxx==1)
34  {
35    maxxx = 2;
36    isIdeal = TRUE;
37  }
38  w1 = new intvec(maxxx+IDELEMS(arg));
39  if (!isIdeal)
40  {
41    for (i=0;i<maxxx;i++)
42    {
43      (*w1)[i] = (**w)[i];
44    }
45  }
46  for (i=maxxx;i<maxxx+IDELEMS(arg);i++)
47  {
48    if (arg->m[i-maxxx]!=NULL)
49    {
[b50590]50      (*w1)[i] = p_FDeg(arg->m[i-maxxx],currRing);
[35aab3]51      if (pGetComp(arg->m[i-maxxx])!=0)
52      {
53        (*w1)[i]+=(**w)[pGetComp(arg->m[i-maxxx])-1];
54      }
55    }
56  }
57  delete (*w);
58  *w = new intvec(IDELEMS(arg)+1);
59  for (i=0;i<IDELEMS(arg);i++)
60  {
61     (**w)[i+1] = (*w1)[i+maxxx];
62  }
63  return w1;
64}
65
66static void syDeleteAbove(ideal up, int k)
67{
68  if (up!=NULL)
69  {
70    for (int i=0;i<IDELEMS(up);i++)
71    {
72      if (up->m[i]!=NULL)
73        pDeleteComp(&(up->m[i]),k+1);
74    }
75  }
76}
77
78/*2
79*minimizes the module mod and cancel superfluous syzygies
80*from syz
81*/
[4e1f7a]82static void syMinStep(ideal mod,ideal &syz,BOOLEAN final=FALSE,ideal up=NULL,
[35aab3]83                      tHomog h=isNotHomog)
84{
85  poly Unit1,Unit2,actWith;
86  int len,i,j,ModComp,m,k,l;
87  BOOLEAN searchUnit,existsUnit;
88
89  if (TEST_OPT_PROT) PrintS("m");
90  if ((final) && (h==isHomog))
91  /*minim is TRUE, we are in the module: maxlength, maxlength <>0*/
92  {
[1be9c5b]93    ideal deg0=id_Jet(syz,0,currRing);
94    id_Delete(&syz,currRing);
[4e1f7a]95    idSkipZeroes(deg0);
[35aab3]96    syz=deg0;
97  }
98/*--cancels empty entees and their related components above--*/
99  j = IDELEMS(syz);
100  while ((j>0) && (!syz->m[j-1])) j--;
101  k = 0;
102  while (k<j)
103  {
104    if (syz->m[k]!=NULL)
105      k++;
106    else
107    {
108      if (TEST_OPT_PROT) PrintS(".");
109      for (l=k;l<j-1;l++) syz->m[l] = syz->m[l+1];
110      syz->m[j-1] = NULL;
111      syDeleteAbove(up,k);
112      j--;
113    }
114  }
115/*--searches for syzygies coming from superfluous elements
116* in the module below--*/
117  searchUnit = TRUE;
[b50590]118  int curr_syz_limit = rGetCurrSyzLimit(currRing);
[86ac3fc]119  BOOLEAN bHasGlobalOrdering=rHasGlobalOrdering(currRing);
120  BOOLEAN bField_has_simple_inverse=rField_has_simple_inverse(currRing);
[35aab3]121  while (searchUnit)
122  {
123    i=0;
124    j=IDELEMS(syz);
125    while ((j>0) && (syz->m[j-1]==NULL)) j--;
126    existsUnit = FALSE;
[86ac3fc]127    if (bHasGlobalOrdering)
[35aab3]128    {
129      while ((i<j) && (!existsUnit))
130      {
131        existsUnit = pVectorHasUnitB(syz->m[i],&ModComp);
132        i++;
133      }
134    }
135    else
136    {
137      int I=0;
138      l = 0;
139      len=0;
140      for (i=0;i<IDELEMS(syz);i++)
141      {
142        if (syz->m[i]!=NULL)
143        {
144          pVectorHasUnit(syz->m[i],&m, &l);
145          if ((len==0) ||((l>0) && (l<len)))
146          {
147            len = l;
148            ModComp = m;
149            I = i;
150          }
151        }
152      }
153//Print("Laenge ist: %d\n",len);
154      if (len>0) existsUnit = TRUE;
155      i = I+1;
156    }
157    if (existsUnit)
158    {
159      i--;
160//--takes out the founded syzygy--
161      if (TEST_OPT_PROT) PrintS("f");
162      actWith = syz->m[i];
[86ac3fc]163      if (!bField_has_simple_inverse) p_Cleardenom(actWith, currRing);
[35aab3]164//Print("actWith: ");pWrite(actWith);
165      syz->m[i] = NULL;
166      for (k=i;k<j-1;k++)  syz->m[k] = syz->m[k+1];
167      syz->m[j-1] = NULL;
168      syDeleteAbove(up,i);
169      j--;
170//--makes Gauss alg. for the column ModComp--
171      Unit1 = pTakeOutComp(&(actWith), ModComp);
172//PrintS("actWith now: ");pWrite(actWith);
173//Print("Unit1: ");pWrite(Unit1);
174      k=0;
175//Print("j= %d",j);
176      while (k<j)
177      {
178        if (syz->m[k]!=NULL)
179        {
180          Unit2 = pTakeOutComp(&(syz->m[k]), ModComp);
181//Print("element %d: ",k);pWrite(syz->m[k]);
182//PrintS("Unit2: ");pWrite(Unit2);
[8049ad]183          syz->m[k] = pMult(pCopy(Unit1),syz->m[k]);
[35aab3]184          syz->m[k] = pSub(syz->m[k],
185            pMult(Unit2,pCopy(actWith)));
186          if (syz->m[k]==NULL)
187          {
188            for (l=k;l<j-1;l++)
189              syz->m[l] = syz->m[l+1];
190            syz->m[j-1] = NULL;
191            j--;
192            syDeleteAbove(up,k);
193            k--;
194          }
195        }
196        k++;
197      }
198      pDelete(&actWith);
199      pDelete(&Unit1);
200//--deletes superfluous elements from the module below---
201      pDelete(&(mod->m[ModComp-1 - curr_syz_limit]));
202      for (k=ModComp-1 - curr_syz_limit;k<IDELEMS(mod)-1;k++)
203        mod->m[k] = mod->m[k+1];
204      mod->m[IDELEMS(mod)-1] = NULL;
205    }
206    else
207      searchUnit = FALSE;
208  }
209  if (TEST_OPT_PROT) PrintLn();
210  idSkipZeroes(mod);
211  idSkipZeroes(syz);
212}
213
214/*2
215* make Gauss with the element elnum in the module component ModComp
216* for the generators from - till
217*/
218void syGaussForOne(ideal syz, int elnum, int ModComp,int from,int till)
219{
[0d1a36]220  int /*k,j,i,*/lu;
[35aab3]221  poly unit1,unit2;
222  poly actWith=syz->m[elnum];
223
224  if (from<0) from = 0;
225  if ((till<=0) || (till>IDELEMS(syz))) till = IDELEMS(syz);
226  syz->m[elnum] = NULL;
[3342df]227  if (!rField_has_simple_inverse(currRing)) p_Cleardenom(actWith, currRing);
[35aab3]228/*--makes Gauss alg. for the column ModComp--*/
229  pTakeOutComp(&(actWith), ModComp, &unit1, &lu);
230  while (from<till)
231  {
[3705d5c]232    poly tmp=syz->m[from];
233    if (/*syz->m[from]*/ tmp!=NULL)
[35aab3]234    {
[3705d5c]235      pTakeOutComp(&(tmp), ModComp, &unit2, &lu);
236      tmp = pMult(pCopy(unit1),tmp);
237      syz->m[from] = pSub(tmp,
[b6201e]238        pMult(unit2,pCopy(actWith)));
[35aab3]239    }
240    from++;
241  }
242  pDelete(&actWith);
243  pDelete(&unit1);
244}
245static void syDeleteAbove1(ideal up, int k)
246{
[6909cfb]247  poly p/*,pp*/;
[35aab3]248  if (up!=NULL)
249  {
250    for (int i=0;i<IDELEMS(up);i++)
251    {
252      p = up->m[i];
253      while ((p!=NULL) && (pGetComp(p)==k))
254      {
[fe3266]255        /*
[35aab3]256        pp = pNext(p);
257        pNext(p) = NULL;
258        pDelete(&p);
259        p = pp;
[fe3266]260        */
261        pLmDelete(&p);
[35aab3]262      }
263      up->m[i] = p;
264      if (p!=NULL)
265      {
266        while (pNext(p)!=NULL)
267        {
268          if (pGetComp(pNext(p))==k)
269          {
[fe3266]270            /*
[35aab3]271            pp = pNext(pNext(p));
272            pNext(pNext(p)) = NULL;
273            pDelete(&pNext(p));
274            pNext(p) = pp;
[fe3266]275            */
276            pLmDelete(&pNext(p));
[35aab3]277          }
278          else
279            pIter(p);
280        }
281      }
282    }
283  }
284}
285/*2
286*minimizes the resolution res
287*assumes homogeneous or local case
288*/
289static void syMinStep1(resolvente res, int length)
290{
[d30a399]291  int i,j,k,index=0;
[35aab3]292  poly p;
293  intvec *have_del=NULL,*to_del=NULL;
294
295  while ((index<length) && (res[index]!=NULL))
296  {
[4f8fd1d]297/*---we take out dependent elements from syz---------------------*/
[35aab3]298    if (res[index+1]!=NULL)
299    {
[4e1f7a]300      ideal deg0 = id_Jet(res[index+1],0,currRing);
301      ideal reddeg0 = kInterRedOld(deg0);
[35aab3]302      idDelete(&deg0);
303      have_del = new intvec(IDELEMS(res[index]));
304      for (i=0;i<IDELEMS(reddeg0);i++)
305      {
306        if (reddeg0->m[i]!=NULL)
307        {
308          j = pGetComp(reddeg0->m[i]);
309          pDelete(&(res[index]->m[j-1]));
[fe3266]310          /*res[index]->m[j-1] = NULL;*/
[35aab3]311          (*have_del)[j-1] = 1;
312        }
313      }
314      idDelete(&reddeg0);
315    }
316    if (index>0)
317    {
318/*--- we search for units and perform Gaussian elimination------*/
319      j = to_del->length();
320      while (j>0)
321      {
322        if ((*to_del)[j-1]==1)
323        {
324          k = 0;
325          while (k<IDELEMS(res[index]))
326          {
327            p = res[index]->m[k];
[0338ae9]328            while ((p!=NULL) && ((!pLmIsConstantComp(p)) || (pGetComp(p)!=j)))
[35aab3]329              pIter(p);
[0338ae9]330            if ((p!=NULL) && (pLmIsConstantComp(p)) && (pGetComp(p)==j)) break;
[35aab3]331            k++;
332          }
333          if (k>=IDELEMS(res[index]))
334          {
335            PrintS("out of range\n");
336          }
337          syGaussForOne(res[index],k,j);
338          if (res[index+1]!=NULL)
339            syDeleteAbove1(res[index+1],k+1);
340          (*to_del)[j-1] = 0;
341        }
342        j--;
343      }
344    }
345    if (to_del!=NULL) delete to_del;
346    to_del = have_del;
347    have_del = NULL;
348    index++;
349  }
350  if (TEST_OPT_PROT) PrintLn();
351  syKillEmptyEntres(res,length);
[4e1f7a]352  if (to_del!=NULL) delete to_del;
[35aab3]353}
354
355void syMinimizeResolvente(resolvente res, int length, int first)
356{
357  int syzIndex=first;
358  intvec *dummy;
359
360  if (syzIndex<1) syzIndex=1;
[86ac3fc]361  if ((syzIndex==1) && (!rIsPluralRing(currRing)) && (idHomModule(res[0],currRing->qideal,&dummy)))
[35aab3]362  {
363    syMinStep1(res,length);
364    delete dummy;
365    return;
366  }
367  while ((syzIndex<length-1) && (res[syzIndex]!=NULL) && (res[syzIndex+1]!=NULL))
368  {
369    syMinStep(res[syzIndex-1],res[syzIndex],FALSE,res[syzIndex+1]);
370    syzIndex++;
371  }
372  if (res[syzIndex]!=NULL)
373    syMinStep(res[syzIndex-1],res[syzIndex]);
374  if (!idIs0(res[0]))
375    idMinEmbedding(res[0],TRUE);
376}
377
378/*2
379* resolution of ideal/module arg, <=maxlength steps, (r[0..maxlength])
380*   no limitation in length if maxlength==0
381* input:arg
382*       minim: TRUE means mres cmd, FALSE nres cmd.
383*       if *len!=0: module weights: weights[0]
384*          (and weights is defined:weights[0..len-1]
385*
386* output:resolvente r[0..length-1],
387*        module weights: weights[0..length-1]
388*/
389resolvente syResolvente(ideal arg, int maxlength, int * length,
390                        intvec *** weights, BOOLEAN minim)
391{
[d30a399]392  BITSET save1;
393  SI_SAVE_OPT1(save1);
[35aab3]394  resolvente newres;
395  tHomog hom=isNotHomog;
396  intvec *w = NULL,**tempW;
[7b25fe]397  int i,k,syzIndex = 0,j,rk_arg=si_max(1,(int)id_RankFreeModule(arg,currRing));
[35aab3]398  int Kstd1_OldDeg=Kstd1_deg;
399  BOOLEAN completeMinim;
400  BOOLEAN oldDegBound=TEST_OPT_DEGBOUND;
401  BOOLEAN setRegularity=TRUE;
402  int wlength=*length;
403
404  if (maxlength!=-1) *length = maxlength+1;
405  else              *length = 5;
406  if ((wlength!=0) && (*length!=wlength))
407  {
408    intvec **wtmp = (intvec**)omAlloc0((*length)*sizeof(intvec*));
409    wtmp[0]=(*weights)[0];
410    omFreeSize((ADDRESS)*weights,wlength*sizeof(intvec*));
411    *weights=wtmp;
412  }
[1be9c5b]413  resolvente res = (resolvente)omAlloc0((*length)*sizeof(ideal));
[35aab3]414
415/*--- initialize the syzygy-ring -----------------------------*/
416  ring origR = currRing;
[7e53d9]417  ring syz_ring = rAssure_SyzComp(origR, TRUE); // will do rChangeCurrRing if needed
[23c7a8]418  rSetSyzComp(rk_arg, syz_ring);
[35aab3]419
420  if (syz_ring != origR)
421  {
[1b420fd]422    rChangeCurrRing(syz_ring);
423    res[0] = idrCopyR_NoSort(arg, origR, syz_ring);
[35aab3]424  }
425  else
426  {
427    res[0] = idCopy(arg);
428  }
429
430/*--- creating weights for the module components ---------------*/
[a63c24]431  if ((weights!=NULL) && (*weights!=NULL)&& ((*weights)[0]!=NULL))
432  {
[ac00e2f]433    if (!idTestHomModule(res[0],currRing->qideal,(*weights)[0]))
[a63c24]434    {
[e6a019]435      WarnS("wrong weights given(1):"); (*weights)[0]->show();PrintLn();
[ac00e2f]436      idHomModule(res[0],currRing->qideal,&w);
[a63c24]437      w->show();PrintLn();
438      *weights=NULL;
439    }
440  }
441
[35aab3]442  if ((weights==NULL) || (*weights==NULL) || ((*weights)[0]==NULL))
443  {
[ac00e2f]444    hom=(tHomog)idHomModule(res[0],currRing->qideal,&w);
[35aab3]445    if (hom==isHomog)
446    {
447      *weights = (intvec**)omAlloc0((*length)*sizeof(intvec*));
448      if (w!=NULL) (*weights)[0] = ivCopy(w);
449    }
450  }
451  else
452  {
453    if ((weights!=NULL) && (*weights!=NULL)&& ((*weights)[0]!=NULL))
454    {
455      w = ivCopy((*weights)[0]);
456      hom = isHomog;
457    }
458  }
[b1dfaf]459
460#ifdef HAVE_PLURAL
[d07546]461  if (rIsPluralRing(currRing) && !rIsSCA(currRing) )
[750e70]462  {
[d07546]463// quick solution; need theory to apply homog GB stuff for G-Algebras
464    hom = isNotHomog;
465  }
[b1dfaf]466#endif // HAVE_PLURAL
467
[35aab3]468  if (hom==isHomog)
469  {
470    intvec *w1 = syPrepareModComp(res[0],&w);
471    if (w!=NULL) { delete w;w=NULL; }
472    w = w1;
473    j = 0;
474    while ((j<IDELEMS(res[0])) && (res[0]->m[j]==NULL)) j++;
475    if (j<IDELEMS(res[0]))
476    {
[b50590]477      if (p_FDeg(res[0]->m[j],currRing)!=pTotaldegree(res[0]->m[j]))
[35aab3]478        setRegularity = FALSE;
479    }
480  }
481  else
482  {
483    setRegularity = FALSE;
484  }
485
486/*--- the main loop --------------------------------------*/
[9e8bfa]487  while ((res[syzIndex]!=NULL) && (!idIs0(res[syzIndex])) &&
[35aab3]488         ((maxlength==-1) || (syzIndex<=maxlength)))
489   // (syzIndex<maxlength+(int)minim)))
490/*--- compute one step more for minimizing-----------------*/
491  {
492    if (Kstd1_deg!=0) Kstd1_deg++;
493    if (syzIndex+1==*length)
494    {
495      newres = (resolvente)omAlloc0((*length+5)*sizeof(ideal));
496      tempW = (intvec**)omAlloc0((*length+5)*sizeof(intvec*));
497      for (j=0;j<*length;j++)
498      {
499        newres[j] = res[j];
500        if (*weights!=NULL) tempW[j] = (*weights)[j];
501        /*else              tempW[j] = NULL;*/
502      }
503      omFreeSize((ADDRESS)res,*length*sizeof(ideal));
504      if (*weights != NULL) omFreeSize((ADDRESS)*weights,*length*sizeof(intvec*));
505      *length += 5;
506      res=newres;
507      *weights = tempW;
508    }
509/*--- interreducing first -----------------------------------*/
510    if (syzIndex>0)
511    {
[7b25fe]512      int rkI=id_RankFreeModule(res[syzIndex],currRing);
[23c7a8]513      rSetSyzComp(rkI, currRing);
[35aab3]514    }
[750e70]515    if(! TEST_OPT_NO_SYZ_MINIM )
[35aab3]516    if (minim || (syzIndex!=0))
517    {
[1be9c5b]518      ideal temp = kInterRedOld(res[syzIndex],currRing->qideal);
[35aab3]519      idDelete(&res[syzIndex]);
520      idSkipZeroes(temp);
521      res[syzIndex] = temp;
522    }
523/*--- computing the syzygy modules --------------------------------*/
[ac00e2f]524    if ((currRing->qideal==NULL)&&(syzIndex==0)&& (!TEST_OPT_DEGBOUND))
[35aab3]525    {
526      res[/*syzIndex+*/1] = idSyzygies(res[0/*syzIndex*/],hom,&w,FALSE,setRegularity,&Kstd1_deg);
[7069f2]527      if ((!TEST_OPT_NOTREGULARITY) && (Kstd1_deg>0)
[2505ad]528      && (!rField_is_Ring(currRing))
[7069f2]529      ) si_opt_1 |= Sy_bit(OPT_DEGBOUND);
[35aab3]530    }
531    else
532    {
533        res[syzIndex+1] = idSyzygies(res[syzIndex],hom,&w,FALSE);
534    }
535    completeMinim=(syzIndex!=maxlength) || (maxlength ==-1) || (hom!=isHomog);
536    syzIndex++;
537    if (TEST_OPT_PROT) Print("[%d]\n",syzIndex);
[750e70]538
539    if(! TEST_OPT_NO_SYZ_MINIM )
[35aab3]540    {
[750e70]541      if ((minim)||(syzIndex>1))
542        syMinStep(res[syzIndex-1],res[syzIndex],!completeMinim,NULL,hom);
543      if (!completeMinim)
544      /*minim is TRUE, we are in the module: maxlength, maxlength <>0*/
545      {
546        idDelete(&res[syzIndex]);
547      }
[35aab3]548    }
549/*---creating the iterated weights for module components ---------*/
[9e8bfa]550    if ((hom == isHomog) && (res[syzIndex]!=NULL) && (!idIs0(res[syzIndex])))
[35aab3]551    {
552//Print("die %d Modulegewichte sind:\n",w1->length());
553//w1->show();
554//PrintLn();
[7b25fe]555      int max_comp = id_RankFreeModule(res[syzIndex],currRing);
[b50590]556      k = max_comp - rGetCurrSyzLimit(currRing);
[35aab3]557      assume(w != NULL);
558      if (w != NULL)
559        w->resize(max_comp+IDELEMS(res[syzIndex]));
560      else
561        w = new intvec(max_comp+IDELEMS(res[syzIndex]));
562      (*weights)[syzIndex] = new intvec(k);
563      for (i=0;i<k;i++)
564      {
565        if (res[syzIndex-1]->m[i]!=NULL) // hs
566        {
[b50590]567          (*w)[i + rGetCurrSyzLimit(currRing)] = p_FDeg(res[syzIndex-1]->m[i],currRing);
[35aab3]568          if (pGetComp(res[syzIndex-1]->m[i])>0)
[b50590]569            (*w)[i + rGetCurrSyzLimit(currRing)]
[35aab3]570              += (*w)[pGetComp(res[syzIndex-1]->m[i])-1];
[b50590]571          (*((*weights)[syzIndex]))[i] = (*w)[i+rGetCurrSyzLimit(currRing)];
[35aab3]572        }
573      }
574      for (i=k;i<k+IDELEMS(res[syzIndex]);i++)
575      {
576        if (res[syzIndex]->m[i-k]!=NULL)
[b50590]577          (*w)[i+rGetCurrSyzLimit(currRing)] = p_FDeg(res[syzIndex]->m[i-k],currRing)
[35aab3]578                    +(*w)[pGetComp(res[syzIndex]->m[i-k])-1];
579      }
580    }
581  }
582/*--- end of the main loop --------------------------------------*/
583/*--- deleting the temporare data structures --------------------*/
584  if ((syzIndex!=0) && (res[syzIndex]!=NULL) && (idIs0(res[syzIndex])))
585    idDelete(&res[syzIndex]);
586  if (w !=NULL) delete w;
587
588  Kstd1_deg=Kstd1_OldDeg;
589  if (!oldDegBound)
[d30a399]590    si_opt_1 &= ~Sy_bit(OPT_DEGBOUND);
[35aab3]591
592  for (i=1; i<=syzIndex; i++)
593  {
[9e8bfa]594    if ((res[i]!=NULL) && ! idIs0(res[i]))
[35aab3]595    {
[741464]596      id_Shift(res[i],-rGetMaxSyzComp(i, currRing),currRing);
[35aab3]597    }
598  }
599/*--- going back to the original ring -------------------------*/
600  if (origR != syz_ring)
601  {
[1b420fd]602    rChangeCurrRing(origR); // should not be needed now?
[35aab3]603    for (i=0; i<=syzIndex; i++)
604    {
[1b420fd]605      res[i] = idrMoveR_NoSort(res[i], syz_ring, origR);
[35aab3]606    }
[5fe834]607    rDelete(syz_ring);
[35aab3]608  }
[d30a399]609  SI_RESTORE_OPT1(save1);
[35aab3]610  return res;
611}
612
613syStrategy syResolution(ideal arg, int maxlength,intvec * w, BOOLEAN minim)
614{
[8b09ef]615
[90adb6a]616#ifdef HAVE_PLURAL
[57bfa2]617  const ideal idSaveCurrRingQuotient = currRing->qideal;
[90adb6a]618  if( rIsSCA(currRing) )
619  {
[57bfa2]620    if( ncExtensions(TESTSYZSCAMASK) )
621    {
[ac00e2f]622      currRing->qideal = SCAQuotient(currRing);
[57bfa2]623    }
[90adb6a]624    const unsigned int m_iFirstAltVar = scaFirstAltVar(currRing);
625    const unsigned int m_iLastAltVar  = scaLastAltVar(currRing);
[e2efe91]626    arg = id_KillSquares(arg, m_iFirstAltVar, m_iLastAltVar, currRing, false); // kill suares in input!
[90adb6a]627  }
628#endif
[750e70]629
[35aab3]630  syStrategy result=(syStrategy)omAlloc0(sizeof(ssyStrategy));
631
[ac00e2f]632  if ((w!=NULL) && (!idTestHomModule(arg,currRing->qideal,w))) // is this right in SCA case???
[a63c24]633  {
[e6a019]634    WarnS("wrong weights given(2):");w->show();PrintLn();
[ac00e2f]635    idHomModule(arg,currRing->qideal,&w);
[a63c24]636    w->show();PrintLn();
637    w=NULL;
638  }
[35aab3]639  if (w!=NULL)
640  {
[835d83]641    result->weights = (intvec**)omAlloc0Bin(char_ptr_bin);
[35aab3]642    (result->weights)[0] = ivCopy(w);
643    result->length = 1;
644  }
[1b8e4b1]645  resolvente fr = syResolvente(arg,maxlength,&(result->length),&(result->weights),minim);
646  resolvente fr1;
[35aab3]647  if (minim)
648  {
649    result->minres = (resolvente)omAlloc0((result->length+1)*sizeof(ideal));
650    fr1 =  result->minres;
651  }
652  else
653  {
654    result->fullres = (resolvente)omAlloc0((result->length+1)*sizeof(ideal));
655    fr1 =  result->fullres;
656  }
657  for (int i=result->length-1;i>=0;i--)
658  {
659    if (fr[i]!=NULL)
660      fr1[i] = fr[i];
661    fr[i] = NULL;
662  }
663  omFreeSize((ADDRESS)fr,(result->length)*sizeof(ideal));
[90adb6a]664
665#ifdef HAVE_PLURAL
666  if( rIsSCA(currRing) )
667  {
[57bfa2]668    if( ncExtensions(TESTSYZSCAMASK) )
669    {
670      currRing->qideal = idSaveCurrRingQuotient;
671    }
[90adb6a]672    id_Delete(&arg, currRing);
673  }
674#endif
675
[35aab3]676  return result;
677}
678
679static poly sypCopyConstant(poly inp)
680{
681  poly outp=NULL,q;
682
683  while (inp!=NULL)
684  {
685    if (pLmIsConstantComp(inp))
686    {
687      if (outp==NULL)
688      {
689        q = outp = pHead(inp);
690      }
691      else
692      {
693        pNext(q) = pHead(inp);
694        pIter(q);
695      }
696    }
697    pIter(inp);
698  }
699  return outp;
700}
701int syDetect(ideal id,int index,BOOLEAN homog,int * degrees,int * tocancel)
702{
[d30a399]703  int i, j, k, subFromRank=0;
[35aab3]704  ideal temp;
705
706  if (idIs0(id)) return 0;
707  temp = idInit(IDELEMS(id),id->rank);
708  for (i=0;i<IDELEMS(id);i++)
709  {
710    temp->m[i] = sypCopyConstant(id->m[i]);
711  }
712  i = IDELEMS(id);
713  while ((i>0) && (temp->m[i-1]==NULL)) i--;
714  if (i==0)
715  {
716    idDelete(&temp);
717    return 0;
718  }
719  j = 0;
720  while ((j<i) && (temp->m[j]==NULL)) j++;
721  while (j<i)
722  {
723    if (homog)
724    {
[b50590]725      if (index==0) k = p_FDeg(temp->m[j],currRing)+degrees[pGetComp(temp->m[j])];
[b130fb]726      else          k = degrees[pGetComp(temp->m[j])];
[35aab3]727      if (k>=index) tocancel[k-index]++;
728      if ((k>=0) && (index==0)) subFromRank++;
729    }
730    else
731    {
732      tocancel[0]--;
733    }
734    syGaussForOne(temp,j,pGetComp(temp->m[j]),j+1,i);
735    j++;
736    while ((j<i) && (temp->m[j]==NULL)) j++;
737  }
738  idDelete(&temp);
739  return subFromRank;
740}
741
742void syDetect(ideal id,int index,int rsmin, BOOLEAN homog,
743              intvec * degrees,intvec * tocancel)
744{
745  int * deg=NULL;
746  int * tocan=(int*) omAlloc0(tocancel->length()*sizeof(int));
747  int i;
748
749  if (homog)
750  {
751    deg = (int*) omAlloc0(degrees->length()*sizeof(int));
752    for (i=degrees->length();i>0;i--)
753      deg[i-1] = (*degrees)[i-1]-rsmin;
754  }
[d30a399]755  syDetect(id,index,homog,deg,tocan);
[35aab3]756  for (i=tocancel->length();i>0;i--)
757    (*tocancel)[i-1] = tocan[i-1];
758  if (homog)
759    omFreeSize((ADDRESS)deg,degrees->length()*sizeof(int));
760  omFreeSize((ADDRESS)tocan,tocancel->length()*sizeof(int));
761}
762
763/*2
764* computes the betti numbers from a given resolution
765* of length 'length' (0..length-1), not necessairily minimal,
766* (if weights are given, they are used)
767* returns the int matrix of betti numbers
768* and the regularity
769*/
770intvec * syBetti(resolvente res,int length, int * regularity,
771                 intvec* weights,BOOLEAN tomin,int * row_shift)
772{
773//#define BETTI_WITH_ZEROS
774  //tomin = FALSE;
775  int i,j=0,k=0,l,rows,cols,mr;
776  int *temp1,*temp2,*temp3;/*used to compute degrees*/
777  int *tocancel; /*(BOOLEAN)tocancel[i]=element is superfluous*/
778  int r0_len;
779
780  /*------ compute size --------------*/
781  *regularity = -1;
782  cols = length;
783  while ((cols>0)
784  && ((res[cols-1]==NULL)
785    || (idIs0(res[cols-1]))))
786  {
787    cols--;
788  }
789  intvec * result;
790  if (idIs0(res[0]))
791  {
792    if (res[0]==NULL)
793      result = new intvec(1,1,1);
794    else
795      result = new intvec(1,1,res[0]->rank);
796    return result;
797  }
798  intvec *w=NULL;
[a63c24]799  if (weights!=NULL)
800  {
[ac00e2f]801    if (!idTestHomModule(res[0],currRing->qideal,weights))
[a63c24]802    {
[e6a019]803      WarnS("wrong weights given(3):");weights->show();PrintLn();
[ac00e2f]804      idHomModule(res[0],currRing->qideal,&w);
[a63c24]805      if (w!=NULL) { w->show();PrintLn();}
806      weights=NULL;
807    }
808  }
[8ebd1a0]809#if 0
[ac00e2f]810  if (idHomModule(res[0],currRing->qideal,&w)!=isHomog)
[8ebd1a0]811  {
[aad4ca4]812    WarnS("betti-command: Input is not homogeneous!");
[8ebd1a0]813    weights=NULL;
814  }
815#endif
[a63c24]816  if (weights==NULL) weights=w;
817  else delete w;
818  r0_len=IDELEMS(res[0]);
819  while ((r0_len>0) && (res[0]->m[r0_len-1]==NULL)) r0_len--;
820  #ifdef SHOW_W
[192b522]821  PrintS("weights:");if (weights!=NULL) weights->show(); else Print("NULL"); PrintLn();
[a63c24]822  #endif
[7b25fe]823  int rkl=l = si_max(id_RankFreeModule(res[0],currRing),res[0]->rank);
[35aab3]824  i = 0;
825  while ((i<length) && (res[i]!=NULL))
826  {
827    if (IDELEMS(res[i])>l) l = IDELEMS(res[i]);
828    i++;
829  }
830  temp1 = (int*)omAlloc0((l+1)*sizeof(int));
831  temp2 = (int*)omAlloc((l+1)*sizeof(int));
832  rows = 1;
833  mr = 1;
834  cols++;
835  for (i=0;i<cols-1;i++)
836  {
[b44f5e]837    if ((i==0) && (weights!=NULL)) p_SetModDeg(weights, currRing);
[35aab3]838    memset(temp2,0,(l+1)*sizeof(int));
839    for (j=0;j<IDELEMS(res[i]);j++)
840    {
841      if (res[i]->m[j]!=NULL)
842      {
[0338ae9]843        if ((pGetComp(res[i]->m[j])>l)
[1c1d660]844        // usual resolutions do not the following, but artifulal built may: (tr. #763)
845        //|| ((i>1) && (res[i-1]->m[pGetComp(res[i]->m[j])-1]==NULL))
846        )
[35aab3]847        {
[df0a9d3]848          WerrorS("input not a resolution");
[35aab3]849          omFreeSize((ADDRESS)temp1,(l+1)*sizeof(int));
850          omFreeSize((ADDRESS)temp2,(l+1)*sizeof(int));
851          return NULL;
852        }
[b50590]853        temp2[j+1] = p_FDeg(res[i]->m[j],currRing)+temp1[pGetComp(res[i]->m[j])];
[35aab3]854        if (temp2[j+1]-i>rows) rows = temp2[j+1]-i;
855        if (temp2[j+1]-i<mr) mr = temp2[j+1]-i;
856      }
857    }
[b44f5e]858    if ((i==0) && (weights!=NULL)) p_SetModDeg(NULL, currRing);
[35aab3]859    temp3 = temp1;
860    temp1 = temp2;
861    temp2 = temp3;
862  }
863  mr--;
[3e1bc9]864  if (weights!=NULL)
865  {
866    for(j=0;j<weights->length();j++)
867    {
868      if (rows <(*weights)[j]+1) rows=(-mr)+(*weights)[j]+1;
869    }
870  }
[35aab3]871  /*------ computation betti numbers --------------*/
872  rows -= mr;
[a63c24]873  result = new intvec(rows+1,cols,0);
874  if (weights!=NULL)
875  {
876    for(j=0;j<weights->length();j++)
877    {
878      IMATELEM((*result),(-mr)+(*weights)[j]+1,1) ++;
[192b522]879      //Print("imat(%d,%d)++ -> %d\n",(-mr)+(*weights)[j]+1, 1, IMATELEM((*result),(-mr)+(*weights)[j]+1,1));
[a63c24]880    }
881  }
882  else
883  {
884    (*result)[(-mr)*cols] = /*idRankFreeModule(res[0])*/ rkl;
885    if ((!idIs0(res[0])) && ((*result)[(-mr)*cols]==0))
886      (*result)[(-mr)*cols] = 1;
887  }
[35aab3]888  tocancel = (int*)omAlloc0((rows+1)*sizeof(int));
889  memset(temp1,0,(l+1)*sizeof(int));
[a63c24]890  if (weights!=NULL)
891  {
892    memset(temp2,0,l*sizeof(int));
[b44f5e]893    p_SetModDeg(weights, currRing);
[a63c24]894  }
895  else
896    memset(temp2,0,l*sizeof(int));
[d30a399]897  syDetect(res[0],0,TRUE,temp2,tocancel);
[b44f5e]898  if (weights!=NULL) p_SetModDeg(NULL, currRing);
[35aab3]899  if (tomin)
[a63c24]900  {
901    //(*result)[(-mr)*cols] -= dummy;
902    for(j=0;j<=rows+mr;j++)
903    {
[192b522]904      //Print("tocancel[%d]=%d imat(%d,%d)=%d\n",j,tocancel[j],(-mr)+j+1,1,IMATELEM((*result),(-mr)+j+1,1));
[a63c24]905      IMATELEM((*result),(-mr)+j+1,1) -= tocancel[j];
906    }
907  }
[35aab3]908  for (i=0;i<cols-1;i++)
909  {
[b44f5e]910    if ((i==0) && (weights!=NULL)) p_SetModDeg(weights, currRing);
[35aab3]911    memset(temp2,0,l*sizeof(int));
912    for (j=0;j<IDELEMS(res[i]);j++)
913    {
914      if (res[i]->m[j]!=NULL)
915      {
[b50590]916        temp2[j+1] = p_FDeg(res[i]->m[j],currRing)+temp1[pGetComp(res[i]->m[j])];
[35aab3]917        //(*result)[i+1+(temp2[j+1]-i-1)*cols]++;
918        //if (temp2[j+1]>i) IMATELEM((*result),temp2[j+1]-i-mr,i+2)++;
919        IMATELEM((*result),temp2[j+1]-i-mr,i+2)++;
920      }
921      else if (i==0)
922      {
923        if (j<r0_len) IMATELEM((*result),-mr,2)++;
924      }
925    }
926  /*------ computation betti numbers, if res not minimal --------------*/
927    if (tomin)
928    {
929      for (j=mr;j<rows+mr;j++)
930      {
931        //(*result)[i+1+j*cols] -= tocancel[j+1];
932        IMATELEM((*result),j+1-mr,i+2) -= tocancel[j+1];
933      }
934      if ((i<length-1) && (res[i+1]!=NULL))
935      {
936        memset(tocancel,0,(rows+1)*sizeof(int));
[d30a399]937        syDetect(res[i+1],i+1,TRUE,temp2,tocancel);
[35aab3]938        for (j=0;j<rows;j++)
939        {
940          //(*result)[i+1+j*cols] -= tocancel[j];
941          IMATELEM((*result),j+1,i+2) -= tocancel[j];
942        }
943      }
944    }
945    temp3 = temp1;
946    temp1 = temp2;
947    temp2 = temp3;
[a63c24]948    for (j=0;j<=rows;j++)
[35aab3]949    {
950    //  if (((*result)[i+1+j*cols]!=0) && (j>*regularity)) *regularity = j;
951      if ((IMATELEM((*result),j+1,i+2)!=0) && (j>*regularity)) *regularity = j;
952    }
[b44f5e]953    if ((i==0) && (weights!=NULL)) p_SetModDeg(NULL, currRing);
[35aab3]954  }
[a63c24]955  // Print("nach minim:\n"); result->show(); PrintLn();
[35aab3]956  /*------ clean up --------------*/
957  omFreeSize((ADDRESS)tocancel,(rows+1)*sizeof(int));
958  omFreeSize((ADDRESS)temp1,(l+1)*sizeof(int));
959  omFreeSize((ADDRESS)temp2,(l+1)*sizeof(int));
960  if ((tomin) && (mr<0))  // deletes the first (zero) line
961  {
[a63c24]962    for (j=1;j<=rows+mr+1;j++)
[35aab3]963    {
964      for (k=1;k<=cols;k++)
965      {
966        IMATELEM((*result),j,k) = IMATELEM((*result),j-mr,k);
967      }
968    }
[a63c24]969    for (j=rows+mr+1;j<=rows+1;j++)
[35aab3]970    {
971      for (k=1;k<=cols;k++)
972      {
973        IMATELEM((*result),j,k) = 0;
974      }
975    }
976  }
977  j = 0;
978  k = 0;
[a63c24]979  for (i=1;i<=result->rows();i++)
[35aab3]980  {
[a63c24]981    for(l=1;l<=result->cols();l++)
982    if (IMATELEM((*result),i,l) != 0)
[35aab3]983    {
[a63c24]984      j = si_max(j, i-1);
985      k = si_max(k, l-1);
[35aab3]986    }
987  }
988  intvec * exactresult=new intvec(j+1,k+1,0);
989  for (i=0;i<exactresult->rows();i++)
990  {
991    for (j=0;j<exactresult->cols();j++)
992    {
993      IMATELEM(*exactresult,i+1,j+1) = IMATELEM(*result,i+1,j+1);
994    }
995  }
996  if (row_shift!=NULL) *row_shift = mr;
997  delete result;
998  return exactresult;
999}
1000
1001/*2
1002* minbare via syzygies
1003*/
1004ideal syMinBase(ideal arg)
1005{
1006  intvec ** weights=NULL;
1007  int leng;
1008  if (idIs0(arg)) return idInit(1,arg->rank);
1009  resolvente res=syResolvente(arg,1,&leng,&weights,TRUE);
1010  ideal result=res[0];
1011  omFreeSize((ADDRESS)res,leng*sizeof(ideal));
1012  if (weights!=NULL)
1013  {
1014    if (*weights!=NULL)
1015    {
1016      delete (*weights);
1017      *weights=NULL;
1018    }
1019    if ((leng>=1) && (*(weights+1)!=NULL))
1020    {
1021      delete *(weights+1);
1022      *(weights+1)=NULL;
1023    }
1024  }
1025  idSkipZeroes(result);
1026  return result;
1027}
1028
[f815bd7]1029#if 0  /* currently used: syBetti */
[35aab3]1030/*2
1031* computes Betti-numbers from a resolvente of
1032* (non-)homogeneous objects
1033* the numbers of entrees !=NULL in res and weights must be equal
1034* and < length
1035*/
1036intvec * syNewBetti(resolvente res, intvec ** weights, int length)
1037{
1038  intvec * result,*tocancel;
1039  int i,j,k,rsmin=0,rsmax=0,rs=0;
1040  BOOLEAN homog=TRUE;
1041
1042  if (weights!=NULL)           //---homogeneous Betti numbers
1043  {
1044/*--------------computes size of the field----------------------*/
1045    for (i=1;i<length;i++)
1046    {
1047      if (weights[i] != NULL)
1048      {
1049        for (j=1;j<(weights[i])->length();j++)
1050        {
1051          if ((*(weights[i]))[j]-i<rsmin) rsmin = (*(weights[i]))[j]-i;
1052          if ((*(weights[i]))[j]-i>rsmax) rsmax = (*(weights[i]))[j]-i;
1053        }
1054      }
1055    }
1056    i = 0;
1057    while (weights[i] != NULL) i++;
1058    i--;
1059    for (j=0;j<IDELEMS(res[i]);j++)
1060    {
1061      if (res[i]->m[j]!=NULL)
1062      {
[b50590]1063        k = p_FDeg(res[i]->m[j],currRing)+(*(weights[i]))[pGetComp(res[i]->m[j])]-i-1;
[35aab3]1064        if (k>rsmax) rsmax = k;
1065        if (k<rsmin) rsmin = k;
1066      }
1067    }
1068    for (j=1;j<(weights[0])->length();j++)
1069    {
1070      if ((*weights[0])[j]>rsmax) rsmax = (*weights[0])[j];
1071      if ((*weights[0])[j]<rsmin) rsmin = (*weights[0])[j];
1072    }
1073//Print("rsmax = %d\n",rsmax);
1074//Print("rsmin = %d\n",rsmin);
1075    rs = rsmax-rsmin+1;
1076    result = new intvec(rs,i+2,0);
1077    tocancel = new intvec(rs);
1078/*-----------enter the Betti numbers-------------------------------*/
1079    if (/*idRankFreeModule(res[0])*/ res[0]->rank==0)
1080    {
1081      IMATELEM(*result,1-rsmin,1)=1;
1082    }
1083    else
1084    {
1085      for (i=1;i<(weights[0])->length();i++)
1086        IMATELEM(*result,(*weights[0])[i]+1-rsmin,1)++;
1087    }
1088    i = 1;
1089    while (weights[i]!=NULL)
1090    {
1091      for (j=1;j<(weights[i])->length();j++)
1092      {
1093        IMATELEM(*result,(*(weights[i]))[j]-i+1-rsmin,i+1)++;
1094      }
1095      i++;
1096    }
1097    i--;
1098    for (j=0;j<IDELEMS(res[i]);j++)
1099    {
1100      if (res[i]->m[j]!=NULL)
1101      {
[b50590]1102        k = p_FDeg(res[i]->m[j],currRing)+(*(weights[i]))[pGetComp(res[i]->m[j])]-i;
[35aab3]1103        IMATELEM(*result,k-rsmin,i+2)++;
1104      }
1105    }
1106  }
1107  else                //-----the non-homgeneous case
1108  {
1109    homog = FALSE;
1110    tocancel = new intvec(1);
1111    k = length;
1112    while ((k>0) && (idIs0(res[k-1]))) k--;
1113    result = new intvec(1,k+1,0);
1114    (*result)[0] = res[0]->rank;
1115    for (i=0;i<length;i++)
1116    {
1117      if (res[i]!=NULL)
1118      {
1119        for (j=0;j<IDELEMS(res[i]);j++)
1120        {
1121          if (res[i]->m[j]!=NULL) (*result)[i+1]++;
1122        }
1123      }
1124    }
1125  }
1126/*--------computes the Betti numbers for the minimized reolvente----*/
1127
1128  i = 1;
1129  while ((res[i]!=NULL) && (weights[i]!=NULL))
1130  {
1131    syDetect(res[i],i,rsmin,homog,weights[i],tocancel);
1132    if (homog)
1133    {
1134      for (j=0;j<rs-1;j++)
1135      {
1136        IMATELEM((*result),j+1,i+1) -= (*tocancel)[j];
1137        IMATELEM((*result),j+1,i+2) -= (*tocancel)[j+1];
1138      }
1139      IMATELEM((*result),rs,i+1) -= (*tocancel)[rs-1];
1140    }
1141    else
1142    {
1143      (*result)[i+1] -= (*tocancel)[0];
1144      (*result)[i+2] -= (*tocancel)[0];
1145    }
1146    i++;
1147  }
1148
1149/*--------print Betti numbers for control---------------------------*/
1150  for(i=rsmin;i<=rsmax;i++)
1151  {
1152    Print("%2d:",i);
1153    for(j=1;j<=result->cols();j++)
1154    {
1155      Print(" %5d",IMATELEM(*result,i-rsmin+1,j));
1156    }
1157    PrintLn();
1158  }
1159  return result;
1160}
[f815bd7]1161#endif
1162
Note: See TracBrowser for help on using the repository browser.