source: git/kernel/syz2.cc @ f3a8c2e

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