source: git/Singular/syz.cc @ b719a3

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