source: git/kernel/syz.cc @ b1dfaf

spielwiese
Last change on this file since b1dfaf was b1dfaf, checked in by Frank Seelisch <seelisch@…>, 14 years ago
patch from Kai (checked for problems under Windows OS: no problems) git-svn-id: file:///usr/local/Singular/svn/trunk@13210 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 29.8 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/* $Id$ */
5
6/*
7* ABSTRACT: resolutions
8*/
9
10
11#include <kernel/mod2.h>
12#include <kernel/options.h>
13#include <omalloc/omalloc.h>
14#include <kernel/polys.h>
15#include <kernel/febase.h>
16#include <kernel/kstd1.h>
17#include <kernel/kutil.h>
18#include <kernel/stairc.h>
19#include <kernel/intvec.h>
20#include <kernel/numbers.h>
21#include <kernel/ideals.h>
22#include <kernel/intvec.h>
23#include <kernel/ring.h>
24#include <kernel/syz.h>
25#include <kernel/prCopy.h>
26
27#include <kernel/sca.h>
28
29void syDeleteRes(resolvente * res,int length)
30{
31  for (int i=0;i<length;i++)
32  {
33    if (!idIs0((*res)[i]))
34      idDelete(&((*res)[i]));
35  }
36  omFreeSize((ADDRESS)res,length*sizeof(ideal));
37  *res = NULL;
38}
39
40static intvec * syPrepareModComp(ideal arg,intvec ** w)
41{
42  intvec *w1 = NULL;
43  int i;
44  BOOLEAN isIdeal=FALSE;
45
46  if ((w==NULL) || (*w==NULL)) return w1;
47  int maxxx = (*w)->length();
48  if (maxxx==1)
49  {
50    maxxx = 2;
51    isIdeal = TRUE;
52  }
53  w1 = new intvec(maxxx+IDELEMS(arg));
54  if (!isIdeal)
55  {
56    for (i=0;i<maxxx;i++)
57    {
58      (*w1)[i] = (**w)[i];
59    }
60  }
61  for (i=maxxx;i<maxxx+IDELEMS(arg);i++)
62  {
63    if (arg->m[i-maxxx]!=NULL)
64    {
65      (*w1)[i] = pFDeg(arg->m[i-maxxx],currRing);
66      if (pGetComp(arg->m[i-maxxx])!=0)
67      {
68        (*w1)[i]+=(**w)[pGetComp(arg->m[i-maxxx])-1];
69      }
70    }
71  }
72  delete (*w);
73  *w = new intvec(IDELEMS(arg)+1);
74  for (i=0;i<IDELEMS(arg);i++)
75  {
76     (**w)[i+1] = (*w1)[i+maxxx];
77  }
78  return w1;
79}
80
81static void syDeleteAbove(ideal up, int k)
82{
83  if (up!=NULL)
84  {
85    for (int i=0;i<IDELEMS(up);i++)
86    {
87      if (up->m[i]!=NULL)
88        pDeleteComp(&(up->m[i]),k+1);
89    }
90  }
91}
92
93/*2
94*minimizes the module mod and cancel superfluous syzygies
95*from syz
96*/
97static void syMinStep(ideal mod,ideal syz,BOOLEAN final=FALSE,ideal up=NULL,
98                      tHomog h=isNotHomog)
99{
100  ideal deg0=NULL;
101  poly Unit1,Unit2,actWith;
102  int len,i,j,ModComp,m,k,l;
103  BOOLEAN searchUnit,existsUnit;
104
105  if (TEST_OPT_PROT) PrintS("m");
106  if ((final) && (h==isHomog))
107  /*minim is TRUE, we are in the module: maxlength, maxlength <>0*/
108  {
109    deg0=idJet(syz,0);
110    idSkipZeroes(deg0);
111    syz=deg0;
112  }
113/*--cancels empty entees and their related components above--*/
114  j = IDELEMS(syz);
115  while ((j>0) && (!syz->m[j-1])) j--;
116  k = 0;
117  while (k<j)
118  {
119    if (syz->m[k]!=NULL)
120      k++;
121    else
122    {
123      if (TEST_OPT_PROT) PrintS(".");
124      for (l=k;l<j-1;l++) syz->m[l] = syz->m[l+1];
125      syz->m[j-1] = NULL;
126      syDeleteAbove(up,k);
127      j--;
128    }
129  }
130/*--searches for syzygies coming from superfluous elements
131* in the module below--*/
132  searchUnit = TRUE;
133  int curr_syz_limit = rGetCurrSyzLimit();
134  while (searchUnit)
135  {
136    i=0;
137    j=IDELEMS(syz);
138    while ((j>0) && (syz->m[j-1]==NULL)) j--;
139    existsUnit = FALSE;
140    if (rHasGlobalOrdering_currRing())
141    {
142      while ((i<j) && (!existsUnit))
143      {
144        existsUnit = pVectorHasUnitB(syz->m[i],&ModComp);
145        i++;
146      }
147    }
148    else
149    {
150      int I=0;
151      l = 0;
152      len=0;
153      for (i=0;i<IDELEMS(syz);i++)
154      {
155        if (syz->m[i]!=NULL)
156        {
157          pVectorHasUnit(syz->m[i],&m, &l);
158          if ((len==0) ||((l>0) && (l<len)))
159          {
160            len = l;
161            ModComp = m;
162            I = i;
163          }
164        }
165      }
166//Print("Laenge ist: %d\n",len);
167      if (len>0) existsUnit = TRUE;
168      i = I+1;
169    }
170    if (existsUnit)
171    {
172      i--;
173//--takes out the founded syzygy--
174      if (TEST_OPT_PROT) PrintS("f");
175      actWith = syz->m[i];
176      if (!rField_has_simple_inverse()) p_Cleardenom(actWith, currRing);
177//Print("actWith: ");pWrite(actWith);
178      syz->m[i] = NULL;
179      for (k=i;k<j-1;k++)  syz->m[k] = syz->m[k+1];
180      syz->m[j-1] = NULL;
181      syDeleteAbove(up,i);
182      j--;
183//--makes Gauss alg. for the column ModComp--
184      Unit1 = pTakeOutComp(&(actWith), ModComp);
185//PrintS("actWith now: ");pWrite(actWith);
186//Print("Unit1: ");pWrite(Unit1);
187      k=0;
188//Print("j= %d",j);
189      while (k<j)
190      {
191        if (syz->m[k]!=NULL)
192        {
193          Unit2 = pTakeOutComp(&(syz->m[k]), ModComp);
194//Print("element %d: ",k);pWrite(syz->m[k]);
195//PrintS("Unit2: ");pWrite(Unit2);
196          syz->m[k] = pMult(pCopy(Unit1),syz->m[k]);
197          syz->m[k] = pSub(syz->m[k],
198            pMult(Unit2,pCopy(actWith)));
199          if (syz->m[k]==NULL)
200          {
201            for (l=k;l<j-1;l++)
202              syz->m[l] = syz->m[l+1];
203            syz->m[j-1] = NULL;
204            j--;
205            syDeleteAbove(up,k);
206            k--;
207          }
208        }
209        k++;
210      }
211      pDelete(&actWith);
212      pDelete(&Unit1);
213//--deletes superfluous elements from the module below---
214      pDelete(&(mod->m[ModComp-1 - curr_syz_limit]));
215      for (k=ModComp-1 - curr_syz_limit;k<IDELEMS(mod)-1;k++)
216        mod->m[k] = mod->m[k+1];
217      mod->m[IDELEMS(mod)-1] = NULL;
218    }
219    else
220      searchUnit = FALSE;
221  }
222  if (TEST_OPT_PROT) PrintLn();
223  idSkipZeroes(mod);
224  idSkipZeroes(syz);
225  if (deg0!=NULL)
226    idDelete(&deg0);
227}
228
229/*2
230* make Gauss with the element elnum in the module component ModComp
231* for the generators from - till
232*/
233void syGaussForOne(ideal syz, int elnum, int ModComp,int from,int till)
234{
235  int k,j,i,lu;
236  poly unit1,unit2;
237  poly actWith=syz->m[elnum];
238
239  if (from<0) from = 0;
240  if ((till<=0) || (till>IDELEMS(syz))) till = IDELEMS(syz);
241  syz->m[elnum] = NULL;
242  if (!rField_has_simple_inverse()) p_Cleardenom(actWith, currRing);
243/*--makes Gauss alg. for the column ModComp--*/
244  pTakeOutComp(&(actWith), ModComp, &unit1, &lu);
245  while (from<till)
246  {
247    poly tmp=syz->m[from];
248    if (/*syz->m[from]*/ tmp!=NULL)
249    {
250      pTakeOutComp(&(tmp), ModComp, &unit2, &lu);
251      tmp = pMult(pCopy(unit1),tmp);
252      syz->m[from] = pSub(tmp,
253        pMult(unit2,pCopy(actWith)));
254    }
255    from++;
256  }
257  pDelete(&actWith);
258  pDelete(&unit1);
259}
260static void syDeleteAbove1(ideal up, int k)
261{
262  poly p,pp;
263  if (up!=NULL)
264  {
265    for (int i=0;i<IDELEMS(up);i++)
266    {
267      p = up->m[i];
268      while ((p!=NULL) && (pGetComp(p)==k))
269      {
270        /*
271        pp = pNext(p);
272        pNext(p) = NULL;
273        pDelete(&p);
274        p = pp;
275        */
276        pLmDelete(&p);
277      }
278      up->m[i] = p;
279      if (p!=NULL)
280      {
281        while (pNext(p)!=NULL)
282        {
283          if (pGetComp(pNext(p))==k)
284          {
285            /*
286            pp = pNext(pNext(p));
287            pNext(pNext(p)) = NULL;
288            pDelete(&pNext(p));
289            pNext(p) = pp;
290            */
291            pLmDelete(&pNext(p));
292          }
293          else
294            pIter(p);
295        }
296      }
297    }
298  }
299}
300/*2
301*minimizes the resolution res
302*assumes homogeneous or local case
303*/
304static void syMinStep1(resolvente res, int length)
305{
306  int i,j,k,l,index=0;
307  poly p;
308  ideal deg0=NULL,reddeg0=NULL;
309  intvec *have_del=NULL,*to_del=NULL;
310
311  while ((index<length) && (res[index]!=NULL))
312  {
313/*---we take out dependend elements from syz---------------------*/
314    if (res[index+1]!=NULL)
315    {
316      deg0 = idJet(res[index+1],0);
317      reddeg0 = kInterRedOld(deg0);
318      idDelete(&deg0);
319      have_del = new intvec(IDELEMS(res[index]));
320      for (i=0;i<IDELEMS(reddeg0);i++)
321      {
322        if (reddeg0->m[i]!=NULL)
323        {
324          j = pGetComp(reddeg0->m[i]);
325          pDelete(&(res[index]->m[j-1]));
326          /*res[index]->m[j-1] = NULL;*/
327          (*have_del)[j-1] = 1;
328        }
329      }
330      idDelete(&reddeg0);
331    }
332    if (index>0)
333    {
334/*--- we search for units and perform Gaussian elimination------*/
335      j = to_del->length();
336      while (j>0)
337      {
338        if ((*to_del)[j-1]==1)
339        {
340          k = 0;
341          while (k<IDELEMS(res[index]))
342          {
343            p = res[index]->m[k];
344            while ((p!=NULL) && ((!pLmIsConstantComp(p)) || (pGetComp(p)!=j)))
345              pIter(p);
346            if ((p!=NULL) && (pLmIsConstantComp(p)) && (pGetComp(p)==j)) break;
347            k++;
348          }
349          if (k>=IDELEMS(res[index]))
350          {
351            PrintS("out of range\n");
352          }
353          syGaussForOne(res[index],k,j);
354          if (res[index+1]!=NULL)
355            syDeleteAbove1(res[index+1],k+1);
356          (*to_del)[j-1] = 0;
357        }
358        j--;
359      }
360    }
361    if (to_del!=NULL) delete to_del;
362    to_del = have_del;
363    have_del = NULL;
364    index++;
365  }
366  if (TEST_OPT_PROT) PrintLn();
367  syKillEmptyEntres(res,length);
368}
369
370void syMinimizeResolvente(resolvente res, int length, int first)
371{
372  int syzIndex=first;
373  intvec *dummy;
374
375  if (syzIndex<1) syzIndex=1;
376  if ((syzIndex==1) && (idHomModule(res[0],currQuotient,&dummy)) && (!rIsPluralRing(currRing)))
377  {
378    syMinStep1(res,length);
379    delete dummy;
380    return;
381  }
382  while ((syzIndex<length-1) && (res[syzIndex]!=NULL) && (res[syzIndex+1]!=NULL))
383  {
384    syMinStep(res[syzIndex-1],res[syzIndex],FALSE,res[syzIndex+1]);
385    syzIndex++;
386  }
387  if (res[syzIndex]!=NULL)
388    syMinStep(res[syzIndex-1],res[syzIndex]);
389  if (!idIs0(res[0]))
390    idMinEmbedding(res[0],TRUE);
391}
392
393/*2
394* resolution of ideal/module arg, <=maxlength steps, (r[0..maxlength])
395*   no limitation in length if maxlength==0
396* input:arg
397*       minim: TRUE means mres cmd, FALSE nres cmd.
398*       if *len!=0: module weights: weights[0]
399*          (and weights is defined:weights[0..len-1]
400*
401* output:resolvente r[0..length-1],
402*        module weights: weights[0..length-1]
403*/
404resolvente syResolvente(ideal arg, int maxlength, int * length,
405                        intvec *** weights, BOOLEAN minim)
406{
407  resolvente res;
408  resolvente newres;
409  tHomog hom=isNotHomog;
410  ideal temp=NULL;
411  intvec *w = NULL,**tempW;
412  int i,k,syzIndex = 0,j,rk_arg=si_max(1,(int)idRankFreeModule(arg));
413  int Kstd1_OldDeg=Kstd1_deg;
414  BOOLEAN completeMinim;
415  BOOLEAN oldDegBound=TEST_OPT_DEGBOUND;
416  BOOLEAN setRegularity=TRUE;
417  int wlength=*length;
418
419  if (maxlength!=-1) *length = maxlength+1;
420  else              *length = 5;
421  if ((wlength!=0) && (*length!=wlength))
422  {
423    intvec **wtmp = (intvec**)omAlloc0((*length)*sizeof(intvec*));
424    wtmp[0]=(*weights)[0];
425    omFreeSize((ADDRESS)*weights,wlength*sizeof(intvec*));
426    *weights=wtmp;
427  }
428  res = (resolvente)omAlloc0((*length)*sizeof(ideal));
429
430/*--- initialize the syzygy-ring -----------------------------*/
431  ring origR = currRing;
432  ring syz_ring = rCurrRingAssure_SyzComp();
433  rSetSyzComp(rk_arg);
434
435  if (syz_ring != origR)
436  {
437    res[0] = idrCopyR_NoSort(arg, origR);
438  }
439  else
440  {
441    res[0] = idCopy(arg);
442  }
443
444/*--- creating weights for the module components ---------------*/
445  if ((weights!=NULL) && (*weights!=NULL)&& ((*weights)[0]!=NULL))
446  {
447    if (!idTestHomModule(res[0],currQuotient,(*weights)[0]))
448    {
449      WarnS("wrong weights given(1):"); (*weights)[0]->show();PrintLn();
450      idHomModule(res[0],currQuotient,&w);
451      w->show();PrintLn();
452      *weights=NULL;
453    }
454  }
455
456  if ((weights==NULL) || (*weights==NULL) || ((*weights)[0]==NULL))
457  {
458    hom=(tHomog)idHomModule(res[0],currQuotient,&w);
459    if (hom==isHomog)
460    {
461      *weights = (intvec**)omAlloc0((*length)*sizeof(intvec*));
462      if (w!=NULL) (*weights)[0] = ivCopy(w);
463    }
464  }
465  else
466  {
467    if ((weights!=NULL) && (*weights!=NULL)&& ((*weights)[0]!=NULL))
468    {
469      w = ivCopy((*weights)[0]);
470      hom = isHomog;
471    }
472  }
473
474#ifdef HAVE_PLURAL
475  if (rIsPluralRing(currRing) && !rIsSCA(currRing) )
476  { 
477// quick solution; need theory to apply homog GB stuff for G-Algebras
478    hom = isNotHomog;
479  }
480#endif // HAVE_PLURAL
481
482  if (hom==isHomog)
483  {
484    intvec *w1 = syPrepareModComp(res[0],&w);
485    if (w!=NULL) { delete w;w=NULL; }
486    w = w1;
487    j = 0;
488    while ((j<IDELEMS(res[0])) && (res[0]->m[j]==NULL)) j++;
489    if (j<IDELEMS(res[0]))
490    {
491      if (pFDeg(res[0]->m[j],currRing)!=pTotaldegree(res[0]->m[j]))
492        setRegularity = FALSE;
493    }
494  }
495  else
496  {
497    setRegularity = FALSE;
498  }
499
500/*--- the main loop --------------------------------------*/
501  while ((!idIs0(res[syzIndex])) &&
502         ((maxlength==-1) || (syzIndex<=maxlength)))
503   // (syzIndex<maxlength+(int)minim)))
504/*--- compute one step more for minimizing-----------------*/
505  {
506    if (Kstd1_deg!=0) Kstd1_deg++;
507    if (syzIndex+1==*length)
508    {
509      newres = (resolvente)omAlloc0((*length+5)*sizeof(ideal));
510      tempW = (intvec**)omAlloc0((*length+5)*sizeof(intvec*));
511      for (j=0;j<*length;j++)
512      {
513        newres[j] = res[j];
514        if (*weights!=NULL) tempW[j] = (*weights)[j];
515        /*else              tempW[j] = NULL;*/
516      }
517      omFreeSize((ADDRESS)res,*length*sizeof(ideal));
518      if (*weights != NULL) omFreeSize((ADDRESS)*weights,*length*sizeof(intvec*));
519      *length += 5;
520      res=newres;
521      *weights = tempW;
522    }
523/*--- interreducing first -----------------------------------*/
524    if (syzIndex>0)
525    {
526      int rkI=idRankFreeModule(res[syzIndex]);
527      rSetSyzComp(rkI);
528    }
529    if (minim || (syzIndex!=0))
530    {
531      temp = kInterRedOld(res[syzIndex],currQuotient);
532      idDelete(&res[syzIndex]);
533      idSkipZeroes(temp);
534      res[syzIndex] = temp;
535    }
536    temp = NULL;
537/*--- computing the syzygy modules --------------------------------*/
538    if ((currQuotient==NULL)&&(syzIndex==0)&& (!TEST_OPT_DEGBOUND))
539    {
540      res[/*syzIndex+*/1] = idSyzygies(res[0/*syzIndex*/],hom,&w,FALSE,setRegularity,&Kstd1_deg);
541      if ((!TEST_OPT_NOTREGULARITY) && (Kstd1_deg>0)) test |= Sy_bit(OPT_DEGBOUND);
542    }
543    else
544    {
545        res[syzIndex+1] = idSyzygies(res[syzIndex],hom,&w,FALSE);
546    }
547    completeMinim=(syzIndex!=maxlength) || (maxlength ==-1) || (hom!=isHomog);
548    syzIndex++;
549    if (TEST_OPT_PROT) Print("[%d]\n",syzIndex);
550    if ((minim)||(syzIndex>1))
551      syMinStep(res[syzIndex-1],res[syzIndex],!completeMinim,NULL,hom);
552    if (!completeMinim)
553    /*minim is TRUE, we are in the module: maxlength, maxlength <>0*/
554    {
555      idDelete(&res[syzIndex]);
556    }
557/*---creating the iterated weights for module components ---------*/
558    if ((hom == isHomog) && (!idIs0(res[syzIndex])))
559    {
560//Print("die %d Modulegewichte sind:\n",w1->length());
561//w1->show();
562//PrintLn();
563      int max_comp = idRankFreeModule(res[syzIndex]);
564      k = max_comp - rGetCurrSyzLimit();
565      assume(w != NULL);
566      if (w != NULL)
567        w->resize(max_comp+IDELEMS(res[syzIndex]));
568      else
569        w = new intvec(max_comp+IDELEMS(res[syzIndex]));
570      (*weights)[syzIndex] = new intvec(k);
571      for (i=0;i<k;i++)
572      {
573        if (res[syzIndex-1]->m[i]!=NULL) // hs
574        {
575          (*w)[i + rGetCurrSyzLimit()] = pFDeg(res[syzIndex-1]->m[i],currRing);
576          if (pGetComp(res[syzIndex-1]->m[i])>0)
577            (*w)[i + rGetCurrSyzLimit()]
578              += (*w)[pGetComp(res[syzIndex-1]->m[i])-1];
579          (*((*weights)[syzIndex]))[i] = (*w)[i+rGetCurrSyzLimit()];
580        }
581      }
582      for (i=k;i<k+IDELEMS(res[syzIndex]);i++)
583      {
584        if (res[syzIndex]->m[i-k]!=NULL)
585          (*w)[i+rGetCurrSyzLimit()] = pFDeg(res[syzIndex]->m[i-k],currRing)
586                    +(*w)[pGetComp(res[syzIndex]->m[i-k])-1];
587      }
588    }
589  }
590/*--- end of the main loop --------------------------------------*/
591/*--- deleting the temporare data structures --------------------*/
592  if ((syzIndex!=0) && (res[syzIndex]!=NULL) && (idIs0(res[syzIndex])))
593    idDelete(&res[syzIndex]);
594  if (w !=NULL) delete w;
595
596  Kstd1_deg=Kstd1_OldDeg;
597  if (!oldDegBound)
598    test &= ~Sy_bit(OPT_DEGBOUND);
599
600  for (i=1; i<=syzIndex; i++)
601  {
602    if (! idIs0(res[i]))
603    {
604      for (j=0; j<IDELEMS(res[i]); j++)
605      {
606        pShift(&res[i]->m[j], -rGetMaxSyzComp(i));
607      }
608    }
609  }
610/*--- going back to the original ring -------------------------*/
611  if (origR != syz_ring)
612  {
613    rChangeCurrRing(origR);
614    for (i=0; i<=syzIndex; i++)
615    {
616      res[i] = idrMoveR_NoSort(res[i], syz_ring);
617    }
618    rKill(syz_ring);
619  }
620  return res;
621}
622
623syStrategy syResolution(ideal arg, int maxlength,intvec * w, BOOLEAN minim)
624{
625
626#ifdef HAVE_PLURAL
627
628  const ideal idSaveCurrQuotient = currQuotient;
629  const ideal idSaveCurrRingQuotient = currRing->qideal;
630
631  if( rIsSCA(currRing) )
632  {
633
634#ifdef RDEBUG
635//    rWrite(currRing);
636//    rDebugPrint(currRing);
637#endif
638
639    if( ncExtensions(TESTSYZSCAMASK) )
640    {
641      currQuotient = SCAQuotient(currRing);
642      currRing->qideal = currQuotient;
643    }
644
645    const unsigned int m_iFirstAltVar = scaFirstAltVar(currRing);
646    const unsigned int m_iLastAltVar  = scaLastAltVar(currRing);
647   
648    arg = id_KillSquares(arg, m_iFirstAltVar, m_iLastAltVar, currRing, false); // kill suares in input!
649  }
650#endif
651 
652  int typ0;
653  syStrategy result=(syStrategy)omAlloc0(sizeof(ssyStrategy));
654
655  if ((w!=NULL) && (!idTestHomModule(arg,currQuotient,w))) // is this right in SCA case???
656  {
657    WarnS("wrong weights given(2):");w->show();PrintLn();
658    idHomModule(arg,currQuotient,&w);
659    w->show();PrintLn();
660    w=NULL;
661  }
662  if (w!=NULL)
663  {
664    result->weights = (intvec**)omAlloc0Bin(char_ptr_bin);
665    (result->weights)[0] = ivCopy(w);
666    result->length = 1;
667  }
668  resolvente fr = syResolvente(arg,maxlength,&(result->length),&(result->weights),minim),fr1;
669  if (minim)
670  {
671    result->minres = (resolvente)omAlloc0((result->length+1)*sizeof(ideal));
672    fr1 =  result->minres;
673  }
674  else
675  {
676    result->fullres = (resolvente)omAlloc0((result->length+1)*sizeof(ideal));
677    fr1 =  result->fullres;
678  }
679  for (int i=result->length-1;i>=0;i--)
680  {
681    if (fr[i]!=NULL)
682      fr1[i] = fr[i];
683    fr[i] = NULL;
684  }
685  omFreeSize((ADDRESS)fr,(result->length)*sizeof(ideal));
686
687
688#ifdef HAVE_PLURAL
689  if( rIsSCA(currRing) )
690  {
691
692    if( ncExtensions(TESTSYZSCAMASK) )
693    {
694      currQuotient     = idSaveCurrQuotient; 
695      currRing->qideal = idSaveCurrRingQuotient;
696    }
697
698    id_Delete(&arg, currRing);
699  }
700#endif
701
702
703  return result;
704}
705
706static poly sypCopyConstant(poly inp)
707{
708  poly outp=NULL,q;
709
710  while (inp!=NULL)
711  {
712    if (pLmIsConstantComp(inp))
713    {
714      if (outp==NULL)
715      {
716        q = outp = pHead(inp);
717      }
718      else
719      {
720        pNext(q) = pHead(inp);
721        pIter(q);
722      }
723    }
724    pIter(inp);
725  }
726  return outp;
727}
728int syDetect(ideal id,int index,BOOLEAN homog,int * degrees,int * tocancel)
729{
730  int i, j, k, ModComp,subFromRank=0, lu;
731  poly p, q, qq, Unit1, Unit2;
732  ideal temp;
733
734  if (idIs0(id)) return 0;
735  temp = idInit(IDELEMS(id),id->rank);
736  for (i=0;i<IDELEMS(id);i++)
737  {
738    temp->m[i] = sypCopyConstant(id->m[i]);
739  }
740  i = IDELEMS(id);
741  while ((i>0) && (temp->m[i-1]==NULL)) i--;
742  if (i==0)
743  {
744    idDelete(&temp);
745    return 0;
746  }
747  j = 0;
748  p = NULL;
749  while ((j<i) && (temp->m[j]==NULL)) j++;
750  while (j<i)
751  {
752    if (homog)
753    {
754      if (index==0) k = pFDeg(temp->m[j],currRing)+degrees[pGetComp(temp->m[j])];
755      else          k = degrees[pGetComp(temp->m[j])];
756      if (k>=index) tocancel[k-index]++;
757      if ((k>=0) && (index==0)) subFromRank++;
758    }
759    else
760    {
761      tocancel[0]--;
762    }
763    syGaussForOne(temp,j,pGetComp(temp->m[j]),j+1,i);
764    j++;
765    while ((j<i) && (temp->m[j]==NULL)) j++;
766  }
767  idDelete(&temp);
768  return subFromRank;
769}
770
771void syDetect(ideal id,int index,int rsmin, BOOLEAN homog,
772              intvec * degrees,intvec * tocancel)
773{
774  int * deg=NULL;
775  int * tocan=(int*) omAlloc0(tocancel->length()*sizeof(int));
776  int i;
777
778  if (homog)
779  {
780    deg = (int*) omAlloc0(degrees->length()*sizeof(int));
781    for (i=degrees->length();i>0;i--)
782      deg[i-1] = (*degrees)[i-1]-rsmin;
783  }
784  int dummy=syDetect(id,index,homog,deg,tocan);
785  for (i=tocancel->length();i>0;i--)
786    (*tocancel)[i-1] = tocan[i-1];
787  if (homog)
788    omFreeSize((ADDRESS)deg,degrees->length()*sizeof(int));
789  omFreeSize((ADDRESS)tocan,tocancel->length()*sizeof(int));
790}
791
792/*2
793* computes the betti numbers from a given resolution
794* of length 'length' (0..length-1), not necessairily minimal,
795* (if weights are given, they are used)
796* returns the int matrix of betti numbers
797* and the regularity
798*/
799intvec * syBetti(resolvente res,int length, int * regularity,
800                 intvec* weights,BOOLEAN tomin,int * row_shift)
801{
802//#define BETTI_WITH_ZEROS
803  //tomin = FALSE;
804  int i,j=0,k=0,l,rows,cols,mr;
805  int *temp1,*temp2,*temp3;/*used to compute degrees*/
806  int *tocancel; /*(BOOLEAN)tocancel[i]=element is superfluous*/
807  int r0_len;
808
809  /*------ compute size --------------*/
810  *regularity = -1;
811  cols = length;
812  while ((cols>0)
813  && ((res[cols-1]==NULL)
814    || (idIs0(res[cols-1]))))
815  {
816    cols--;
817  }
818  intvec * result;
819  if (idIs0(res[0]))
820  {
821    if (res[0]==NULL)
822      result = new intvec(1,1,1);
823    else
824      result = new intvec(1,1,res[0]->rank);
825    return result;
826  }
827  intvec *w=NULL;
828  if (weights!=NULL)
829  {
830    if (!idTestHomModule(res[0],currQuotient,weights))
831    {
832      WarnS("wrong weights given(3):");weights->show();PrintLn();
833      idHomModule(res[0],currQuotient,&w);
834      if (w!=NULL) { w->show();PrintLn();}
835      weights=NULL;
836    }
837  }
838#if 0
839  if (idHomModule(res[0],currQuotient,&w)!=isHomog)
840  {
841    Warn("betti-command: Input is not homogeneous!");
842    weights=NULL;
843  }
844#endif
845  if (weights==NULL) weights=w;
846  else delete w;
847  r0_len=IDELEMS(res[0]);
848  while ((r0_len>0) && (res[0]->m[r0_len-1]==NULL)) r0_len--;
849  #ifdef SHOW_W
850  PrintS("weights:");if (weights!=NULL) weights->show(); else Print("NULL"); PrintLn();
851  #endif
852  int rkl=l = si_max(idRankFreeModule(res[0]),res[0]->rank);
853  i = 0;
854  while ((i<length) && (res[i]!=NULL))
855  {
856    if (IDELEMS(res[i])>l) l = IDELEMS(res[i]);
857    i++;
858  }
859  temp1 = (int*)omAlloc0((l+1)*sizeof(int));
860  temp2 = (int*)omAlloc((l+1)*sizeof(int));
861  rows = 1;
862  mr = 1;
863  cols++;
864  for (i=0;i<cols-1;i++)
865  {
866    if ((i==0) && (weights!=NULL)) pSetModDeg(weights);
867    memset(temp2,0,(l+1)*sizeof(int));
868    for (j=0;j<IDELEMS(res[i]);j++)
869    {
870      if (res[i]->m[j]!=NULL)
871      {
872        if ((pGetComp(res[i]->m[j])>l)
873        || ((i>1) && (res[i-1]->m[pGetComp(res[i]->m[j])-1]==NULL)))
874        {
875          WerrorS("input not a resolvent");
876          omFreeSize((ADDRESS)temp1,(l+1)*sizeof(int));
877          omFreeSize((ADDRESS)temp2,(l+1)*sizeof(int));
878          return NULL;
879        }
880        temp2[j+1] = pFDeg(res[i]->m[j],currRing)+temp1[pGetComp(res[i]->m[j])];
881        if (temp2[j+1]-i>rows) rows = temp2[j+1]-i;
882        if (temp2[j+1]-i<mr) mr = temp2[j+1]-i;
883      }
884    }
885    if ((i==0) && (weights!=NULL)) pSetModDeg(NULL);
886    temp3 = temp1;
887    temp1 = temp2;
888    temp2 = temp3;
889  }
890  mr--;
891  if (weights!=NULL)
892  {
893    for(j=0;j<weights->length();j++)
894    {
895      if (rows <(*weights)[j]+1) rows=(-mr)+(*weights)[j]+1;
896    }
897  }
898  /*------ computation betti numbers --------------*/
899  rows -= mr;
900  result = new intvec(rows+1,cols,0);
901  if (weights!=NULL)
902  {
903    for(j=0;j<weights->length();j++)
904    {
905      IMATELEM((*result),(-mr)+(*weights)[j]+1,1) ++;
906      //Print("imat(%d,%d)++ -> %d\n",(-mr)+(*weights)[j]+1, 1, IMATELEM((*result),(-mr)+(*weights)[j]+1,1));
907    }
908  }
909  else
910  {
911    (*result)[(-mr)*cols] = /*idRankFreeModule(res[0])*/ rkl;
912    if ((!idIs0(res[0])) && ((*result)[(-mr)*cols]==0))
913      (*result)[(-mr)*cols] = 1;
914  }
915  tocancel = (int*)omAlloc0((rows+1)*sizeof(int));
916  memset(temp1,0,(l+1)*sizeof(int));
917  if (weights!=NULL)
918  {
919    memset(temp2,0,l*sizeof(int));
920    pSetModDeg(weights);
921  }
922  else
923    memset(temp2,0,l*sizeof(int));
924  int dummy = syDetect(res[0],0,TRUE,temp2,tocancel);
925  if (weights!=NULL) pSetModDeg(NULL);
926  if (tomin)
927  {
928    //(*result)[(-mr)*cols] -= dummy;
929    for(j=0;j<=rows+mr;j++)
930    {
931      //Print("tocancel[%d]=%d imat(%d,%d)=%d\n",j,tocancel[j],(-mr)+j+1,1,IMATELEM((*result),(-mr)+j+1,1));
932      IMATELEM((*result),(-mr)+j+1,1) -= tocancel[j];
933    }
934  }
935  for (i=0;i<cols-1;i++)
936  {
937    if ((i==0) && (weights!=NULL)) pSetModDeg(weights);
938    memset(temp2,0,l*sizeof(int));
939    for (j=0;j<IDELEMS(res[i]);j++)
940    {
941      if (res[i]->m[j]!=NULL)
942      {
943        temp2[j+1] = pFDeg(res[i]->m[j],currRing)+temp1[pGetComp(res[i]->m[j])];
944        //(*result)[i+1+(temp2[j+1]-i-1)*cols]++;
945        //if (temp2[j+1]>i) IMATELEM((*result),temp2[j+1]-i-mr,i+2)++;
946        IMATELEM((*result),temp2[j+1]-i-mr,i+2)++;
947      }
948      else if (i==0)
949      {
950        if (j<r0_len) IMATELEM((*result),-mr,2)++;
951      }
952    }
953  /*------ computation betti numbers, if res not minimal --------------*/
954    if (tomin)
955    {
956      for (j=mr;j<rows+mr;j++)
957      {
958        //(*result)[i+1+j*cols] -= tocancel[j+1];
959        IMATELEM((*result),j+1-mr,i+2) -= tocancel[j+1];
960      }
961      if ((i<length-1) && (res[i+1]!=NULL))
962      {
963        memset(tocancel,0,(rows+1)*sizeof(int));
964        dummy = syDetect(res[i+1],i+1,TRUE,temp2,tocancel);
965        for (j=0;j<rows;j++)
966        {
967          //(*result)[i+1+j*cols] -= tocancel[j];
968          IMATELEM((*result),j+1,i+2) -= tocancel[j];
969        }
970      }
971    }
972    temp3 = temp1;
973    temp1 = temp2;
974    temp2 = temp3;
975    for (j=0;j<=rows;j++)
976    {
977    //  if (((*result)[i+1+j*cols]!=0) && (j>*regularity)) *regularity = j;
978      if ((IMATELEM((*result),j+1,i+2)!=0) && (j>*regularity)) *regularity = j;
979    }
980    if ((i==0) && (weights!=NULL)) pSetModDeg(NULL);
981  }
982  // Print("nach minim:\n"); result->show(); PrintLn();
983  /*------ clean up --------------*/
984  omFreeSize((ADDRESS)tocancel,(rows+1)*sizeof(int));
985  omFreeSize((ADDRESS)temp1,(l+1)*sizeof(int));
986  omFreeSize((ADDRESS)temp2,(l+1)*sizeof(int));
987  if ((tomin) && (mr<0))  // deletes the first (zero) line
988  {
989    for (j=1;j<=rows+mr+1;j++)
990    {
991      for (k=1;k<=cols;k++)
992      {
993        IMATELEM((*result),j,k) = IMATELEM((*result),j-mr,k);
994      }
995    }
996    for (j=rows+mr+1;j<=rows+1;j++)
997    {
998      for (k=1;k<=cols;k++)
999      {
1000        IMATELEM((*result),j,k) = 0;
1001      }
1002    }
1003  }
1004  j = 0;
1005  k = 0;
1006  for (i=1;i<=result->rows();i++)
1007  {
1008    for(l=1;l<=result->cols();l++)
1009    if (IMATELEM((*result),i,l) != 0)
1010    {
1011      j = si_max(j, i-1);
1012      k = si_max(k, l-1);
1013    }
1014  }
1015  intvec * exactresult=new intvec(j+1,k+1,0);
1016  for (i=0;i<exactresult->rows();i++)
1017  {
1018    for (j=0;j<exactresult->cols();j++)
1019    {
1020      IMATELEM(*exactresult,i+1,j+1) = IMATELEM(*result,i+1,j+1);
1021    }
1022  }
1023  if (row_shift!=NULL) *row_shift = mr;
1024  delete result;
1025  return exactresult;
1026}
1027
1028/*2
1029* minbare via syzygies
1030*/
1031ideal syMinBase(ideal arg)
1032{
1033  intvec ** weights=NULL;
1034  int leng;
1035  if (idIs0(arg)) return idInit(1,arg->rank);
1036  resolvente res=syResolvente(arg,1,&leng,&weights,TRUE);
1037  ideal result=res[0];
1038  omFreeSize((ADDRESS)res,leng*sizeof(ideal));
1039  if (weights!=NULL)
1040  {
1041    if (*weights!=NULL)
1042    {
1043      delete (*weights);
1044      *weights=NULL;
1045    }
1046    if ((leng>=1) && (*(weights+1)!=NULL))
1047    {
1048      delete *(weights+1);
1049      *(weights+1)=NULL;
1050    }
1051  }
1052  idSkipZeroes(result);
1053  return result;
1054}
1055
1056/*2
1057* computes Betti-numbers from a resolvente of
1058* (non-)homogeneous objects
1059* the numbers of entrees !=NULL in res and weights must be equal
1060* and < length
1061*/
1062intvec * syNewBetti(resolvente res, intvec ** weights, int length)
1063{
1064  intvec * result,*tocancel;
1065  int i,j,k,rsmin=0,rsmax=0,rs=0;
1066  BOOLEAN homog=TRUE;
1067
1068  if (weights!=NULL)           //---homogeneous Betti numbers
1069  {
1070/*--------------computes size of the field----------------------*/
1071    for (i=1;i<length;i++)
1072    {
1073      if (weights[i] != NULL)
1074      {
1075        for (j=1;j<(weights[i])->length();j++)
1076        {
1077          if ((*(weights[i]))[j]-i<rsmin) rsmin = (*(weights[i]))[j]-i;
1078          if ((*(weights[i]))[j]-i>rsmax) rsmax = (*(weights[i]))[j]-i;
1079        }
1080      }
1081    }
1082    i = 0;
1083    while (weights[i] != NULL) i++;
1084    i--;
1085    for (j=0;j<IDELEMS(res[i]);j++)
1086    {
1087      if (res[i]->m[j]!=NULL)
1088      {
1089        k = pFDeg(res[i]->m[j],currRing)+(*(weights[i]))[pGetComp(res[i]->m[j])]-i-1;
1090        if (k>rsmax) rsmax = k;
1091        if (k<rsmin) rsmin = k;
1092      }
1093    }
1094    for (j=1;j<(weights[0])->length();j++)
1095    {
1096      if ((*weights[0])[j]>rsmax) rsmax = (*weights[0])[j];
1097      if ((*weights[0])[j]<rsmin) rsmin = (*weights[0])[j];
1098    }
1099//Print("rsmax = %d\n",rsmax);
1100//Print("rsmin = %d\n",rsmin);
1101    rs = rsmax-rsmin+1;
1102    result = new intvec(rs,i+2,0);
1103    tocancel = new intvec(rs);
1104/*-----------enter the Betti numbers-------------------------------*/
1105    if (/*idRankFreeModule(res[0])*/ res[0]->rank==0)
1106    {
1107      IMATELEM(*result,1-rsmin,1)=1;
1108    }
1109    else
1110    {
1111      for (i=1;i<(weights[0])->length();i++)
1112        IMATELEM(*result,(*weights[0])[i]+1-rsmin,1)++;
1113    }
1114    i = 1;
1115    while (weights[i]!=NULL)
1116    {
1117      for (j=1;j<(weights[i])->length();j++)
1118      {
1119        IMATELEM(*result,(*(weights[i]))[j]-i+1-rsmin,i+1)++;
1120      }
1121      i++;
1122    }
1123    i--;
1124    for (j=0;j<IDELEMS(res[i]);j++)
1125    {
1126      if (res[i]->m[j]!=NULL)
1127      {
1128        k = pFDeg(res[i]->m[j],currRing)+(*(weights[i]))[pGetComp(res[i]->m[j])]-i;
1129        IMATELEM(*result,k-rsmin,i+2)++;
1130      }
1131    }
1132  }
1133  else                //-----the non-homgeneous case
1134  {
1135    homog = FALSE;
1136    tocancel = new intvec(1);
1137    k = length;
1138    while ((k>0) && (idIs0(res[k-1]))) k--;
1139    result = new intvec(1,k+1,0);
1140    (*result)[0] = res[0]->rank;
1141    for (i=0;i<length;i++)
1142    {
1143      if (res[i]!=NULL)
1144      {
1145        for (j=0;j<IDELEMS(res[i]);j++)
1146        {
1147          if (res[i]->m[j]!=NULL) (*result)[i+1]++;
1148        }
1149      }
1150    }
1151  }
1152/*--------computes the Betti numbers for the minimized reolvente----*/
1153
1154  i = 1;
1155  while ((res[i]!=NULL) && (weights[i]!=NULL))
1156  {
1157    syDetect(res[i],i,rsmin,homog,weights[i],tocancel);
1158    if (homog)
1159    {
1160      for (j=0;j<rs-1;j++)
1161      {
1162        IMATELEM((*result),j+1,i+1) -= (*tocancel)[j];
1163        IMATELEM((*result),j+1,i+2) -= (*tocancel)[j+1];
1164      }
1165      IMATELEM((*result),rs,i+1) -= (*tocancel)[rs-1];
1166    }
1167    else
1168    {
1169      (*result)[i+1] -= (*tocancel)[0];
1170      (*result)[i+2] -= (*tocancel)[0];
1171    }
1172    i++;
1173  }
1174
1175/*--------print Betti numbers for control---------------------------*/
1176  for(i=rsmin;i<=rsmax;i++)
1177  {
1178    Print("%2d:",i);
1179    for(j=1;j<=result->cols();j++)
1180    {
1181      Print(" %5d",IMATELEM(*result,i-rsmin+1,j));
1182    }
1183    PrintLn();
1184  }
1185  return result;
1186}
1187
1188/*2
1189* is looking for the minimal minimized module of a resolvente
1190* i.e. returns 0 if res comes from a mres-command and 1 in the
1191* case of res-commands
1192*/
1193int syIsMinimizedFrom(resolvente res,int length)
1194{
1195  poly p;
1196  int i,j=length;
1197
1198  while ((j>0) && (res[j-1]==NULL)) j--;
1199  while (j>0)
1200  {
1201    for (i=0;i<IDELEMS(res[j-1]);i++)
1202    {
1203      p = res[j-1]->m[i];
1204      while (p!=NULL)
1205      {
1206        if (pLmIsConstantComp(p)) return j;
1207        p = pNext(p);
1208      }
1209    }
1210    j--;
1211  }
1212  return j;
1213}
Note: See TracBrowser for help on using the repository browser.