source: git/Singular/syz2.cc @ 512a2b

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