source: git/Singular/syz.cc @ 82716e

spielwiese
Last change on this file since 82716e was 25003c, checked in by Thomas Siebert <siebert@…>, 26 years ago
*** empty log message *** git-svn-id: file:///usr/local/Singular/svn/trunk@1532 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 24.9 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/* $Id: syz.cc,v 1.9 1998-04-29 07:05:30 siebert Exp $ */
5
6/*
7* ABSTRACT: resolutions
8*/
9
10
11#include "mod2.h"
12#include "tok.h"
13#include "mmemory.h"
14#include "polys.h"
15#include "febase.h"
16#include "kstd1.h"
17#include "kutil.h"
18#include "spolys.h"
19#include "stairc.h"
20#include "ipid.h"
21#include "cntrlc.h"
22#include "ipid.h"
23#include "intvec.h"
24#include "ipshell.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=NULL;
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 (deg0!=NULL)
216    idDelete(&deg0);
217}
218
219void syMinimizeResolvente(resolvente res, int length, int first)
220{
221  int syzIndex=first;
222
223  if (syzIndex<1) syzIndex=1;
224  while ((syzIndex<length-1) && (res[syzIndex]!=NULL) && (res[syzIndex+1]!=NULL))
225  {
226    syMinStep(res[syzIndex-1],res[syzIndex],FALSE,res[syzIndex+1]);
227    syzIndex++;
228  }
229  if (res[syzIndex]!=NULL)
230    syMinStep(res[syzIndex-1],res[syzIndex]);
231}
232
233/*2
234* resolution of ideal/module arg, <=maxlength steps, (r[0..maxlength])
235*   no limitation in length if maxlength==0
236* input:arg
237*       minim: TRUE means mres cmd, FALSE res cmd.
238*       if *len!=0: modul weights: weights[0]
239*          (and weights is defined:weights[0..len-1]
240*
241* output:resolvente r[0..length-1],
242*        modul weights: weights[0..length-1]
243*/
244resolvente syResolvente(ideal arg, int maxlength, int * length,
245                        intvec *** weights, BOOLEAN minim)
246{
247  resolvente res;
248  resolvente newres;
249  tHomog hom=isNotHomog;
250  ideal temp=NULL;
251  intvec *w,*w1=NULL,**tempW;
252  int i,k,syzIndex = 0,j;
253  int Kstd1_OldDeg=Kstd1_deg;
254  BOOLEAN completeMinim;
255  BOOLEAN oldDegBound=TEST_OPT_DEGBOUND;
256  int wlength=*length;
257
258  if (maxlength!=-1) *length = maxlength+1;
259  else              *length = 5;
260  if ((wlength!=0) && (*length!=wlength))
261  {
262    intvec **wtmp = (intvec**)Alloc0((*length)*sizeof(intvec*));
263    wtmp[0]=(*weights)[0];
264    Free((ADDRESS)*weights,wlength*sizeof(intvec*));
265    *weights=wtmp;
266  }
267  res = (resolvente)Alloc0((*length)*sizeof(ideal));
268  res[0] = idCopy(arg);
269  if ((weights==NULL) || (*weights==NULL) || ((*weights)[0]==NULL))
270  {
271    hom=(tHomog)idHomModule(res[0],currQuotient,&w);
272    if (hom==isHomog)
273    {
274      *weights = (intvec**)Alloc0((*length)*sizeof(intvec*));
275      if (w!=NULL) (*weights)[0] = ivCopy(w);
276//Print("die %d Modulegewichte sind:\n",w->length());
277//w->show();
278//Print("\n");
279    }
280  }
281  else
282  {
283    if ((weights!=NULL) && (*weights!=NULL)&& ((*weights)[0]!=NULL))
284    {
285      w = ivCopy((*weights)[0]);
286      hom = isHomog;
287    }
288  }
289  if (hom==isHomog)
290  {
291    w1 = syPrepareModComp(res[0],&w);
292    if (w!=NULL) { delete w;w=NULL; }
293  }
294  while ((!idIs0(res[syzIndex])) &&
295         ((maxlength==-1) || (syzIndex<=maxlength)))
296   // (syzIndex<maxlength+(int)minim)))
297  /*compute one step more for minimizing*/
298  {
299    if (Kstd1_deg!=0) Kstd1_deg++;
300    if (syzIndex+1==*length)
301    {
302      newres = (resolvente)Alloc0((*length+5)*sizeof(ideal));
303      tempW = (intvec**)Alloc0((*length+5)*sizeof(intvec*));
304      for (j=0;j<*length;j++)
305      {
306        newres[j] = res[j];
307        if (*weights!=NULL) tempW[j] = (*weights)[j];
308        /*else              tempW[j] = NULL;*/
309      }
310      Free((ADDRESS)res,*length*sizeof(ideal));
311      Free((ADDRESS)*weights,*length*sizeof(intvec*));
312      *length += 5;
313      res=newres;
314      *weights = tempW;
315    }
316    if (minim || (syzIndex!=0))
317    {
318      temp = kInterRed(res[syzIndex],currQuotient);
319      idDelete(&res[syzIndex]);
320      idSkipZeroes(temp);
321      res[syzIndex] = temp;
322    }
323    temp = NULL;
324    if ((currQuotient==NULL)&&(syzIndex==0)&& (!TEST_OPT_DEGBOUND))
325    {
326      res[/*syzIndex+*/1] = idSyzygies(res[0/*syzIndex*/],
327        NULL/*currQuotient*/,hom,&w1,TRUE,Kstd1_deg);
328      if ((!TEST_OPT_NOTREGULARITY) && (Kstd1_deg>0)) test |= Sy_bit(OPT_DEGBOUND);
329    }
330    else
331      res[syzIndex+1] = idSyzygies(res[syzIndex],currQuotient,hom,&w1);
332    completeMinim=(syzIndex!=maxlength) || (maxlength ==-1) || (hom!=isHomog);
333    syzIndex++;
334    if (TEST_OPT_PROT) Print("[%d]\n",syzIndex);
335    if ((minim)||(syzIndex>1))
336      syMinStep(res[syzIndex-1],res[syzIndex],!completeMinim,NULL,hom);
337    if (!completeMinim)
338    /*minim is TRUE, we are in the module: maxlength, maxlength <>0*/
339    {
340      idDelete(&res[syzIndex]);
341    }
342    if ((hom == isHomog) && (!idIs0(res[syzIndex])))
343    {
344//Print("die %d Modulegewichte sind:\n",w1->length());
345//w1->show();
346//Print("\n");
347      k = idRankFreeModule(res[syzIndex]);
348      w = new intvec(k+IDELEMS(res[syzIndex]));
349      (*weights)[syzIndex] = new intvec(k);
350      for (i=0;i<k;i++)
351      {
352        if (res[syzIndex-1]->m[i]!=NULL) // hs
353        {
354          (*w)[i] = pFDeg(res[syzIndex-1]->m[i]);
355          if (pGetComp(res[syzIndex-1]->m[i])>0)
356            (*w)[i] += (*w1)[pGetComp(res[syzIndex-1]->m[i])-1];
357          (*((*weights)[syzIndex]))[i] = (*w)[i];
358        }
359      }
360      for (i=k;i<k+IDELEMS(res[syzIndex]);i++)
361      {
362        if (res[syzIndex]->m[i-k]!=NULL)
363          (*w)[i] = pFDeg(res[syzIndex]->m[i-k])
364                    +(*w)[pGetComp(res[syzIndex]->m[i-k])-1];
365      }
366      delete w1;
367      w1 = w;
368      w = NULL;
369    }
370  }
371  if ((syzIndex!=0) && (res[syzIndex]!=NULL) && (idIs0(res[syzIndex])))
372    idDelete(&res[syzIndex]);
373  if (w1!=NULL) delete w1;
374  //if ((maxlength!=-1) && (!minim))
375  //{
376  //  while (syzIndex<=maxlength)
377  //  {
378  //    res[syzIndex]=idInit(1,1);
379  //    syzIndex++;
380  //    if (syzIndex<=maxlength)
381  //    {
382  //      res[syzIndex]=idFreeModule(1);
383  //      syzIndex++;
384  //    }
385  //  }
386  //}
387  Kstd1_deg=Kstd1_OldDeg;
388  if (!oldDegBound)
389    test &= ~Sy_bit(OPT_DEGBOUND);
390  return res;
391}
392
393syStrategy syResolution(ideal arg, int maxlength,intvec * w, BOOLEAN minim)
394{
395  int typ0;
396  syStrategy result=(syStrategy)Alloc0(sizeof(ssyStrategy));
397
398  if (w!=NULL)
399  {
400    result->weights = (intvec**)Alloc0(sizeof(intvec*));
401    (result->weights)[0] = ivCopy(w);
402    result->length = 1;
403  }
404  resolvente fr = syResolvente(arg,maxlength,&(result->length),&(result->weights),minim),fr1;
405  if (minim)
406  {
407    result->minres = (resolvente)Alloc0((result->length+1)*sizeof(ideal));
408    fr1 =  result->minres;
409  }
410  else
411  {
412    result->fullres = (resolvente)Alloc0((result->length+1)*sizeof(ideal));
413    fr1 =  result->fullres;
414  }
415  for (int i=result->length-1;i>=0;i--)
416  {
417    if (fr[i]!=NULL)
418      fr1[i] = fr[i];
419    fr[i] = NULL;
420  }
421  Free((ADDRESS)fr,(result->length)*sizeof(ideal));
422  return result;
423}
424
425
426resolvente syMinRes(ideal arg, int maxlength, int * length, BOOLEAN minim)
427{
428  resolvente res;
429  resolvente newres;
430  tHomog hom;
431  ideal temp=NULL;
432  intvec *w,*w1=NULL;
433  int i,k,syzIndex = 0,j;
434  int Kstd1_OldDeg=Kstd1_deg;
435  BOOLEAN isIdeal=FALSE,oldDegBound=TEST_OPT_DEGBOUND;
436
437  if (maxlength!=-1) *length = maxlength+1;
438  else              *length = 4;
439  res = (resolvente)Alloc0((*length)*sizeof(ideal));
440  res[0] = idCopy(arg);
441  hom=(tHomog)idHomModule(res[0],currQuotient,&w);
442//Print("die %d Modulegewichte sind:\n",w->length());
443//w->show();
444//PrintLn();
445  if (hom==isHomog)
446  {
447    //if (w!=NULL)
448    //{
449    //  w->show();
450    //  PrintLn();
451    //  if (res[0]!=NULL) jjPRINT_MA0(idModule2Matrix(idCopy(res[0])),NULL);
452    //}
453    int maxxx = 0;
454    if (w!=NULL) maxxx = w->length();
455    if (maxxx==1)
456    {
457      maxxx = 2;
458      isIdeal = TRUE;
459    }
460    w1 = new intvec(maxxx+IDELEMS(arg));
461    if (!isIdeal)
462    {
463      for (i=0;i<maxxx;i++)
464      {
465        (*w1)[i] = (*w)[i];
466      }
467    }
468    for (i=maxxx;i<maxxx+IDELEMS(arg);i++)
469    {
470      if (res[0]->m[i-maxxx]!=0)
471        (*w1)[i] = pFDeg(res[0]->m[i-maxxx])+(*w)[pGetComp(res[0]->m[i-maxxx])];
472    }
473    if (w!=NULL) {delete w;w=NULL;}
474  }
475  while ((!idIs0(res[syzIndex]))
476  && ((maxlength==-1)
477    || (syzIndex<maxlength)))
478   // (syzIndex<maxlength+(int)minim)))
479  /*compute one step more for minimizing*/
480  {
481    if (Kstd1_deg!=0) Kstd1_deg++;
482    if (syzIndex+1==*length)
483    {
484      newres = (resolvente)Alloc0((*length+4)*sizeof(ideal));
485      for (j=0;j<*length;j++) newres[j] = res[j];
486      Free((ADDRESS)res,*length*sizeof(ideal));
487      *length += 4;
488      res=newres;
489    }
490    if (minim || (syzIndex!=0))
491    {
492      temp = kInterRed(res[syzIndex],currQuotient);
493      idDelete(&res[syzIndex]);
494      res[syzIndex] = temp;
495    }
496    temp = NULL;
497    if ((hom==isHomog) && (currRing->OrdSgn==1))
498    {
499//PrintS("Input: ");
500//for (i=0;i<IDELEMS(res[syzIndex]);i++)
501//{
502//Print("Element %d: ",i);pWrite(res[syzIndex]->m[i]);
503//}
504      if ((currQuotient==NULL)&&(syzIndex==0)&& (!TEST_OPT_DEGBOUND))
505      {
506        res[/*syzIndex+*/1] = idSyzMin(res[0/*syzIndex*/],
507          NULL/*currQuotient*/,hom,&w1,TRUE,Kstd1_deg);
508        if ((!TEST_OPT_NOTREGULARITY) && (Kstd1_deg>0)) test |= Sy_bit(OPT_DEGBOUND);
509//Print("Regularity is %d \n",Kstd1_deg);
510      }
511      else
512        res[syzIndex+1] = idSyzMin(res[syzIndex],
513          currQuotient,hom,&w1,FALSE,Kstd1_deg);
514//PrintS("Input: ");
515//for (i=0;i<IDELEMS(res[syzIndex+1]);i++)
516//{
517//Print("Element %d: ",i);pWrite(res[syzIndex+1]->m[i]);
518//}
519    }
520    else
521      res[syzIndex+1] = idSyzygies(res[syzIndex],currQuotient,hom,&w1);
522//Print("im %d-ten Syzygienmodul: \n",syzIndex+1);
523//for (i=0;i<IDELEMS(res[syzIndex+1]);i++)
524//{
525//Print("El %d: ",i);pWrite(res[syzIndex+1]->m[i]);
526//}
527    syzIndex++;
528    if ((maxlength!=-1) && (syzIndex==maxlength)) idDelete(&res[syzIndex]);
529    if (TEST_OPT_PROT) Print("[%d]\n",syzIndex);
530    if  ((hom==isHomog) && (res[syzIndex]!=NULL))
531    {
532//Print("die %d Modulegewichte sind:\n",w1->length());
533//w1->show();
534//PrintLn();
535      k = idRankFreeModule(res[syzIndex]);
536      w = new intvec(k+IDELEMS(res[syzIndex])+1);
537      for (i=1;i<=k;i++)
538      {
539        (*w)[i] = pFDeg(res[syzIndex-1]->m[i-1])
540                  +(*w1)[pGetComp(res[syzIndex-1]->m[i-1])];
541      }
542      for (i=k+1;i<k+1+IDELEMS(res[syzIndex]);i++)
543      {
544        if (res[syzIndex]->m[i-k-1]!=NULL)
545          (*w)[i] = pFDeg(res[syzIndex]->m[i-k-1])
546                    +(*w)[pGetComp(res[syzIndex]->m[i-k-1])];
547      }
548      delete w1;
549      w1 = w;
550      w = NULL;
551    }
552  }
553  if (w1!=NULL) delete w1;
554  if (maxlength!=-1)
555  {
556    if ((syzIndex<maxlength-1) && (res[syzIndex]!=NULL)) syzIndex++;
557    while (syzIndex<maxlength)
558    {
559      res[syzIndex]=idInit(1,1);
560      syzIndex++;
561    }
562  }
563  Kstd1_deg=Kstd1_OldDeg;
564  if (!oldDegBound)
565    test &= ~Sy_bit(OPT_DEGBOUND);
566  return res;
567}
568
569void syDetect(ideal id,int index,BOOLEAN homog,int * degrees,int * tocancel)
570{
571  int i, j, k, ModComp;
572  poly p, q, qq, Unit1, Unit2;
573  ideal temp;
574
575  if (idIs0(id)) return;
576  temp = idInit(IDELEMS(id),id->rank);
577  for (i=0;i<IDELEMS(id);i++)
578  {
579    p = id->m[i];
580    while (p!=NULL)
581    {
582      if (pIsConstantComp(p))
583      {
584        if (temp->m[i]==NULL)
585        {
586          temp->m[i] = pHead(p);
587          q = temp->m[i];
588        }
589        else
590        {
591          pNext(q) = pHead(p);
592          pIter(q);
593        }
594      }
595      pIter(p);
596    }
597  }
598  i = IDELEMS(id);
599  while ((i>0) && (temp->m[i-1]==NULL)) i--;
600  if (i==0)
601  {
602    idDelete(&temp);
603    return;
604  }
605  j = 0;
606  while ((j<i) && (temp->m[j]==NULL)) j++;
607  if (j<i)
608  {
609    p = temp->m[j];
610    temp->m[j] = NULL;
611  }
612  else
613  {
614    p = NULL;
615  }
616  while (p!=NULL)
617  {
618    if (homog)
619    {
620      k = pFDeg(id->m[j])+degrees[pGetComp(id->m[j])];
621      if (k>=index) tocancel[k-index]++;
622    }
623    else
624    {
625      tocancel[0]--;
626    }
627    ModComp = pGetComp(p);
628//Print("Element %d: ",j);pWrite(p);
629    Unit1 = pTakeOutComp1(&p,ModComp);
630    pSetComp(Unit1,0);
631    for (k=j+1;k<i;k++)
632    {
633      if (temp->m[k]!=NULL)
634      {
635//Print("erst %d: ",k);pWrite(temp->m[k]);
636        Unit2 = pTakeOutComp1(&(temp->m[k]),ModComp);
637        if (Unit2!=NULL)
638        {
639          number tmp=nCopy(pGetCoeff(Unit2));
640          tmp=nNeg(tmp);
641          qq = pCopy(p);
642          pMultN(qq,tmp);
643          nDelete(&tmp);
644          pMultN(temp->m[k],pGetCoeff(Unit1));
645          temp->m[k] = pAdd(temp->m[k],qq);
646//Print("dann %d: ",k);pWrite(temp->m[k]);
647          pDelete(&Unit2);
648        }
649      }
650    }
651//PrintLn();
652    pDelete(&Unit1);
653    pDelete(&p);
654    j++;
655    while ((j<i) && (temp->m[j]==NULL)) j++;
656    if (j>=i)
657    {
658      idDelete(&temp);
659      return;
660    }
661    p = temp->m[j];
662    temp->m[j] = NULL;
663  }
664  idDelete(&temp);
665}
666
667void syDetect(ideal id,int index,int rsmin, BOOLEAN homog,
668              intvec * degrees,intvec * tocancel)
669{
670  int * deg=NULL;
671  int * tocan=(int*) Alloc0(tocancel->length()*sizeof(int));
672  int i;
673
674  if (homog)
675  {
676    deg = (int*) Alloc0(degrees->length()*sizeof(int));
677    for (i=degrees->length();i>0;i--)
678      deg[i-1] = (*degrees)[i-1]-rsmin;
679  }
680  syDetect(id,index,homog,deg,tocan);
681  for (i=tocancel->length();i>0;i--)
682    (*tocancel)[i-1] = tocan[i-1];
683  if (homog)
684    Free((ADDRESS)deg,degrees->length()*sizeof(int));
685  Free((ADDRESS)tocan,tocancel->length()*sizeof(int));
686}
687
688/*2
689* computes the betti numbers from a given resolution
690* of length 'length' (0..length-1), not necessairily minimal,
691* (if weights are given, they are used)
692* returns the int matrix of betti numbers
693* and the regularity
694*/
695intvec * syBetti(resolvente res,int length, int * regularity,
696                 intvec* weights)
697{
698  int i,j=0,k=0,l,rows,cols;
699  int *temp1,*temp2,*temp3;/*used to compute degrees*/
700  int *tocancel; /*(BOOLEAN)tocancel[i]=element is superfluous*/
701
702  /*------ compute size --------------*/
703  *regularity = -1;
704  cols = length;
705  while ((cols>0)
706  && ((res[cols-1]==NULL)
707    || (idIs0(res[cols-1]))))
708  {
709    cols--;
710  }
711  intvec * result;
712  if (idIs0(res[0]))
713  {
714    if (res[0]==NULL)
715      result = new intvec(1,1,1);
716    else
717      result = new intvec(1,1,res[0]->rank);
718    return result;
719  }
720  l = max(idRankFreeModule(res[0]),res[0]->rank);
721  i = 0;
722  while ((i<length) && (res[i]!=NULL))
723  {
724    if (IDELEMS(res[i])>l) l = IDELEMS(res[i]);
725    i++;
726  }
727  temp1 = (int*)Alloc0((l+1)*sizeof(int));
728  temp2 = (int*)Alloc((l+1)*sizeof(int));
729  rows = 1;
730  cols++;
731  for (i=0;i<cols-1;i++)
732  {
733    if ((i==0) && (weights!=NULL)) pSetModDeg(weights);
734    memset(temp2,0,(l+1)*sizeof(int));
735    for (j=0;j<IDELEMS(res[i]);j++)
736    {
737      if (res[i]->m[j]!=NULL)
738     {
739        if ((pGetComp(res[i]->m[j])>l)
740        || ((i>1) && (res[i-1]->m[pGetComp(res[i]->m[j])-1]==NULL)))
741        {
742          WerrorS("input not a resolvent");
743          Free((ADDRESS)temp1,(l+1)*sizeof(int));
744          Free((ADDRESS)temp2,(l+1)*sizeof(int));
745          return NULL;
746        }
747        temp2[j+1] = pFDeg(res[i]->m[j])+temp1[pGetComp(res[i]->m[j])];
748        if (temp2[j+1]-i>rows) rows = temp2[j+1]-i;
749      }
750    }
751    if ((i==0) && (weights!=NULL)) pSetModDeg(NULL);
752    temp3 = temp1;
753    temp1 = temp2;
754    temp2 = temp3;
755  }
756  /*------ computation betti numbers --------------*/
757  result = new intvec(rows,cols,0);
758  (*result)[0] = /*idRankFreeModule(res[0])*/ res[0]->rank;
759  if ((!idIs0(res[0])) && ((*result)[0]==0)) (*result)[0] = 1;
760  tocancel = (int*)Alloc0((rows+1)*sizeof(int));
761  memset(temp2,0,l*sizeof(int));
762  memset(temp1,0,(l+1)*sizeof(int));
763  syDetect(res[0],0,TRUE,temp2,tocancel);
764  (*result)[0] -= tocancel[0];
765  for (i=0;i<cols-1;i++)
766  {
767    if ((i==0) && (weights!=NULL)) pSetModDeg(weights);
768    memset(temp2,0,l*sizeof(int));
769    for (j=0;j<IDELEMS(res[i]);j++)
770    {
771      if (res[i]->m[j]!=NULL)
772      {
773        temp2[j+1] = pFDeg(res[i]->m[j])+temp1[pGetComp(res[i]->m[j])];
774        //(*result)[i+1+(temp2[j+1]-i-1)*cols]++;
775        if (temp2[j+1]>i) IMATELEM((*result),temp2[j+1]-i,i+2)++;
776      }
777    }
778  /*------ computation betti numbers, if res not minimal --------------*/
779    for (j=0;j<rows-1;j++)
780    {
781      //(*result)[i+1+j*cols] -= tocancel[j+1];
782      IMATELEM((*result),j+1,i+2) -= tocancel[j+1];
783    }
784    if ((i<length-1) && (res[i+1]!=NULL))
785    {
786      memset(tocancel,0,(rows+1)*sizeof(int));
787      syDetect(res[i+1],i+1,TRUE,temp2,tocancel);
788      for (j=0;j<rows;j++)
789      {
790        //(*result)[i+1+j*cols] -= tocancel[j];
791        IMATELEM((*result),j+1,i+2) -= tocancel[j];
792      }
793    }
794    temp3 = temp1;
795    temp1 = temp2;
796    temp2 = temp3;
797    for (j=0;j<rows;j++)
798    {
799    //  if (((*result)[i+1+j*cols]!=0) && (j>*regularity)) *regularity = j;
800      if ((IMATELEM((*result),j+1,i+2)!=0) && (j>*regularity)) *regularity = j;
801    }
802    if ((i==0) && (weights!=NULL)) pSetModDeg(NULL);
803  }
804  /*------ clean up --------------*/
805  Free((ADDRESS)tocancel,(rows+1)*sizeof(int));
806  Free((ADDRESS)temp1,(l+1)*sizeof(int));
807  Free((ADDRESS)temp2,(l+1)*sizeof(int));
808  j = 0;
809  k = 0;
810  for (i=0;i<rows*cols;i++)
811  {
812    if ((*result)[i] != 0)
813    {
814      j = i/cols;
815      if (i>k+j*cols) k = i-j*cols;
816    }
817  }
818  intvec * exactresult=new intvec(j+1,k+1,0);
819  for (i=0;i<exactresult->rows();i++)
820  {
821    for (j=0;j<exactresult->cols();j++)
822    {
823      IMATELEM(*exactresult,i+1,j+1) = IMATELEM(*result,i+1,j+1);
824    }
825  }
826  delete result;
827  return exactresult;
828}
829
830/*2
831* minbare via syzygies
832*/
833ideal syMinBase(ideal arg)
834{
835  intvec ** weights=NULL;
836  int leng;
837  if (idIs0(arg)) return idInit(1,arg->rank);
838  resolvente res=syResolvente(arg,1,&leng,&weights,TRUE);
839  ideal result=res[0];
840  Free((ADDRESS)res,leng*sizeof(ideal));
841  if (weights!=NULL)
842  {
843    if (*weights!=NULL)
844    {
845      delete (*weights);
846      *weights=NULL;
847    }
848    if ((leng>=1) && (*(weights+1)!=NULL))
849    {
850      delete *(weights+1);
851      *(weights+1)=NULL;
852    }
853  }
854  idSkipZeroes(result);
855  return result;
856}
857
858/*2
859* computes Betti-numbers from a resolvente of
860* (non-)homogeneous objects
861* the numbers of entrees !=NULL in res and weights must be equal
862* and < length
863*/
864intvec * syNewBetti(resolvente res, intvec ** weights, int length)
865{
866  intvec * result,*tocancel;
867  int i,j,k,rsmin=0,rsmax=0,rs=0;
868  BOOLEAN homog=TRUE;
869
870  if (weights!=NULL)           //---homogeneous Betti numbers
871  {
872/*--------------computes size of the field----------------------*/
873    for (i=1;i<length;i++)
874    {
875      if (weights[i] != NULL)
876      {
877        for (j=1;j<(weights[i])->length();j++)
878        {
879          if ((*(weights[i]))[j]-i<rsmin) rsmin = (*(weights[i]))[j]-i;
880          if ((*(weights[i]))[j]-i>rsmax) rsmax = (*(weights[i]))[j]-i;
881        }
882      }
883    }
884    i = 0;
885    while (weights[i] != NULL) i++;
886    i--;
887    for (j=0;j<IDELEMS(res[i]);j++)
888    {
889      if (res[i]->m[j]!=NULL)
890      {
891        k = pFDeg(res[i]->m[j])+(*(weights[i]))[pGetComp(res[i]->m[j])]-i-1;
892        if (k>rsmax) rsmax = k;
893        if (k<rsmin) rsmin = k;
894      }
895    }
896    for (j=1;j<(weights[0])->length();j++)
897    {
898      if ((*weights[0])[j]>rsmax) rsmax = (*weights[0])[j];
899      if ((*weights[0])[j]<rsmin) rsmin = (*weights[0])[j];
900    }
901//Print("rsmax = %d\n",rsmax);
902//Print("rsmin = %d\n",rsmin);
903    rs = rsmax-rsmin+1;
904    result = new intvec(rs,i+2,0);
905    tocancel = new intvec(rs);
906/*-----------enter the Betti numbers-------------------------------*/
907    if (/*idRankFreeModule(res[0])*/ res[0]->rank==0)
908    {
909      IMATELEM(*result,1-rsmin,1)=1;
910    }
911    else
912    {
913      for (i=1;i<(weights[0])->length();i++)
914        IMATELEM(*result,(*weights[0])[i]+1-rsmin,1)++;
915    }
916    i = 1;
917    while (weights[i]!=NULL)
918    {
919      for (j=1;j<(weights[i])->length();j++)
920      {
921        IMATELEM(*result,(*(weights[i]))[j]-i+1-rsmin,i+1)++;
922      }
923      i++;
924    }
925    i--;
926    for (j=0;j<IDELEMS(res[i]);j++)
927    {
928      if (res[i]->m[j]!=NULL)
929      {
930        k = pFDeg(res[i]->m[j])+(*(weights[i]))[pGetComp(res[i]->m[j])]-i;
931        IMATELEM(*result,k-rsmin,i+2)++;
932      }
933    }
934  }
935  else                //-----the non-homgeneous case
936  {
937    homog = FALSE;
938    tocancel = new intvec(1);
939    k = length;
940    while ((k>0) && (idIs0(res[k-1]))) k--;
941    result = new intvec (1,k+1,0);
942    (*result)[0] = res[0]->rank;
943    for (i=0;i<length;i++)
944    {
945      if (res[i]!=NULL)
946      {
947        for (j=0;j<IDELEMS(res[i]);j++)
948        {
949          if (res[i]->m[j]!=NULL) (*result)[i+1]++;
950        }
951      }
952    }
953  }
954/*--------computes the Betti numbers for the minimized reolvente----*/
955
956  i = 1;
957  while ((res[i]!=NULL) && (weights[i]!=NULL))
958  {
959    syDetect(res[i],i,rsmin,homog,weights[i],tocancel);
960    if (homog)
961    {
962      for (j=0;j<rs-1;j++)
963      {
964        IMATELEM((*result),j+1,i+1) -= (*tocancel)[j];
965        IMATELEM((*result),j+1,i+2) -= (*tocancel)[j+1];
966      }
967      IMATELEM((*result),rs,i+1) -= (*tocancel)[rs-1];
968    }
969    else
970    {
971      (*result)[i+1] -= (*tocancel)[0];
972      (*result)[i+2] -= (*tocancel)[0];
973    }
974    i++;
975  }
976
977/*--------print Betti numbers for control---------------------------*/
978  for(i=rsmin;i<=rsmax;i++)
979  {
980    Print("%2d:",i);
981    for(j=1;j<=result->cols();j++)
982    {
983      Print(" %5d",IMATELEM(*result,i-rsmin+1,j));
984    }
985    PrintLn();
986  }
987  return result;
988}
989
990/*2
991* is looking for the minimal minimized module of a resolvente
992* i.e. returns 0 if res comes from a mres-command and 1 in the
993* case of res-commands
994*/
995int syIsMinimizedFrom(resolvente res,int length)
996{
997  poly p;
998  int i,j=length;
999
1000  while ((j>0) && (res[j-1]==NULL)) j--;
1001  while (j>0)
1002  {
1003    for (i=0;i<IDELEMS(res[j-1]);i++)
1004    {
1005      p = res[j-1]->m[i];
1006      while (p!=NULL)
1007      {
1008        if (pIsConstantComp(p)) return j;
1009        p = pNext(p);
1010      }
1011    }
1012    j--;
1013  }
1014  return j;
1015}
Note: See TracBrowser for help on using the repository browser.