source: git/kernel/syz.cc @ 76e501

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