source: git/kernel/syz2.cc @ 16f511

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