source: git/kernel/GBEngine/syz2.cc @ 9e8bfa

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