source: git/Singular/syz.cc @ 38e02cf

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