source: git/kernel/GBEngine/syz.cc @ 9f7665

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