source: git/kernel/syz2.cc @ abe5c8

fieker-DuValspielwiese
Last change on this file since abe5c8 was 6ce030f, checked in by Oleksandr Motsak <motsak@…>, 12 years ago
removal of the $Id$ svn tag from everywhere NOTE: the git SHA1 may be used instead (only on special places) NOTE: the libraries Singular/LIB/*.lib still contain the marker due to our current use of svn
  • Property mode set to 100644
File size: 29.0 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/*
5* ABSTRACT: resolutions
6*/
7
8#include "config.h"
9#include "mod2.h"
10
11#include <omalloc/omalloc.h>
12
13#include <misc/mylimits.h>
14#include <misc/options.h>
15#include <misc/intvec.h>
16//#include "cntrlc.h"
17
18#include <coeffs/coeffs.h>
19#include <coeffs/numbers.h>
20
21#include <polys/monomials/ring.h>
22#include <polys/kbuckets.h>
23#include <polys/prCopy.h>
24
25#include <kernel/polys.h>
26
27// #include <kernel/modulop.h>
28
29#include <kernel/febase.h>
30#include <kernel/kstd1.h>
31#include <kernel/kutil.h>
32#include <kernel/stairc.h>
33
34#include <kernel/syz.h>
35#include <kernel/ideals.h>
36
37
38
39
40
41//#define SHOW_PROT
42//#define SHOW_RED
43//#define SHOW_HILB
44//#define SHOW_RESULT
45//#define INVERT_PAIRS
46//#define SHOW_CRIT
47//#define SHOW_SPRFL
48#define USE_CHAINCRIT
49#define poly_write(p) wrp(p);PrintLn()
50/*--- some heuristics to optimize the pair sets---*/
51/*--- chose only one (!!!) at the same time ------*/
52//#define USE_HEURISTIC1
53#define USE_HEURISTIC2
54
55#ifdef SHOW_CRIT
56static int crit;
57static int crit1;
58static int spfl;
59static int cons_pairs;
60static int crit_fails;
61#endif
62typedef struct sopen_pairs open_pairs;
63typedef open_pairs* crit_pairs;
64struct sopen_pairs
65{
66  crit_pairs next;
67  int first_poly;
68  int second_poly;
69};
70/*3
71* computes pairs from the new elements (beginning with the element newEl)
72* in the module index
73*/
74static void syCreateNewPairs_Hilb(syStrategy syzstr, int index,
75            int actdeg)
76{
77  SSet temp;
78  SObject tso;
79  poly toHandle,tsyz=NULL,p,pp;
80  int r1,r2=0,rr,l=(*syzstr->Tl)[index];
81  int i,j,r=0,ti;
82  BOOLEAN toComp=FALSE;
83#ifdef USE_CHAINCRIT
84  crit_pairs cp=NULL,tcp;
85#endif
86  actdeg += index;
87  long * ShiftedComponents = syzstr->ShiftedComponents[index-1];
88  int* Components = syzstr->truecomponents[index-1];
89
90  while ((l>0) && ((syzstr->resPairs[index])[l-1].lcm==NULL)) l--;
91  rr = l-1;
92  while ((rr>=0) && (((syzstr->resPairs[index])[rr].p==NULL) || 
93        ((syzstr->resPairs[index])[rr].order>actdeg))) rr--;
94  r2 = rr+1; 
95  while ((rr>=0) && ((syzstr->resPairs[index])[rr].order==actdeg)
96         && ((syzstr->resPairs[index])[rr].syzind<0))
97  {
98    rr--;
99    r++;
100  }
101  if (r==0) return;
102  ideal nP=idInit(l,syzstr->res[index]->rank);
103#ifdef INVERT_PAIRS
104  r1 = rr+1;
105#else
106  r1 = r2-1;
107#endif
108/*---------- there are new pairs ------------------------------*/
109  loop
110  {
111/*--- chose first new elements --------------------------------*/
112    toComp = FALSE;
113    toHandle = (syzstr->resPairs[index])[r1].p;
114    if (toHandle!=NULL)
115    {
116      int tc=pGetComp(toHandle);
117      (syzstr->resPairs[index])[r1].syzind = 0;
118      for (i=0; i<r1;i++)
119      {
120        if (((syzstr->resPairs[index])[i].p!=NULL) &&
121            (pGetComp((syzstr->resPairs[index])[i].p)==tc))
122        {
123#ifdef USE_CHAINCRIT
124          tcp = cp;
125          if (tcp!=NULL)
126          {
127            while ((tcp!=NULL) && 
128              ((tcp->first_poly!=i)||(tcp->second_poly!=r1))) tcp = tcp->next;
129          }
130          if (tcp==NULL)
131          {
132#endif
133            p = pOne();
134            pLcm((syzstr->resPairs[index])[i].p,toHandle,p);
135            pSetm(p);
136            j = 0;
137            while (j<i)
138            { 
139              if (nP->m[j]!=NULL)
140              {
141                if (pLmDivisibleByNoComp(nP->m[j],p))
142                {
143                  pDelete(&p);
144                  /* p = NULL;*/
145                  break;
146                }
147                else if (pLmDivisibleByNoComp(p,nP->m[j]))
148                {
149                  pDelete(&(nP->m[j]));
150                  /* nP->m[j] = NULL;*/
151                }
152#ifdef USE_CHAINCRIT
153                else
154                {
155                  poly p1,p2;
156                  int ip= currRing->N;
157                  p1 = pDivide(p,(syzstr->resPairs[index])[i].p);
158                  p2 = pDivide(nP->m[j],(syzstr->resPairs[index])[j].p);
159                  while ((ip>0) && (pGetExp(p1,ip)*pGetExp(p2,ip)==0)) ip--;
160                  if (ip==0)
161                  {
162#ifdef SHOW_SPRFL
163Print("Hier: %d, %d\n",j,i);
164#endif
165                    if (i>rr)
166                    {
167                      tcp=(crit_pairs)omAlloc0(sizeof(sopen_pairs));
168                      tcp->next = cp;
169                      tcp->first_poly = j;
170                      tcp->second_poly = i;
171                      cp = tcp;
172                      tcp = NULL;
173                    }
174                    else
175                    {
176                      ti=0;
177                      while ((ti<l) && (((syzstr->resPairs[index])[ti].ind1!=j)||
178                             ((syzstr->resPairs[index])[ti].ind2!=i))) ti++;
179                      if (ti<l) 
180                      {
181#ifdef SHOW_SPRFL
182Print("gefunden in Mod %d: ",index); poly_write((syzstr->resPairs[index])[ti].lcm);
183#endif
184                        syDeletePair(&(syzstr->resPairs[index])[ti]);
185#ifdef SHOW_CRIT
186                        crit1++;
187#endif
188                        toComp = TRUE;
189                      }
190                    }
191                  }
192                  pLmDelete(&p1);
193                  pLmDelete(&p2);
194                }
195#endif
196              }
197              j++;
198            }
199            if (p!=NULL)
200            {
201              nP->m[i] = p;
202            }
203#ifdef USE_CHAINCRIT
204          }
205          else
206          {
207#ifdef SHOW_CRIT
208            crit1++;
209#endif
210          }
211#endif
212        }
213      }
214      if (toComp) syCompactify1(syzstr->resPairs[index],&l,r1);
215      for (i=0;i<r1;i++)
216      {
217        if (nP->m[i]!=NULL)
218        {
219          tso.lcm = p = nP->m[i];
220          nP->m[i] = NULL;
221          tso.order = pTotaldegree(p);
222          if ((syzstr->cw!=NULL) && (index>0) && (pGetComp(p)>0))
223          {
224            int ii=index-1,jj=pGetComp(p);
225            while (ii>0)
226            {
227              jj = pGetComp(syzstr->res[ii]->m[jj-1]);
228              ii--;
229            }
230            tso.order += (*syzstr->cw)[jj-1];
231          }
232          tso.p1 = (syzstr->resPairs[index])[i].p;
233          tso.p2 = toHandle;
234          tso.ind1 = i;
235          tso.ind2 = r1;
236          tso.syzind = -1;
237          tso.isNotMinimal = (poly)1;
238          tso.p = NULL;
239          tso.length = -1;
240          number coefgcd =
241            nGcd(pGetCoeff(tso.p1),pGetCoeff(tso.p2),currRing);
242          tso.syz = pCopy((syzstr->resPairs[index])[i].syz);
243          poly tt = pDivide(tso.lcm,tso.p1);
244          pSetCoeff(tt,nDiv(pGetCoeff(tso.p1),coefgcd));
245          tso.syz = pMult_mm(tso.syz,tt);
246          pLmDelete(&tt);
247          coefgcd = nNeg(coefgcd);
248          pp = pCopy((syzstr->resPairs[index])[r1].syz);
249          tt = pDivide(tso.lcm,tso.p2);
250          pSetCoeff(tt,nDiv(pGetCoeff(tso.p2),coefgcd));
251          pp = pMult_mm(pp,tt);
252          pLmDelete(&tt);
253          tso.syz = pAdd(pp,tso.syz);
254          nDelete(&coefgcd);
255          pSetComp(tso.lcm,pGetComp((syzstr->resPairs[index])[r1].syz));
256#ifdef SHOW_PROT
257Print("erzeuge Paar im Modul %d,%d mit: \n",index,tso.order);
258PrintS("poly1: ");poly_write(tso.p1);
259PrintS("poly2: ");poly_write(tso.p2);
260PrintS("syz: ");poly_write(tso.syz);
261PrintS("sPoly: ");poly_write(tso.p);
262PrintLn();
263#endif
264          syEnterPair(syzstr,&tso,&l,index);
265        }
266      }
267    }
268#ifdef INVERT_PAIRS
269    r1++;
270    if (r1>=r2) break;
271#else
272    r1--;
273    if (r1<=rr) break;
274#endif
275  }
276  idDelete(&nP);
277#ifdef USE_CHAINCRIT
278  while (cp!=NULL)
279  {
280    tcp = cp;
281    cp = cp->next;
282    omFreeSize((ADDRESS)tcp,sizeof(sopen_pairs));
283  }
284#endif
285}
286
287/*3
288* determines the place of a polynomial in the right ordered resolution
289* set the vectors of truecomponents
290*/
291static void syOrder_Hilb(poly p,syStrategy syzstr,int index)
292{
293  int i=IDELEMS(syzstr->orderedRes[index]);
294
295  while ((i>0) && (syzstr->orderedRes[index]->m[i-1]==NULL)) i--;
296  syzstr->orderedRes[index]->m[i] = p;
297}
298
299static void syHalfPair(poly syz, int newEl, syStrategy syzstr, int index)
300{
301  SObject tso;
302  memset(&tso,0,sizeof(tso));
303  int l=(*syzstr->Tl)[index];
304
305  while ((l>0) && ((syzstr->resPairs[index])[l-1].syz==NULL)) l--;
306  if ((syzstr->cw!=NULL) && (index>0) && (pGetComp(syz)>0))
307  {
308    int ii=index-1,jj=pGetComp(syz);
309    while (ii>0)
310    {
311      jj = pGetComp(syzstr->res[ii]->m[jj-1]);
312      ii--;
313    }
314    tso.order += (*syzstr->cw)[jj-1];
315  }
316  tso.p1 = NULL;
317  tso.p2 = NULL;
318  tso.ind1 = 0;
319  tso.ind2 = 0;
320  tso.syzind = -1;
321  tso.isNotMinimal = NULL;
322  tso.p = syz;
323  tso.order = pTotaldegree(syz);
324  tso.syz = pHead(syz);
325  pSetComp(tso.syz,newEl+1);
326  pSetm(tso.syz);
327  tso.lcm = pHead(tso.syz);
328  tso.length = pLength(syz);
329  syOrder_Hilb(syz,syzstr,index);
330#ifdef SHOW_PROT
331Print("erzeuge Halbpaar im Module %d,%d mit: \n",index,tso.order);
332PrintS("syz: ");poly_write(tso.syz);
333PrintS("sPoly: ");poly_write(tso.p);
334PrintLn();
335#endif
336  syEnterPair(syzstr,&tso,&l,index);
337}
338/*3
339* computes the order of pairs of same degree
340* must be monoton
341*/
342static intvec* syLinStrat2(SSet nextPairs, syStrategy syzstr,
343                          int howmuch, int index,intvec ** secondpairs)
344{
345  ideal o_r=syzstr->res[index+1];
346  int i=0,i1=0,i2=0,l,ll=IDELEMS(o_r);
347  intvec *result=new intvec(howmuch+1);
348  BOOLEAN isDivisible;
349  SObject tso;
350
351#ifndef USE_HEURISTIC2
352  while (i1<howmuch)
353  {
354    (*result)[i1] = i1+1;
355    i1++;
356  }
357  return result;
358#else
359  while ((ll>0) && (o_r->m[ll-1]==NULL)) ll--;
360  while (i<howmuch)
361  {
362    tso = nextPairs[i];
363    isDivisible = FALSE;
364    l = 0;
365    while ((l<ll) && (!isDivisible))
366    {
367      if (o_r->m[l]!=NULL)
368      {
369        isDivisible = isDivisible ||
370          pLmDivisibleBy(o_r->m[l],tso.lcm);
371      }
372      l++;
373    }
374    if (isDivisible)
375    {
376#ifdef SHOW_PROT
377Print("streiche Paar im Modul %d,%d mit: \n",index,nextPairs[i].order);
378PrintS("poly1: ");poly_write(nextPairs[i].p1);
379PrintS("poly2: ");poly_write(nextPairs[i].p2);
380PrintS("syz: ");poly_write(nextPairs[i].syz);
381PrintS("sPoly: ");poly_write(nextPairs[i].p);
382PrintLn();
383#endif
384      //syDeletePair(&nextPairs[i]);
385      if (*secondpairs==NULL) *secondpairs = new intvec(howmuch);
386      (**secondpairs)[i2] = i+1;
387      i2++;
388#ifdef SHOW_CRIT
389      crit++;
390#endif
391    }
392    else
393    {
394//      nextPairs[i].p = sySPoly(tso.p1, tso.p2,tso.lcm);
395      (*result)[i1] = i+1;
396      i1++;
397    }
398    i++;
399  }
400  return result;
401#endif
402}
403
404inline void sySPRedSyz(syStrategy syzstr,sSObject redWith,poly q=NULL)
405{
406  poly p=pDivide(q,redWith.p);
407  pSetCoeff(p,nDiv(pGetCoeff(q),pGetCoeff(redWith.p)));
408  int il=-1;
409  kBucket_Minus_m_Mult_p(syzstr->syz_bucket,p,redWith.syz,&il,NULL);
410  pLmDelete(&p);
411}
412
413static poly syRed_Hilb(poly toRed,syStrategy syzstr,int index)
414{
415  ideal redWith=syzstr->res[index];
416  if (redWith==NULL) return toRed;
417  int j=IDELEMS(redWith),i;
418  poly q,result=NULL,resultp;
419
420  while ((j>0) && (redWith->m[j-1]==NULL)) j--;
421  if ((toRed==NULL) || (j==0)) return toRed;
422  kBucketInit(syzstr->bucket,toRed,-1);
423  q = kBucketGetLm(syzstr->bucket);
424  loop
425  {
426    if (q==NULL)
427    {
428      break;
429    }
430    i = 0;
431    loop
432    {
433      if (pLmDivisibleBy(redWith->m[i],q))
434      {
435        number up = kBucketPolyRed(syzstr->bucket,redWith->m[i],
436                         pLength(redWith->m[i]), NULL);
437        nDelete(&up);
438        q = kBucketGetLm(syzstr->bucket);
439        if (toRed==NULL) break;
440        i = 0;
441      }
442      else
443      {
444        i++;
445      }
446      if ((i>=j) || (q==NULL)) break;
447    }
448    if (q!=NULL)
449    {
450      if (result==NULL)
451      {
452        resultp = result = kBucketExtractLm(syzstr->bucket);
453      }
454      else
455      {
456        pNext(resultp) = kBucketExtractLm(syzstr->bucket);
457        pIter(resultp);
458      }
459      q = kBucketGetLm(syzstr->bucket);
460    }
461  }
462  kBucketClear(syzstr->bucket,&q,&i);
463  if (q!=NULL) PrintS("Hier ist was schief gelaufen!\n");
464  return result;
465}
466
467#ifdef USE_HEURISTIC1
468intvec *ivStrip(intvec* arg)
469{
470  int l=arg->rows()*arg->cols(),i=0,ii=0;
471  intvec *tempV=new intvec(l);
472
473  while (i+ii<l)
474  {
475    if ((*arg)[i+ii]==0)
476    {
477      ii++;
478    }
479    else
480    {
481      (*tempV)[i] = (*arg)[i+ii];
482      i++;
483    }
484  }
485  if (i==0)
486  {
487    delete tempV;
488    return NULL;
489  }
490  intvec * result=new intvec(i+1);
491  for (ii=0;ii<i;ii++)
492   (*result)[ii] = (*tempV)[ii];
493  delete tempV;
494  return result;
495}
496#endif
497
498/*3
499* reduces all pairs of degree deg in the module index
500* put the reduced generators to the resolvente which contains
501* the truncated kStd
502*/
503static void syRedNextPairs_Hilb(SSet nextPairs, syStrategy syzstr,
504               int howmuch, int index,int actord,int* toSub,
505               int *maxindex,int *maxdeg)
506{
507  int i,j,k=IDELEMS(syzstr->res[index]);
508  int ks=IDELEMS(syzstr->res[index+1]),kk,l,ll;
509  int ks1=IDELEMS(syzstr->orderedRes[index+1]);
510  int kres=(*syzstr->Tl)[index];
511  int toGo=0;
512  int il;
513  number coefgcd,n;
514  SSet redset=syzstr->resPairs[index];
515  poly p=NULL,q,tp;
516  intvec *spl1;
517  SObject tso;
518  intvec *spl3=NULL;
519#ifdef USE_HEURISTIC1
520  intvec *spl2=new intvec(howmuch+1,howmuch+1,0);
521  int there_are_superfluous=0;
522  int step=1,jj,j1,j2;
523#endif
524  long * ShiftedComponents = syzstr->ShiftedComponents[index];
525  int* Components = syzstr->truecomponents[index];
526  assume(Components != NULL && ShiftedComponents != NULL);
527  BOOLEAN need_reset;
528
529  actord += index;
530  if ((nextPairs==NULL) || (howmuch==0)) return;
531  while ((k>0) && (syzstr->res[index]->m[k-1]==NULL)) k--;
532  while ((ks>0) && (syzstr->res[index+1]->m[ks-1]==NULL)) ks--;
533  while ((ks1>0) && (syzstr->orderedRes[index+1]->m[ks1-1]==NULL)) ks1--;
534  while ((kres>0) &&
535        ((redset[kres-1].p==NULL) || (redset[kres-1].order>actord))) kres--;
536  while ((kres<(*syzstr->Tl)[index]) &&
537         (redset[kres-1].order!=0) && (redset[kres-1].order<=actord)) kres++;
538  spl1 = syLinStrat2(nextPairs,syzstr,howmuch,index,&spl3);
539#ifdef SHOW_PROT
540PrintS("spl1 ist hier: ");spl1->show(0,0);
541#endif
542  i=0;
543  kk = (*spl1)[i]-1;
544  if (index==1)
545  {
546    intvec * temp1_hilb = hHstdSeries(syzstr->res[index],NULL,NULL,NULL);
547    if (actord<temp1_hilb->length())
548    {
549      toGo = (*temp1_hilb)[actord];
550#ifdef SHOW_HILB
551Print("\nStze toGo im Modul %d und Grad %d auf: %d\n",1,actord-1,toGo);
552#endif
553    }
554    delete temp1_hilb;
555  }
556  else
557  {
558    if (actord<=(syzstr->hilb_coeffs[index])->length())
559    {
560      toGo = (*syzstr->hilb_coeffs[index])[actord-1];
561#ifdef SHOW_HILB
562Print("\nStze toGo im Modul %d und Grad %d auf: %d\n",index,actord-1,toGo);
563#endif
564    }
565  }
566  if ((syzstr->hilb_coeffs[index+1]!=NULL) &&
567      (actord<=(syzstr->hilb_coeffs[index+1])->length()))
568  {
569    toGo += (*syzstr->hilb_coeffs[index+1])[actord-1];
570#ifdef SHOW_HILB
571Print("\nAddiere zu toGo aus Modul %d und Grad %d: %d\n",index+1,actord-1,(*syzstr->hilb_coeffs[index+1])[actord-1]);
572#endif
573  }
574#ifdef SHOW_HILB
575Print("<H%d>",toGo);
576#endif
577  while (kk>=0)
578  {
579    if (toGo==0) 
580    {
581      while (kk>=0)
582      {
583        pDelete(&nextPairs[kk].p);
584        pDelete(&nextPairs[kk].syz);
585        syDeletePair(&nextPairs[kk]);
586        nextPairs[kk].p = nextPairs[kk].syz = nextPairs[kk].lcm = NULL;
587        i++;
588        kk = (*spl1)[i]-1;
589#ifdef USE_HEURISTIC2
590        if (kk<0)
591        {
592          i = 0;
593          delete spl1;
594          spl1 = spl3;
595          spl3 = NULL;
596          if (spl1!=NULL)
597            kk = (*spl1)[i]-1;
598        }
599#endif
600      }
601      if (spl1!=NULL) delete spl1;
602      break;
603    }
604    tso = nextPairs[kk];
605    if ((tso.p1!=NULL) && (tso.p2!=NULL))
606    {
607#ifdef SHOW_CRIT
608      cons_pairs++;
609#endif
610      //tso.p = sySPoly(tso.p1, tso.p2,tso.lcm);
611      tso.p = ksOldCreateSpoly(tso.p2, tso.p1);
612#ifdef SHOW_PROT
613PrintS("reduziere Paar mit: \n");
614PrintS("poly1: ");poly_write(tso.p1);
615PrintS("poly2: ");poly_write(tso.p2);
616PrintS("syz: ");poly_write(tso.syz);
617PrintS("sPoly: ");poly_write(tso.p);
618#endif
619      if (tso.p != NULL)
620      {
621        kBucketInit(syzstr->bucket,tso.p,-1);
622        kBucketInit(syzstr->syz_bucket,tso.syz,-1);
623        q = kBucketGetLm(syzstr->bucket);
624        j = 0;
625        while (j<kres) 
626        {
627          if ((redset[j].p!=NULL) && (pLmDivisibleBy(redset[j].p,q)) 
628              && ((redset[j].ind1!=tso.ind1) || (redset[j].ind2!=tso.ind2)))
629          {
630#ifdef SHOW_RED
631kBucketClear(syzstr->bucket,&tso.p,&tso.length);
632kBucketClear(syzstr->syz_bucket,&tso.syz,&il);
633PrintS("reduziere: ");poly_write(tso.p);
634PrintS("syz: ");poly_write(tso.syz);
635PrintS("mit: ");poly_write(redset[j].p);
636PrintS("syz: ");poly_write(redset[j].syz);
637kBucketInit(syzstr->bucket,tso.p,tso.length);
638kBucketInit(syzstr->syz_bucket,tso.syz,il);
639#endif
640            sySPRedSyz(syzstr,redset[j],q);
641            number up = kBucketPolyRed(syzstr->bucket,redset[j].p,
642                         redset[j].length, NULL);
643            nDelete(&up);
644            q = kBucketGetLm(syzstr->bucket);
645#ifdef SHOW_RED
646kBucketClear(syzstr->bucket,&tso.p,&tso.length);
647kBucketClear(syzstr->syz_bucket,&tso.syz,&il);
648PrintS("zu: ");poly_write(tso.p);
649PrintS("syz: ");poly_write(tso.syz);
650kBucketInit(syzstr->bucket,tso.p,tso.length);
651kBucketInit(syzstr->syz_bucket,tso.syz,il);
652PrintLn();
653#endif
654            if (q==NULL) break;
655            j = 0;
656          }
657          else
658          {
659            j++;
660          }
661        }
662        kBucketClear(syzstr->bucket,&tso.p,&tso.length);
663        kBucketClear(syzstr->syz_bucket,&tso.syz,&il);
664      }
665#ifdef SHOW_PROT
666PrintS("erhalte Paar mit: \n");
667PrintS("syz: ");poly_write(tso.syz);
668PrintS("sPoly: ");poly_write(tso.p);
669PrintLn();
670#endif
671#ifdef SHOW_SPRFL
672//PrintLn();
673wrp(tso.lcm);
674Print(" mit index %d, %d ",tso.ind1,tso.ind2);
675#endif
676      if (tso.p != NULL)
677      {
678        if (TEST_OPT_PROT) PrintS("g");
679        (*toSub)++;
680        toGo--;
681        if (!nIsOne(pGetCoeff(tso.p)))
682        {
683          number n=nInvers(pGetCoeff(tso.p));
684          pNorm(tso.p);
685          pMult_nn(tso.syz,n);
686          nDelete(&n);
687        }
688        if (k==IDELEMS((syzstr->res)[index]))
689          syEnlargeFields(syzstr,index);
690        syzstr->res[index]->m[k] = tso.p;
691        k++;
692      }
693      else
694      {
695        if (ks==IDELEMS(syzstr->res[index+1]))
696          syEnlargeFields(syzstr,index+1);
697        syzstr->res[index+1]->m[ks] = syRed_Hilb(tso.syz,syzstr,index+1);
698        if (syzstr->res[index+1]->m[ks]!=NULL)
699        {
700          if (TEST_OPT_PROT) PrintS("s");
701          toGo--;
702          pNorm(syzstr->res[index+1]->m[ks]);
703          syHalfPair(syzstr->res[index+1]->m[ks],ks1,syzstr,index+1);
704          ks++;
705          ks1++;
706          if (index+1>*maxindex) *maxindex = index+1;
707          if (actord-index>*maxdeg) *maxdeg = actord-index;
708        }
709        else
710        {
711          if (TEST_OPT_PROT) PrintS("-");
712#ifdef SHOW_CRIT
713          spfl++;
714#endif
715#ifdef USE_HEURISTIC1
716          if (there_are_superfluous>=0)
717          {
718            j = i+1;
719            jj = (*spl1)[j]-1;
720            j1 = 1;
721            while (jj>=0)
722            {
723              if (tso.ind2==nextPairs[jj].ind2)
724              {
725                IMATELEM(*spl2,j1,step) = jj+1;
726                j1++;
727                for (j2=j;j2<spl1->length()-1;j2++)
728                {
729                  (*spl1)[j2] = (*spl1)[j2+1];
730                }
731              }
732              else
733              {
734                j++;
735              }
736              jj = (*spl1)[j]-1;
737            }
738            step++;
739            if (there_are_superfluous==0) there_are_superfluous = 1;
740          }
741#endif
742#ifdef SHOW_SPRFL
743Print("ist ueberfluessig in Mod %d",index);
744//Print("\n ueberfluessig in Mod %d:",index);
745//wrp(tso.lcm);
746//PrintLn();
747#endif
748        }
749        tso.syz = NULL;
750        syDeletePair(&tso);
751        tso.p = tso.syz = tso.lcm = NULL;
752      }
753      nextPairs[kk] = tso;
754    }
755#ifdef SHOW_SPRFL
756PrintLn();
757#endif
758    i++;
759#ifdef SHOW_PROT
760PrintS("spl1 ist hier: ");spl1->show(0,0);
761Print("naechstes i ist: %d",i);
762#endif
763    kk = (*spl1)[i]-1;
764#ifdef USE_HEURISTIC1
765    if ((kk<0) && (there_are_superfluous>0))
766    {
767      i = 0;
768      delete spl1;
769      spl1 = ivStrip(spl2); 
770      delete spl2;
771      if (spl1!=NULL)
772      {
773        there_are_superfluous = -1;
774        kk = (*spl1)[i]-1;
775      }
776    } 
777#endif
778#ifdef USE_HEURISTIC2
779    if ((kk<0) && (toGo>0))
780    {
781#ifdef SHOW_CRIT
782      crit_fails++;
783#endif
784      i = 0;
785      delete spl1;
786      spl1 = spl3;
787      spl3 = NULL;
788      if (spl1!=NULL)
789        kk = (*spl1)[i]-1;
790    }
791#endif
792  }
793  delete spl1;
794  if (spl3!=NULL) delete spl3;
795}
796
797void sySetNewHilb(syStrategy syzstr, int toSub,int index,int actord)
798{
799  int i;
800  actord += index;
801  intvec * temp_hilb = hHstdSeries(syzstr->res[index+1],NULL,NULL,NULL);
802  intvec * cont_hilb = hHstdSeries(syzstr->res[index],NULL,NULL,NULL);
803  if ((index+1<syzstr->length) && (syzstr->hilb_coeffs[index+1]==NULL))
804  {
805    syzstr->hilb_coeffs[index+1] = new intvec(16*((actord/16)+1));
806  }
807  else if (actord>=syzstr->hilb_coeffs[index+1]->length())
808  {
809    intvec * ttt=new intvec(16*((actord/16)+1));
810    for (i=syzstr->hilb_coeffs[index+1]->length()-1;i>=0;i--)
811    {
812      (*ttt)[i] = (*(syzstr->hilb_coeffs[index+1]))[i];
813    }
814    delete syzstr->hilb_coeffs[index+1];
815    syzstr->hilb_coeffs[index+1] = ttt;
816  }
817  if (actord+1<temp_hilb->length())
818  {
819#ifdef SHOW_HILB
820Print("\nSetze fuer Modul %d im Grad %d die Wert: \n",index+1,actord);
821(temp_hilb)->show(0,0);
822#endif
823    int k=si_min(temp_hilb->length()-1,(syzstr->hilb_coeffs[index+1])->length());
824    for (int j=k;j>actord;j--)
825      (*(syzstr->hilb_coeffs[index+1]))[j-1] = (*temp_hilb)[j];
826  }
827  else
828  {
829    (*(syzstr->hilb_coeffs[index+1]))[actord] = 0;
830  }
831  delete temp_hilb;
832  if ((index>1) && (actord<=syzstr->hilb_coeffs[index]->length()))
833  {
834#ifdef SHOW_HILB
835Print("\nSubtrahiere im Modul %d im Grad %d den Wert: %d\n",index,actord-1,toSub);
836#endif
837    (*syzstr->hilb_coeffs[index])[actord-1]-=toSub;
838  }
839  if (syzstr->hilb_coeffs[index]!=NULL)
840  {
841    if (cont_hilb->length()>syzstr->hilb_coeffs[index]->length())
842      syzstr->hilb_coeffs[index]->resize(cont_hilb->length());
843    for (int j=cont_hilb->length()-1;j>actord;j--)
844      (*(syzstr->hilb_coeffs[index]))[j-1] = (*cont_hilb)[j];
845  }
846  delete cont_hilb;
847#ifdef SHOW_HILB
848Print("<h,%d>",(*(syzstr->hilb_coeffs[index+1]))[actord]);
849#endif
850}
851
852/*3
853* reduces the generators of the module index in degree deg
854* (which are actual syzygies of the module index-1)
855* wrt. the ideal generated by elements of lower degrees
856*/
857static void syRedGenerOfCurrDeg_Hilb(syStrategy syzstr, int deg,int *maxindex,int *maxdeg)
858{
859  ideal res=syzstr->res[1];
860  int i=0,j,k=IDELEMS(res),k1=IDELEMS(syzstr->orderedRes[1]);
861  SSet sPairs1=syzstr->resPairs[1];
862  SSet sPairs=syzstr->resPairs[0];
863
864  while ((k>0) && (res->m[k-1]==NULL)) k--;
865  while ((k1>0) && (syzstr->orderedRes[1]->m[k1-1]==NULL)) k1--;
866  while ((i<(*syzstr->Tl)[0]) && (((sPairs)[i].syz==NULL) ||
867          ((sPairs)[i].order<deg)))
868    i++;
869  if ((i>=(*syzstr->Tl)[0]) || ((sPairs)[i].order>deg)) return;
870  while ((i<(*syzstr->Tl)[0]) && (((sPairs)[i].syz==NULL) ||
871         ((sPairs)[i].order==deg)))
872  {
873    if ((sPairs)[i].syz!=NULL)
874    {
875#ifdef SHOW_PROT
876PrintS("reduziere Erzeuger: \n");
877PrintS("syz: ");poly_write((sPairs)[i].syz);
878#endif
879      (sPairs)[i].syz = syRed_Hilb((sPairs)[i].syz,syzstr,1);
880#ifdef SHOW_PROT
881PrintS("erhalte Erzeuger: \n");
882PrintS("syz: ");poly_write((sPairs)[i].syz);
883PrintLn();
884#endif
885      if ((sPairs)[i].syz != NULL)
886      {
887        if (k==IDELEMS(res))
888        {
889          syEnlargeFields(syzstr,1);
890          res=syzstr->res[1];
891        }
892        if (TEST_OPT_DEBUG)
893        {
894          if ((sPairs)[i].isNotMinimal==NULL)
895          {
896            PrintS("\nminimal generator: ");
897            pWrite((syzstr->resPairs[0])[i].syz);
898            PrintS("comes from: ");pWrite((syzstr->resPairs[0])[i].p1);
899            PrintS("and: ");pWrite((syzstr->resPairs[0])[i].p2);
900          }
901        }
902        res->m[k] = (sPairs)[i].syz;
903        pNorm(res->m[k]);
904        syHalfPair(res->m[k],k1,syzstr,1);
905        k1++;
906        k++;
907        if (1>*maxindex) *maxindex = 1;
908        if (deg-1>*maxdeg) *maxdeg = deg-1;
909      }
910    }
911    i++;
912  }
913}
914
915/*3
916* reorders the result (stored in orderedRes) according
917* to the seqence given by res
918*/
919static void syReOrdResult_Hilb(syStrategy syzstr,int maxindex,int maxdeg)
920{
921  ideal reor,toreor;
922  int i,j,k,l,m,togo;
923  syzstr->betti = new intvec(maxdeg,maxindex+1,0);
924  (*syzstr->betti)[0] = 1;
925  for (i=1;i<=syzstr->length;i++)
926  {
927    if (!idIs0(syzstr->orderedRes[i])) 
928    {
929      toreor = syzstr->orderedRes[i];
930      k = IDELEMS(toreor);
931      while ((k>0) && (toreor->m[k-1]==NULL)) k--;
932      reor = idInit(k,toreor->rank);
933      togo = IDELEMS(syzstr->res[i]);
934      for (j=0;j<k;j++)
935      {
936        if (toreor->m[j]!=NULL) (IMATELEM(*syzstr->betti,p_FDeg(toreor->m[j],currRing)-i+1,i+1))++;
937        reor->m[j] = toreor->m[j];
938        toreor->m[j] = NULL;
939      }
940      m = 0; 
941      for (j=0;j<togo;j++)
942      {
943        if (syzstr->res[i]->m[j]!=NULL)
944        {
945          l = 0;
946          while ((l<k) && (syzstr->res[i]->m[j]!=reor->m[l])) l++;
947          if (l<k)
948          {
949            toreor->m[m] = reor->m[l];
950            reor->m[l] = NULL;
951            m++;
952          }
953        }
954      }
955      idDelete(&reor);
956    }
957  }
958}
959
960/*2
961* the CoCoA-algorithm for free resolutions, using a formula
962* for remaining pairs based on Hilbert-functions
963*/
964syStrategy syHilb(ideal arg,int * length)
965{
966  int i,j,actdeg=32000,index=0,reg=-1;
967  int startdeg,howmuch,toSub=0;
968  int maxindex=0,maxdeg=0;
969  ideal temp=NULL;
970  SSet nextPairs;
971  ring origR = currRing;
972  syStrategy syzstr=(syStrategy)omAlloc0(sizeof(ssyStrategy));
973
974  if ((idIs0(arg)) || (id_RankFreeModule(arg,currRing)>0))
975  {
976    syzstr->minres = (resolvente)omAllocBin(sip_sideal_bin);
977    syzstr->length = 1;
978    syzstr->minres[0] = idInit(1,arg->rank);
979    return syzstr;
980  }
981 
982  // Creare dp,S ring and change to it
983  syzstr->syRing = rAssure_dp_C(origR);
984  rChangeCurrRing(syzstr->syRing);
985
986  // set initial ShiftedComps
987  currcomponents = (int*)omAlloc0((arg->rank+1)*sizeof(int));
988  currShiftedComponents = (long*)omAlloc0((arg->rank+1)*sizeof(long));
989
990/*--- initializes the data structures---------------*/
991#ifdef SHOW_CRIT
992  crit = 0;
993  crit1 = 0;
994  spfl = 0;
995  cons_pairs = 0;
996  crit_fails = 0;
997#endif
998  syzstr->length = *length = currRing->N+2;
999  syzstr->Tl = new intvec(*length+1);
1000  temp = idInit(IDELEMS(arg),arg->rank);
1001  for (i=0;i<IDELEMS(arg);i++)
1002  {
1003    if (origR != syzstr->syRing)
1004      temp->m[i] = prCopyR( arg->m[i], origR, syzstr->syRing);
1005    else
1006      temp->m[i] = pCopy( arg->m[i]);
1007    if (temp->m[i]!=NULL)
1008    {
1009      j = pTotaldegree(temp->m[i]);
1010      if (j<actdeg) actdeg = j;
1011    }
1012  }
1013  idTest(temp);
1014  idSkipZeroes(temp);
1015  syzstr->resPairs = syInitRes(temp,length,syzstr->Tl,syzstr->cw);
1016  omFreeSize((ADDRESS)currcomponents,(arg->rank+1)*sizeof(int));
1017  omFreeSize((ADDRESS)currShiftedComponents,(arg->rank+1)*sizeof(int));
1018  syzstr->res = (resolvente)omAlloc0((*length+1)*sizeof(ideal));
1019  syzstr->orderedRes = (resolvente)omAlloc0((*length+1)*sizeof(ideal));
1020  syzstr->elemLength = (int**)omAlloc0((*length+1)*sizeof(int*));
1021  syzstr->truecomponents = (int**)omAlloc0((*length+1)*sizeof(int*));
1022  syzstr->backcomponents = (int**)omAlloc0((*length+1)*sizeof(int*));
1023  syzstr->ShiftedComponents = (long**)omAlloc0((*length+1)*sizeof(long*));
1024  syzstr->Howmuch = (int**)omAlloc0((*length+1)*sizeof(int*));
1025  syzstr->Firstelem = (int**)omAlloc0((*length+1)*sizeof(int*));
1026  syzstr->hilb_coeffs = (intvec**)omAlloc0((*length+1)*sizeof(intvec*));
1027  syzstr->sev = (unsigned long **)omAlloc0((*length+1)*sizeof(unsigned long*));
1028  syzstr->bucket = kBucketCreate(currRing);
1029  syzstr->syz_bucket = kBucketCreate(currRing);
1030  startdeg = actdeg;
1031  nextPairs = syChosePairs(syzstr,&index,&howmuch,&actdeg);
1032/*--- computes the resolution ----------------------*/
1033  while (nextPairs!=NULL)
1034  {
1035#ifdef SHOW_PROT
1036Print("compute %d Paare im Module %d im Grad %d \n",howmuch,index,actdeg+index);
1037#endif
1038    if (TEST_OPT_PROT) Print("%d",actdeg);
1039    if (TEST_OPT_PROT) Print("(m%d)",index);
1040    if (index==0)
1041      i = syInitSyzMod(syzstr,index,id_RankFreeModule(arg, origR)+1);
1042    else
1043      i = syInitSyzMod(syzstr,index);
1044    j = syInitSyzMod(syzstr,index+1);
1045    if (index>0)
1046    {
1047      syRedNextPairs_Hilb(nextPairs,syzstr,howmuch,index,actdeg,&toSub,&maxindex,&maxdeg);
1048      sySetNewHilb(syzstr,toSub,index,actdeg);
1049      toSub = 0;
1050      syCompactifyPairSet(syzstr->resPairs[index],(*syzstr->Tl)[index],0);
1051    }
1052    else
1053      syRedGenerOfCurrDeg_Hilb(syzstr,actdeg,&maxindex,&maxdeg);
1054/*--- creates new pairs -----------------------------*/
1055#ifdef SHOW_PROT
1056Print("Bilde neue Paare in Modul %d!\n",index);
1057#endif
1058    syCreateNewPairs_Hilb(syzstr,index,actdeg);
1059    if (index<(*length)-1)
1060    {
1061#ifdef SHOW_PROT
1062Print("Bilde neue Paare in Modul %d!\n",index+1);
1063#endif
1064      syCreateNewPairs_Hilb(syzstr,index+1,actdeg-1);
1065    }
1066    index++;
1067    nextPairs = syChosePairs(syzstr,&index,&howmuch,&actdeg);
1068  }
1069  syReOrdResult_Hilb(syzstr,maxindex,maxdeg);
1070#ifdef SHOW_RESULT
1071PrintS("minimal resolution:\n");
1072for (int ti=1;ti<=*length;ti++)
1073{
1074  if (!idIs0(syzstr->orderedRes[ti])) idPrint(syzstr->orderedRes[ti]);
1075}
1076PrintS("full resolution:\n");
1077for (int ti=1;ti<=*length;ti++)
1078{
1079  if (!idIs0(syzstr->res[ti])) idPrint(syzstr->res[ti]);
1080}
1081#endif
1082#ifdef SHOW_CRIT
1083Print("Criterion %d times applied\n",crit);
1084Print("Criterion1 %d times applied\n",crit1);
1085Print("%d superfluous pairs\n",spfl);
1086Print("%d pairs considered\n",cons_pairs);
1087Print("Criterion fails %d times\n",crit_fails);
1088crit = 0;
1089crit1 = 0;
1090spfl = 0;
1091cons_pairs = 0;
1092crit_fails = 0;
1093#endif
1094  if (temp!=NULL) idDelete(&temp);
1095  kBucketDestroy(&(syzstr->bucket));
1096  kBucketDestroy(&(syzstr->syz_bucket));
1097  if (origR != syzstr->syRing) 
1098    rChangeCurrRing(origR);
1099  else
1100    currRing =  origR;
1101  if (TEST_OPT_PROT) PrintLn();
1102  return syzstr;
1103}
Note: See TracBrowser for help on using the repository browser.