source: git/kernel/GBEngine/syz3.cc @ 5aba1a

fieker-DuValspielwiese
Last change on this file since 5aba1a was 5aba1a, checked in by Hans Schoenemann <hannes@…>, 8 years ago
fix: compilation debug/opt. version
  • Property mode set to 100644
File size: 61.1 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/*
5* ABSTRACT: resolutions
6*/
7
8#include <kernel/mod2.h>
9#include <omalloc/omalloc.h>
10
11#include <misc/mylimits.h>
12#include <misc/options.h>
13#include <misc/intvec.h>
14
15#include <coeffs/numbers.h>
16
17#include <polys/monomials/ring.h>
18#include <polys/kbuckets.h>
19#include <polys/prCopy.h>
20#include <polys/matpol.h>
21
22#include <kernel/combinatorics/stairc.h>
23#include <kernel/combinatorics/hilb.h>
24
25#include <kernel/GBEngine/kstd1.h>
26#include <kernel/GBEngine/kutil.h>
27#include <kernel/GBEngine/syz.h>
28
29#include <kernel/ideals.h>
30#include <kernel/polys.h>
31
32//#define SHOW_PROT
33//#define SHOW_RED
34//#define SHOW_Kosz
35//#define SHOW_RESULT
36//#define INVERT_PAIRS
37//#define ONLY_STD
38//#define EXPERIMENT1    //Hier stimmt was mit der Anzahl der Erzeuger in xyz11 nicht!!
39#define EXPERIMENT2
40#define EXPERIMENT3
41#define WITH_BUCKET     //Use of buckets in EXPERIMENT3 (Product criterion)
42#define WITH_SCHREYER_ORD
43#define USE_CHAINCRIT
44#define USE_CHAINCRIT0
45#define USE_PROD_CRIT
46#define USE_REGULARITY
47#define WITH_SORT
48//#define FULL_TOTAKE
49int discard_pairs;
50int short_pairs;
51
52/*3
53* assumes the ideals old_ideal and new_ideal to be homogeneous
54* tests wether the new_ideal is a regular extension of the old_ideal
55*/
56static BOOLEAN syIsRegular(ideal old_ideal,ideal new_ideal,int deg)
57{
58  intvec * old_hilbs=hHstdSeries(old_ideal,NULL,NULL,NULL);
59  intvec * new_hilbs=hHstdSeries(new_ideal,NULL,NULL,NULL);
60  int biggest_length=si_max(old_hilbs->length()+deg,new_hilbs->length());
61  intvec * shifted_old_hilbs=new intvec(biggest_length);
62  intvec * old_hilb1=new intvec(biggest_length);
63  intvec * new_hilb1=new intvec(biggest_length);
64  int i;
65  BOOLEAN isRegular=TRUE;
66
67  for (i=old_hilbs->length()+deg-1;i>=deg;i--)
68    (*shifted_old_hilbs)[i] = (*old_hilbs)[i-deg];
69  for (i=old_hilbs->length()-1;i>=0;i--)
70    (*old_hilb1)[i] = (*old_hilbs)[i]-(*shifted_old_hilbs)[i];
71  for (i=old_hilbs->length()+deg-1;i>=old_hilbs->length();i--)
72    (*old_hilb1)[i] = -(*shifted_old_hilbs)[i];
73  for (i=new_hilbs->length()-1;i>=0;i--)
74    (*new_hilb1)[i] = (*new_hilbs)[i];
75  i = 0;
76  while ((i<biggest_length) && isRegular)
77  {
78    isRegular = isRegular && ((*old_hilb1)[i] == (*new_hilb1)[i]);
79    i++;
80  }
81  delete old_hilbs;
82  delete new_hilbs;
83  delete old_hilb1;
84  delete new_hilb1;
85  delete shifted_old_hilbs;
86  return isRegular;
87}
88
89/*3
90* shows the resolution stored in syzstr->orderedRes
91*/
92#if 0 /* unused*/
93static void syShowRes(syStrategy syzstr)
94{
95  int i=0;
96
97  while ((i<syzstr->length) && (!idIs0(syzstr->res[i])))
98  {
99    Print("aktueller hoechster index ist: %d\n",(*syzstr->Tl)[i]);
100    Print("der %d-te modul ist:\n",i);
101    idPrint(syzstr->res[i]);
102    PrintS("Seine Darstellung:\n");
103    idPrint(syzstr->orderedRes[i]);
104    i++;
105  }
106}
107#endif
108
109/*3
110* produces the next subresolution for a regular extension
111*/
112static void syCreateRegularExtension(syStrategy syzstr,ideal old_ideal,
113            ideal old_repr,int old_tl, poly next_generator,resolvente totake)
114{
115  int index=syzstr->length-1,i,j,start,start_ttk/*,new_tl*/;
116  poly gen=pCopy(next_generator),p;
117  poly neg_gen=pCopy(next_generator);
118  ideal current_ideal,current_repr;
119  int current_tl;
120  poly w_gen=pHead(next_generator);
121  pSetComp(w_gen,0);
122  pSetmComp(w_gen);
123
124  //syShowRes(syzstr);
125  neg_gen = pNeg(neg_gen);
126  if (pGetComp(gen)>0)
127  {
128    p_Shift(&gen,-1,currRing);
129    p_Shift(&neg_gen,-1,currRing);
130  }
131  while (index>0)
132  {
133    if (index%2==0)
134      p = gen;
135    else
136      p = neg_gen;
137    if (index>1)
138    {
139      current_ideal = syzstr->res[index-1];
140      current_repr = syzstr->orderedRes[index-1];
141      current_tl = (*syzstr->Tl)[index-1];
142    }
143    else
144    {
145      current_ideal = old_ideal;
146      current_repr = old_repr;
147      current_tl = old_tl;
148    }
149    if (!idIs0(current_ideal))
150    {
151      if (idIs0(syzstr->res[index]))
152      {
153        syzstr->res[index] = idInit(IDELEMS(current_ideal),
154          current_ideal->rank+current_tl);
155        syzstr->orderedRes[index] = idInit(IDELEMS(current_ideal),
156          current_ideal->rank);
157        start = 0;
158      }
159      else
160      {
161        start = IDELEMS(syzstr->res[index]);
162        while ((start>0) && (syzstr->res[index]->m[start-1]==NULL)) start--;
163        if (IDELEMS(syzstr->res[index])<start+IDELEMS(current_ideal))
164        {
165          pEnlargeSet(&syzstr->res[index]->m,IDELEMS(syzstr->res[index]),
166                   IDELEMS(current_ideal));
167          IDELEMS(syzstr->res[index]) += IDELEMS(current_ideal);
168          pEnlargeSet(&syzstr->orderedRes[index]->m,IDELEMS(syzstr->orderedRes[index]),
169                   IDELEMS(current_ideal));
170          IDELEMS(syzstr->orderedRes[index]) += IDELEMS(current_ideal);
171        }
172      }
173      if (idIs0(totake[index]))
174      {
175        totake[index] = idInit(IDELEMS(current_ideal),
176          current_ideal->rank+current_tl);
177        start_ttk = 0;
178      }
179      else
180      {
181        start_ttk = IDELEMS(totake[index]);
182        while ((start_ttk>0) && (totake[index]->m[start_ttk-1]==NULL)) start_ttk--;
183        if (IDELEMS(totake[index])<start_ttk+IDELEMS(current_ideal))
184        {
185          pEnlargeSet(&totake[index]->m,IDELEMS(totake[index]),
186                   IDELEMS(current_ideal));
187          for (j=IDELEMS(totake[index]);j<IDELEMS(totake[index])+
188                                  IDELEMS(current_ideal);j++)
189            totake[index]->m[j] = NULL;
190          IDELEMS(totake[index]) += IDELEMS(current_ideal);
191        }
192      }
193      for (i=0;i<IDELEMS(current_ideal);i++)
194      {
195        if (current_ideal->m[i]!=NULL)
196        {
197          syzstr->res[index]->m[i+start] = pCopy(current_ideal->m[i]);
198          syzstr->res[index]->m[i+start] = pMult_mm(syzstr->res[index]->m[i+start],w_gen);
199          p_Shift(&syzstr->res[index]->m[i+start],current_tl,currRing);
200          syzstr->res[index]->m[i+start] = pAdd(syzstr->res[index]->m[i+start],
201            ppMult_qq(current_repr->m[i],p));
202          syzstr->orderedRes[index]->m[i+start] = pCopy(current_repr->m[i]);
203          syzstr->orderedRes[index]->m[i+start] =
204            pMult_mm(syzstr->orderedRes[index]->m[i+start],w_gen);
205          if ((*syzstr->Tl)[index]!=0)
206            p_Shift(&syzstr->orderedRes[index]->m[i+start],(*syzstr->Tl)[index],currRing);
207        }
208      }
209      for (i=0;i<IDELEMS(totake[index-1]);i++)
210      {
211        if (totake[index-1]->m[i]!=NULL)
212        {
213          if ((index==1) && ((i==IDELEMS(current_ideal) ||
214               (totake[index-1]->m[i+1]==NULL)))) break;
215          totake[index]->m[i+start_ttk] =
216            pMult_mm(pCopy(totake[index-1]->m[i]),w_gen);
217          p_Shift(&totake[index]->m[i+start_ttk],current_tl,currRing);
218#ifdef FULL_TOTAKE
219          poly pp=pCopy(p);
220          p_Shift(&pp,i+1,currRing);
221          totake[index]->m[i+start_ttk] = pAdd(totake[index]->m[i+start_ttk],pp);
222#endif
223        }
224      }
225      (*syzstr->Tl)[index] += current_tl;
226    }
227    index--;
228  }
229  pDelete(&gen);
230  pDelete(&neg_gen);
231  pDelete(&w_gen);
232  //syShowRes(syzstr);
233}
234
235/*3
236* proves the consistence of the pairset resPairs with the corresponding
237* set of generators;
238* only for tests
239*/
240static void syTestPairs(SSet resPairs,int length,ideal old_generators)
241{
242  int i=0;
243
244  while (i<length)
245  {
246    if (resPairs[i].lcm!=NULL)
247    {
248      if (resPairs[i].p1!=NULL)
249        assume(resPairs[i].p1==old_generators->m[resPairs[i].ind1]);
250      if (resPairs[i].p2!=NULL)
251        assume(resPairs[i].p2==old_generators->m[resPairs[i].ind2]);
252    }
253    i++;
254  }
255}
256
257/*3
258* cancels the weight monomials given by the leading terms of totake
259* from the resolution res;
260* works in place on res, but reads only from totake
261*/
262void syReorder_Kosz(syStrategy syzstr)
263{
264  int length=syzstr->length;
265  int syzIndex=length-1,i,j;
266  resolvente res=syzstr->fullres;
267  poly p;
268
269  while ((syzIndex!=0) && (res[syzIndex]==NULL)) syzIndex--;
270  while (syzIndex>0)
271  {
272    for(i=0;i<IDELEMS(res[syzIndex]);i++)
273    {
274#ifdef USE_REGULARITY
275      if ((syzstr->regularity>0) && (res[syzIndex]->m[i]!=NULL))
276      {
277        if (p_FDeg(res[syzIndex]->m[i],currRing)>=syzstr->regularity+syzIndex)
278          pDelete(&res[syzIndex]->m[i]);
279      }
280#endif
281      p = res[syzIndex]->m[i];
282      while (p!=NULL)
283      {
284        if (res[syzIndex-1]->m[pGetComp(p)-1]!=NULL)
285        {
286          for(j=1;j<=(currRing->N);j++)
287          {
288            pSetExp(p,j,pGetExp(p,j)
289                        -pGetExp(res[syzIndex-1]->m[pGetComp(p)-1],j));
290          }
291        }
292        else
293          PrintS("error in the resolvent\n");
294        pSetm(p);
295        pIter(p);
296      }
297    }
298    syzIndex--;
299  }
300}
301
302/*3
303* updates the pairset resPairs by generating all pairs including the
304* new_generators in the 0-th modul;
305* the new_generators are inserted in the old_generators;
306* new_generators is empty after the procedure;
307*/
308static void updatePairs(SSet *resPairs,int *l_pairs,syStrategy syzstr,
309       int index,ideal new_generators,ideal new_repr,int crit_comp)
310{
311  if (idIs0(new_generators)) return;
312  ideal old_generators=syzstr->res[index];
313  ideal old_repr=syzstr->orderedRes[index];
314  int i=0,j,k,kk,og_elem=0,og_idel=IDELEMS(old_generators),l=*l_pairs,jj,ll,j1;
315  int og_ini=0;
316  ideal pairs=idInit(og_idel+IDELEMS(new_generators),old_generators->rank);
317  polyset prs=pairs->m;
318  poly p=NULL;
319  SObject tso;
320
321  syInitializePair(&tso);
322  while ((og_elem<og_idel) && (old_generators->m[og_elem]!=NULL))
323  {
324    if ((index>0) && (pGetComp(old_generators->m[og_elem])<=crit_comp))
325      og_ini = og_elem;
326    og_elem++;
327  }
328  while ((l>0) && ((*resPairs)[l-1].lcm==NULL)) l--;
329  while ((i<IDELEMS(new_generators)) && (new_generators->m[i]!=NULL))
330  {
331    syTestPairs(*resPairs,*l_pairs,old_generators);
332    if (IDELEMS(old_generators)==og_elem)
333    {
334      pEnlargeSet(&old_generators->m,IDELEMS(old_generators),16);
335      IDELEMS(old_generators) += 16;
336      pEnlargeSet(&old_repr->m,IDELEMS(old_repr),16);
337      IDELEMS(old_repr) += 16;
338    }
339    k = p_FDeg(new_generators->m[i],currRing);
340    kk = pGetComp(new_generators->m[i]);
341    j = og_ini;
342    while ((j<og_elem) && (old_generators->m[j]!=NULL) &&
343           (pGetComp(old_generators->m[j])<kk)) j++;
344    while ((j<og_elem) && (old_generators->m[j]!=NULL) &&
345           (p_FDeg(old_generators->m[j],currRing)<=k)) j++;
346    for (jj=og_elem;jj>j;jj--)
347    {
348      old_generators->m[jj] = old_generators->m[jj-1];
349      old_repr->m[jj] = old_repr->m[jj-1];
350    }
351    old_generators->m[j] = new_generators->m[i];
352    new_generators->m[i] = NULL;
353    old_repr->m[j] = new_repr->m[i];
354    new_repr->m[i] = NULL;
355    og_elem++;
356    for (jj=0;jj<*l_pairs;jj++)
357    {
358      if ((*resPairs)[jj].lcm!=NULL)
359      {
360        if ((*resPairs)[jj].ind1>=j) (*resPairs)[jj].ind1++;
361        if ((*resPairs)[jj].ind2>=j) (*resPairs)[jj].ind2++;
362      }
363    }
364    syTestPairs(*resPairs,*l_pairs,old_generators);
365    for (jj=og_ini;jj<og_elem;jj++)
366    {
367      if ((j!=jj) && (pGetComp(old_generators->m[jj])==pGetComp(old_generators->m[j])))
368      {
369        p = pOne();
370        pLcm(old_generators->m[jj],old_generators->m[j],p);
371        pSetComp(p,j+1);
372        pSetm(p);
373        j1 = 0;
374        while (j1<jj)
375        {
376          if (prs[j1]!=NULL)
377          {
378            if (pLmDivisibleByNoComp(prs[j1],p))
379            {
380              pDelete(&p);
381              break;
382            }
383            else if (pLmDivisibleByNoComp(p,prs[j1]))
384            {
385              pDelete(&(prs[j1]));
386            }
387#ifdef USE_CHAINCRIT0
388            else
389            {
390              poly p1,p2;
391              int ip=(currRing->N);
392              p1 = pDivide(p,old_generators->m[jj]);
393              p2 = pDivide(prs[j1],old_generators->m[j1]);
394              while ((ip>0) && (pGetExp(p1,ip)*pGetExp(p2,ip)==0)) ip--;
395              if (ip==0)
396              {
397                int ti=0;
398                while ((ti<l) && (((*resPairs)[ti].ind1!=j1)|| ((*resPairs)[ti].ind2!=jj))) ti++;
399                if (ti<l)
400                {
401                  if (TEST_OPT_PROT) PrintS("cc");
402                  syDeletePair(&(*resPairs)[ti]);
403                  syCompactifyPairSet(*resPairs,*l_pairs,ti);
404                  l--;
405                }
406              }
407              pDelete(&p1);
408              pDelete(&p2);
409            }
410#endif
411          }
412          j1++;
413        }
414        if (p!=NULL)
415          prs[jj] = p;
416      }
417    }
418    for (jj=og_ini;jj<og_elem;jj++)
419    {
420      if (prs[jj] !=NULL)
421      {
422        if (l>=*l_pairs)
423        {
424          SSet temp = (SSet)omAlloc0((*l_pairs+16)*sizeof(SObject));
425          for (ll=0;ll<*l_pairs;ll++)
426          {
427            temp[ll].p = (*resPairs)[ll].p;
428            temp[ll].p1 = (*resPairs)[ll].p1;
429            temp[ll].p2 = (*resPairs)[ll].p2;
430            temp[ll].syz = (*resPairs)[ll].syz;
431            temp[ll].lcm = (*resPairs)[ll].lcm;
432            temp[ll].ind1 = (*resPairs)[ll].ind1;
433            temp[ll].ind2 = (*resPairs)[ll].ind2;
434            temp[ll].syzind = (*resPairs)[ll].syzind;
435            temp[ll].order = (*resPairs)[ll].order;
436            temp[ll].isNotMinimal = (*resPairs)[ll].isNotMinimal;
437          }
438          omFreeSize((ADDRESS)(*resPairs),*l_pairs*sizeof(SObject));
439          *l_pairs += 16;
440          (*resPairs) = temp;
441        }
442        tso.lcm = prs[jj];
443        prs[jj] = NULL;
444        tso.order = p_FDeg(tso.lcm,currRing);
445        tso.p1 = old_generators->m[jj];
446        tso.p2 = old_generators->m[j];
447        tso.ind1 = jj;
448        tso.ind2 = j;
449        tso.syzind = -1;
450        tso.isNotMinimal = NULL;
451        tso.p = NULL;
452        tso.syz = NULL;
453        SSet rP=*resPairs;
454#ifdef SHOW_PROT
455Print("erzeuge Paar im Modul %d,%d mit: \n",index,tso.order);
456PrintS("poly1: ");pWrite(tso.p1);
457PrintS("poly2: ");pWrite(tso.p2);
458PrintS("syz: ");pWrite(tso.syz);
459PrintS("sPoly: ");pWrite(tso.p);
460PrintLn();
461#endif
462        syEnterPair(rP,&tso,&l,index);
463        syInitializePair(&tso);
464      }
465    }
466    i++;
467  }
468  idDelete(&pairs);
469}
470
471/*3
472* performs the modification of a single reduction on the syzygy-level
473*/
474inline void sySPRedSyz_Kosz(syStrategy syzstr,poly redWith,poly syz,poly q=NULL,int l_syz=-1)
475{
476  poly p=pDivide(q,redWith);
477  pSetCoeff(p,nDiv(pGetCoeff(q),pGetCoeff(redWith)));
478  kBucket_Minus_m_Mult_p(syzstr->syz_bucket,p,syz,&l_syz,NULL);
479  pDelete(&p);
480}
481
482/*3
483* normalizes the poly bucket by the ideal;
484* stops the reduction whenever the leading component is less than the
485* crit_comp;
486* returns the changing status
487*/
488static BOOLEAN syRedSyz(kBucket_pt bucket,ideal red,int crit_comp,int* g_l)
489{
490  poly p = kBucketGetLm(bucket);
491  int j = 0,i=IDELEMS(red)-1;
492  number n;
493  BOOLEAN isChanged=FALSE;
494
495  loop
496  {
497    if ((j>=i) || (p==NULL) || (pGetComp(p)<=crit_comp)) break;
498    if ((red->m[j]!=NULL) && (pDivisibleBy(red->m[j],p)))
499    {
500      n = kBucketPolyRed(bucket,red->m[j], g_l[j], NULL);
501      nDelete(&n);
502      p = kBucketGetLm(bucket);
503      isChanged = TRUE;
504      j = 0;
505    }
506    else
507      j++;
508  }
509  return isChanged;
510}
511
512/*3
513* a tail reduction for the syzygies yielding new generators
514*/
515static poly syRedTailSyz(poly tored,ideal red,ideal sec_red,int crit_comp,syStrategy syzstr,
516            int * gen_length,int * secgen_length,int * tored_length)
517{
518  int i=IDELEMS(red)-1,num_mon,num_tail;
519  poly h,hn;
520  // BOOLEAN dummy;
521
522  while ((i>0) && (red->m[i-1]==NULL)) i--;
523  i--;
524  h = tored;
525  if ((h!=NULL) && (pGetComp(h)>crit_comp))
526  {
527    num_mon = 1;
528    hn = pNext(h);
529    num_tail = *tored_length-1;
530    while (hn!=NULL)
531    {
532      kBucketInit(syzstr->syz_bucket,hn,num_tail);
533      /*dummy =*/ (void) syRedSyz(syzstr->syz_bucket,red,crit_comp,gen_length);
534      kBucketClear(syzstr->syz_bucket,&hn,&num_tail);
535      pNext(h) = hn;
536      if ((hn==NULL) || (pGetComp(hn)<=crit_comp))
537        break;
538      else
539      {
540        pIter(h);
541        pIter(hn);
542        num_mon++;
543        num_tail--;
544      }
545    }
546    if (sec_red!=NULL)
547    {
548      while (hn!=NULL)
549      {
550        kBucketInit(syzstr->syz_bucket,hn,num_tail);
551        /*dummy =*/ (void) syRedSyz(syzstr->syz_bucket,sec_red,crit_comp,secgen_length);
552        kBucketClear(syzstr->syz_bucket,&hn,&num_tail);
553        pNext(h) = hn;
554        if (hn==NULL)
555          break;
556        else
557        {
558          pIter(h);
559          pIter(hn);
560          num_mon++;
561          num_tail--;
562        }
563      }
564    }
565    *tored_length = num_mon+num_tail;
566  }
567  assume(pLength(tored)==*tored_length);
568  return tored;
569}
570
571/*3
572* the complete reduction of a single pair which is just stored
573* in bucket and syz_bucket
574*/
575static BOOLEAN syRedSyzPair(syStrategy syzstr,int index,int* g_l,int* orp_l)
576{
577  kBucket_pt bucket=syzstr->bucket;
578  poly p = kBucketGetLm(bucket);
579  ideal red=syzstr->res[index],repr=syzstr->orderedRes[index];
580  int j = 0,i=IDELEMS(red)-1;
581  number n;
582  BOOLEAN isChanged=FALSE;
583
584  loop
585  {
586    if ((j>=i) || (p==NULL)) break;
587    if ((red->m[j]!=NULL) && (pDivisibleBy(red->m[j],p)))
588    {
589      sySPRedSyz_Kosz(syzstr,red->m[j],repr->m[j],p,orp_l[j]);
590      n = kBucketPolyRed(bucket,red->m[j], g_l[j], NULL);
591      nDelete(&n);
592      p = kBucketGetLm(bucket);
593      isChanged = TRUE;
594      j = 0;
595    }
596    else
597      j++;
598  }
599  return isChanged;
600}
601
602/*3
603* the tailreduction for generators (which includes the correction of
604* the corresponding representation)
605*/
606#if 0 /*unused*/
607static void syRedTailSyzPair(SObject tso,syStrategy syzstr,int index,
608            int * gen_length,int* orp_l,int * tored_l,int * syzred_l)
609{
610  int num_mon,num_tail,syz_l;
611  poly h,hn;
612  BOOLEAN dummy;
613
614  h = tso.p;
615  kBucketInit(syzstr->syz_bucket,tso.syz,*syzred_l);
616  if (h!=NULL)
617  {
618    num_mon = 1;
619    hn = pNext(h);
620    num_tail = *tored_l-1;
621    while (hn!=NULL)
622    {
623      kBucketInit(syzstr->bucket,hn,num_tail);
624      dummy = syRedSyzPair(syzstr,index,gen_length,orp_l);
625      kBucketClear(syzstr->bucket,&hn,&num_tail);
626      pNext(h) = hn;
627      if (hn==NULL)
628        break;
629      else
630      {
631        pIter(h);
632        pIter(hn);
633        num_mon++;
634        num_tail--;
635      }
636    }
637    *tored_l = num_mon+num_tail;
638  }
639  kBucketClear(syzstr->syz_bucket,&tso.syz,&syz_l);
640  assume(pLength(tso.syz)==syz_l);
641  assume(pLength(tso.p)==*tored_l);
642}
643#endif
644
645/*3
646* the reduction of a pair in the 0-th module
647*/
648static void redOnePair(SSet resPairs,int itso,int l, ideal syzygies,
649            int crit_comp, syStrategy syzstr,int index,ideal new_generators,
650            ideal new_repr,int * ogm_l,int * orp_l)
651{
652  SObject tso = resPairs[itso];
653  assume (tso.lcm!=NULL);
654  ideal old_generators=syzstr->res[index];
655  ideal old_repr=syzstr->orderedRes[index];
656  int og_idel=IDELEMS(old_generators),ng_place=IDELEMS(new_generators);
657  int toReplace=0;
658  int i,j,syz_l;
659  number /*coefgcd,*/n;
660  polyset ogm=old_generators->m;
661  poly p;
662  BOOLEAN deleteP=FALSE;
663#ifdef EXPERIMENT1
664  poly syzp;
665#endif
666  int syz_place=IDELEMS(syzygies);
667
668  while ((syz_place>0) && (syzygies->m[syz_place-1]==NULL)) syz_place--;
669  while ((ng_place>0) && (new_generators->m[ng_place-1]==NULL)) ng_place--;
670  while ((og_idel>0) && (old_generators->m[og_idel-1]==NULL)) og_idel--;
671  assume (tso.ind1<og_idel);
672  assume (tso.ind2<og_idel);
673  assume (tso.ind1!=tso.ind2);
674  assume (tso.p1 == old_generators->m[tso.ind1]);
675  assume (tso.p2 == old_generators->m[tso.ind2]);
676  tso.p1 = old_generators->m[tso.ind1];
677  tso.p2 = old_generators->m[tso.ind2];
678  if ((tso.p1!=NULL) && (tso.p2!=NULL))
679  {
680    if (TEST_OPT_PROT)
681      PrintS(".");
682    if (index==0)
683    {
684/*--- tests wether a generator must be replaced (lt(f1)|lt(f2)!)--*/
685      if (p_FDeg(tso.p1,currRing)==p_FDeg(tso.lcm,currRing))
686        toReplace = tso.ind1+1;
687      else if (p_FDeg(tso.p2,currRing)==p_FDeg(tso.lcm,currRing))
688        toReplace = tso.ind2+1;
689    }
690#ifdef EXPERIMENT3
691/*--- tests wether the product criterion applies --------------*/
692    if ((index==0) && (old_generators->rank==1) &&
693        (p_FDeg(tso.p1,currRing)+p_FDeg(tso.p2,currRing)==tso.order))
694    {
695      tso.p = NULL;
696      p = pCopy(tso.p1);
697      p_Shift(&p,-1,currRing);
698#ifdef WITH_BUCKET
699      poly pp;
700      pp = pMult_mm(pCopy(old_repr->m[tso.ind2]),p);
701      kBucketInit(syzstr->syz_bucket,pp,-1);
702      pLmDelete(&p);
703      p = pNeg(p);
704      pp = pCopy(old_repr->m[tso.ind2]);
705      int il=-1;
706      while (p!=NULL)
707      {
708        kBucket_Minus_m_Mult_p(syzstr->syz_bucket,p,pp,&il,NULL);
709        pLmDelete(&p);
710      }
711      pDelete(&pp);
712      p = pCopy(tso.p2);
713      p_Shift(&p,-1,currRing);
714      pp = pCopy(old_repr->m[tso.ind1]);
715      il=-1;
716      while (p!=NULL)
717      {
718        kBucket_Minus_m_Mult_p(syzstr->syz_bucket,p,pp,&il,NULL);
719        pLmDelete(&p);
720      }
721      pDelete(&pp);
722      kBucketClear(syzstr->syz_bucket,&tso.syz,&j);
723#else
724      tso.syz = pMult(p,pCopy(old_repr->m[tso.ind2]));
725      p = pCopy(tso.p2);
726      p_Shift(&p,-1,currRing);
727      tso.syz = pSub(tso.syz,pMult(p,pCopy(old_repr->m[tso.ind1])));
728#endif
729    }
730    else
731#endif
732/*--- the product criterion does not apply --------------------*/
733    {
734      tso.p = ksOldCreateSpoly(tso.p2,tso.p1);
735      number coefgcd = n_Gcd(pGetCoeff(tso.p1),pGetCoeff(tso.p2),currRing->cf);
736      assume (old_repr->m[tso.ind1]!=NULL);
737      tso.syz = pCopy(old_repr->m[tso.ind1]);
738      poly tt = pDivide(tso.lcm,tso.p1);
739      pSetComp(tt,0);
740      pSetmComp(tt);
741      pSetCoeff(tt,nDiv(pGetCoeff(tso.p1),coefgcd));
742      tso.syz = pMult_mm(tso.syz,tt);
743      pDelete(&tt);
744      coefgcd = nInpNeg(coefgcd);
745      assume (old_repr->m[tso.ind2]!=NULL);
746      p = pCopy(old_repr->m[tso.ind2]);
747      tt = pDivide(tso.lcm,tso.p2);
748      pSetComp(tt,0);
749      pSetmComp(tt);
750      pSetCoeff(tt,nDiv(pGetCoeff(tso.p2),coefgcd));
751      p = pMult_mm(p,tt);
752      pDelete(&tt);
753      tso.syz = pAdd(p,tso.syz);
754#ifdef EXPERIMENT2
755      if ((tso.syz!=NULL) && (pGetComp(tso.syz)<=crit_comp))
756      {
757/*--- breaks when the leading component is less than crit_comp ------*/
758        deleteP = TRUE;
759        discard_pairs++;
760      }
761#endif
762      nDelete(&coefgcd);
763    }                             //End of the else-part of EXPERIMENT3
764#ifdef SHOW_PROT
765Print("reduziere Paar im Module %d mit: \n",index);
766PrintS("poly1: ");pWrite(tso.p1);
767PrintS("poly2: ");pWrite(tso.p2);
768PrintS("syz: ");pWrite(tso.syz);
769PrintS("sPoly: ");pWrite(tso.p);
770#endif
771    assume(tso.syz!=NULL);
772    kBucketInit(syzstr->syz_bucket,tso.syz,-1);
773    if ((tso.p!=NULL) && (!deleteP))
774    {
775      kBucketInit(syzstr->bucket,tso.p,-1);
776      p = kBucketGetLm(syzstr->bucket);
777      j = 0;
778      loop
779      {
780        if (j>=og_idel)
781        {
782/*--- reduction with generators computed in this procedure ---*/
783          j = 0;
784          while ((j<ng_place) && (!pDivisibleBy(new_generators->m[j],p))) j++;
785          if (j>=ng_place) break;
786          assume (new_repr->m[j]!=NULL);
787          sySPRedSyz_Kosz(syzstr,new_generators->m[j],new_repr->m[j],p);
788          n = kBucketPolyRed(syzstr->bucket,new_generators->m[j],pLength(new_generators->m[j]), NULL);
789          p = kBucketGetLm(syzstr->bucket);
790#ifdef EXPERIMENT1
791          syzp = kBucketGetLm(syzstr->syz_bucket);
792          if ((syzp!=NULL) && (pGetComp(syzp)<=crit_comp))
793          {
794            deleteP =TRUE;
795            break;
796          }
797          //if (syzp==NULL)
798            //assume(p==NULL);
799          //else
800            //if (pGetComp(syzp)<=crit_comp) short_pairs++;
801#endif
802          if (p==NULL) break;
803          j = 0;
804        }
805        if (pDivisibleBy(ogm[j],p))
806        {
807/*--- reduction with general old generators ---------------------*/
808          assume (old_repr->m[j]!=NULL);
809          sySPRedSyz_Kosz(syzstr,ogm[j],old_repr->m[j],p,orp_l[j]);
810          n = kBucketPolyRed(syzstr->bucket,ogm[j],ogm_l[j], NULL);
811          p = kBucketGetLm(syzstr->bucket);
812#ifdef EXPERIMENT1
813          syzp = kBucketGetLm(syzstr->syz_bucket);
814          if ((syzp!=NULL) && (pGetComp(syzp)<=crit_comp))
815          {
816            break;
817            deleteP =TRUE;
818          }
819          //if (syzp==NULL)
820            //assume(p==NULL);
821          //else
822            //if ((pGetComp(syzp)<=crit_comp) && (p!=NULL)) short_pairs++;
823#endif
824          if (p==NULL) break;
825          j = 0;
826        }
827        else
828          j++;
829      }
830      kBucketClear(syzstr->bucket,&tso.p,&tso.length);
831    }
832    kBucketClear(syzstr->syz_bucket,&tso.syz,&syz_l);
833    if (deleteP)
834    {
835      pDelete(&tso.p);
836      pDelete(&tso.syz);
837    }
838  }
839  else
840  {
841    PrintS("Shit happens!\n");
842  }
843#ifdef SHOW_PROT
844Print("erhalte Paar im Module %d mit: \n",index);
845PrintS("syz: ");pWrite(tso.syz);
846PrintS("sPoly: ");pWrite(tso.p);
847PrintLn();
848#endif
849  if (toReplace)
850  {
851/*-- replaces the generator if neccesary ------------------*/
852    pDelete(&old_generators->m[toReplace-1]);
853    pDelete(&old_repr->m[toReplace-1]);
854    for (i=toReplace-1;i<og_idel-1;i++)
855    {
856      old_generators->m[i] = old_generators->m[i+1];
857      old_repr->m[i] = old_repr->m[i+1];
858    }
859    old_generators->m[og_idel-1] = NULL;
860    old_repr->m[og_idel-1] = NULL;
861    for (i=itso+1;i<l;i++)
862    {
863      if (resPairs[i].lcm!=NULL)
864      {
865        if ((resPairs[i].ind1==toReplace-1)||(resPairs[i].ind2==toReplace-1))
866          syDeletePair(&resPairs[i]);
867        else
868        {
869          if (resPairs[i].ind1>=toReplace)
870            (resPairs[i].ind1)--;
871          if (resPairs[i].ind2>=toReplace)
872            (resPairs[i].ind2)--;
873        }
874      }
875    }
876    syCompactifyPairSet(resPairs,l,itso+1);
877  }
878  if (tso.p!=NULL)
879  {
880/*-- stores the new generator ---------------------------------*/
881    //syRedTailSyzPair(tso,syzstr,index,ogm_l,orp_l,&tso.length,&syz_l);
882    if (ng_place>=IDELEMS(new_generators))
883    {
884      pEnlargeSet(&new_generators->m,IDELEMS(new_generators),16);
885      IDELEMS(new_generators) += 16;
886      pEnlargeSet(&new_repr->m,IDELEMS(new_repr),16);
887      IDELEMS(new_repr) += 16;
888    }
889    if (!nIsOne(pGetCoeff(tso.p)))
890    {
891      n=nInvers(pGetCoeff(tso.p));
892      pNorm(tso.p);
893      pMult_nn(tso.syz,n);
894      nDelete(&n);
895    }
896    new_generators->m[ng_place] = tso.p;
897    tso.p = NULL;
898    new_repr->m[ng_place] = tso.syz;
899    tso.syz = NULL;
900  }
901  else
902  {
903/*--- takes the syzygy as new generator of the next module ---*/
904    if (tso.syz==NULL)
905    {
906#ifndef EXPERIMENT2
907#ifdef EXPERIMENT3
908      short_pairs++;
909#endif
910#endif
911    }
912    else if (pGetComp(tso.syz)<=crit_comp)
913    {
914      pDelete(&tso.syz);
915    }
916    else
917    {
918      if (syz_place>=IDELEMS(syzygies))
919      {
920        pEnlargeSet(&syzygies->m,IDELEMS(syzygies),16);
921        IDELEMS(syzygies) += 16;
922      }
923      syzygies->m[syz_place] = tso.syz;
924      tso.syz = NULL;
925      pNorm(syzygies->m[syz_place]);
926    }
927  }
928  resPairs[itso] = tso;
929  syDeletePair(&resPairs[itso]);
930  syTestPairs(resPairs,l,old_generators);
931}
932
933/*3
934* reduction of all pairs of a fixed degree of the 0-th module
935*/
936static BOOLEAN redPairs(SSet resPairs,int l_pairs, ideal syzygies,
937  ideal new_generators,ideal new_repr, int crit_comp,syStrategy syzstr,
938  int index)
939{
940  if (resPairs[0].lcm==NULL) return TRUE;
941  int i,j,actdeg=resPairs[0].order;
942  int * ogm_l=(int*)omAlloc0(IDELEMS(syzstr->res[index])*sizeof(int));
943  int * orp_l=(int*)omAlloc0(IDELEMS(syzstr->orderedRes[index])*sizeof(int));
944  // int t1=IDELEMS(syzstr->res[index]),t2=IDELEMS(syzstr->orderedRes[index]);
945
946  for (j=IDELEMS(syzstr->res[index])-1;j>=0;j--)
947  {
948    if (syzstr->res[index]->m[j]!=NULL)
949      ogm_l[j] = pLength(syzstr->res[index]->m[j]);
950  }
951  for (j=IDELEMS(syzstr->orderedRes[index])-1;j>=0;j--)
952  {
953    if (syzstr->orderedRes[index]->m[j]!=NULL)
954      orp_l[j] = pLength(syzstr->orderedRes[index]->m[j]);
955  }
956  loop
957  {
958    i = 0;
959    if (TEST_OPT_PROT)
960      Print("(%d,%d)",index,resPairs[0].order);
961    while (resPairs[i].order==actdeg)
962    {
963      syTestPairs(resPairs,l_pairs,syzstr->res[index]);
964      redOnePair(resPairs,i,l_pairs,syzygies,crit_comp,syzstr,index,
965                 new_generators, new_repr,ogm_l,orp_l);
966      i++;
967      syTestPairs(resPairs,l_pairs,syzstr->res[index]);
968    }
969    syTestPairs(resPairs,l_pairs,syzstr->res[index]);
970    syCompactifyPairSet(resPairs,l_pairs,0);
971    syTestPairs(resPairs,l_pairs,syzstr->res[index]);
972    if (!idIs0(new_generators))
973      break;
974    else if (resPairs[0].lcm==NULL)  //there are no pairs left and no new_gens
975    {
976      omFreeSize((ADDRESS)ogm_l,IDELEMS(syzstr->res[index])*sizeof(int));
977      omFreeSize((ADDRESS)orp_l,IDELEMS(syzstr->orderedRes[index])*sizeof(int));
978      return TRUE;
979    }
980    else
981      actdeg = resPairs[0].order;
982  }
983  syTestPairs(resPairs,l_pairs,syzstr->res[index]);
984  omFreeSize((ADDRESS)ogm_l,IDELEMS(syzstr->res[index])*sizeof(int));
985  omFreeSize((ADDRESS)orp_l,IDELEMS(syzstr->orderedRes[index])*sizeof(int));
986  return FALSE;
987}
988
989/*3
990* extends the standard basis old_generators with new_generators;
991* returns the syzygies which involve the new elements;
992* assumes that the components of the new_generators are sperated
993* from those of old_generators, i.e. whenever the leading term
994* of a syzygy lies in the part of the old_generators, the syzygy
995* lie just in the module old_generators
996* assumes that the new_generators are reduced w.r.t. old_generators
997*/
998static ideal kosz_std(ideal new_generators,ideal new_repr,syStrategy syzstr,
999                      int index,int next_comp)
1000{
1001  int og_idel=IDELEMS(syzstr->res[index]);
1002  int l_pairs=2*og_idel;
1003  ideal syzygies=idInit(16,syzstr->res[index]->rank+1);
1004  if ((idIs0(new_generators)) || (new_generators->m[0]==NULL))
1005  {
1006    Werror("Hier ist was faul!\n");
1007    return NULL;
1008  }
1009  SSet resPairs=(SSet)omAlloc0(l_pairs*sizeof(SObject));
1010  loop
1011  {
1012    updatePairs(&resPairs,&l_pairs,syzstr,index,
1013                new_generators,new_repr,next_comp);
1014    if (redPairs(resPairs,l_pairs,syzygies, new_generators,new_repr,
1015                 next_comp,syzstr,index)) break;
1016  }
1017  omFreeSize((SSet)resPairs,l_pairs*sizeof(SObject));
1018  return syzygies;
1019}
1020
1021/*3
1022* normalizes the incoming generators
1023*/
1024static poly normalize(poly next_p,ideal add_generators, syStrategy syzstr,
1025                      int * g_l,int * p_l,int crit_comp)
1026{
1027  int j=0,i=IDELEMS(add_generators);
1028  kBucketInit(syzstr->bucket,next_p,pLength(next_p));
1029  poly p = kBucketGetLm(syzstr->bucket),result;
1030  number n;
1031
1032  loop
1033  {
1034    if ((j>=i) || (p==NULL) || (pGetComp(p)<=crit_comp)) break;
1035    if ((add_generators->m[j]!=NULL) && (pDivisibleBy(add_generators->m[j],p)))
1036    {
1037      n = kBucketPolyRed(syzstr->bucket,add_generators->m[j], g_l[j], NULL);
1038      nDelete(&n);
1039      p = kBucketGetLm(syzstr->bucket);
1040      j = 0;
1041    }
1042    else
1043      j++;
1044  }
1045  kBucketClear(syzstr->bucket,&result,p_l);
1046  return result;
1047}
1048
1049/*3
1050* updates the pairs inthe higher modules
1051*/
1052static void updatePairsHIndex(SSet *resPairs,int *l_pairs,syStrategy /*syzstr*/,
1053       int index,ideal add_generators,ideal /*add_repr*/,ideal /*new_generators*/,
1054       ideal /*new_repr*/,int /*crit_comp*/,int* first_new)
1055{
1056  int i=*first_new,l=*l_pairs,j,ll,j1,add_idel=IDELEMS(add_generators);
1057  ideal pairs=idInit(add_idel,add_generators->rank);
1058  polyset prs=pairs->m;
1059  poly p=NULL;
1060  SObject tso;
1061
1062  syInitializePair(&tso);
1063  while ((l>0) && ((*resPairs)[l-1].lcm==NULL)) l--;
1064  while ((i<add_idel) && (add_generators->m[i]!=NULL))
1065  {
1066    for (j=0;j<i;j++)
1067    {
1068      if (pGetComp(add_generators->m[j]) == pGetComp(add_generators->m[i]))
1069      {
1070        p = pOne();
1071        pLcm(add_generators->m[j],add_generators->m[i],p);
1072        pSetComp(p,i+1);
1073        pSetm(p);
1074        j1 = 0;
1075        while (j1<j)
1076        {
1077          if (prs[j1]!=NULL)
1078          {
1079            if (pLmDivisibleByNoComp(prs[j1],p))
1080            {
1081              pDelete(&p);
1082              break;
1083            }
1084            else if (pLmDivisibleByNoComp(p,prs[j1]))
1085            {
1086              pDelete(&(prs[j1]));
1087            }
1088#ifdef USE_CHAINCRIT
1089            else
1090            {
1091              poly p1,p2;
1092              int ip=(currRing->N);
1093              p1 = pDivide(p,add_generators->m[j]);
1094              p2 = pDivide(prs[j1],add_generators->m[j1]);
1095              while ((ip>0) && (pGetExp(p1,ip)*pGetExp(p2,ip)==0)) ip--;
1096              if (ip==0)
1097              {
1098                int ti=0;
1099                while ((ti<l) && (((*resPairs)[ti].ind1!=j1)|| ((*resPairs)[ti].ind2!=j))) ti++;
1100                if (ti<l)
1101                {
1102                  if (TEST_OPT_PROT) PrintS("cc");
1103                  syDeletePair(&(*resPairs)[ti]);
1104                  syCompactifyPairSet(*resPairs,*l_pairs,ti);
1105                  l--;
1106                }
1107              }
1108              pDelete(&p1);
1109              pDelete(&p2);
1110            }
1111#endif
1112          }
1113          j1++;
1114        }
1115        if (p!=NULL)
1116          prs[j] = p;
1117      }
1118    }
1119    for (j=0;j<i;j++)
1120    {
1121      if (prs[j] !=NULL)
1122      {
1123        if (l>=*l_pairs)
1124        {
1125          SSet temp = (SSet)omAlloc0((*l_pairs+16)*sizeof(SObject));
1126          for (ll=0;ll<*l_pairs;ll++)
1127          {
1128            temp[ll].p = (*resPairs)[ll].p;
1129            temp[ll].p1 = (*resPairs)[ll].p1;
1130            temp[ll].p2 = (*resPairs)[ll].p2;
1131            temp[ll].syz = (*resPairs)[ll].syz;
1132            temp[ll].lcm = (*resPairs)[ll].lcm;
1133            temp[ll].ind1 = (*resPairs)[ll].ind1;
1134            temp[ll].ind2 = (*resPairs)[ll].ind2;
1135            temp[ll].syzind = (*resPairs)[ll].syzind;
1136            temp[ll].order = (*resPairs)[ll].order;
1137            temp[ll].isNotMinimal = (*resPairs)[ll].isNotMinimal;
1138          }
1139          omFreeSize((ADDRESS)(*resPairs),*l_pairs*sizeof(SObject));
1140          *l_pairs += 16;
1141          (*resPairs) = temp;
1142        }
1143        tso.lcm = prs[j];
1144        prs[j] = NULL;
1145        tso.order = p_FDeg(tso.lcm,currRing);
1146        tso.p1 = add_generators->m[j];
1147        tso.p2 = add_generators->m[i];
1148        tso.ind1 = j;
1149        tso.ind2 = i;
1150        tso.syzind = -1;
1151        tso.isNotMinimal = NULL;
1152        tso.p = NULL;
1153        tso.syz = NULL;
1154        SSet rP=*resPairs;
1155#ifdef SHOW_PROT
1156Print("erzeuge Paar im Modul %d,%d mit: \n",index,tso.order);
1157PrintS("poly1: ");pWrite(tso.p1);
1158PrintS("poly2: ");pWrite(tso.p2);
1159PrintS("syz: ");pWrite(tso.syz);
1160PrintS("sPoly: ");pWrite(tso.p);
1161PrintLn();
1162#endif
1163        syEnterPair(rP,&tso,&l,index);
1164        syInitializePair(&tso);
1165      }
1166    }
1167    i++;
1168  }
1169  *first_new = i;
1170  idDelete(&pairs);
1171}
1172
1173/*3
1174* reduction of a single pair in the higher moduls
1175*/
1176#ifdef SHOW_PROT
1177static void redOnePairHIndex(SSet resPairs,int itso, int crit_comp,
1178            syStrategy syzstr,int index,ideal add_generators, ideal add_repr,
1179            ideal new_generators, ideal new_repr,int * next_place_add,int ** g_l,
1180            poly deg_soc)
1181#else
1182static void redOnePairHIndex(SSet resPairs,int itso, int crit_comp,
1183            syStrategy syzstr,int /*index*/,ideal add_generators, ideal add_repr,
1184            ideal new_generators, ideal new_repr,int * next_place_add,int ** g_l,
1185            poly deg_soc)
1186#endif
1187{
1188  SObject tso = resPairs[itso];
1189  assume (tso.lcm!=NULL);
1190  int ng_place=IDELEMS(new_generators);
1191  int i,j;
1192  number n;
1193  poly p;
1194#ifdef EXPERIMENT1
1195  poly syzp;
1196#endif
1197
1198  assume (tso.ind1<*next_place_add);
1199  assume (tso.ind2<*next_place_add);
1200  assume (tso.ind1!=tso.ind2);
1201  assume (tso.p1 == add_generators->m[tso.ind1]);
1202  assume (tso.p2 == add_generators->m[tso.ind2]);
1203  tso.p1 = add_generators->m[tso.ind1];
1204  tso.p2 = add_generators->m[tso.ind2];
1205  if ((tso.p1!=NULL) && (tso.p2!=NULL))
1206  {
1207    if (TEST_OPT_PROT)
1208      PrintS(".");
1209#ifdef USE_PROD_CRIT
1210    if (p_FDeg(tso.p1,currRing)+p_FDeg(tso.p2,currRing)==tso.order+p_FDeg(deg_soc,currRing))
1211    {
1212      if (TEST_OPT_PROT) PrintS("pc");
1213      int ac=pGetComp(tso.p1);
1214      assume(ac=pGetComp(tso.p2));
1215      poly p1=pCopy(tso.p1);
1216      poly p2=pCopy(tso.p2);
1217      poly pp1,pp2,tp1,tp2;
1218      poly sp1=pCopy(add_repr->m[tso.ind1]),sp2=pCopy(add_repr->m[tso.ind2]);
1219      pp1 = p1;
1220      pp2 = p2;
1221      loop
1222      {
1223        assume(pp1!=NULL);
1224        for(i=(int)(currRing->N); i; i--)
1225          pSetExp(pp1,i, pGetExp(pp1,i)- pGetExp(deg_soc,i));
1226        pSetComp(pp1, 0);
1227        pSetm(pp1);
1228        if ((pNext(pp1)!=NULL) && (pGetComp(pNext(pp1))!=ac))  break;
1229        pIter(pp1);
1230      }
1231      loop
1232      {
1233        assume(pp2!=NULL);
1234        for(i=(int)(currRing->N); i; i--)
1235          pSetExp(pp2,i, pGetExp(pp2,i)- pGetExp(deg_soc,i));
1236        pSetComp(pp2, 0);
1237        pSetm(pp2);
1238        if ((pNext(pp2)!=NULL) && (pGetComp(pNext(pp2))!=ac)) break;
1239        pIter(pp2);
1240      }
1241      tp1 = pNext(pp1);
1242      tp2 = pNext(pp2);
1243      pNext(pp1) = NULL;
1244      pNext(pp2) = NULL;
1245      //p_Shift(&p1,-ac,currRing);
1246      //p_Shift(&p2,-ac,currRing);
1247      tp1 = pMult(tp1,pCopy(p2));
1248      tp2 = pMult(tp2,pCopy(p1));
1249      sp1 = pMult(p2,sp1);
1250      sp2 = pMult(p1,sp2);
1251      tso.p = pSub(tp1,tp2);
1252      tso.syz = pSub(sp1,sp2);
1253    }
1254    else
1255#endif
1256    {
1257      tso.p = ksOldCreateSpoly(tso.p2,tso.p1);
1258      number coefgcd = n_Gcd(pGetCoeff(tso.p1),pGetCoeff(tso.p2),currRing->cf);
1259      assume (add_repr->m[tso.ind1]!=NULL);
1260      tso.syz = pCopy(add_repr->m[tso.ind1]);
1261      poly tt = pDivide(tso.lcm,tso.p1);
1262      pSetComp(tt,0);
1263      pSetmComp(tt);
1264      pSetCoeff(tt,nDiv(pGetCoeff(tso.p1),coefgcd));
1265      tso.syz = pMult_mm(tso.syz,tt);
1266      pDelete(&tt);
1267      coefgcd = nInpNeg(coefgcd);
1268      assume (add_repr->m[tso.ind2]!=NULL);
1269      p = pCopy(add_repr->m[tso.ind2]);
1270      tt = pDivide(tso.lcm,tso.p2);
1271      pSetComp(tt,0);
1272      pSetmComp(tt);
1273      pSetCoeff(tt,nDiv(pGetCoeff(tso.p2),coefgcd));
1274      p = pMult_mm(p,tt);
1275      pDelete(&tt);
1276      tso.syz = pAdd(p,tso.syz);
1277      nDelete(&coefgcd);
1278    }
1279#ifdef SHOW_PROT
1280Print("reduziere Paar im Module %d mit: \n",index);
1281PrintS("poly1: ");pWrite(tso.p1);
1282PrintS("poly2: ");pWrite(tso.p2);
1283PrintS("syz: ");pWrite(tso.syz);
1284PrintS("sPoly: ");pWrite(tso.p);
1285#endif
1286    assume(tso.syz!=NULL);
1287    kBucketInit(syzstr->syz_bucket,tso.syz,-1);
1288    if (tso.p!=NULL)
1289    {
1290      kBucketInit(syzstr->bucket,tso.p,-1);
1291      p = kBucketGetLm(syzstr->bucket);
1292      j = 0;
1293      loop
1294      {
1295        if (j>=*next_place_add) break;
1296        if (pDivisibleBy(add_generators->m[j],p))
1297        {
1298          assume (add_repr->m[j]!=NULL);
1299          sySPRedSyz_Kosz(syzstr,add_generators->m[j],add_repr->m[j],p);
1300          n = kBucketPolyRed(syzstr->bucket,add_generators->m[j],
1301                   pLength(add_generators->m[j]), NULL);
1302          p = kBucketGetLm(syzstr->bucket);
1303          if ((p==NULL) || (pGetComp(p)<=crit_comp)) break;
1304          j = 0;
1305        }
1306        else
1307          j++;
1308      }
1309      kBucketClear(syzstr->bucket,&tso.p,&tso.length);
1310    }
1311    kBucketClear(syzstr->syz_bucket,&tso.syz,&j);
1312  }
1313  else
1314  {
1315    PrintS("Shit happens!\n");
1316  }
1317#ifdef SHOW_PROT
1318Print("erhalte Paar im Module %d mit: \n",index);
1319PrintS("syz: ");pWrite(tso.syz);
1320PrintS("sPoly: ");pWrite(tso.p);
1321PrintLn();
1322#endif
1323  if (tso.p!=NULL)
1324  {
1325    if (!nIsOne(pGetCoeff(tso.p)))
1326    {
1327      n=nInvers(pGetCoeff(tso.p));
1328      pNorm(tso.p);
1329      pMult_nn(tso.syz,n);
1330      nDelete(&n);
1331    }
1332  }
1333  if ((TEST_OPT_PROT) && (tso.syz==NULL)) PrintS("null");
1334  if ((tso.p!=NULL) && (pGetComp(tso.p)>crit_comp))
1335  {
1336    if (*next_place_add>=IDELEMS(add_generators))
1337    {
1338      pEnlargeSet(&add_generators->m,IDELEMS(add_generators),16);
1339      pEnlargeSet(&add_repr->m,IDELEMS(add_repr),16);
1340      *g_l = (int*)omRealloc0Size((ADDRESS)*g_l, IDELEMS(add_generators)*sizeof(int),
1341                            (IDELEMS(add_generators)+16)*sizeof(int));
1342      IDELEMS(add_generators) += 16;
1343      IDELEMS(add_repr) += 16;
1344    }
1345    assume(add_repr->m[*next_place_add]==NULL);
1346    add_generators->m[*next_place_add] = tso.p;
1347    add_repr->m[*next_place_add] = tso.syz;
1348    (*g_l)[*next_place_add] = tso.length;
1349    (*next_place_add)++;
1350  }
1351  else
1352  {
1353    while ((ng_place>0) && (new_generators->m[ng_place-1]==NULL) &&
1354          (new_repr->m[ng_place-1]==NULL)) ng_place--;
1355    if (ng_place>=IDELEMS(new_generators))
1356    {
1357      pEnlargeSet(&new_generators->m,IDELEMS(new_generators),16);
1358      IDELEMS(new_generators) += 16;
1359      pEnlargeSet(&new_repr->m,IDELEMS(new_repr),16);
1360      IDELEMS(new_repr) += 16;
1361    }
1362    new_generators->m[ng_place] = tso.p;
1363    new_repr->m[ng_place] = tso.syz;
1364  }
1365  tso.p = NULL;
1366  tso.syz = NULL;
1367  resPairs[itso] = tso;
1368  syDeletePair(&resPairs[itso]);
1369}
1370
1371/*3
1372* reduction of all pairs of a fixed degree of a fixed module
1373*/
1374static BOOLEAN reducePairsHIndex(SSet resPairs,int l_pairs,syStrategy syzstr,
1375       int index,ideal add_generators,ideal add_repr,ideal new_generators,
1376       ideal new_repr,int crit_comp,int * red_deg,int * next_place_add,int **g_l,
1377       resolvente totake)
1378{
1379  if (resPairs[0].lcm==NULL) return FALSE;
1380  int i=0;
1381  poly deg_soc;
1382
1383  if (TEST_OPT_PROT)
1384    Print("(%d,%d)",index,resPairs[0].order);
1385  while ((i<l_pairs) && (resPairs[i].order==*red_deg))
1386  {
1387    assume(totake[index-1]!=NULL);
1388    assume(pGetComp(resPairs[i].p1)<=IDELEMS(totake[index-1]));
1389    assume(totake[index-1]->m[pGetComp(resPairs[i].p1)-1]!=NULL);
1390    deg_soc = totake[index-1]->m[pGetComp(resPairs[i].p1)-1];
1391    redOnePairHIndex(resPairs,i,crit_comp,syzstr,index, add_generators,add_repr,
1392                     new_generators, new_repr,next_place_add,g_l,deg_soc);
1393    i++;
1394  }
1395  syCompactifyPairSet(resPairs,l_pairs,0);
1396  if (resPairs[0].lcm==NULL)  //there are no pairs left and no new_gens
1397    return FALSE;
1398  else
1399    *red_deg = resPairs[0].order;
1400  return TRUE;
1401}
1402
1403/*3
1404* we proceed the generators of the next module;
1405* they are stored in add_generators and add_repr;
1406* if the normal form of a new genrators w.r.t. add_generators has
1407* pGetComp<crit_comp it is skipped from the reduction;
1408* new_generators and new_repr (which are empty) stores the result of the
1409* reduction which is normalized afterwards
1410*/
1411static void procedeNextGenerators(ideal temp_generators,ideal /*temp_repr*/,
1412      ideal new_generators, ideal new_repr, ideal add_generators,
1413      ideal add_repr, syStrategy syzstr,int index, int crit_comp,
1414      resolvente totake)
1415{
1416  int i=0,j,next_new_el;
1417  int idel_temp=IDELEMS(temp_generators);
1418  int next_place_add;
1419  int p_length,red_deg,l_pairs=IDELEMS(add_generators);
1420  poly next_p;
1421  int * gen_length=(int*)omAlloc0(IDELEMS(add_generators)*sizeof(int));
1422  int * secgen_length=(int*)omAlloc0(IDELEMS(syzstr->res[index])*sizeof(int));
1423  BOOLEAN pairs_left;
1424  SSet resPairs=(SSet)omAlloc0(l_pairs*sizeof(SObject));
1425
1426  for (j=IDELEMS(syzstr->res[index])-1;j>=0;j--)
1427  {
1428    if (syzstr->res[index]->m[j]!=NULL)
1429      secgen_length[j] = pLength(syzstr->res[index]->m[j]);
1430  }
1431  assume(idIs0(new_generators));
1432  next_place_add = IDELEMS(add_generators);
1433  while ((next_place_add>0) && (add_generators->m[next_place_add-1]==NULL))
1434    next_place_add--;
1435  int next_deg = p_FDeg(temp_generators->m[i],currRing);
1436  next_new_el = next_place_add;
1437/*--- loop about all all elements-----------------------------------*/
1438  while ((i<idel_temp) && (temp_generators->m[i]!=NULL))
1439  {
1440/*--- separates elements of equal degree----------------------------*/
1441#ifdef USE_REGULARITY
1442    if (syzstr->regularity>0)
1443    {
1444      if (next_deg >= syzstr->regularity+index)
1445      {
1446        while ((i<idel_temp) && (temp_generators->m[i]!=NULL))
1447        {
1448          pDelete(&temp_generators->m[i]);
1449          i++;
1450        }
1451        break;
1452      }
1453    }
1454#endif
1455    while ((i<idel_temp) && (p_FDeg(temp_generators->m[i],currRing)==next_deg))
1456    {
1457      next_p = temp_generators->m[i];
1458      temp_generators->m[i] = NULL;
1459      next_p = normalize(next_p,add_generators,syzstr,gen_length,&p_length,
1460                crit_comp);
1461      if (next_p!=NULL)
1462      {
1463        if (pGetComp(next_p)<=crit_comp)
1464        {
1465          pDelete(&next_p);
1466          //if (TEST_OPT_PROT) Print("u(%d)",index);
1467        }
1468        else
1469        {
1470          next_p = syRedTailSyz(next_p,add_generators,syzstr->res[index],crit_comp,syzstr,
1471            gen_length,secgen_length,&p_length);
1472          if (!nIsOne(pGetCoeff(next_p)))
1473            pNorm(next_p);
1474          if (next_place_add>=IDELEMS(add_generators))
1475          {
1476            pEnlargeSet(&add_generators->m,IDELEMS(add_generators),16);
1477            pEnlargeSet(&add_repr->m,IDELEMS(add_repr),16);
1478            gen_length = (int*)omRealloc0Size((ADDRESS)gen_length, IDELEMS(add_generators)*sizeof(int),
1479                                        (IDELEMS(add_generators)+16)*sizeof(int));
1480            IDELEMS(add_generators) += 16;
1481            IDELEMS(add_repr) += 16;
1482          }
1483          add_generators->m[next_place_add] = next_p;
1484          if (totake[index]==NULL)
1485            totake[index] = idInit(16,new_generators->rank);
1486          if ((*syzstr->Tl)[index]==IDELEMS(totake[index]))
1487          {
1488            pEnlargeSet(&totake[index]->m,IDELEMS(totake[index]),
1489                        (*syzstr->Tl)[index]+16-IDELEMS(totake[index]));
1490            for (j=IDELEMS(totake[index]);j<(*syzstr->Tl)[index]+16;j++)
1491              totake[index]->m[j] = NULL;
1492            IDELEMS(totake[index]) = (*syzstr->Tl)[index]+16;
1493          }
1494#ifdef FULL_TOTAKE
1495          totake[index]->m[(*syzstr->Tl)[index]] = pCopy(next_p);
1496#else
1497          totake[index]->m[(*syzstr->Tl)[index]] = pHead(next_p);
1498#endif
1499          assume(add_repr->m[next_place_add]==NULL);
1500#ifdef WITH_SCHREYER_ORD
1501          add_repr->m[next_place_add] = pHead(add_generators->m[next_place_add]);
1502#else
1503          add_repr->m[next_place_add] = pOne();
1504#endif
1505          ((*syzstr->Tl)[index])++;
1506          pSetComp(add_repr->m[next_place_add],(*syzstr->Tl)[index]);
1507          pSetmComp(add_repr->m[next_place_add]);
1508          gen_length[next_place_add] = p_length;
1509          next_place_add++;
1510        }
1511      }
1512      i++;
1513    }                        //end inner loop
1514    red_deg = next_deg;
1515    if (i<idel_temp)
1516      next_deg = p_FDeg(temp_generators->m[i],currRing);
1517    else
1518      next_deg = -1;
1519    if ((next_place_add>next_new_el) || (next_deg<0))  //there are new generators or pairs
1520    {
1521/*-reducing and generating pairs untill the degree of the next generators-*/
1522      pairs_left = TRUE;
1523      while (pairs_left && ((next_deg<0) || (red_deg<= next_deg)))
1524      {
1525        updatePairsHIndex(&resPairs,&l_pairs,syzstr,index,add_generators,
1526          add_repr,new_generators,new_repr,crit_comp,&next_new_el);
1527        pairs_left = reducePairsHIndex(resPairs,l_pairs,syzstr,index,add_generators,
1528           add_repr,new_generators,new_repr,crit_comp,&red_deg,&next_place_add,&gen_length,
1529           totake);
1530      }
1531    }
1532  }
1533  omFreeSize((SSet)resPairs,l_pairs*sizeof(SObject));
1534  omFreeSize((ADDRESS)gen_length,IDELEMS(add_generators)*sizeof(int));
1535  omFreeSize((ADDRESS)secgen_length,IDELEMS(syzstr->res[index])*sizeof(int));
1536}
1537
1538/*3
1539* normalizes the part of the next reduction lying within the block
1540* of former generators (old_generators);
1541*/
1542static ideal normalizeOldPart(ideal new_generators,ideal new_repr,
1543                      syStrategy syzstr,int index,int /*crit_comp*/)
1544{
1545  ideal old_generators= syzstr->res[index];
1546  ideal old_repr= syzstr->orderedRes[index];
1547  int i,j=0,ii=IDELEMS(old_generators)-1,dummy;
1548  poly p;
1549  number n;
1550  int * g_l=(int*)omAlloc0(IDELEMS(old_generators)*sizeof(int));
1551
1552  for (i=0;i<IDELEMS(old_generators);i++)
1553  {
1554    if (old_generators->m[i]!=NULL)
1555    {
1556      g_l[i] = pLength(old_generators->m[i]);
1557    }
1558  }
1559  for (i=IDELEMS(new_generators)-1;i>=0;i--)
1560  {
1561    if (new_generators->m[i]!=NULL)
1562    {
1563      kBucketInit(syzstr->bucket,new_generators->m[i],
1564                   pLength(new_generators->m[i]));
1565      kBucketInit(syzstr->syz_bucket,new_repr->m[i],
1566                   pLength(new_repr->m[i]));
1567      p = kBucketGetLm(syzstr->bucket);
1568      loop
1569      {
1570        if ((j>=ii) || (p==NULL)) break;
1571        if ((old_generators->m[j]!=NULL) &&
1572            (pDivisibleBy(old_generators->m[j],p)))
1573        {
1574          sySPRedSyz_Kosz(syzstr,old_generators->m[j],old_repr->m[j],p);
1575          n = kBucketPolyRed(syzstr->bucket,old_generators->m[j], g_l[j], NULL);
1576          nDelete(&n);
1577          p = kBucketGetLm(syzstr->bucket);
1578          j = 0;
1579        }
1580        else
1581          j++;
1582      }
1583      assume (p==NULL);
1584      kBucketClear(syzstr->bucket,&new_generators->m[i],&dummy);
1585      kBucketClear(syzstr->syz_bucket,&new_repr->m[i],&dummy);
1586    }
1587  }
1588  ideal result=idInit(IDELEMS(new_repr),new_repr->rank);
1589  for (j=IDELEMS(new_repr)-1;j>=0;j--)
1590  {
1591    result->m[j] = new_repr->m[j];
1592    if ((result->m[j]!=NULL) && (!nIsOne(pGetCoeff(result->m[j]))))
1593      pNorm(result->m[j]);
1594    new_repr->m[j] = NULL;
1595  }
1596  omFreeSize((ADDRESS)g_l,IDELEMS(old_generators)*sizeof(int));
1597  return result;
1598}
1599
1600/*3
1601* constructs the new subresolution for a nonregular extension
1602*/
1603static ideal kosz_ext(ideal new_generators,ideal new_repr,syStrategy syzstr,
1604                      int index,int next_comp,resolvente totake)
1605{
1606  ideal temp_generators =idInit(IDELEMS(new_generators),new_generators->rank);
1607  ideal temp_repr=idInit(IDELEMS(new_repr),new_repr->rank);
1608  ideal add_generators =idInit(IDELEMS(new_generators),new_generators->rank);
1609  ideal add_repr=idInit(IDELEMS(new_repr),new_repr->rank);
1610  int min_deg=-1;
1611  int j,jj,k,deg_p,idel_temp=IDELEMS(temp_generators);
1612  poly p;
1613/*--reorder w.r.t. the degree----------------------------------------*/
1614  for (j=IDELEMS(new_generators)-1;j>=0;j--)
1615  {
1616    if (new_generators->m[j]!=NULL)
1617    {
1618      p = new_generators->m[j];
1619      new_generators->m[j] = NULL;
1620      deg_p = p_FDeg(p,currRing);
1621      if (min_deg<0)
1622      {
1623        min_deg = deg_p;
1624      }
1625      else
1626      {
1627        if (deg_p<min_deg) min_deg = deg_p;
1628      }
1629      k = 0;
1630      while ((k<idel_temp) && (temp_generators->m[k]!=NULL) &&
1631             (p_FDeg(temp_generators->m[k],currRing)<=deg_p)) k++;
1632      for (jj=idel_temp-1;jj>k;jj--)
1633      {
1634        temp_generators->m[jj] = temp_generators->m[jj-1];
1635      }
1636      temp_generators->m[k] = p;
1637    }
1638  }
1639/*--- computing the standard basis in the resolution of the extension -*/
1640  procedeNextGenerators(temp_generators,temp_repr,new_generators,new_repr,
1641    add_generators,add_repr,syzstr,index,next_comp,totake);
1642  j = IDELEMS(syzstr->res[index]);
1643  while ((j>0) && (syzstr->res[index]->m[j-1]==NULL)) j--;
1644  jj = IDELEMS(add_generators);
1645  while ((jj>0) && (add_generators->m[jj-1]==NULL)) jj--;
1646  if (j+jj>=IDELEMS(syzstr->res[index]))
1647  {
1648    pEnlargeSet(&syzstr->res[index]->m,IDELEMS(syzstr->res[index]),
1649                j+jj+1-IDELEMS(syzstr->res[index]));
1650    IDELEMS(syzstr->res[index]) = j+jj+1;
1651    pEnlargeSet(&syzstr->orderedRes[index]->m,IDELEMS(syzstr->orderedRes[index]),
1652                j+jj+1-IDELEMS(syzstr->orderedRes[index]));
1653    IDELEMS(syzstr->orderedRes[index]) = j+jj+1;
1654  }
1655  for (k=0;k<jj;k++)
1656  {
1657    syzstr->res[index]->m[j+k] = add_generators->m[k];
1658    syzstr->orderedRes[index]->m[j+k] = add_repr->m[k];
1659    add_generators->m[k] = NULL;
1660    add_repr->m[k] = NULL;
1661  }
1662  assume(idIs0(add_generators));
1663  assume(idIs0(add_repr));
1664  idDelete(&add_generators);
1665  idDelete(&add_repr);
1666  idDelete(&temp_generators);
1667  idDelete(&temp_repr);
1668/*--- normalizing the rest to get the syzygies ------------------------*/
1669  return normalizeOldPart(new_generators,new_repr,syzstr,index,next_comp);
1670}
1671
1672/*
1673* this procedure assumes that the first order is C !!!
1674* INPUT: old_generators - the generators of the actual module
1675*                         computed so far (they are mixed vectors)
1676*        old_repr       - the representations of the old generators
1677*        new_generators - generators coming from reductions below,
1678*                         they must have leading terms in new components
1679*                         (they live only in the module part)
1680*  (*syzstr->Tl)[index] - the last used component in the syzygy
1681* OUTPUT: old_generators is updated
1682*         new_generators is empty
1683*         the return value is a set of new generators for the syzygies,
1684*/
1685static ideal syAppendSyz(ideal new_generators, syStrategy syzstr,int index,int crit_comp,
1686                         resolvente totake)
1687{
1688  int i,j;
1689  ideal result;
1690  int rk_new_gens = id_RankFreeModule(new_generators,currRing);
1691  if (syzstr->res[index]==NULL)
1692  {
1693    syzstr->res[index] = idInit(1,si_max(rk_new_gens,1));
1694    syzstr->orderedRes[index] = idInit(1,si_max(rk_new_gens,1));
1695  }
1696  int ng_idel=IDELEMS(new_generators);
1697  ideal new_repr =idInit(ng_idel, crit_comp+ng_idel);
1698
1699  if (index==0)
1700  {
1701    //int * og_l=(int*)omAlloc0(IDELEMS(syzstr->res[0])*sizeof(int));
1702    //for (i=IDELEMS(syzstr->res[0])-1;i>=0;i--)
1703    //{
1704      //if (syzstr->res[0]->m[i]!=NULL)
1705        //og_l[i] = pLength(syzstr->res[0]->m[i]);
1706    //}
1707    for (i=0;i<ng_idel;i++)
1708    {
1709      if (new_generators->m[i]!=NULL)
1710      {
1711        //int ng_l=pLength(new_generators->m[i]);
1712        //new_generators->m[i] = syRedTailSyz(new_generators->m[i],syzstr->res[0],NULL,0,syzstr,
1713            //og_l,NULL,&ng_l);
1714        if (totake[index]==NULL)
1715          totake[index] = idInit(16,new_generators->rank);
1716        if ((*syzstr->Tl)[index]>=IDELEMS(totake[index]))
1717        {
1718          pEnlargeSet(&totake[index]->m,IDELEMS(totake[index]),
1719                      (*syzstr->Tl)[index]+16-IDELEMS(totake[index]));
1720          for (j=IDELEMS(totake[index]);j<(*syzstr->Tl)[index]+16;j++)
1721            totake[index]->m[j] = NULL;
1722          IDELEMS(totake[index]) = (*syzstr->Tl)[index]+16;
1723        }
1724#ifdef FULL_TOTAKE
1725        totake[index]->m[(*syzstr->Tl)[index]] = pCopy(new_generators->m[i]);
1726#else
1727        totake[index]->m[(*syzstr->Tl)[index]] = pHead(new_generators->m[i]);
1728#endif
1729#ifdef WITH_SCHREYER_ORD
1730        new_repr->m[i] = pHead(new_generators->m[i]);
1731#else
1732        new_repr->m[i] = pOne();
1733#endif
1734        ((*syzstr->Tl)[index])++;
1735        pSetComp(new_repr->m[i],(*syzstr->Tl)[index]);
1736        pSetmComp(new_repr->m[i]);
1737      }
1738    }
1739    //omFreeSize((ADDRESS)og_l,IDELEMS(syzstr->res[0])*sizeof(int));
1740#ifdef SHOW_PROT
1741PrintS("Add new generators:\n");
1742idPrint(new_generators);
1743PrintS("with representaions:\n");
1744idPrint(new_repr);
1745#endif
1746    result = kosz_std(new_generators,new_repr,syzstr,index,crit_comp);
1747  }
1748  else
1749  {
1750    result = kosz_ext(new_generators,new_repr,syzstr,index,crit_comp,totake);
1751  }
1752  idSkipZeroes(result);
1753  assume(idIs0(new_repr));
1754  idDelete(&new_repr);
1755  return result;
1756}
1757
1758/*
1759* main call of the extended Koszul-resolution
1760*/
1761syStrategy syKosz(ideal arg,int * length)
1762{
1763  int i,j,jj,k=0,index=0,rk_arg/*,next_syz=0*/;
1764  int crit_comp,t_comp,next_deg,old_tl;
1765  ideal temp=NULL,old_ideal,old_repr;
1766  ring origR = currRing;
1767  poly next_gen;
1768  BOOLEAN isRegular;
1769
1770  discard_pairs = 0;
1771  short_pairs = 0;
1772  if (idIs0(arg)) return NULL;
1773  rk_arg = id_RankFreeModule(arg,currRing);
1774  syStrategy syzstr=(syStrategy)omAlloc0(sizeof(ssyStrategy));
1775/*--- changes to a Cdp-ring ----------------------------*/
1776  syzstr->syRing = rAssure_C_dp(origR); rChangeCurrRing(syzstr->syRing);
1777/*--- initializes the data structures---------------*/
1778  syzstr->length = *length = (syzstr->syRing->N)+2;
1779  syzstr->regularity = -1;
1780  if (origR!=syzstr->syRing)
1781    temp = idrCopyR(arg, origR, syzstr->syRing);
1782  else
1783    temp = idCopy(arg);
1784  if (rk_arg==0)
1785  {
1786    id_Shift(temp,1,currRing);
1787  }
1788  idSkipZeroes(temp);
1789#ifdef WITH_SORT
1790  if (temp->m[0]!=NULL)
1791  {
1792    int md;
1793    int maxdeg=p_FDeg(temp->m[IDELEMS(temp)-1],currRing);
1794    ideal temp1=idInit(IDELEMS(temp),temp->rank);
1795    for (j=IDELEMS(temp)-2;j>=0;j--)
1796    {
1797      jj = p_FDeg(temp->m[j],currRing);
1798      if (jj>maxdeg) maxdeg = jj;
1799    }
1800    while (!idIs0(temp))
1801    {
1802      md = maxdeg;
1803      for (j=IDELEMS(temp)-1;j>=0;j--)
1804      {
1805        if (temp->m[j]!=NULL)
1806        {
1807          jj = p_FDeg(temp->m[j],currRing);
1808          if (jj<md) md = jj;
1809        }
1810      }
1811      for (j=0;j<IDELEMS(temp);j++)
1812      {
1813        if ((temp->m[j]!=NULL) && (p_FDeg(temp->m[j],currRing)==md))
1814        {
1815          temp1->m[k] = temp->m[j];
1816          temp->m[j] = NULL;
1817          k++;
1818        }
1819      }
1820    }
1821    idDelete(&temp);
1822    temp = temp1;
1823    temp1 = NULL;
1824  }
1825#endif
1826#ifdef USE_REGULARITY
1827  int last_generator=IDELEMS(temp)-1;
1828  while ((last_generator>=0) && (temp->m[last_generator]==NULL))
1829    last_generator--;
1830#endif
1831  syzstr->res = (resolvente)omAlloc0((*length+1)*sizeof(ideal));
1832  syzstr->orderedRes = (resolvente)omAlloc0((*length+1)*sizeof(ideal));
1833  resolvente totake=(resolvente)omAlloc0((*length+1)*sizeof(ideal));
1834  syzstr->Tl = new intvec(*length+1);
1835  syzstr->bucket = kBucketCreate(currRing);
1836  syzstr->syz_bucket = kBucketCreate(currRing);
1837  ideal new_generators=idInit(1,si_max(rk_arg,1));
1838  ideal temp_gens,old_std;
1839  syzstr->res[0] = idInit(1,1);
1840  if (rk_arg>1) syzstr->res[0]->rank = rk_arg;
1841  syzstr->orderedRes[0] = idInit(1,1);
1842/*--- computes the resolution ----------------------*/
1843  i = 0;
1844  while (i<IDELEMS(temp))
1845  {
1846    if (temp->m[i]!=NULL)
1847    {
1848      new_generators->m[0] = kNF(syzstr->res[0],currRing->qideal,temp->m[i]);
1849      if (!nIsOne(pGetCoeff(new_generators->m[0])))
1850        pNorm(new_generators->m[0]);
1851      next_deg = p_FDeg(new_generators->m[0],currRing);
1852      next_gen = pCopy(new_generators->m[0]);
1853    }
1854    if (!idIs0(new_generators))
1855    {
1856      index = 0;
1857      while (index<=*length)
1858      {
1859        if (index==0)
1860        {
1861          old_ideal = idCopy(syzstr->res[0]);
1862          old_repr = idCopy(syzstr->orderedRes[0]);
1863          old_tl = (*syzstr->Tl)[0];
1864          old_std = id_Head(syzstr->res[0],currRing);
1865        }
1866        t_comp = (*syzstr->Tl)[index];
1867        if (index==0) crit_comp = t_comp;
1868        temp_gens = syAppendSyz(new_generators,syzstr, index,crit_comp,totake);
1869        crit_comp = t_comp;
1870        if (index==0)
1871        {
1872          isRegular = syIsRegular(old_std,syzstr->res[0],next_deg);
1873#ifndef ONLY_STD
1874          if (isRegular)
1875            syCreateRegularExtension(syzstr,old_ideal,old_repr,old_tl,next_gen,
1876                                     totake);
1877#ifdef USE_REGULARITY
1878        if ((index==0) && (!isRegular) && (i==last_generator))
1879        {
1880/*----------- we are computing the regularity -----------------------*/
1881          ideal initial=id_Head(syzstr->res[0],currRing);
1882          int len=0,reg=0;
1883          intvec *w=NULL;
1884          ring dp_C_ring = rAssure_dp_C(currRing); rChangeCurrRing(dp_C_ring);
1885          initial = idrMoveR_NoSort(initial, syzstr->syRing, dp_C_ring);
1886          resolvente res = sySchreyerResolvente(initial,-1,&len,TRUE, TRUE);
1887          intvec * dummy = syBetti(res,len,&reg, w);
1888          syzstr->regularity = reg+2;
1889          delete dummy;
1890          delete w;
1891          for (j=0;j<len;j++)
1892          {
1893            if (res[j]!=NULL) idDelete(&(res[j]));
1894          }
1895          omFreeSize((ADDRESS)res,len*sizeof(ideal));
1896          idDelete(&initial);
1897          rChangeCurrRing(syzstr->syRing);
1898          rDelete(dp_C_ring);
1899        }
1900#endif
1901#endif
1902          idDelete(&old_ideal);
1903          idDelete(&old_repr);
1904          idDelete(&old_std);
1905          if (TEST_OPT_PROT)
1906          {
1907            if (isRegular)
1908              PrintS("\n regular\n");
1909            else
1910              PrintS("\n not regular\n");
1911          }
1912          if (next_gen!=NULL)
1913            pDelete(&next_gen);
1914          if (isRegular)
1915          {
1916            idDelete(&temp_gens);
1917            break;
1918          }
1919        }
1920        idDelete(&new_generators);
1921        new_generators = temp_gens;
1922#ifdef ONLY_STD
1923        break;
1924#endif
1925        if (idIs0(new_generators)) break;
1926        index++;
1927      }
1928      if (!idIs0(new_generators))
1929      {
1930        for (j=0;j<IDELEMS(new_generators);j++)
1931        {
1932          if (new_generators->m[j]!=NULL)
1933          {
1934            pDelete(&new_generators->m[j]);
1935            new_generators->m[j] = NULL;
1936          }
1937        }
1938      }
1939    }
1940    i++;
1941  }
1942  if (idIs0(new_generators) && new_generators!=NULL) idDelete(&new_generators);
1943  if (temp!=NULL) idDelete(&temp);
1944  kBucketDestroy(&(syzstr->bucket));
1945  kBucketDestroy(&(syzstr->syz_bucket));
1946  index = 0;
1947  syzstr->fullres = syzstr->res;
1948  syzstr->res = NULL;
1949  index = 0;
1950  while ((index<=*length) && (syzstr->fullres[index]!=NULL))
1951  {
1952#ifdef SHOW_RESULT
1953    Print("The %d-th syzygy-module is now:\n",index);
1954    ideal ttt=id_Head(syzstr->fullres[index],currRing);
1955    idShow(ttt);
1956    idDelete(&ttt);
1957    //if (index>0)
1958    //{
1959      //Print("The related module is: \n");
1960      //idPrint(totake[index-1]);
1961    //}
1962    //Print("The %d-th module of the minimal resolution is:\n",index);
1963    if (!idIs0(totake[index]))
1964      idShow(totake[index]);
1965    //Print("with standard basis:\n");
1966    //idPrint(syzstr->fullres[index]);
1967    //if ((index<*length) && (totake[index+1]!=NULL))
1968    //{
1969      //Print("The %d-th syzygy-module is now:\n",index+1);
1970      //idPrint(totake[index+1]);
1971      //matrix m1=idModule2Matrix(totake[index]);
1972      //matrix m2=idModule2Matrix(totake[index+1]);
1973      //matrix m3=mpMult(m1,m2);
1974      //idPrint((ideal)m3);
1975    //}
1976#endif
1977    if (!idIs0(totake[index]))
1978    {
1979      for(i=0;i<IDELEMS(totake[index]);i++)
1980      {
1981        if (totake[index]->m[i]!=NULL)
1982        {
1983          j=0;
1984          while ((j<IDELEMS(syzstr->fullres[index])) &&
1985            ((syzstr->fullres[index]->m[j]==NULL) ||
1986            (!pLmEqual(syzstr->fullres[index]->m[j],totake[index]->m[i])))) j++;
1987          if (j<IDELEMS(syzstr->fullres[index]))
1988          {
1989            pDelete(&totake[index]->m[i]);
1990            totake[index]->m[i] = syzstr->fullres[index]->m[j];
1991            syzstr->fullres[index]->m[j] = NULL;
1992          }
1993          else
1994          {
1995            PrintS("Da ist was faul!!!\n");
1996            Print("Aber: Regularitaet %d, Grad %ld\n",
1997                   syzstr->regularity,p_FDeg(totake[index]->m[i],currRing));
1998          }
1999        }
2000      }
2001      idDelete(&syzstr->fullres[index]);
2002      syzstr->fullres[index] = totake[index];
2003    }
2004#ifdef SHOW_RESULT
2005    idShow(syzstr->fullres[index]);
2006#endif
2007    index++;
2008  }
2009  syReorder_Kosz(syzstr);
2010  index = 0;
2011  while ((index<=*length) && (syzstr->orderedRes[index]!=NULL))
2012  {
2013    idDelete(&(syzstr->orderedRes[index]));
2014    index++;
2015  }
2016  if (origR!=syzstr->syRing)
2017  {
2018    rChangeCurrRing(origR);
2019    index = 0;
2020    while ((index<=*length) && (syzstr->fullres[index]!=NULL))
2021    {
2022      syzstr->fullres[index] = idrMoveR(syzstr->fullres[index],syzstr->syRing, origR);
2023      index++;
2024    }
2025  }
2026  delete syzstr->Tl;
2027  syzstr->Tl = NULL;
2028  rDelete(syzstr->syRing);
2029  syzstr->syRing = NULL;
2030  omFreeSize((ADDRESS)totake,(*length+1)*sizeof(ideal));
2031  omFreeSize((ADDRESS)syzstr->orderedRes,(*length+1)*sizeof(ideal));
2032//Print("Pairs to discard: %d\n",discard_pairs);
2033//Print("Pairs shorter reduced: %d\n",short_pairs);
2034//discard_pairs = 0;
2035//short_pairs = 0;
2036  return syzstr;
2037}
2038
Note: See TracBrowser for help on using the repository browser.