source: git/kernel/syz.cc @ 8c5988

spielwiese
Last change on this file since 8c5988 was 192b522, checked in by Hans Schönemann <hannes@…>, 19 years ago
*hannes: idHomModule test fails git-svn-id: file:///usr/local/Singular/svn/trunk@7786 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 28.6 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/* $Id: syz.cc,v 1.8 2005-03-17 14:14:17 Singular 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],currRing);
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(pCopy(unit1),syz->m[from]);
250      syz->m[from] = pSub(syz->m[from],
251        pMult(unit2,pCopy(actWith)));
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        /*
269        pp = pNext(p);
270        pNext(p) = NULL;
271        pDelete(&p);
272        p = pp;
273        */
274        pLmDelete(&p);
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            /*
284            pp = pNext(pNext(p));
285            pNext(pNext(p)) = NULL;
286            pDelete(&pNext(p));
287            pNext(p) = pp;
288            */
289            pLmDelete(&pNext(p));
290          }
291          else
292            pIter(p);
293        }
294      }
295    }
296  }
297}
298/*2
299*minimizes the resolution res
300*assumes homogeneous or local case
301*/
302static void syMinStep1(resolvente res, int length)
303{
304  int i,j,k,l,index=0;
305  poly p;
306  ideal deg0=NULL,reddeg0=NULL;
307  intvec *have_del=NULL,*to_del=NULL;
308
309  while ((index<length) && (res[index]!=NULL))
310  {
311/*---we take out dependend elements from syz---------------------*/
312    if (res[index+1]!=NULL)
313    {
314      deg0 = idJet(res[index+1],0);
315      reddeg0 = kInterRed(deg0);
316      idDelete(&deg0);
317      have_del = new intvec(IDELEMS(res[index]));
318      for (i=0;i<IDELEMS(reddeg0);i++)
319      {
320        if (reddeg0->m[i]!=NULL)
321        {
322          j = pGetComp(reddeg0->m[i]);
323          pDelete(&(res[index]->m[j-1]));
324          /*res[index]->m[j-1] = NULL;*/
325          (*have_del)[j-1] = 1;
326        }
327      }
328      idDelete(&reddeg0);
329    }
330    if (index>0)
331    {
332/*--- we search for units and perform Gaussian elimination------*/
333      j = to_del->length();
334      while (j>0)
335      {
336        if ((*to_del)[j-1]==1)
337        {
338          k = 0;
339          while (k<IDELEMS(res[index]))
340          {
341            p = res[index]->m[k];
342            while ((p!=NULL) && ((!pLmIsConstantComp(p)) || (pGetComp(p)!=j)))
343              pIter(p);
344            if ((p!=NULL) && (pLmIsConstantComp(p)) && (pGetComp(p)==j)) break;
345            k++;
346          }
347          if (k>=IDELEMS(res[index]))
348          {
349            PrintS("out of range\n");
350          }
351          syGaussForOne(res[index],k,j);
352          if (res[index+1]!=NULL)
353            syDeleteAbove1(res[index+1],k+1);
354          (*to_del)[j-1] = 0;
355        }
356        j--;
357      }
358    }
359    if (to_del!=NULL) delete to_del;
360    to_del = have_del;
361    have_del = NULL;
362    index++;
363  }
364  if (TEST_OPT_PROT) PrintLn();
365  syKillEmptyEntres(res,length);
366}
367
368void syMinimizeResolvente(resolvente res, int length, int first)
369{
370  int syzIndex=first;
371  intvec *dummy;
372
373  if (syzIndex<1) syzIndex=1;
374  if ((syzIndex==1) && (idHomModule(res[0],currQuotient,&dummy)) && (!rIsPluralRing(currRing)))
375  {
376    syMinStep1(res,length);
377    delete dummy;
378    return;
379  }
380  while ((syzIndex<length-1) && (res[syzIndex]!=NULL) && (res[syzIndex+1]!=NULL))
381  {
382    syMinStep(res[syzIndex-1],res[syzIndex],FALSE,res[syzIndex+1]);
383    syzIndex++;
384  }
385  if (res[syzIndex]!=NULL)
386    syMinStep(res[syzIndex-1],res[syzIndex]);
387  if (!idIs0(res[0]))
388    idMinEmbedding(res[0],TRUE);
389}
390
391/*2
392* resolution of ideal/module arg, <=maxlength steps, (r[0..maxlength])
393*   no limitation in length if maxlength==0
394* input:arg
395*       minim: TRUE means mres cmd, FALSE nres cmd.
396*       if *len!=0: module weights: weights[0]
397*          (and weights is defined:weights[0..len-1]
398*
399* output:resolvente r[0..length-1],
400*        module weights: weights[0..length-1]
401*/
402resolvente syResolvente(ideal arg, int maxlength, int * length,
403                        intvec *** weights, BOOLEAN minim)
404{
405  resolvente res;
406  resolvente newres;
407  tHomog hom=isNotHomog;
408  ideal temp=NULL;
409  intvec *w = NULL,**tempW;
410  int i,k,syzIndex = 0,j,rk_arg=si_max(1,(int)idRankFreeModule(arg));
411  int Kstd1_OldDeg=Kstd1_deg;
412  BOOLEAN completeMinim;
413  BOOLEAN oldDegBound=TEST_OPT_DEGBOUND;
414  BOOLEAN setRegularity=TRUE;
415  int wlength=*length;
416
417  if (maxlength!=-1) *length = maxlength+1;
418  else              *length = 5;
419  if ((wlength!=0) && (*length!=wlength))
420  {
421    intvec **wtmp = (intvec**)omAlloc0((*length)*sizeof(intvec*));
422    wtmp[0]=(*weights)[0];
423    omFreeSize((ADDRESS)*weights,wlength*sizeof(intvec*));
424    *weights=wtmp;
425  }
426  res = (resolvente)omAlloc0((*length)*sizeof(ideal));
427
428/*--- initialize the syzygy-ring -----------------------------*/
429  ring origR = currRing;
430  ring syz_ring = rCurrRingAssure_SyzComp();
431  rSetSyzComp(rk_arg);
432
433  if (syz_ring != origR)
434  {
435    res[0] = idrCopyR_NoSort(arg, origR);
436  }
437  else
438  {
439    res[0] = idCopy(arg);
440  }
441
442/*--- creating weights for the module components ---------------*/
443  if ((weights!=NULL) && (*weights!=NULL)&& ((*weights)[0]!=NULL))
444  {
445    if (!idTestHomModule(res[0],currQuotient,(*weights)[0]))
446    {
447      WarnS("1wrong weights given:"); (*weights)[0]->show();PrintLn();
448      idHomModule(res[0],currQuotient,&w);
449      w->show();PrintLn();
450      *weights=NULL;
451    }
452  }
453
454  if ((weights==NULL) || (*weights==NULL) || ((*weights)[0]==NULL))
455  {
456    hom=(tHomog)idHomModule(res[0],currQuotient,&w);
457    if (hom==isHomog)
458    {
459      *weights = (intvec**)omAlloc0((*length)*sizeof(intvec*));
460      if (w!=NULL) (*weights)[0] = ivCopy(w);
461    }
462  }
463  else
464  {
465    if ((weights!=NULL) && (*weights!=NULL)&& ((*weights)[0]!=NULL))
466    {
467      w = ivCopy((*weights)[0]);
468      hom = isHomog;
469    }
470  }
471  if (hom==isHomog)
472  {
473    intvec *w1 = syPrepareModComp(res[0],&w);
474    if (w!=NULL) { delete w;w=NULL; }
475    w = w1;
476    j = 0;
477    while ((j<IDELEMS(res[0])) && (res[0]->m[j]==NULL)) j++;
478    if (j<IDELEMS(res[0]))
479    {
480      if (pFDeg(res[0]->m[j],currRing)!=pTotaldegree(res[0]->m[j]))
481        setRegularity = FALSE;
482    }
483  }
484  else
485  {
486    setRegularity = FALSE;
487  }
488
489/*--- the main loop --------------------------------------*/
490  while ((!idIs0(res[syzIndex])) &&
491         ((maxlength==-1) || (syzIndex<=maxlength)))
492   // (syzIndex<maxlength+(int)minim)))
493/*--- compute one step more for minimizing-----------------*/
494  {
495    if (Kstd1_deg!=0) Kstd1_deg++;
496    if (syzIndex+1==*length)
497    {
498      newres = (resolvente)omAlloc0((*length+5)*sizeof(ideal));
499      tempW = (intvec**)omAlloc0((*length+5)*sizeof(intvec*));
500      for (j=0;j<*length;j++)
501      {
502        newres[j] = res[j];
503        if (*weights!=NULL) tempW[j] = (*weights)[j];
504        /*else              tempW[j] = NULL;*/
505      }
506      omFreeSize((ADDRESS)res,*length*sizeof(ideal));
507      if (*weights != NULL) omFreeSize((ADDRESS)*weights,*length*sizeof(intvec*));
508      *length += 5;
509      res=newres;
510      *weights = tempW;
511    }
512/*--- interreducing first -----------------------------------*/
513    if (syzIndex>0)
514    {
515      int rkI=idRankFreeModule(res[syzIndex]);
516      rSetSyzComp(rkI);
517    }
518    if (minim || (syzIndex!=0))
519    {
520      temp = kInterRed(res[syzIndex],currQuotient);
521      idDelete(&res[syzIndex]);
522      idSkipZeroes(temp);
523      res[syzIndex] = temp;
524    }
525    temp = NULL;
526/*--- computing the syzygy modules --------------------------------*/
527    if ((currQuotient==NULL)&&(syzIndex==0)&& (!TEST_OPT_DEGBOUND))
528    {
529      res[/*syzIndex+*/1] = idSyzygies(res[0/*syzIndex*/],hom,&w,FALSE,setRegularity,&Kstd1_deg);
530      if ((!TEST_OPT_NOTREGULARITY) && (Kstd1_deg>0)) test |= Sy_bit(OPT_DEGBOUND);
531    }
532    else
533    {
534        res[syzIndex+1] = idSyzygies(res[syzIndex],hom,&w,FALSE);
535    }
536    completeMinim=(syzIndex!=maxlength) || (maxlength ==-1) || (hom!=isHomog);
537    syzIndex++;
538    if (TEST_OPT_PROT) Print("[%d]\n",syzIndex);
539    if ((minim)||(syzIndex>1))
540      syMinStep(res[syzIndex-1],res[syzIndex],!completeMinim,NULL,hom);
541    if (!completeMinim)
542    /*minim is TRUE, we are in the module: maxlength, maxlength <>0*/
543    {
544      idDelete(&res[syzIndex]);
545    }
546/*---creating the iterated weights for module components ---------*/
547    if ((hom == isHomog) && (!idIs0(res[syzIndex])))
548    {
549//Print("die %d Modulegewichte sind:\n",w1->length());
550//w1->show();
551//PrintLn();
552      int max_comp = idRankFreeModule(res[syzIndex]);
553      k = max_comp - rGetCurrSyzLimit();
554      assume(w != NULL);
555      if (w != NULL)
556        w->resize(max_comp+IDELEMS(res[syzIndex]));
557      else
558        w = new intvec(max_comp+IDELEMS(res[syzIndex]));
559      (*weights)[syzIndex] = new intvec(k);
560      for (i=0;i<k;i++)
561      {
562        if (res[syzIndex-1]->m[i]!=NULL) // hs
563        {
564          (*w)[i + rGetCurrSyzLimit()] = pFDeg(res[syzIndex-1]->m[i],currRing);
565          if (pGetComp(res[syzIndex-1]->m[i])>0)
566            (*w)[i + rGetCurrSyzLimit()]
567              += (*w)[pGetComp(res[syzIndex-1]->m[i])-1];
568          (*((*weights)[syzIndex]))[i] = (*w)[i+rGetCurrSyzLimit()];
569        }
570      }
571      for (i=k;i<k+IDELEMS(res[syzIndex]);i++)
572      {
573        if (res[syzIndex]->m[i-k]!=NULL)
574          (*w)[i+rGetCurrSyzLimit()] = pFDeg(res[syzIndex]->m[i-k],currRing)
575                    +(*w)[pGetComp(res[syzIndex]->m[i-k])-1];
576      }
577    }
578  }
579/*--- end of the main loop --------------------------------------*/
580/*--- deleting the temporare data structures --------------------*/
581  if ((syzIndex!=0) && (res[syzIndex]!=NULL) && (idIs0(res[syzIndex])))
582    idDelete(&res[syzIndex]);
583  if (w !=NULL) delete w;
584
585  Kstd1_deg=Kstd1_OldDeg;
586  if (!oldDegBound)
587    test &= ~Sy_bit(OPT_DEGBOUND);
588
589  for (i=1; i<=syzIndex; i++)
590  {
591    if (! idIs0(res[i]))
592    {
593      for (j=0; j<IDELEMS(res[i]); j++)
594      {
595        pShift(&res[i]->m[j], -rGetMaxSyzComp(i));
596      }
597    }
598  }
599/*--- going back to the original ring -------------------------*/
600  if (origR != syz_ring)
601  {
602    rChangeCurrRing(origR);
603    for (i=0; i<=syzIndex; i++)
604    {
605      res[i] = idrMoveR_NoSort(res[i], syz_ring);
606    }
607    rKill(syz_ring);
608  }
609  return res;
610}
611
612syStrategy syResolution(ideal arg, int maxlength,intvec * w, BOOLEAN minim)
613{
614  int typ0;
615  syStrategy result=(syStrategy)omAlloc0(sizeof(ssyStrategy));
616
617  if ((w!=NULL) && (!idTestHomModule(arg,currQuotient,w)))
618  {
619    WarnS("wrong weights given");w->show();PrintLn();
620    idHomModule(arg,currQuotient,&w);
621    w->show();PrintLn();
622    w=NULL;
623  }
624  if (w!=NULL)
625  {
626    result->weights = (intvec**)omAlloc0Bin(void_ptr_bin);
627    (result->weights)[0] = ivCopy(w);
628    result->length = 1;
629  }
630  resolvente fr = syResolvente(arg,maxlength,&(result->length),&(result->weights),minim),fr1;
631  if (minim)
632  {
633    result->minres = (resolvente)omAlloc0((result->length+1)*sizeof(ideal));
634    fr1 =  result->minres;
635  }
636  else
637  {
638    result->fullres = (resolvente)omAlloc0((result->length+1)*sizeof(ideal));
639    fr1 =  result->fullres;
640  }
641  for (int i=result->length-1;i>=0;i--)
642  {
643    if (fr[i]!=NULL)
644      fr1[i] = fr[i];
645    fr[i] = NULL;
646  }
647  omFreeSize((ADDRESS)fr,(result->length)*sizeof(ideal));
648  return result;
649}
650
651static poly sypCopyConstant(poly inp)
652{
653  poly outp=NULL,q;
654
655  while (inp!=NULL)
656  {
657    if (pLmIsConstantComp(inp))
658    {
659      if (outp==NULL)
660      {
661        q = outp = pHead(inp);
662      }
663      else
664      {
665        pNext(q) = pHead(inp);
666        pIter(q);
667      }
668    }
669    pIter(inp);
670  }
671  return outp;
672}
673int syDetect(ideal id,int index,BOOLEAN homog,int * degrees,int * tocancel)
674{
675  int i, j, k, ModComp,subFromRank=0, lu;
676  poly p, q, qq, Unit1, Unit2;
677  ideal temp;
678
679  if (idIs0(id)) return 0;
680  temp = idInit(IDELEMS(id),id->rank);
681  for (i=0;i<IDELEMS(id);i++)
682  {
683    temp->m[i] = sypCopyConstant(id->m[i]);
684  }
685  i = IDELEMS(id);
686  while ((i>0) && (temp->m[i-1]==NULL)) i--;
687  if (i==0)
688  {
689    idDelete(&temp);
690    return 0;
691  }
692  j = 0;
693  p = NULL;
694  while ((j<i) && (temp->m[j]==NULL)) j++;
695  while (j<i)
696  {
697    if (homog)
698    {
699      if (index==0) k = pFDeg(temp->m[j],currRing)+degrees[pGetComp(temp->m[j])];
700      else          k = degrees[pGetComp(temp->m[j])];
701      if (k>=index) tocancel[k-index]++;
702      if ((k>=0) && (index==0)) subFromRank++;
703    }
704    else
705    {
706      tocancel[0]--;
707    }
708    syGaussForOne(temp,j,pGetComp(temp->m[j]),j+1,i);
709    j++;
710    while ((j<i) && (temp->m[j]==NULL)) j++;
711  }
712  idDelete(&temp);
713  return subFromRank;
714}
715
716void syDetect(ideal id,int index,int rsmin, BOOLEAN homog,
717              intvec * degrees,intvec * tocancel)
718{
719  int * deg=NULL;
720  int * tocan=(int*) omAlloc0(tocancel->length()*sizeof(int));
721  int i;
722
723  if (homog)
724  {
725    deg = (int*) omAlloc0(degrees->length()*sizeof(int));
726    for (i=degrees->length();i>0;i--)
727      deg[i-1] = (*degrees)[i-1]-rsmin;
728  }
729  int dummy=syDetect(id,index,homog,deg,tocan);
730  for (i=tocancel->length();i>0;i--)
731    (*tocancel)[i-1] = tocan[i-1];
732  if (homog)
733    omFreeSize((ADDRESS)deg,degrees->length()*sizeof(int));
734  omFreeSize((ADDRESS)tocan,tocancel->length()*sizeof(int));
735}
736
737/*2
738* computes the betti numbers from a given resolution
739* of length 'length' (0..length-1), not necessairily minimal,
740* (if weights are given, they are used)
741* returns the int matrix of betti numbers
742* and the regularity
743*/
744intvec * syBetti(resolvente res,int length, int * regularity,
745                 intvec* weights,BOOLEAN tomin,int * row_shift)
746{
747//#define BETTI_WITH_ZEROS
748  //tomin = FALSE;
749  int i,j=0,k=0,l,rows,cols,mr;
750  int *temp1,*temp2,*temp3;/*used to compute degrees*/
751  int *tocancel; /*(BOOLEAN)tocancel[i]=element is superfluous*/
752  int r0_len;
753
754  /*------ compute size --------------*/
755  *regularity = -1;
756  cols = length;
757  while ((cols>0)
758  && ((res[cols-1]==NULL)
759    || (idIs0(res[cols-1]))))
760  {
761    cols--;
762  }
763  intvec * result;
764  if (idIs0(res[0]))
765  {
766    if (res[0]==NULL)
767      result = new intvec(1,1,1);
768    else
769      result = new intvec(1,1,res[0]->rank);
770    return result;
771  }
772  intvec *w=NULL;
773  if (weights!=NULL)
774  {
775    if (!idTestHomModule(res[0],currQuotient,weights))
776    {
777      WarnS("wrong weights given");weights->show();PrintLn();
778      idHomModule(res[0],currQuotient,&w);
779      if (w!=NULL) { w->show();PrintLn();}
780      weights=NULL;
781    }
782  }
783  //if (idHomModule(res[0],currQuotient,&w)!=isHomog)
784  //{
785  //  Warn("betti-command: Input is not homogeneous!");
786  //  weights=NULL;
787  //}
788  if (weights==NULL) weights=w;
789  else delete w;
790  r0_len=IDELEMS(res[0]);
791  while ((r0_len>0) && (res[0]->m[r0_len-1]==NULL)) r0_len--;
792  #ifdef SHOW_W
793  PrintS("weights:");if (weights!=NULL) weights->show(); else Print("NULL"); PrintLn();
794  #endif
795  int rkl=l = si_max(idRankFreeModule(res[0]),res[0]->rank);
796  i = 0;
797  while ((i<length) && (res[i]!=NULL))
798  {
799    if (IDELEMS(res[i])>l) l = IDELEMS(res[i]);
800    i++;
801  }
802  temp1 = (int*)omAlloc0((l+1)*sizeof(int));
803  temp2 = (int*)omAlloc((l+1)*sizeof(int));
804  rows = 1;
805  mr = 1;
806  cols++;
807  for (i=0;i<cols-1;i++)
808  {
809    if ((i==0) && (weights!=NULL)) pSetModDeg(weights);
810    memset(temp2,0,(l+1)*sizeof(int));
811    for (j=0;j<IDELEMS(res[i]);j++)
812    {
813      if (res[i]->m[j]!=NULL)
814      {
815        if ((pGetComp(res[i]->m[j])>l)
816        || ((i>1) && (res[i-1]->m[pGetComp(res[i]->m[j])-1]==NULL)))
817        {
818          WerrorS("input not a resolvent");
819          omFreeSize((ADDRESS)temp1,(l+1)*sizeof(int));
820          omFreeSize((ADDRESS)temp2,(l+1)*sizeof(int));
821          return NULL;
822        }
823        temp2[j+1] = pFDeg(res[i]->m[j],currRing)+temp1[pGetComp(res[i]->m[j])];
824        if (temp2[j+1]-i>rows) rows = temp2[j+1]-i;
825        if (temp2[j+1]-i<mr) mr = temp2[j+1]-i;
826      }
827    }
828    if ((i==0) && (weights!=NULL)) pSetModDeg(NULL);
829    temp3 = temp1;
830    temp1 = temp2;
831    temp2 = temp3;
832  }
833  mr--;
834  if (weights!=NULL)
835  {
836    for(j=0;j<weights->length();j++)
837    {
838      if (rows <(*weights)[j]+1) rows=(-mr)+(*weights)[j]+1;
839    }
840  }
841  /*------ computation betti numbers --------------*/
842  rows -= mr;
843  result = new intvec(rows+1,cols,0);
844  if (weights!=NULL)
845  {
846    for(j=0;j<weights->length();j++)
847    {
848      IMATELEM((*result),(-mr)+(*weights)[j]+1,1) ++;
849      //Print("imat(%d,%d)++ -> %d\n",(-mr)+(*weights)[j]+1, 1, IMATELEM((*result),(-mr)+(*weights)[j]+1,1));
850    }
851  }
852  else
853  {
854    (*result)[(-mr)*cols] = /*idRankFreeModule(res[0])*/ rkl;
855    if ((!idIs0(res[0])) && ((*result)[(-mr)*cols]==0))
856      (*result)[(-mr)*cols] = 1;
857  }
858  tocancel = (int*)omAlloc0((rows+1)*sizeof(int));
859  memset(temp1,0,(l+1)*sizeof(int));
860  if (weights!=NULL)
861  {
862    memset(temp2,0,l*sizeof(int));
863    pSetModDeg(weights);
864  }
865  else
866    memset(temp2,0,l*sizeof(int));
867  int dummy = syDetect(res[0],0,TRUE,temp2,tocancel);
868  if (weights!=NULL) pSetModDeg(NULL);
869  if (tomin)
870  {
871    //(*result)[(-mr)*cols] -= dummy;
872    for(j=0;j<=rows+mr;j++)
873    {
874      //Print("tocancel[%d]=%d imat(%d,%d)=%d\n",j,tocancel[j],(-mr)+j+1,1,IMATELEM((*result),(-mr)+j+1,1));
875      IMATELEM((*result),(-mr)+j+1,1) -= tocancel[j];
876    }
877  }
878  for (i=0;i<cols-1;i++)
879  {
880    if ((i==0) && (weights!=NULL)) pSetModDeg(weights);
881    memset(temp2,0,l*sizeof(int));
882    for (j=0;j<IDELEMS(res[i]);j++)
883    {
884      if (res[i]->m[j]!=NULL)
885      {
886        temp2[j+1] = pFDeg(res[i]->m[j],currRing)+temp1[pGetComp(res[i]->m[j])];
887        //(*result)[i+1+(temp2[j+1]-i-1)*cols]++;
888        //if (temp2[j+1]>i) IMATELEM((*result),temp2[j+1]-i-mr,i+2)++;
889        IMATELEM((*result),temp2[j+1]-i-mr,i+2)++;
890      }
891      else if (i==0)
892      {
893        if (j<r0_len) IMATELEM((*result),-mr,2)++;
894      }
895    }
896  /*------ computation betti numbers, if res not minimal --------------*/
897    if (tomin)
898    {
899      for (j=mr;j<rows+mr;j++)
900      {
901        //(*result)[i+1+j*cols] -= tocancel[j+1];
902        IMATELEM((*result),j+1-mr,i+2) -= tocancel[j+1];
903      }
904      if ((i<length-1) && (res[i+1]!=NULL))
905      {
906        memset(tocancel,0,(rows+1)*sizeof(int));
907        dummy = syDetect(res[i+1],i+1,TRUE,temp2,tocancel);
908        for (j=0;j<rows;j++)
909        {
910          //(*result)[i+1+j*cols] -= tocancel[j];
911          IMATELEM((*result),j+1,i+2) -= tocancel[j];
912        }
913      }
914    }
915    temp3 = temp1;
916    temp1 = temp2;
917    temp2 = temp3;
918    for (j=0;j<=rows;j++)
919    {
920    //  if (((*result)[i+1+j*cols]!=0) && (j>*regularity)) *regularity = j;
921      if ((IMATELEM((*result),j+1,i+2)!=0) && (j>*regularity)) *regularity = j;
922    }
923    if ((i==0) && (weights!=NULL)) pSetModDeg(NULL);
924  }
925  // Print("nach minim:\n"); result->show(); PrintLn();
926  /*------ clean up --------------*/
927  omFreeSize((ADDRESS)tocancel,(rows+1)*sizeof(int));
928  omFreeSize((ADDRESS)temp1,(l+1)*sizeof(int));
929  omFreeSize((ADDRESS)temp2,(l+1)*sizeof(int));
930  if ((tomin) && (mr<0))  // deletes the first (zero) line
931  {
932    for (j=1;j<=rows+mr+1;j++)
933    {
934      for (k=1;k<=cols;k++)
935      {
936        IMATELEM((*result),j,k) = IMATELEM((*result),j-mr,k);
937      }
938    }
939    for (j=rows+mr+1;j<=rows+1;j++)
940    {
941      for (k=1;k<=cols;k++)
942      {
943        IMATELEM((*result),j,k) = 0;
944      }
945    }
946  }
947  j = 0;
948  k = 0;
949  for (i=1;i<=result->rows();i++)
950  {
951    for(l=1;l<=result->cols();l++)
952    if (IMATELEM((*result),i,l) != 0)
953    {
954      j = si_max(j, i-1);
955      k = si_max(k, l-1);
956    }
957  }
958  intvec * exactresult=new intvec(j+1,k+1,0);
959  for (i=0;i<exactresult->rows();i++)
960  {
961    for (j=0;j<exactresult->cols();j++)
962    {
963      IMATELEM(*exactresult,i+1,j+1) = IMATELEM(*result,i+1,j+1);
964    }
965  }
966  if (row_shift!=NULL) *row_shift = mr;
967  delete result;
968  return exactresult;
969}
970
971/*2
972* minbare via syzygies
973*/
974ideal syMinBase(ideal arg)
975{
976  intvec ** weights=NULL;
977  int leng;
978  if (idIs0(arg)) return idInit(1,arg->rank);
979  resolvente res=syResolvente(arg,1,&leng,&weights,TRUE);
980  ideal result=res[0];
981  omFreeSize((ADDRESS)res,leng*sizeof(ideal));
982  if (weights!=NULL)
983  {
984    if (*weights!=NULL)
985    {
986      delete (*weights);
987      *weights=NULL;
988    }
989    if ((leng>=1) && (*(weights+1)!=NULL))
990    {
991      delete *(weights+1);
992      *(weights+1)=NULL;
993    }
994  }
995  idSkipZeroes(result);
996  return result;
997}
998
999/*2
1000* computes Betti-numbers from a resolvente of
1001* (non-)homogeneous objects
1002* the numbers of entrees !=NULL in res and weights must be equal
1003* and < length
1004*/
1005intvec * syNewBetti(resolvente res, intvec ** weights, int length)
1006{
1007  intvec * result,*tocancel;
1008  int i,j,k,rsmin=0,rsmax=0,rs=0;
1009  BOOLEAN homog=TRUE;
1010
1011  if (weights!=NULL)           //---homogeneous Betti numbers
1012  {
1013/*--------------computes size of the field----------------------*/
1014    for (i=1;i<length;i++)
1015    {
1016      if (weights[i] != NULL)
1017      {
1018        for (j=1;j<(weights[i])->length();j++)
1019        {
1020          if ((*(weights[i]))[j]-i<rsmin) rsmin = (*(weights[i]))[j]-i;
1021          if ((*(weights[i]))[j]-i>rsmax) rsmax = (*(weights[i]))[j]-i;
1022        }
1023      }
1024    }
1025    i = 0;
1026    while (weights[i] != NULL) i++;
1027    i--;
1028    for (j=0;j<IDELEMS(res[i]);j++)
1029    {
1030      if (res[i]->m[j]!=NULL)
1031      {
1032        k = pFDeg(res[i]->m[j],currRing)+(*(weights[i]))[pGetComp(res[i]->m[j])]-i-1;
1033        if (k>rsmax) rsmax = k;
1034        if (k<rsmin) rsmin = k;
1035      }
1036    }
1037    for (j=1;j<(weights[0])->length();j++)
1038    {
1039      if ((*weights[0])[j]>rsmax) rsmax = (*weights[0])[j];
1040      if ((*weights[0])[j]<rsmin) rsmin = (*weights[0])[j];
1041    }
1042//Print("rsmax = %d\n",rsmax);
1043//Print("rsmin = %d\n",rsmin);
1044    rs = rsmax-rsmin+1;
1045    result = new intvec(rs,i+2,0);
1046    tocancel = new intvec(rs);
1047/*-----------enter the Betti numbers-------------------------------*/
1048    if (/*idRankFreeModule(res[0])*/ res[0]->rank==0)
1049    {
1050      IMATELEM(*result,1-rsmin,1)=1;
1051    }
1052    else
1053    {
1054      for (i=1;i<(weights[0])->length();i++)
1055        IMATELEM(*result,(*weights[0])[i]+1-rsmin,1)++;
1056    }
1057    i = 1;
1058    while (weights[i]!=NULL)
1059    {
1060      for (j=1;j<(weights[i])->length();j++)
1061      {
1062        IMATELEM(*result,(*(weights[i]))[j]-i+1-rsmin,i+1)++;
1063      }
1064      i++;
1065    }
1066    i--;
1067    for (j=0;j<IDELEMS(res[i]);j++)
1068    {
1069      if (res[i]->m[j]!=NULL)
1070      {
1071        k = pFDeg(res[i]->m[j],currRing)+(*(weights[i]))[pGetComp(res[i]->m[j])]-i;
1072        IMATELEM(*result,k-rsmin,i+2)++;
1073      }
1074    }
1075  }
1076  else                //-----the non-homgeneous case
1077  {
1078    homog = FALSE;
1079    tocancel = new intvec(1);
1080    k = length;
1081    while ((k>0) && (idIs0(res[k-1]))) k--;
1082    result = new intvec(1,k+1,0);
1083    (*result)[0] = res[0]->rank;
1084    for (i=0;i<length;i++)
1085    {
1086      if (res[i]!=NULL)
1087      {
1088        for (j=0;j<IDELEMS(res[i]);j++)
1089        {
1090          if (res[i]->m[j]!=NULL) (*result)[i+1]++;
1091        }
1092      }
1093    }
1094  }
1095/*--------computes the Betti numbers for the minimized reolvente----*/
1096
1097  i = 1;
1098  while ((res[i]!=NULL) && (weights[i]!=NULL))
1099  {
1100    syDetect(res[i],i,rsmin,homog,weights[i],tocancel);
1101    if (homog)
1102    {
1103      for (j=0;j<rs-1;j++)
1104      {
1105        IMATELEM((*result),j+1,i+1) -= (*tocancel)[j];
1106        IMATELEM((*result),j+1,i+2) -= (*tocancel)[j+1];
1107      }
1108      IMATELEM((*result),rs,i+1) -= (*tocancel)[rs-1];
1109    }
1110    else
1111    {
1112      (*result)[i+1] -= (*tocancel)[0];
1113      (*result)[i+2] -= (*tocancel)[0];
1114    }
1115    i++;
1116  }
1117
1118/*--------print Betti numbers for control---------------------------*/
1119  for(i=rsmin;i<=rsmax;i++)
1120  {
1121    Print("%2d:",i);
1122    for(j=1;j<=result->cols();j++)
1123    {
1124      Print(" %5d",IMATELEM(*result,i-rsmin+1,j));
1125    }
1126    PrintLn();
1127  }
1128  return result;
1129}
1130
1131/*2
1132* is looking for the minimal minimized module of a resolvente
1133* i.e. returns 0 if res comes from a mres-command and 1 in the
1134* case of res-commands
1135*/
1136int syIsMinimizedFrom(resolvente res,int length)
1137{
1138  poly p;
1139  int i,j=length;
1140
1141  while ((j>0) && (res[j-1]==NULL)) j--;
1142  while (j>0)
1143  {
1144    for (i=0;i<IDELEMS(res[j-1]);i++)
1145    {
1146      p = res[j-1]->m[i];
1147      while (p!=NULL)
1148      {
1149        if (pLmIsConstantComp(p)) return j;
1150        p = pNext(p);
1151      }
1152    }
1153    j--;
1154  }
1155  return j;
1156}
Note: See TracBrowser for help on using the repository browser.