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

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