source: git/Singular/syz.cc @ 4a81ec

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