source: git/kernel/GBEngine/syz.cc

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