source: git/kernel/syz.cc @ 62dd9b

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