source: git/kernel/syz2.cc @ a5b80a

spielwiese
Last change on this file since a5b80a was 76cfef, checked in by Oleksandr Motsak <motsak@…>, 13 years ago
FIX: fixed #includes in the kernel/ sources ADD: dummy headers in some strange cases
  • Property mode set to 100644
File size: 28.9 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/* $Id$ */
5/*
6* ABSTRACT: resolutions
7*/
8
9#include <kernel/mod2.h>
10#include <omalloc/mylimits.h>
11#include <misc/options.h>
12#include <omalloc/omalloc.h>
13#include <kernel/syz.h>
14#include <polys/polys.h>
15#include <kernel/febase.h>
16#include <kernel/kstd1.h>
17#include <kernel/kutil.h>
18#include <kernel/stairc.h>
19//#include "cntrlc.h"
20#include <misc/intvec.h>
21#include <coeffs/numbers.h>
22#include <kernel/modulop.h>
23#include <kernel/ideals.h>
24#include <misc/intvec.h>
25#include <polys/monomials/ring.h>
26#include <polys/kbuckets.h>
27#include <polys/prCopy.h>
28
29//#define SHOW_PROT
30//#define SHOW_RED
31//#define SHOW_HILB
32//#define SHOW_RESULT
33//#define INVERT_PAIRS
34//#define SHOW_CRIT
35//#define SHOW_SPRFL
36#define USE_CHAINCRIT
37#define poly_write(p) wrp(p);PrintLn()
38/*--- some heuristics to optimize the pair sets---*/
39/*--- chose only one (!!!) at the same time ------*/
40//#define USE_HEURISTIC1
41#define USE_HEURISTIC2
42
43#ifdef SHOW_CRIT
44static int crit;
45static int crit1;
46static int spfl;
47static int cons_pairs;
48static int crit_fails;
49#endif
50typedef struct sopen_pairs open_pairs;
51typedef open_pairs* crit_pairs;
52struct sopen_pairs
53{
54  crit_pairs next;
55  int first_poly;
56  int second_poly;
57};
58/*3
59* computes pairs from the new elements (beginning with the element newEl)
60* in the module index
61*/
62static void syCreateNewPairs_Hilb(syStrategy syzstr, int index,
63            int actdeg)
64{
65  SSet temp;
66  SObject tso;
67  poly toHandle,tsyz=NULL,p,pp;
68  int r1,r2=0,rr,l=(*syzstr->Tl)[index];
69  int i,j,r=0,ti;
70  BOOLEAN toComp=FALSE;
71#ifdef USE_CHAINCRIT
72  crit_pairs cp=NULL,tcp;
73#endif
74  actdeg += index;
75  long * ShiftedComponents = syzstr->ShiftedComponents[index-1];
76  int* Components = syzstr->truecomponents[index-1];
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=pVariables;
145                  p1 = pDivide(p,(syzstr->resPairs[index])[i].p);
146                  p2 = pDivide(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                  pLmDelete(&p1);
181                  pLmDelete(&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            nGcd(pGetCoeff(tso.p1),pGetCoeff(tso.p2),currRing);
230          tso.syz = pCopy((syzstr->resPairs[index])[i].syz);
231          poly tt = pDivide(tso.lcm,tso.p1);
232          pSetCoeff(tt,nDiv(pGetCoeff(tso.p1),coefgcd));
233          tso.syz = pMult_mm(tso.syz,tt);
234          pLmDelete(&tt);
235          coefgcd = nNeg(coefgcd);
236          pp = pCopy((syzstr->resPairs[index])[r1].syz);
237          tt = pDivide(tso.lcm,tso.p2);
238          pSetCoeff(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=pDivide(q,redWith.p);
395  pSetCoeff(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,l,ll;
497  int ks1=IDELEMS(syzstr->orderedRes[index+1]);
498  int kres=(*syzstr->Tl)[index];
499  int toGo=0;
500  int il;
501  number coefgcd,n;
502  SSet redset=syzstr->resPairs[index];
503  poly p=NULL,q,tp;
504  intvec *spl1;
505  SObject tso;
506  intvec *spl3=NULL;
507#ifdef USE_HEURISTIC1
508  intvec *spl2=new intvec(howmuch+1,howmuch+1,0);
509  int there_are_superfluous=0;
510  int step=1,jj,j1,j2;
511#endif
512  long * ShiftedComponents = syzstr->ShiftedComponents[index];
513  int* Components = syzstr->truecomponents[index];
514  assume(Components != NULL && ShiftedComponents != NULL);
515  BOOLEAN need_reset;
516
517  actord += index;
518  if ((nextPairs==NULL) || (howmuch==0)) return;
519  while ((k>0) && (syzstr->res[index]->m[k-1]==NULL)) k--;
520  while ((ks>0) && (syzstr->res[index+1]->m[ks-1]==NULL)) ks--;
521  while ((ks1>0) && (syzstr->orderedRes[index+1]->m[ks1-1]==NULL)) ks1--;
522  while ((kres>0) &&
523        ((redset[kres-1].p==NULL) || (redset[kres-1].order>actord))) kres--;
524  while ((kres<(*syzstr->Tl)[index]) &&
525         (redset[kres-1].order!=0) && (redset[kres-1].order<=actord)) kres++;
526  spl1 = syLinStrat2(nextPairs,syzstr,howmuch,index,&spl3);
527#ifdef SHOW_PROT
528PrintS("spl1 ist hier: ");spl1->show(0,0);
529#endif
530  i=0;
531  kk = (*spl1)[i]-1;
532  if (index==1)
533  {
534    intvec * temp1_hilb = hHstdSeries(syzstr->res[index],NULL,NULL,NULL);
535    if (actord<temp1_hilb->length())
536    {
537      toGo = (*temp1_hilb)[actord];
538#ifdef SHOW_HILB
539Print("\nStze toGo im Modul %d und Grad %d auf: %d\n",1,actord-1,toGo);
540#endif
541    }
542    delete temp1_hilb;
543  }
544  else
545  {
546    if (actord<=(syzstr->hilb_coeffs[index])->length())
547    {
548      toGo = (*syzstr->hilb_coeffs[index])[actord-1];
549#ifdef SHOW_HILB
550Print("\nStze toGo im Modul %d und Grad %d auf: %d\n",index,actord-1,toGo);
551#endif
552    }
553  }
554  if ((syzstr->hilb_coeffs[index+1]!=NULL) &&
555      (actord<=(syzstr->hilb_coeffs[index+1])->length()))
556  {
557    toGo += (*syzstr->hilb_coeffs[index+1])[actord-1];
558#ifdef SHOW_HILB
559Print("\nAddiere zu toGo aus Modul %d und Grad %d: %d\n",index+1,actord-1,(*syzstr->hilb_coeffs[index+1])[actord-1]);
560#endif
561  }
562#ifdef SHOW_HILB
563Print("<H%d>",toGo);
564#endif
565  while (kk>=0)
566  {
567    if (toGo==0) 
568    {
569      while (kk>=0)
570      {
571        pDelete(&nextPairs[kk].p);
572        pDelete(&nextPairs[kk].syz);
573        syDeletePair(&nextPairs[kk]);
574        nextPairs[kk].p = nextPairs[kk].syz = nextPairs[kk].lcm = NULL;
575        i++;
576        kk = (*spl1)[i]-1;
577#ifdef USE_HEURISTIC2
578        if (kk<0)
579        {
580          i = 0;
581          delete spl1;
582          spl1 = spl3;
583          spl3 = NULL;
584          if (spl1!=NULL)
585            kk = (*spl1)[i]-1;
586        }
587#endif
588      }
589      if (spl1!=NULL) delete spl1;
590      break;
591    }
592    tso = nextPairs[kk];
593    if ((tso.p1!=NULL) && (tso.p2!=NULL))
594    {
595#ifdef SHOW_CRIT
596      cons_pairs++;
597#endif
598      //tso.p = sySPoly(tso.p1, tso.p2,tso.lcm);
599      tso.p = ksOldCreateSpoly(tso.p2, tso.p1);
600#ifdef SHOW_PROT
601PrintS("reduziere Paar mit: \n");
602PrintS("poly1: ");poly_write(tso.p1);
603PrintS("poly2: ");poly_write(tso.p2);
604PrintS("syz: ");poly_write(tso.syz);
605PrintS("sPoly: ");poly_write(tso.p);
606#endif
607      if (tso.p != NULL)
608      {
609        kBucketInit(syzstr->bucket,tso.p,-1);
610        kBucketInit(syzstr->syz_bucket,tso.syz,-1);
611        q = kBucketGetLm(syzstr->bucket);
612        j = 0;
613        while (j<kres) 
614        {
615          if ((redset[j].p!=NULL) && (pLmDivisibleBy(redset[j].p,q)) 
616              && ((redset[j].ind1!=tso.ind1) || (redset[j].ind2!=tso.ind2)))
617          {
618#ifdef SHOW_RED
619kBucketClear(syzstr->bucket,&tso.p,&tso.length);
620kBucketClear(syzstr->syz_bucket,&tso.syz,&il);
621PrintS("reduziere: ");poly_write(tso.p);
622PrintS("syz: ");poly_write(tso.syz);
623PrintS("mit: ");poly_write(redset[j].p);
624PrintS("syz: ");poly_write(redset[j].syz);
625kBucketInit(syzstr->bucket,tso.p,tso.length);
626kBucketInit(syzstr->syz_bucket,tso.syz,il);
627#endif
628            sySPRedSyz(syzstr,redset[j],q);
629            number up = kBucketPolyRed(syzstr->bucket,redset[j].p,
630                         redset[j].length, NULL);
631            nDelete(&up);
632            q = kBucketGetLm(syzstr->bucket);
633#ifdef SHOW_RED
634kBucketClear(syzstr->bucket,&tso.p,&tso.length);
635kBucketClear(syzstr->syz_bucket,&tso.syz,&il);
636PrintS("zu: ");poly_write(tso.p);
637PrintS("syz: ");poly_write(tso.syz);
638kBucketInit(syzstr->bucket,tso.p,tso.length);
639kBucketInit(syzstr->syz_bucket,tso.syz,il);
640PrintLn();
641#endif
642            if (q==NULL) break;
643            j = 0;
644          }
645          else
646          {
647            j++;
648          }
649        }
650        kBucketClear(syzstr->bucket,&tso.p,&tso.length);
651        kBucketClear(syzstr->syz_bucket,&tso.syz,&il);
652      }
653#ifdef SHOW_PROT
654PrintS("erhalte Paar mit: \n");
655PrintS("syz: ");poly_write(tso.syz);
656PrintS("sPoly: ");poly_write(tso.p);
657PrintLn();
658#endif
659#ifdef SHOW_SPRFL
660//PrintLn();
661wrp(tso.lcm);
662Print(" mit index %d, %d ",tso.ind1,tso.ind2);
663#endif
664      if (tso.p != NULL)
665      {
666        if (TEST_OPT_PROT) PrintS("g");
667        (*toSub)++;
668        toGo--;
669        if (!nIsOne(pGetCoeff(tso.p)))
670        {
671          number n=nInvers(pGetCoeff(tso.p));
672          pNorm(tso.p);
673          pMult_nn(tso.syz,n);
674          nDelete(&n);
675        }
676        if (k==IDELEMS((syzstr->res)[index]))
677          syEnlargeFields(syzstr,index);
678        syzstr->res[index]->m[k] = tso.p;
679        k++;
680      }
681      else
682      {
683        if (ks==IDELEMS(syzstr->res[index+1]))
684          syEnlargeFields(syzstr,index+1);
685        syzstr->res[index+1]->m[ks] = syRed_Hilb(tso.syz,syzstr,index+1);
686        if (syzstr->res[index+1]->m[ks]!=NULL)
687        {
688          if (TEST_OPT_PROT) PrintS("s");
689          toGo--;
690          pNorm(syzstr->res[index+1]->m[ks]);
691          syHalfPair(syzstr->res[index+1]->m[ks],ks1,syzstr,index+1);
692          ks++;
693          ks1++;
694          if (index+1>*maxindex) *maxindex = index+1;
695          if (actord-index>*maxdeg) *maxdeg = actord-index;
696        }
697        else
698        {
699          if (TEST_OPT_PROT) PrintS("-");
700#ifdef SHOW_CRIT
701          spfl++;
702#endif
703#ifdef USE_HEURISTIC1
704          if (there_are_superfluous>=0)
705          {
706            j = i+1;
707            jj = (*spl1)[j]-1;
708            j1 = 1;
709            while (jj>=0)
710            {
711              if (tso.ind2==nextPairs[jj].ind2)
712              {
713                IMATELEM(*spl2,j1,step) = jj+1;
714                j1++;
715                for (j2=j;j2<spl1->length()-1;j2++)
716                {
717                  (*spl1)[j2] = (*spl1)[j2+1];
718                }
719              }
720              else
721              {
722                j++;
723              }
724              jj = (*spl1)[j]-1;
725            }
726            step++;
727            if (there_are_superfluous==0) there_are_superfluous = 1;
728          }
729#endif
730#ifdef SHOW_SPRFL
731Print("ist ueberfluessig in Mod %d",index);
732//Print("\n ueberfluessig in Mod %d:",index);
733//wrp(tso.lcm);
734//PrintLn();
735#endif
736        }
737        tso.syz = NULL;
738        syDeletePair(&tso);
739        tso.p = tso.syz = tso.lcm = NULL;
740      }
741      nextPairs[kk] = tso;
742    }
743#ifdef SHOW_SPRFL
744PrintLn();
745#endif
746    i++;
747#ifdef SHOW_PROT
748PrintS("spl1 ist hier: ");spl1->show(0,0);
749Print("naechstes i ist: %d",i);
750#endif
751    kk = (*spl1)[i]-1;
752#ifdef USE_HEURISTIC1
753    if ((kk<0) && (there_are_superfluous>0))
754    {
755      i = 0;
756      delete spl1;
757      spl1 = ivStrip(spl2); 
758      delete spl2;
759      if (spl1!=NULL)
760      {
761        there_are_superfluous = -1;
762        kk = (*spl1)[i]-1;
763      }
764    } 
765#endif
766#ifdef USE_HEURISTIC2
767    if ((kk<0) && (toGo>0))
768    {
769#ifdef SHOW_CRIT
770      crit_fails++;
771#endif
772      i = 0;
773      delete spl1;
774      spl1 = spl3;
775      spl3 = NULL;
776      if (spl1!=NULL)
777        kk = (*spl1)[i]-1;
778    }
779#endif
780  }
781  delete spl1;
782  if (spl3!=NULL) delete spl3;
783}
784
785void sySetNewHilb(syStrategy syzstr, int toSub,int index,int actord)
786{
787  int i;
788  actord += index;
789  intvec * temp_hilb = hHstdSeries(syzstr->res[index+1],NULL,NULL,NULL);
790  intvec * cont_hilb = hHstdSeries(syzstr->res[index],NULL,NULL,NULL);
791  if ((index+1<syzstr->length) && (syzstr->hilb_coeffs[index+1]==NULL))
792  {
793    syzstr->hilb_coeffs[index+1] = new intvec(16*((actord/16)+1));
794  }
795  else if (actord>=syzstr->hilb_coeffs[index+1]->length())
796  {
797    intvec * ttt=new intvec(16*((actord/16)+1));
798    for (i=syzstr->hilb_coeffs[index+1]->length()-1;i>=0;i--)
799    {
800      (*ttt)[i] = (*(syzstr->hilb_coeffs[index+1]))[i];
801    }
802    delete syzstr->hilb_coeffs[index+1];
803    syzstr->hilb_coeffs[index+1] = ttt;
804  }
805  if (actord+1<temp_hilb->length())
806  {
807#ifdef SHOW_HILB
808Print("\nSetze fuer Modul %d im Grad %d die Wert: \n",index+1,actord);
809(temp_hilb)->show(0,0);
810#endif
811    int k=si_min(temp_hilb->length()-1,(syzstr->hilb_coeffs[index+1])->length());
812    for (int j=k;j>actord;j--)
813      (*(syzstr->hilb_coeffs[index+1]))[j-1] = (*temp_hilb)[j];
814  }
815  else
816  {
817    (*(syzstr->hilb_coeffs[index+1]))[actord] = 0;
818  }
819  delete temp_hilb;
820  if ((index>1) && (actord<=syzstr->hilb_coeffs[index]->length()))
821  {
822#ifdef SHOW_HILB
823Print("\nSubtrahiere im Modul %d im Grad %d den Wert: %d\n",index,actord-1,toSub);
824#endif
825    (*syzstr->hilb_coeffs[index])[actord-1]-=toSub;
826  }
827  if (syzstr->hilb_coeffs[index]!=NULL)
828  {
829    if (cont_hilb->length()>syzstr->hilb_coeffs[index]->length())
830      syzstr->hilb_coeffs[index]->resize(cont_hilb->length());
831    for (int j=cont_hilb->length()-1;j>actord;j--)
832      (*(syzstr->hilb_coeffs[index]))[j-1] = (*cont_hilb)[j];
833  }
834  delete cont_hilb;
835#ifdef SHOW_HILB
836Print("<h,%d>",(*(syzstr->hilb_coeffs[index+1]))[actord]);
837#endif
838}
839
840/*3
841* reduces the generators of the module index in degree deg
842* (which are actual syzygies of the module index-1)
843* wrt. the ideal generated by elements of lower degrees
844*/
845static void syRedGenerOfCurrDeg_Hilb(syStrategy syzstr, int deg,int *maxindex,int *maxdeg)
846{
847  ideal res=syzstr->res[1];
848  int i=0,j,k=IDELEMS(res),k1=IDELEMS(syzstr->orderedRes[1]);
849  SSet sPairs1=syzstr->resPairs[1];
850  SSet sPairs=syzstr->resPairs[0];
851
852  while ((k>0) && (res->m[k-1]==NULL)) k--;
853  while ((k1>0) && (syzstr->orderedRes[1]->m[k1-1]==NULL)) k1--;
854  while ((i<(*syzstr->Tl)[0]) && (((sPairs)[i].syz==NULL) ||
855          ((sPairs)[i].order<deg)))
856    i++;
857  if ((i>=(*syzstr->Tl)[0]) || ((sPairs)[i].order>deg)) return;
858  while ((i<(*syzstr->Tl)[0]) && (((sPairs)[i].syz==NULL) ||
859         ((sPairs)[i].order==deg)))
860  {
861    if ((sPairs)[i].syz!=NULL)
862    {
863#ifdef SHOW_PROT
864PrintS("reduziere Erzeuger: \n");
865PrintS("syz: ");poly_write((sPairs)[i].syz);
866#endif
867      (sPairs)[i].syz = syRed_Hilb((sPairs)[i].syz,syzstr,1);
868#ifdef SHOW_PROT
869PrintS("erhalte Erzeuger: \n");
870PrintS("syz: ");poly_write((sPairs)[i].syz);
871PrintLn();
872#endif
873      if ((sPairs)[i].syz != NULL)
874      {
875        if (k==IDELEMS(res))
876        {
877          syEnlargeFields(syzstr,1);
878          res=syzstr->res[1];
879        }
880        if (TEST_OPT_DEBUG)
881        {
882          if ((sPairs)[i].isNotMinimal==NULL)
883          {
884            PrintS("\nminimal generator: ");
885            pWrite((syzstr->resPairs[0])[i].syz);
886            PrintS("comes from: ");pWrite((syzstr->resPairs[0])[i].p1);
887            PrintS("and: ");pWrite((syzstr->resPairs[0])[i].p2);
888          }
889        }
890        res->m[k] = (sPairs)[i].syz;
891        pNorm(res->m[k]);
892        syHalfPair(res->m[k],k1,syzstr,1);
893        k1++;
894        k++;
895        if (1>*maxindex) *maxindex = 1;
896        if (deg-1>*maxdeg) *maxdeg = deg-1;
897      }
898    }
899    i++;
900  }
901}
902
903/*3
904* reorders the result (stored in orderedRes) according
905* to the seqence given by res
906*/
907static void syReOrdResult_Hilb(syStrategy syzstr,int maxindex,int maxdeg)
908{
909  ideal reor,toreor;
910  int i,j,k,l,m,togo;
911  syzstr->betti = new intvec(maxdeg,maxindex+1,0);
912  (*syzstr->betti)[0] = 1;
913  for (i=1;i<=syzstr->length;i++)
914  {
915    if (!idIs0(syzstr->orderedRes[i])) 
916    {
917      toreor = syzstr->orderedRes[i];
918      k = IDELEMS(toreor);
919      while ((k>0) && (toreor->m[k-1]==NULL)) k--;
920      reor = idInit(k,toreor->rank);
921      togo = IDELEMS(syzstr->res[i]);
922      for (j=0;j<k;j++)
923      {
924        if (toreor->m[j]!=NULL) (IMATELEM(*syzstr->betti,pFDeg(toreor->m[j],currRing)-i+1,i+1))++;
925        reor->m[j] = toreor->m[j];
926        toreor->m[j] = NULL;
927      }
928      m = 0; 
929      for (j=0;j<togo;j++)
930      {
931        if (syzstr->res[i]->m[j]!=NULL)
932        {
933          l = 0;
934          while ((l<k) && (syzstr->res[i]->m[j]!=reor->m[l])) l++;
935          if (l<k)
936          {
937            toreor->m[m] = reor->m[l];
938            reor->m[l] = NULL;
939            m++;
940          }
941        }
942      }
943      idDelete(&reor);
944    }
945  }
946}
947
948/*2
949* the CoCoA-algorithm for free resolutions, using a formula
950* for remaining pairs based on Hilbert-functions
951*/
952syStrategy syHilb(ideal arg,int * length)
953{
954  int i,j,actdeg=32000,index=0,reg=-1;
955  int startdeg,howmuch,toSub=0;
956  int maxindex=0,maxdeg=0;
957  ideal temp=NULL;
958  SSet nextPairs;
959  ring origR = currRing;
960  syStrategy syzstr=(syStrategy)omAlloc0(sizeof(ssyStrategy));
961
962  if ((idIs0(arg)) || (idRankFreeModule(arg)>0))
963  {
964    syzstr->minres = (resolvente)omAllocBin(sip_sideal_bin);
965    syzstr->length = 1;
966    syzstr->minres[0] = idInit(1,arg->rank);
967    return syzstr;
968  }
969 
970  // Creare dp,S ring and change to it
971  syzstr->syRing = rCurrRingAssure_dp_C();
972
973  // set initial ShiftedComps
974  currcomponents = (int*)omAlloc0((arg->rank+1)*sizeof(int));
975  currShiftedComponents = (long*)omAlloc0((arg->rank+1)*sizeof(long));
976
977/*--- initializes the data structures---------------*/
978#ifdef SHOW_CRIT
979  crit = 0;
980  crit1 = 0;
981  spfl = 0;
982  cons_pairs = 0;
983  crit_fails = 0;
984#endif
985  syzstr->length = *length = pVariables+2;
986  syzstr->Tl = new intvec(*length+1);
987  temp = idInit(IDELEMS(arg),arg->rank);
988  for (i=0;i<IDELEMS(arg);i++)
989  {
990    if (origR != syzstr->syRing)
991      temp->m[i] = prCopyR( arg->m[i], origR);
992    else
993      temp->m[i] = pCopy( arg->m[i]);
994    if (temp->m[i]!=NULL)
995    {
996      j = pTotaldegree(temp->m[i]);
997      if (j<actdeg) actdeg = j;
998    }
999  }
1000  idTest(temp);
1001  idSkipZeroes(temp);
1002  syzstr->resPairs = syInitRes(temp,length,syzstr->Tl,syzstr->cw);
1003  omFreeSize((ADDRESS)currcomponents,(arg->rank+1)*sizeof(int));
1004  omFreeSize((ADDRESS)currShiftedComponents,(arg->rank+1)*sizeof(int));
1005  syzstr->res = (resolvente)omAlloc0((*length+1)*sizeof(ideal));
1006  syzstr->orderedRes = (resolvente)omAlloc0((*length+1)*sizeof(ideal));
1007  syzstr->elemLength = (int**)omAlloc0((*length+1)*sizeof(int*));
1008  syzstr->truecomponents = (int**)omAlloc0((*length+1)*sizeof(int*));
1009  syzstr->backcomponents = (int**)omAlloc0((*length+1)*sizeof(int*));
1010  syzstr->ShiftedComponents = (long**)omAlloc0((*length+1)*sizeof(long*));
1011  syzstr->Howmuch = (int**)omAlloc0((*length+1)*sizeof(int*));
1012  syzstr->Firstelem = (int**)omAlloc0((*length+1)*sizeof(int*));
1013  syzstr->hilb_coeffs = (intvec**)omAlloc0((*length+1)*sizeof(intvec*));
1014  syzstr->sev = (unsigned long **)omAlloc0((*length+1)*sizeof(unsigned long*));
1015  syzstr->bucket = kBucketCreate();
1016  syzstr->syz_bucket = kBucketCreate();
1017  startdeg = actdeg;
1018  nextPairs = syChosePairs(syzstr,&index,&howmuch,&actdeg);
1019/*--- computes the resolution ----------------------*/
1020  while (nextPairs!=NULL)
1021  {
1022#ifdef SHOW_PROT
1023Print("compute %d Paare im Module %d im Grad %d \n",howmuch,index,actdeg+index);
1024#endif
1025    if (TEST_OPT_PROT) Print("%d",actdeg);
1026    if (TEST_OPT_PROT) Print("(m%d)",index);
1027    if (index==0)
1028      i = syInitSyzMod(syzstr,index,idRankFreeModule(arg, origR)+1);
1029    else
1030      i = syInitSyzMod(syzstr,index);
1031    j = syInitSyzMod(syzstr,index+1);
1032    if (index>0)
1033    {
1034      syRedNextPairs_Hilb(nextPairs,syzstr,howmuch,index,actdeg,&toSub,&maxindex,&maxdeg);
1035      sySetNewHilb(syzstr,toSub,index,actdeg);
1036      toSub = 0;
1037      syCompactifyPairSet(syzstr->resPairs[index],(*syzstr->Tl)[index],0);
1038    }
1039    else
1040      syRedGenerOfCurrDeg_Hilb(syzstr,actdeg,&maxindex,&maxdeg);
1041/*--- creates new pairs -----------------------------*/
1042#ifdef SHOW_PROT
1043Print("Bilde neue Paare in Modul %d!\n",index);
1044#endif
1045    syCreateNewPairs_Hilb(syzstr,index,actdeg);
1046    if (index<(*length)-1)
1047    {
1048#ifdef SHOW_PROT
1049Print("Bilde neue Paare in Modul %d!\n",index+1);
1050#endif
1051      syCreateNewPairs_Hilb(syzstr,index+1,actdeg-1);
1052    }
1053    index++;
1054    nextPairs = syChosePairs(syzstr,&index,&howmuch,&actdeg);
1055  }
1056  syReOrdResult_Hilb(syzstr,maxindex,maxdeg);
1057#ifdef SHOW_RESULT
1058PrintS("minimal resolution:\n");
1059for (int ti=1;ti<=*length;ti++)
1060{
1061  if (!idIs0(syzstr->orderedRes[ti])) idPrint(syzstr->orderedRes[ti]);
1062}
1063PrintS("full resolution:\n");
1064for (int ti=1;ti<=*length;ti++)
1065{
1066  if (!idIs0(syzstr->res[ti])) idPrint(syzstr->res[ti]);
1067}
1068#endif
1069#ifdef SHOW_CRIT
1070Print("Criterion %d times applied\n",crit);
1071Print("Criterion1 %d times applied\n",crit1);
1072Print("%d superfluous pairs\n",spfl);
1073Print("%d pairs considered\n",cons_pairs);
1074Print("Criterion fails %d times\n",crit_fails);
1075crit = 0;
1076crit1 = 0;
1077spfl = 0;
1078cons_pairs = 0;
1079crit_fails = 0;
1080#endif
1081  if (temp!=NULL) idDelete(&temp);
1082  kBucketDestroy(&(syzstr->bucket));
1083  kBucketDestroy(&(syzstr->syz_bucket));
1084  if (origR != syzstr->syRing) 
1085    rChangeCurrRing(origR);
1086  else
1087    currRing =  origR;
1088  if (TEST_OPT_PROT) PrintLn();
1089  return syzstr;
1090}
Note: See TracBrowser for help on using the repository browser.