source: git/kernel/syz2.cc @ fbc7cb

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