source: git/kernel/syz.cc @ b6201e

spielwiese
Last change on this file since b6201e was b6201e, checked in by Viktor Levandovskyy <levandov@…>, 19 years ago
*levandov: for prune - in syzGaussOne right multiplication changed to the left one git-svn-id: file:///usr/local/Singular/svn/trunk@7594 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 28.5 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/* $Id: syz.cc,v 1.6 2004-10-26 15:40:28 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(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])!=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]);
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])
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])+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"); P
794rintLn();
795  #endif
796  int rkl=l = si_max(idRankFreeModule(res[0]),res[0]->rank);
797  i = 0;
798  while ((i<length) && (res[i]!=NULL))
799  {
800    if (IDELEMS(res[i])>l) l = IDELEMS(res[i]);
801    i++;
802  }
803  temp1 = (int*)omAlloc0((l+1)*sizeof(int));
804  temp2 = (int*)omAlloc((l+1)*sizeof(int));
805  rows = 1;
806  mr = 1;
807  cols++;
808  for (i=0;i<cols-1;i++)
809  {
810    if ((i==0) && (weights!=NULL)) pSetModDeg(weights);
811    memset(temp2,0,(l+1)*sizeof(int));
812    for (j=0;j<IDELEMS(res[i]);j++)
813    {
814      if (res[i]->m[j]!=NULL)
815      {
816        if ((pGetComp(res[i]->m[j])>l)
817        || ((i>1) && (res[i-1]->m[pGetComp(res[i]->m[j])-1]==NULL)))
818        {
819          WerrorS("input not a resolvent");
820          omFreeSize((ADDRESS)temp1,(l+1)*sizeof(int));
821          omFreeSize((ADDRESS)temp2,(l+1)*sizeof(int));
822          return NULL;
823        }
824        temp2[j+1] = pFDeg(res[i]->m[j])+temp1[pGetComp(res[i]->m[j])];
825        if (temp2[j+1]-i>rows) rows = temp2[j+1]-i;
826        if (temp2[j+1]-i<mr) mr = temp2[j+1]-i;
827      }
828    }
829    if ((i==0) && (weights!=NULL)) pSetModDeg(NULL);
830    temp3 = temp1;
831    temp1 = temp2;
832    temp2 = temp3;
833  }
834  mr--;
835  if (weights!=NULL)
836  {
837    for(j=0;j<weights->length();j++)
838    {
839      if (rows <(*weights)[j]+1) rows=(-mr)+(*weights)[j]+1;
840    }
841  }
842  /*------ computation betti numbers --------------*/
843  rows -= mr;
844  result = new intvec(rows+1,cols,0);
845  if (weights!=NULL)
846  {
847    for(j=0;j<weights->length();j++)
848    {
849      IMATELEM((*result),(-mr)+(*weights)[j]+1,1) ++;
850      // Print("imat(%d,%d)++ -> %d\n",(-mr)+(*weights)[j]+1, 1, IMATELEM((*result),(-mr)+(*weights)[j]+1,1));
851    }
852  }
853  else
854  {
855    (*result)[(-mr)*cols] = /*idRankFreeModule(res[0])*/ rkl;
856    if ((!idIs0(res[0])) && ((*result)[(-mr)*cols]==0))
857      (*result)[(-mr)*cols] = 1;
858  }
859  tocancel = (int*)omAlloc0((rows+1)*sizeof(int));
860  memset(temp1,0,(l+1)*sizeof(int));
861  if (weights!=NULL)
862  {
863    memset(temp2,0,l*sizeof(int));
864    pSetModDeg(weights);
865  }
866  else
867    memset(temp2,0,l*sizeof(int));
868  int dummy = syDetect(res[0],0,TRUE,temp2,tocancel);
869  if (weights!=NULL) pSetModDeg(NULL);
870  if (tomin)
871  {
872    //(*result)[(-mr)*cols] -= dummy;
873    for(j=0;j<=rows+mr;j++)
874    {
875      // Print("tocancel[%d]=%d imat(%d,%d)=%d\n",j,tocancel[j],(-mr)+j+1,1,IMATELEM((*result),(-mr)+j+1,1));
876      IMATELEM((*result),(-mr)+j+1,1) -= tocancel[j];
877    }
878  }
879  for (i=0;i<cols-1;i++)
880  {
881    if ((i==0) && (weights!=NULL)) pSetModDeg(weights);
882    memset(temp2,0,l*sizeof(int));
883    for (j=0;j<IDELEMS(res[i]);j++)
884    {
885      if (res[i]->m[j]!=NULL)
886      {
887        temp2[j+1] = pFDeg(res[i]->m[j])+temp1[pGetComp(res[i]->m[j])];
888        //(*result)[i+1+(temp2[j+1]-i-1)*cols]++;
889        //if (temp2[j+1]>i) IMATELEM((*result),temp2[j+1]-i-mr,i+2)++;
890        IMATELEM((*result),temp2[j+1]-i-mr,i+2)++;
891      }
892      else if (i==0)
893      {
894        if (j<r0_len) IMATELEM((*result),-mr,2)++;
895      }
896    }
897  /*------ computation betti numbers, if res not minimal --------------*/
898    if (tomin)
899    {
900      for (j=mr;j<rows+mr;j++)
901      {
902        //(*result)[i+1+j*cols] -= tocancel[j+1];
903        IMATELEM((*result),j+1-mr,i+2) -= tocancel[j+1];
904      }
905      if ((i<length-1) && (res[i+1]!=NULL))
906      {
907        memset(tocancel,0,(rows+1)*sizeof(int));
908        dummy = syDetect(res[i+1],i+1,TRUE,temp2,tocancel);
909        for (j=0;j<rows;j++)
910        {
911          //(*result)[i+1+j*cols] -= tocancel[j];
912          IMATELEM((*result),j+1,i+2) -= tocancel[j];
913        }
914      }
915    }
916    temp3 = temp1;
917    temp1 = temp2;
918    temp2 = temp3;
919    for (j=0;j<=rows;j++)
920    {
921    //  if (((*result)[i+1+j*cols]!=0) && (j>*regularity)) *regularity = j;
922      if ((IMATELEM((*result),j+1,i+2)!=0) && (j>*regularity)) *regularity = j;
923    }
924    if ((i==0) && (weights!=NULL)) pSetModDeg(NULL);
925  }
926  // Print("nach minim:\n"); result->show(); PrintLn();
927  /*------ clean up --------------*/
928  omFreeSize((ADDRESS)tocancel,(rows+1)*sizeof(int));
929  omFreeSize((ADDRESS)temp1,(l+1)*sizeof(int));
930  omFreeSize((ADDRESS)temp2,(l+1)*sizeof(int));
931  if ((tomin) && (mr<0))  // deletes the first (zero) line
932  {
933    for (j=1;j<=rows+mr+1;j++)
934    {
935      for (k=1;k<=cols;k++)
936      {
937        IMATELEM((*result),j,k) = IMATELEM((*result),j-mr,k);
938      }
939    }
940    for (j=rows+mr+1;j<=rows+1;j++)
941    {
942      for (k=1;k<=cols;k++)
943      {
944        IMATELEM((*result),j,k) = 0;
945      }
946    }
947  }
948  j = 0;
949  k = 0;
950  for (i=1;i<=result->rows();i++)
951  {
952    for(l=1;l<=result->cols();l++)
953    if (IMATELEM((*result),i,l) != 0)
954    {
955      j = si_max(j, i-1);
956      k = si_max(k, l-1);
957    }
958  }
959  intvec * exactresult=new intvec(j+1,k+1,0);
960  for (i=0;i<exactresult->rows();i++)
961  {
962    for (j=0;j<exactresult->cols();j++)
963    {
964      IMATELEM(*exactresult,i+1,j+1) = IMATELEM(*result,i+1,j+1);
965    }
966  }
967  if (row_shift!=NULL) *row_shift = mr;
968  delete result;
969  return exactresult;
970}
971
972/*2
973* minbare via syzygies
974*/
975ideal syMinBase(ideal arg)
976{
977  intvec ** weights=NULL;
978  int leng;
979  if (idIs0(arg)) return idInit(1,arg->rank);
980  resolvente res=syResolvente(arg,1,&leng,&weights,TRUE);
981  ideal result=res[0];
982  omFreeSize((ADDRESS)res,leng*sizeof(ideal));
983  if (weights!=NULL)
984  {
985    if (*weights!=NULL)
986    {
987      delete (*weights);
988      *weights=NULL;
989    }
990    if ((leng>=1) && (*(weights+1)!=NULL))
991    {
992      delete *(weights+1);
993      *(weights+1)=NULL;
994    }
995  }
996  idSkipZeroes(result);
997  return result;
998}
999
1000/*2
1001* computes Betti-numbers from a resolvente of
1002* (non-)homogeneous objects
1003* the numbers of entrees !=NULL in res and weights must be equal
1004* and < length
1005*/
1006intvec * syNewBetti(resolvente res, intvec ** weights, int length)
1007{
1008  intvec * result,*tocancel;
1009  int i,j,k,rsmin=0,rsmax=0,rs=0;
1010  BOOLEAN homog=TRUE;
1011
1012  if (weights!=NULL)           //---homogeneous Betti numbers
1013  {
1014/*--------------computes size of the field----------------------*/
1015    for (i=1;i<length;i++)
1016    {
1017      if (weights[i] != NULL)
1018      {
1019        for (j=1;j<(weights[i])->length();j++)
1020        {
1021          if ((*(weights[i]))[j]-i<rsmin) rsmin = (*(weights[i]))[j]-i;
1022          if ((*(weights[i]))[j]-i>rsmax) rsmax = (*(weights[i]))[j]-i;
1023        }
1024      }
1025    }
1026    i = 0;
1027    while (weights[i] != NULL) i++;
1028    i--;
1029    for (j=0;j<IDELEMS(res[i]);j++)
1030    {
1031      if (res[i]->m[j]!=NULL)
1032      {
1033        k = pFDeg(res[i]->m[j])+(*(weights[i]))[pGetComp(res[i]->m[j])]-i-1;
1034        if (k>rsmax) rsmax = k;
1035        if (k<rsmin) rsmin = k;
1036      }
1037    }
1038    for (j=1;j<(weights[0])->length();j++)
1039    {
1040      if ((*weights[0])[j]>rsmax) rsmax = (*weights[0])[j];
1041      if ((*weights[0])[j]<rsmin) rsmin = (*weights[0])[j];
1042    }
1043//Print("rsmax = %d\n",rsmax);
1044//Print("rsmin = %d\n",rsmin);
1045    rs = rsmax-rsmin+1;
1046    result = new intvec(rs,i+2,0);
1047    tocancel = new intvec(rs);
1048/*-----------enter the Betti numbers-------------------------------*/
1049    if (/*idRankFreeModule(res[0])*/ res[0]->rank==0)
1050    {
1051      IMATELEM(*result,1-rsmin,1)=1;
1052    }
1053    else
1054    {
1055      for (i=1;i<(weights[0])->length();i++)
1056        IMATELEM(*result,(*weights[0])[i]+1-rsmin,1)++;
1057    }
1058    i = 1;
1059    while (weights[i]!=NULL)
1060    {
1061      for (j=1;j<(weights[i])->length();j++)
1062      {
1063        IMATELEM(*result,(*(weights[i]))[j]-i+1-rsmin,i+1)++;
1064      }
1065      i++;
1066    }
1067    i--;
1068    for (j=0;j<IDELEMS(res[i]);j++)
1069    {
1070      if (res[i]->m[j]!=NULL)
1071      {
1072        k = pFDeg(res[i]->m[j])+(*(weights[i]))[pGetComp(res[i]->m[j])]-i;
1073        IMATELEM(*result,k-rsmin,i+2)++;
1074      }
1075    }
1076  }
1077  else                //-----the non-homgeneous case
1078  {
1079    homog = FALSE;
1080    tocancel = new intvec(1);
1081    k = length;
1082    while ((k>0) && (idIs0(res[k-1]))) k--;
1083    result = new intvec(1,k+1,0);
1084    (*result)[0] = res[0]->rank;
1085    for (i=0;i<length;i++)
1086    {
1087      if (res[i]!=NULL)
1088      {
1089        for (j=0;j<IDELEMS(res[i]);j++)
1090        {
1091          if (res[i]->m[j]!=NULL) (*result)[i+1]++;
1092        }
1093      }
1094    }
1095  }
1096/*--------computes the Betti numbers for the minimized reolvente----*/
1097
1098  i = 1;
1099  while ((res[i]!=NULL) && (weights[i]!=NULL))
1100  {
1101    syDetect(res[i],i,rsmin,homog,weights[i],tocancel);
1102    if (homog)
1103    {
1104      for (j=0;j<rs-1;j++)
1105      {
1106        IMATELEM((*result),j+1,i+1) -= (*tocancel)[j];
1107        IMATELEM((*result),j+1,i+2) -= (*tocancel)[j+1];
1108      }
1109      IMATELEM((*result),rs,i+1) -= (*tocancel)[rs-1];
1110    }
1111    else
1112    {
1113      (*result)[i+1] -= (*tocancel)[0];
1114      (*result)[i+2] -= (*tocancel)[0];
1115    }
1116    i++;
1117  }
1118
1119/*--------print Betti numbers for control---------------------------*/
1120  for(i=rsmin;i<=rsmax;i++)
1121  {
1122    Print("%2d:",i);
1123    for(j=1;j<=result->cols();j++)
1124    {
1125      Print(" %5d",IMATELEM(*result,i-rsmin+1,j));
1126    }
1127    PrintLn();
1128  }
1129  return result;
1130}
1131
1132/*2
1133* is looking for the minimal minimized module of a resolvente
1134* i.e. returns 0 if res comes from a mres-command and 1 in the
1135* case of res-commands
1136*/
1137int syIsMinimizedFrom(resolvente res,int length)
1138{
1139  poly p;
1140  int i,j=length;
1141
1142  while ((j>0) && (res[j-1]==NULL)) j--;
1143  while (j>0)
1144  {
1145    for (i=0;i<IDELEMS(res[j-1]);i++)
1146    {
1147      p = res[j-1]->m[i];
1148      while (p!=NULL)
1149      {
1150        if (pLmIsConstantComp(p)) return j;
1151        p = pNext(p);
1152      }
1153    }
1154    j--;
1155  }
1156  return j;
1157}
Note: See TracBrowser for help on using the repository browser.