source: git/kernel/GBEngine/syz2.cc @ c796c8

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