source: git/Singular/syz2.cc @ fdc537

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