source: git/kernel/syz.cc @ 737a68

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