source: git/kernel/syz.cc @ fbc7cb

jengelh-datetimespielwiese
Last change on this file since fbc7cb was ba5e9e, checked in by Oleksandr Motsak <motsak@…>, 9 years ago
Changed configure-scripts to generate individual public config files for each package: resources, libpolys, singular (main) fix: sources should include correct corresponding config headers.
  • 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#ifdef HAVE_CONFIG_H
11#include "singularconfig.h"
12#endif /* HAVE_CONFIG_H */
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/kstd1.h>
19#include <kernel/kutil.h>
20#include <kernel/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/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)) si_opt_1 |= Sy_bit(OPT_DEGBOUND);
537    }
538    else
539    {
540        res[syzIndex+1] = idSyzygies(res[syzIndex],hom,&w,FALSE);
541    }
542    completeMinim=(syzIndex!=maxlength) || (maxlength ==-1) || (hom!=isHomog);
543    syzIndex++;
544    if (TEST_OPT_PROT) Print("[%d]\n",syzIndex);
545
546    if(! TEST_OPT_NO_SYZ_MINIM )
547    {
548      if ((minim)||(syzIndex>1))
549        syMinStep(res[syzIndex-1],res[syzIndex],!completeMinim,NULL,hom);
550      if (!completeMinim)
551      /*minim is TRUE, we are in the module: maxlength, maxlength <>0*/
552      {
553        idDelete(&res[syzIndex]);
554      }
555    }
556/*---creating the iterated weights for module components ---------*/
557    if ((hom == isHomog) && (!idIs0(res[syzIndex])))
558    {
559//Print("die %d Modulegewichte sind:\n",w1->length());
560//w1->show();
561//PrintLn();
562      int max_comp = id_RankFreeModule(res[syzIndex],currRing);
563      k = max_comp - rGetCurrSyzLimit(currRing);
564      assume(w != NULL);
565      if (w != NULL)
566        w->resize(max_comp+IDELEMS(res[syzIndex]));
567      else
568        w = new intvec(max_comp+IDELEMS(res[syzIndex]));
569      (*weights)[syzIndex] = new intvec(k);
570      for (i=0;i<k;i++)
571      {
572        if (res[syzIndex-1]->m[i]!=NULL) // hs
573        {
574          (*w)[i + rGetCurrSyzLimit(currRing)] = p_FDeg(res[syzIndex-1]->m[i],currRing);
575          if (pGetComp(res[syzIndex-1]->m[i])>0)
576            (*w)[i + rGetCurrSyzLimit(currRing)]
577              += (*w)[pGetComp(res[syzIndex-1]->m[i])-1];
578          (*((*weights)[syzIndex]))[i] = (*w)[i+rGetCurrSyzLimit(currRing)];
579        }
580      }
581      for (i=k;i<k+IDELEMS(res[syzIndex]);i++)
582      {
583        if (res[syzIndex]->m[i-k]!=NULL)
584          (*w)[i+rGetCurrSyzLimit(currRing)] = p_FDeg(res[syzIndex]->m[i-k],currRing)
585                    +(*w)[pGetComp(res[syzIndex]->m[i-k])-1];
586      }
587    }
588  }
589/*--- end of the main loop --------------------------------------*/
590/*--- deleting the temporare data structures --------------------*/
591  if ((syzIndex!=0) && (res[syzIndex]!=NULL) && (idIs0(res[syzIndex])))
592    idDelete(&res[syzIndex]);
593  if (w !=NULL) delete w;
594
595  Kstd1_deg=Kstd1_OldDeg;
596  if (!oldDegBound)
597    si_opt_1 &= ~Sy_bit(OPT_DEGBOUND);
598
599  for (i=1; i<=syzIndex; i++)
600  {
601    if (! idIs0(res[i]))
602    {
603      for (j=0; j<IDELEMS(res[i]); j++)
604      {
605        p_Shift(&res[i]->m[j], -rGetMaxSyzComp(i, currRing),currRing);
606      }
607    }
608  }
609/*--- going back to the original ring -------------------------*/
610  if (origR != syz_ring)
611  {
612    rChangeCurrRing(origR); // should not be needed now?
613    for (i=0; i<=syzIndex; i++)
614    {
615      res[i] = idrMoveR_NoSort(res[i], syz_ring, origR);
616    }
617    rDelete(syz_ring);
618  }
619  SI_RESTORE_OPT1(save1);
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  syStrategy result=(syStrategy)omAlloc0(sizeof(ssyStrategy));
653
654  if ((w!=NULL) && (!idTestHomModule(arg,currQuotient,w))) // is this right in SCA case???
655  {
656    WarnS("wrong weights given(2):");w->show();PrintLn();
657    idHomModule(arg,currQuotient,&w);
658    w->show();PrintLn();
659    w=NULL;
660  }
661  if (w!=NULL)
662  {
663    result->weights = (intvec**)omAlloc0Bin(char_ptr_bin);
664    (result->weights)[0] = ivCopy(w);
665    result->length = 1;
666  }
667  resolvente fr = syResolvente(arg,maxlength,&(result->length),&(result->weights),minim),fr1;
668  if (minim)
669  {
670    result->minres = (resolvente)omAlloc0((result->length+1)*sizeof(ideal));
671    fr1 =  result->minres;
672  }
673  else
674  {
675    result->fullres = (resolvente)omAlloc0((result->length+1)*sizeof(ideal));
676    fr1 =  result->fullres;
677  }
678  for (int i=result->length-1;i>=0;i--)
679  {
680    if (fr[i]!=NULL)
681      fr1[i] = fr[i];
682    fr[i] = NULL;
683  }
684  omFreeSize((ADDRESS)fr,(result->length)*sizeof(ideal));
685
686
687#ifdef HAVE_PLURAL
688  if( rIsSCA(currRing) )
689  {
690
691    if( ncExtensions(TESTSYZSCAMASK) )
692    {
693      currQuotient     = idSaveCurrQuotient;
694      currRing->qideal = idSaveCurrRingQuotient;
695    }
696
697    id_Delete(&arg, currRing);
698  }
699#endif
700
701
702  return result;
703}
704
705static poly sypCopyConstant(poly inp)
706{
707  poly outp=NULL,q;
708
709  while (inp!=NULL)
710  {
711    if (pLmIsConstantComp(inp))
712    {
713      if (outp==NULL)
714      {
715        q = outp = pHead(inp);
716      }
717      else
718      {
719        pNext(q) = pHead(inp);
720        pIter(q);
721      }
722    }
723    pIter(inp);
724  }
725  return outp;
726}
727int syDetect(ideal id,int index,BOOLEAN homog,int * degrees,int * tocancel)
728{
729  int i, j, k, subFromRank=0;
730  ideal temp;
731
732  if (idIs0(id)) return 0;
733  temp = idInit(IDELEMS(id),id->rank);
734  for (i=0;i<IDELEMS(id);i++)
735  {
736    temp->m[i] = sypCopyConstant(id->m[i]);
737  }
738  i = IDELEMS(id);
739  while ((i>0) && (temp->m[i-1]==NULL)) i--;
740  if (i==0)
741  {
742    idDelete(&temp);
743    return 0;
744  }
745  j = 0;
746  while ((j<i) && (temp->m[j]==NULL)) j++;
747  while (j<i)
748  {
749    if (homog)
750    {
751      if (index==0) k = p_FDeg(temp->m[j],currRing)+degrees[pGetComp(temp->m[j])];
752      else          k = degrees[pGetComp(temp->m[j])];
753      if (k>=index) tocancel[k-index]++;
754      if ((k>=0) && (index==0)) subFromRank++;
755    }
756    else
757    {
758      tocancel[0]--;
759    }
760    syGaussForOne(temp,j,pGetComp(temp->m[j]),j+1,i);
761    j++;
762    while ((j<i) && (temp->m[j]==NULL)) j++;
763  }
764  idDelete(&temp);
765  return subFromRank;
766}
767
768void syDetect(ideal id,int index,int rsmin, BOOLEAN homog,
769              intvec * degrees,intvec * tocancel)
770{
771  int * deg=NULL;
772  int * tocan=(int*) omAlloc0(tocancel->length()*sizeof(int));
773  int i;
774
775  if (homog)
776  {
777    deg = (int*) omAlloc0(degrees->length()*sizeof(int));
778    for (i=degrees->length();i>0;i--)
779      deg[i-1] = (*degrees)[i-1]-rsmin;
780  }
781  syDetect(id,index,homog,deg,tocan);
782  for (i=tocancel->length();i>0;i--)
783    (*tocancel)[i-1] = tocan[i-1];
784  if (homog)
785    omFreeSize((ADDRESS)deg,degrees->length()*sizeof(int));
786  omFreeSize((ADDRESS)tocan,tocancel->length()*sizeof(int));
787}
788
789/*2
790* computes the betti numbers from a given resolution
791* of length 'length' (0..length-1), not necessairily minimal,
792* (if weights are given, they are used)
793* returns the int matrix of betti numbers
794* and the regularity
795*/
796intvec * syBetti(resolvente res,int length, int * regularity,
797                 intvec* weights,BOOLEAN tomin,int * row_shift)
798{
799//#define BETTI_WITH_ZEROS
800  //tomin = FALSE;
801  int i,j=0,k=0,l,rows,cols,mr;
802  int *temp1,*temp2,*temp3;/*used to compute degrees*/
803  int *tocancel; /*(BOOLEAN)tocancel[i]=element is superfluous*/
804  int r0_len;
805
806  /*------ compute size --------------*/
807  *regularity = -1;
808  cols = length;
809  while ((cols>0)
810  && ((res[cols-1]==NULL)
811    || (idIs0(res[cols-1]))))
812  {
813    cols--;
814  }
815  intvec * result;
816  if (idIs0(res[0]))
817  {
818    if (res[0]==NULL)
819      result = new intvec(1,1,1);
820    else
821      result = new intvec(1,1,res[0]->rank);
822    return result;
823  }
824  intvec *w=NULL;
825  if (weights!=NULL)
826  {
827    if (!idTestHomModule(res[0],currQuotient,weights))
828    {
829      WarnS("wrong weights given(3):");weights->show();PrintLn();
830      idHomModule(res[0],currQuotient,&w);
831      if (w!=NULL) { w->show();PrintLn();}
832      weights=NULL;
833    }
834  }
835#if 0
836  if (idHomModule(res[0],currQuotient,&w)!=isHomog)
837  {
838    Warn("betti-command: Input is not homogeneous!");
839    weights=NULL;
840  }
841#endif
842  if (weights==NULL) weights=w;
843  else delete w;
844  r0_len=IDELEMS(res[0]);
845  while ((r0_len>0) && (res[0]->m[r0_len-1]==NULL)) r0_len--;
846  #ifdef SHOW_W
847  PrintS("weights:");if (weights!=NULL) weights->show(); else Print("NULL"); PrintLn();
848  #endif
849  int rkl=l = si_max(id_RankFreeModule(res[0],currRing),res[0]->rank);
850  i = 0;
851  while ((i<length) && (res[i]!=NULL))
852  {
853    if (IDELEMS(res[i])>l) l = IDELEMS(res[i]);
854    i++;
855  }
856  temp1 = (int*)omAlloc0((l+1)*sizeof(int));
857  temp2 = (int*)omAlloc((l+1)*sizeof(int));
858  rows = 1;
859  mr = 1;
860  cols++;
861  for (i=0;i<cols-1;i++)
862  {
863    if ((i==0) && (weights!=NULL)) p_SetModDeg(weights, currRing);
864    memset(temp2,0,(l+1)*sizeof(int));
865    for (j=0;j<IDELEMS(res[i]);j++)
866    {
867      if (res[i]->m[j]!=NULL)
868      {
869        if ((pGetComp(res[i]->m[j])>(unsigned)l)
870        || ((i>1) && (res[i-1]->m[pGetComp(res[i]->m[j])-1]==NULL)))
871        {
872          WerrorS("input not a resolvent");
873          omFreeSize((ADDRESS)temp1,(l+1)*sizeof(int));
874          omFreeSize((ADDRESS)temp2,(l+1)*sizeof(int));
875          return NULL;
876        }
877        temp2[j+1] = p_FDeg(res[i]->m[j],currRing)+temp1[pGetComp(res[i]->m[j])];
878        if (temp2[j+1]-i>rows) rows = temp2[j+1]-i;
879        if (temp2[j+1]-i<mr) mr = temp2[j+1]-i;
880      }
881    }
882    if ((i==0) && (weights!=NULL)) p_SetModDeg(NULL, currRing);
883    temp3 = temp1;
884    temp1 = temp2;
885    temp2 = temp3;
886  }
887  mr--;
888  if (weights!=NULL)
889  {
890    for(j=0;j<weights->length();j++)
891    {
892      if (rows <(*weights)[j]+1) rows=(-mr)+(*weights)[j]+1;
893    }
894  }
895  /*------ computation betti numbers --------------*/
896  rows -= mr;
897  result = new intvec(rows+1,cols,0);
898  if (weights!=NULL)
899  {
900    for(j=0;j<weights->length();j++)
901    {
902      IMATELEM((*result),(-mr)+(*weights)[j]+1,1) ++;
903      //Print("imat(%d,%d)++ -> %d\n",(-mr)+(*weights)[j]+1, 1, IMATELEM((*result),(-mr)+(*weights)[j]+1,1));
904    }
905  }
906  else
907  {
908    (*result)[(-mr)*cols] = /*idRankFreeModule(res[0])*/ rkl;
909    if ((!idIs0(res[0])) && ((*result)[(-mr)*cols]==0))
910      (*result)[(-mr)*cols] = 1;
911  }
912  tocancel = (int*)omAlloc0((rows+1)*sizeof(int));
913  memset(temp1,0,(l+1)*sizeof(int));
914  if (weights!=NULL)
915  {
916    memset(temp2,0,l*sizeof(int));
917    p_SetModDeg(weights, currRing);
918  }
919  else
920    memset(temp2,0,l*sizeof(int));
921  syDetect(res[0],0,TRUE,temp2,tocancel);
922  if (weights!=NULL) p_SetModDeg(NULL, currRing);
923  if (tomin)
924  {
925    //(*result)[(-mr)*cols] -= dummy;
926    for(j=0;j<=rows+mr;j++)
927    {
928      //Print("tocancel[%d]=%d imat(%d,%d)=%d\n",j,tocancel[j],(-mr)+j+1,1,IMATELEM((*result),(-mr)+j+1,1));
929      IMATELEM((*result),(-mr)+j+1,1) -= tocancel[j];
930    }
931  }
932  for (i=0;i<cols-1;i++)
933  {
934    if ((i==0) && (weights!=NULL)) p_SetModDeg(weights, currRing);
935    memset(temp2,0,l*sizeof(int));
936    for (j=0;j<IDELEMS(res[i]);j++)
937    {
938      if (res[i]->m[j]!=NULL)
939      {
940        temp2[j+1] = p_FDeg(res[i]->m[j],currRing)+temp1[pGetComp(res[i]->m[j])];
941        //(*result)[i+1+(temp2[j+1]-i-1)*cols]++;
942        //if (temp2[j+1]>i) IMATELEM((*result),temp2[j+1]-i-mr,i+2)++;
943        IMATELEM((*result),temp2[j+1]-i-mr,i+2)++;
944      }
945      else if (i==0)
946      {
947        if (j<r0_len) IMATELEM((*result),-mr,2)++;
948      }
949    }
950  /*------ computation betti numbers, if res not minimal --------------*/
951    if (tomin)
952    {
953      for (j=mr;j<rows+mr;j++)
954      {
955        //(*result)[i+1+j*cols] -= tocancel[j+1];
956        IMATELEM((*result),j+1-mr,i+2) -= tocancel[j+1];
957      }
958      if ((i<length-1) && (res[i+1]!=NULL))
959      {
960        memset(tocancel,0,(rows+1)*sizeof(int));
961        syDetect(res[i+1],i+1,TRUE,temp2,tocancel);
962        for (j=0;j<rows;j++)
963        {
964          //(*result)[i+1+j*cols] -= tocancel[j];
965          IMATELEM((*result),j+1,i+2) -= tocancel[j];
966        }
967      }
968    }
969    temp3 = temp1;
970    temp1 = temp2;
971    temp2 = temp3;
972    for (j=0;j<=rows;j++)
973    {
974    //  if (((*result)[i+1+j*cols]!=0) && (j>*regularity)) *regularity = j;
975      if ((IMATELEM((*result),j+1,i+2)!=0) && (j>*regularity)) *regularity = j;
976    }
977    if ((i==0) && (weights!=NULL)) p_SetModDeg(NULL, currRing);
978  }
979  // Print("nach minim:\n"); result->show(); PrintLn();
980  /*------ clean up --------------*/
981  omFreeSize((ADDRESS)tocancel,(rows+1)*sizeof(int));
982  omFreeSize((ADDRESS)temp1,(l+1)*sizeof(int));
983  omFreeSize((ADDRESS)temp2,(l+1)*sizeof(int));
984  if ((tomin) && (mr<0))  // deletes the first (zero) line
985  {
986    for (j=1;j<=rows+mr+1;j++)
987    {
988      for (k=1;k<=cols;k++)
989      {
990        IMATELEM((*result),j,k) = IMATELEM((*result),j-mr,k);
991      }
992    }
993    for (j=rows+mr+1;j<=rows+1;j++)
994    {
995      for (k=1;k<=cols;k++)
996      {
997        IMATELEM((*result),j,k) = 0;
998      }
999    }
1000  }
1001  j = 0;
1002  k = 0;
1003  for (i=1;i<=result->rows();i++)
1004  {
1005    for(l=1;l<=result->cols();l++)
1006    if (IMATELEM((*result),i,l) != 0)
1007    {
1008      j = si_max(j, i-1);
1009      k = si_max(k, l-1);
1010    }
1011  }
1012  intvec * exactresult=new intvec(j+1,k+1,0);
1013  for (i=0;i<exactresult->rows();i++)
1014  {
1015    for (j=0;j<exactresult->cols();j++)
1016    {
1017      IMATELEM(*exactresult,i+1,j+1) = IMATELEM(*result,i+1,j+1);
1018    }
1019  }
1020  if (row_shift!=NULL) *row_shift = mr;
1021  delete result;
1022  return exactresult;
1023}
1024
1025/*2
1026* minbare via syzygies
1027*/
1028ideal syMinBase(ideal arg)
1029{
1030  intvec ** weights=NULL;
1031  int leng;
1032  if (idIs0(arg)) return idInit(1,arg->rank);
1033  resolvente res=syResolvente(arg,1,&leng,&weights,TRUE);
1034  ideal result=res[0];
1035  omFreeSize((ADDRESS)res,leng*sizeof(ideal));
1036  if (weights!=NULL)
1037  {
1038    if (*weights!=NULL)
1039    {
1040      delete (*weights);
1041      *weights=NULL;
1042    }
1043    if ((leng>=1) && (*(weights+1)!=NULL))
1044    {
1045      delete *(weights+1);
1046      *(weights+1)=NULL;
1047    }
1048  }
1049  idSkipZeroes(result);
1050  return result;
1051}
1052
1053#if 0  /* currently used: syBetti */
1054/*2
1055* computes Betti-numbers from a resolvente of
1056* (non-)homogeneous objects
1057* the numbers of entrees !=NULL in res and weights must be equal
1058* and < length
1059*/
1060intvec * syNewBetti(resolvente res, intvec ** weights, int length)
1061{
1062  intvec * result,*tocancel;
1063  int i,j,k,rsmin=0,rsmax=0,rs=0;
1064  BOOLEAN homog=TRUE;
1065
1066  if (weights!=NULL)           //---homogeneous Betti numbers
1067  {
1068/*--------------computes size of the field----------------------*/
1069    for (i=1;i<length;i++)
1070    {
1071      if (weights[i] != NULL)
1072      {
1073        for (j=1;j<(weights[i])->length();j++)
1074        {
1075          if ((*(weights[i]))[j]-i<rsmin) rsmin = (*(weights[i]))[j]-i;
1076          if ((*(weights[i]))[j]-i>rsmax) rsmax = (*(weights[i]))[j]-i;
1077        }
1078      }
1079    }
1080    i = 0;
1081    while (weights[i] != NULL) i++;
1082    i--;
1083    for (j=0;j<IDELEMS(res[i]);j++)
1084    {
1085      if (res[i]->m[j]!=NULL)
1086      {
1087        k = p_FDeg(res[i]->m[j],currRing)+(*(weights[i]))[pGetComp(res[i]->m[j])]-i-1;
1088        if (k>rsmax) rsmax = k;
1089        if (k<rsmin) rsmin = k;
1090      }
1091    }
1092    for (j=1;j<(weights[0])->length();j++)
1093    {
1094      if ((*weights[0])[j]>rsmax) rsmax = (*weights[0])[j];
1095      if ((*weights[0])[j]<rsmin) rsmin = (*weights[0])[j];
1096    }
1097//Print("rsmax = %d\n",rsmax);
1098//Print("rsmin = %d\n",rsmin);
1099    rs = rsmax-rsmin+1;
1100    result = new intvec(rs,i+2,0);
1101    tocancel = new intvec(rs);
1102/*-----------enter the Betti numbers-------------------------------*/
1103    if (/*idRankFreeModule(res[0])*/ res[0]->rank==0)
1104    {
1105      IMATELEM(*result,1-rsmin,1)=1;
1106    }
1107    else
1108    {
1109      for (i=1;i<(weights[0])->length();i++)
1110        IMATELEM(*result,(*weights[0])[i]+1-rsmin,1)++;
1111    }
1112    i = 1;
1113    while (weights[i]!=NULL)
1114    {
1115      for (j=1;j<(weights[i])->length();j++)
1116      {
1117        IMATELEM(*result,(*(weights[i]))[j]-i+1-rsmin,i+1)++;
1118      }
1119      i++;
1120    }
1121    i--;
1122    for (j=0;j<IDELEMS(res[i]);j++)
1123    {
1124      if (res[i]->m[j]!=NULL)
1125      {
1126        k = p_FDeg(res[i]->m[j],currRing)+(*(weights[i]))[pGetComp(res[i]->m[j])]-i;
1127        IMATELEM(*result,k-rsmin,i+2)++;
1128      }
1129    }
1130  }
1131  else                //-----the non-homgeneous case
1132  {
1133    homog = FALSE;
1134    tocancel = new intvec(1);
1135    k = length;
1136    while ((k>0) && (idIs0(res[k-1]))) k--;
1137    result = new intvec(1,k+1,0);
1138    (*result)[0] = res[0]->rank;
1139    for (i=0;i<length;i++)
1140    {
1141      if (res[i]!=NULL)
1142      {
1143        for (j=0;j<IDELEMS(res[i]);j++)
1144        {
1145          if (res[i]->m[j]!=NULL) (*result)[i+1]++;
1146        }
1147      }
1148    }
1149  }
1150/*--------computes the Betti numbers for the minimized reolvente----*/
1151
1152  i = 1;
1153  while ((res[i]!=NULL) && (weights[i]!=NULL))
1154  {
1155    syDetect(res[i],i,rsmin,homog,weights[i],tocancel);
1156    if (homog)
1157    {
1158      for (j=0;j<rs-1;j++)
1159      {
1160        IMATELEM((*result),j+1,i+1) -= (*tocancel)[j];
1161        IMATELEM((*result),j+1,i+2) -= (*tocancel)[j+1];
1162      }
1163      IMATELEM((*result),rs,i+1) -= (*tocancel)[rs-1];
1164    }
1165    else
1166    {
1167      (*result)[i+1] -= (*tocancel)[0];
1168      (*result)[i+2] -= (*tocancel)[0];
1169    }
1170    i++;
1171  }
1172
1173/*--------print Betti numbers for control---------------------------*/
1174  for(i=rsmin;i<=rsmax;i++)
1175  {
1176    Print("%2d:",i);
1177    for(j=1;j<=result->cols();j++)
1178    {
1179      Print(" %5d",IMATELEM(*result,i-rsmin+1,j));
1180    }
1181    PrintLn();
1182  }
1183  return result;
1184}
1185#endif
1186
Note: See TracBrowser for help on using the repository browser.