source: git/kernel/GBEngine/syz.cc @ 8204ed

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