source: git/Singular/syz.cc @ c5f4b9

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