source: git/Singular/tgb.cc @ 7724a3

spielwiese
Last change on this file since 7724a3 was 7724a3, checked in by Michael Brickenstein <bricken@…>, 21 years ago
*bricken: algorithmic changes git-svn-id: file:///usr/local/Singular/svn/trunk@6289 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 28.5 KB
Line 
1#include "mod2.h"
2#include <omalloc.h>
3#include "p_polys.h"
4#include "prCopy.h"
5#include "ideals.h"
6#include "ring.h"
7#include "febase.h"
8#include "structs.h"
9#include "polys.h"
10#include "stdlib.h"
11
12
13#include "kutil.h"
14#include "kInline.cc"
15#include "kstd1.h"
16#include "kbuckets.h"
17
18#define FULLREDUCTIONS
19#define HALFREDUCTIONS
20#define HEAD_BIN
21//#define HOMOGENEOUS_EXAMPLE
22#define REDTAIL_S
23//#define REDTAIL_PROT
24//#define QUICK_SPOLY_TEST
25//#define DIAGONAL_GOING
26//#define RANDOM_WALK
27struct int_pair_node{
28  int_pair_node* next;
29  int a;
30  int b;
31};
32#ifdef RANDOM_WALK
33int my_rand(int n){
34  //RANDOM integer beetween 0,..,n-1
35  int erg=(double(rand())/RAND_MAX) *n;
36  return ((erg>=n)?erg-1:erg);
37
38}
39#endif
40enum calc_state
41  {
42    UNCALCULATED,
43    HASTREP,
44    UNIMPORTANT
45  };
46struct calc_dat
47{
48  int* rep;
49  char** states;
50  ideal S;
51  ring r;
52  int* lengths;
53  long* short_Exps;
54  kStrategy strat;
55  int** deg;
56  int* T_deg;
57  int* misses;
58  int_pair_node* soon_free;
59#ifdef HEAD_BIN
60  struct omBin_s*   HeadBin;
61#endif
62  int max_misses;
63  int found_i;
64  int found_j;
65  int continue_i;
66  int continue_j;
67  int n;
68  int skipped_i;
69  int normal_forms;
70  int skipped_pairs;
71  int current_degree;
72  int misses_counter;
73  int misses_series;
74};
75bool find_next_pair(calc_dat* c);
76void shorten_tails(calc_dat* c, poly monom);
77void replace_pair(int & i, int & j, calc_dat* c);
78void initial_data(calc_dat* c);
79void add_to_basis(poly h, calc_dat* c);
80void do_this_spoly_stuff(int i,int j,calc_dat* c);
81ideal t_rep_gb(ring r,ideal arg_I);
82bool has_t_rep(const int & arg_i, const int & arg_j, calc_dat* state);
83int* make_connections(int from, poly bound, calc_dat* c);
84int* make_connections(int from, int to, poly bound, calc_dat* c);
85void now_t_rep(int arg_i, int arg_j, calc_dat* c);
86int pLcmDeg(poly a, poly b);
87int simple_posInS (kStrategy strat, poly p,int len);
88
89bool find_next_pair(calc_dat* c)
90{
91  int start_i,start_j,i,j;
92  start_i=0;
93  start_j=-1;
94#ifndef HOMOGENEOUS_EXAMPLE
95  if (c->misses_series>80){
96    c->current_degree++;
97    c->misses_series=0;
98  }
99#endif
100  start_i=c->continue_i;
101  start_j=c->continue_j;
102#ifndef DIAGONAL_GOING
103  for (int i=start_i;i<c->n;i++){
104    if (c->T_deg[i]>c->current_degree)
105      {
106        c->skipped_pairs++;
107        continue;
108      }
109    for(int j=(i==start_i)?start_j+1:0;j<i;j++){
110      // printf("searching at %d,%d",i,j);
111      if (c->misses_counter>=2) {
112        c->skipped_pairs++;
113        break;
114      }
115      if (c->states[i][j]==UNCALCULATED){
116        if(c->deg[i][j]<=c->current_degree)
117          {
118            c->continue_i=c->found_i=i;
119            c->continue_j=c->found_j=j;
120            return true;
121          }
122        else
123          {
124            ++(c->skipped_pairs);
125          }
126      }
127    }
128    c->misses_counter=0;
129  }
130
131
132  if (!((start_i==1) &&(start_j==0))){
133    c->continue_i=1;
134    c->continue_j=0;
135    return find_next_pair(c);
136  }
137
138#else
139  int z=0;
140  i=start_i;
141  j=start_j;
142  z=start_i+start_j;
143  while(z<=(2*c->n-3)){
144
145    while(i>j){
146      if(i>=c->n) {
147        if (c->skipped_i<0) c->skipped_i=i;
148      } else{
149        if(c->states[i][j]==UNCALCULATED){
150          if (c->deg[i][j]<=c->current_degree)
151            {
152              c->continue_i=c->found_i=i;
153              c->continue_j=c->found_j=j;
154              return true;
155            }
156          else
157            {
158              ++(c->skipped_pairs);
159            }
160        }
161      }
162      --i;
163      ++j;
164    }
165    ++z;
166    i=z;
167    j=0;
168  }
169#endif
170  if (c->skipped_pairs>0){
171    ++(c->current_degree);
172    c->skipped_pairs=0;
173    c->continue_i=0;
174    c->continue_j=0;
175    return find_next_pair(c);
176  }
177  return false;
178}
179void move_forward_in_S(int old_pos, int new_pos,kStrategy strat)
180{
181  poly p=strat->S[old_pos];
182  int ecart=strat->ecartS[old_pos];
183  long sev=strat->sevS[old_pos];
184  int s_2_r=strat->S_2_R[old_pos];
185  int length=strat->lenS[old_pos];
186  int i;
187  for (i=old_pos; i>new_pos; i--)
188    {
189      strat->S[i] = strat->S[i-1];
190      strat->ecartS[i] = strat->ecartS[i-1];
191      strat->sevS[i] = strat->sevS[i-1];
192      strat->S_2_R[i] = strat->S_2_R[i-1];
193    }
194  if (strat->lenS!=NULL)
195    for (i=old_pos; i>new_pos; i--)
196      strat->lenS[i] = strat->lenS[i-1];
197
198  strat->S[new_pos]=p;
199  strat->ecartS[new_pos]=ecart;
200  strat->sevS[new_pos]=sev;
201  strat->S_2_R[new_pos]=s_2_r;
202  strat->lenS[new_pos]=length;
203}
204void replace_pair(int & i, int & j, calc_dat* c)
205{
206  c->soon_free=NULL;
207  int curr_deg;
208  poly lm=pOne();
209
210  pLcm(c->S->m[i], c->S->m[j], lm);
211  pSetm(lm);
212  int deciding_deg= pFDeg(lm);
213  int* i_con =make_connections(i,j,lm,c);
214  int z=0;
215
216  for (int n=0;((n<c->n) && (i_con[n]>=0));n++){
217    if (i_con[n]==j){
218      //       curr_deg=pFDeg(lm);
219      //       for(int z1=0;((z1<c->n) && (i_con[z1]>=0));z1++)
220      //         for (int z2=z1+1;((z2<c->n)&&(i_con[z2]>=0));z2++)
221      //         {
222      //           pLcm(c->S->m[i_con[z1]], c->S->m[i_con[z2]], lm);
223      //           pSetm(lm);
224      //           if (pFDeg(lm)==curr_deg)
225      //             now_t_rep(i_con[z1],i_con[z2],c);
226      //         }
227      now_t_rep(i,j,c);
228      omfree(i_con);
229      p_Delete(&lm,c->r);
230      return;
231    }
232  }
233
234  int* j_con =make_connections(j,lm,c);
235  i= i_con[0];
236  j=j_con[0];
237  if(c->n>1){
238    if (i_con[1]>=0)
239      i=i_con[1];
240    else {
241      if (j_con[1]>=0)
242        j=j_con[1];
243    }
244  }
245  pLcm(c->S->m[i], c->S->m[j], lm);
246  pSetm(lm);
247  poly short_s;
248  curr_deg=pFDeg(lm);
249  int_pair_node* last=NULL;
250
251  for (int n=0;((n<c->n) && (j_con[n]>=0));n++){
252    for (int m=0;((m<c->n) && (i_con[m]>=0));m++){
253      pLcm(c->S->m[i_con[m]], c->S->m[j_con[n]], lm);
254      pSetm(lm);
255      if (pFDeg(lm)>=deciding_deg)
256        {
257          int_pair_node* h= (int_pair_node*)omalloc(sizeof(int_pair_node));
258          if (last!=NULL)
259            last->next=h;
260          else
261            c->soon_free=h;
262          h->next=NULL;
263          h->a=i_con[m];
264          h->b=j_con[n];
265          last=h;
266        }
267      //      if ((comp_deg<curr_deg)
268      //  ||
269      //  ((comp_deg==curr_deg) &&
270      short_s=ksCreateShortSpoly(c->S->m[i_con[m]],c->S->m[j_con[n]],c->r);
271      if (short_s==NULL) {
272        i=i_con[m];
273        j=j_con[n];
274        now_t_rep(i_con[m],j_con[n],c);
275        p_Delete(&lm,c->r);
276        omfree(i_con);
277        omfree(j_con);
278
279        return;
280
281      }
282#ifdef QUICK_SPOLY_TEST
283
284      for (int dz=0;dz<=c->n;dz++){
285        if (dz==c->n) {
286          //have found not head reducing pair
287          i=i_con[m];
288          j=j_con[n];
289          p_Delete(&short_s,c->r);
290          p_Delete(&lm,c->r);
291          omfree(i_con);
292          omfree(j_con);
293
294          return;
295        }
296        if (p_LmDivisibleBy(c->S->m[dz],short_s,c->r)) break;
297      }
298#endif
299      int comp_deg(pFDeg(short_s));
300      p_Delete(&short_s,c->r);
301      if ((comp_deg<curr_deg)
302          ||
303          ((comp_deg==curr_deg) &&
304           (c->misses[i]+c->misses[j]
305            <=
306            c->misses[i_con[m]]+c->misses[j_con[n]])))
307        //       if ((comp_deg<curr_deg)
308        //           ||
309        //           ((comp_deg==curr_deg) &&
310        //            (c->lengths[i]+c->lengths[j]
311        //             <=
312        //             c->lengths[i_con[m]]+c->lengths[j_con[n]])))
313
314        {
315          curr_deg=comp_deg;
316          i=i_con[m];
317          j=j_con[n];
318        }
319    }
320  }
321  p_Delete(&lm,c->r);
322  omfree(i_con);
323  omfree(j_con);
324  return;
325}
326int* make_connections(int from, poly bound, calc_dat* c)
327{
328  ideal I=c->S;
329  int s=pFDeg(bound);
330  int* cans=(int*) omalloc(c->n*sizeof(int));
331  int* connected=(int*) omalloc(c->n*sizeof(int));
332  int cans_length=0;
333  connected[0]=from;
334  int connected_length=1;
335  long neg_bounds_short= ~p_GetShortExpVector(bound,c->r);
336  for (int i=0;i<c->n;i++){
337    if (c->T_deg[i]>s) continue;
338    if (i!=from){
339      if(p_LmShortDivisibleBy(I->m[i],c->short_Exps[i],bound,neg_bounds_short,c->r)){
340        cans[cans_length]=i;
341        cans_length++;
342      }
343    }
344  }
345  int not_yet_found=cans_length;
346  int con_checked=0;
347  int pos;
348  while((not_yet_found>0) && (con_checked<connected_length)){
349    pos=connected[con_checked];
350    for(int i=0;i<cans_length;i++){
351      if (cans[i]<0) continue;
352      if (has_t_rep(pos,cans[i],c)){
353
354        connected[connected_length]=cans[i];
355        connected_length++;
356        cans[i]=-1;
357        --not_yet_found;
358
359      }
360    }
361    con_checked++;
362  }
363  if (connected_length<c->n){
364    connected[connected_length]=-1;
365  }
366  omfree(cans);
367  return connected;
368}
369int* make_connections(int from, int to, poly bound, calc_dat* c)
370{
371  ideal I=c->S;
372  int s=pFDeg(bound);
373  int* cans=(int*) omalloc(c->n*sizeof(int));
374  int* connected=(int*) omalloc(c->n*sizeof(int));
375  cans[0]=to;
376  int cans_length=1;
377  connected[0]=from;
378  int last_cans_pos=-1;
379  int connected_length=1;
380  long neg_bounds_short= ~p_GetShortExpVector(bound,c->r);
381  // for (int i=0;i<c->n;i++){
382  //     if (c->T_deg[i]>s) continue;
383  //     if (i!=from){
384  //       if(p_LmShortDivisibleBy(I->m[i],c->short_Exps[i],bound,neg_bounds_short,c->r)){
385  //         cans[cans_length]=i;
386  //         cans_length++;
387  //       }
388  //     }
389  //   }
390  int not_yet_found=cans_length;
391  int con_checked=0;
392  int pos;
393  bool can_find_more=true;
394  while(((not_yet_found>0) && (con_checked<connected_length))||can_find_more){
395    if ((con_checked<connected_length)&& (not_yet_found>0)){
396      pos=connected[con_checked];
397      for(int i=0;i<cans_length;i++){
398        if (cans[i]<0) continue;
399        if (has_t_rep(pos,cans[i],c)){
400
401          connected[connected_length]=cans[i];
402          connected_length++;
403          cans[i]=-1;
404          --not_yet_found;
405
406          if (connected[connected_length-1]==to){
407            if (connected_length<c->n){
408              connected[connected_length]=-1;
409            }
410
411            omfree(cans);
412            return connected;
413          }
414        }
415      }
416      con_checked++;
417    } else {
418
419      for(last_cans_pos++;last_cans_pos<=c->n;last_cans_pos++){
420        if (last_cans_pos==c->n){
421          if (connected_length<c->n){
422            connected[connected_length]=-1;
423          }
424          omfree(cans);
425          return connected;
426        }
427        if ((last_cans_pos==from)||(last_cans_pos==to))
428          continue;
429        if(p_LmShortDivisibleBy(I->m[last_cans_pos],c->short_Exps[last_cans_pos],bound,neg_bounds_short,c->r)){
430          cans[cans_length]=last_cans_pos;
431          cans_length++;
432          break;
433        }
434      }
435      not_yet_found++;
436      for (int i=0;i<con_checked;i++){
437        if (has_t_rep(connected[i],last_cans_pos,c)){
438
439          connected[connected_length]=last_cans_pos;
440          connected_length++;
441          cans[cans_length-1]=-1;
442
443          --not_yet_found;
444          if (connected[connected_length-1]==to){
445            if (connected_length<c->n){
446              connected[connected_length]=-1;
447            }
448
449            omfree(cans);
450            return connected;
451          }
452          break;
453        }
454
455      }
456    }
457  }
458  if (connected_length<c->n){
459    connected[connected_length]=-1;
460  }
461
462  omfree(cans);
463  return connected;
464}
465#ifdef HEAD_BIN
466static inline poly p_MoveHead(poly p, omBin b)
467{
468  poly np;
469  omTypeAllocBin(poly, np, b);
470  memmove(np, p, b->sizeW*sizeof(long));
471  omFreeBinAddr(p);
472  return np;
473}
474#endif
475void initial_data(calc_dat* c){
476  void* h;
477  int i,j;
478  c->misses_series=0;
479  c->soon_free=NULL;
480  c->misses_counter=0;
481  c->max_misses=0;
482  c->normal_forms=0;
483  c->current_degree=1;
484  c->skipped_pairs=0;
485  int n=c->S->idelems();
486  c->n=n;
487  c->T_deg=(int*) omalloc(n*sizeof(int));
488  c->continue_i=0;
489  c->continue_j=0;
490  c->skipped_i=-1;
491#ifdef HEAD_BIN
492  c->HeadBin=omGetSpecBin(POLYSIZE + (currRing->ExpL_Size)*sizeof(long));
493#endif
494  /* omUnGetSpecBin(&(c->HeadBin)); */
495  h=omalloc(n*sizeof(char*));
496  if (h!=NULL)
497    c->states=(char**) h;
498  else
499    exit(1);
500  c->misses=(int*) omalloc(n*sizeof(int));
501  c->deg=(int **) omalloc(n*sizeof(int*));
502  h=omalloc(n*sizeof(int));
503  if (h!=NULL)
504    c->lengths=(int*) h;
505  else
506    exit(1);
507
508  h=omalloc(n*sizeof(int));
509  if (h!=NULL)
510    c->rep=(int*) h;
511  else
512    exit(1);
513  c->short_Exps=(long*) omalloc(n*sizeof(long));
514  for (i=0;i<n;i++)
515    {
516#ifdef HEAD_BIN
517      c->S->m[i]=p_MoveHead(c->S->m[i],c->HeadBin);
518#endif
519      c->T_deg[i]=pFDeg(c->S->m[i]);
520      c->lengths[i]=pLength(c->S->m[i]);
521      c->misses[i]=0;
522      c->states[i]=(char *)omAlloc((i+1)*sizeof(char));
523      c->deg[i]=(int*) omAlloc((i+1)*sizeof(int));
524      for (j=0;j<i;j++){
525        //check product criterion
526        if (pHasNotCF(c->S->m[i],c->S->m[j])){
527          c->states[i][j]=HASTREP;
528        } else {
529          c->states[i][j]=UNCALCULATED;
530        }
531        c->deg[i][j]=pLcmDeg(c->S->m[i],c->S->m[j]);
532      }
533      if ((c->lengths[i]==1) && (c->lengths[j]==1))
534        c->states[i][j]=HASTREP;
535      c->rep[i]=i;
536      c->short_Exps[i]=p_GetShortExpVector(c->S->m[i],c->r);
537
538    }
539
540
541  c->strat=new skStrategy;
542  c->strat->syzComp = 0;
543  initBuchMoraCrit(c->strat);
544  initBuchMoraPos(c->strat);
545  c->strat->initEcart = initEcartBBA;
546  c->strat->enterS = enterSBba;
547  c->strat->sl = -1;
548  /* initS(c->S,NULL,c->strat); */
549/* intS start: */
550 i=((i+IDELEMS(c->S)+15)/16)*16;
551 c->strat->ecartS=(intset)omAlloc(i*sizeof(int)); /*initec(i);*/
552 c->strat->sevS=(unsigned long*)omAlloc0(i*sizeof(unsigned long));
553 /*initsevS(i);*/
554 c->strat->S_2_R=(int*)omAlloc0(i*sizeof(int));/*initS_2_R(i);*/
555 c->strat->fromQ=NULL;
556 c->strat->Shdl=idInit(i,1);
557 c->strat->S=c->strat->Shdl->m;
558 c->strat->lenS=(int*)omAlloc0(i*sizeof(int));
559 for (i=0; i<IDELEMS(c->S); i++)
560   {
561     int pos;
562     assume (c->S->m[i]!=NULL);
563     LObject h;
564     h.p = c->S->m[i];
565     h.pNorm();
566     c->strat->initEcart(&h);
567     assume(c->lengths[i]==pLength(h.p));
568     if (c->strat->sl==-1) pos=0;
569     else pos = simple_posInS(c->strat,h.p,c->lengths[i]);
570     h.sev = pGetShortExpVector(h.p);
571     c->strat->enterS(h,pos,c->strat);
572     c->strat->lenS[pos]=c->lengths[i];
573   }
574 //c->strat->lenS=(int*)omAlloc0(IDELEMS(c->strat->Shdl)*sizeof(int));
575 //for (i=c->strat->sl;i>=0;i--)
576 //{
577 //  pNorm(c->strat->S[i]);
578 //  c->strat->lenS[i]=pLength(c->strat->S[i]);
579 //}
580 /* initS end */
581}
582int simple_posInS (kStrategy strat, poly p,int len)
583{
584  if(strat->sl==-1) return 0;
585  polyset set=strat->S;
586  intset setL=strat->lenS;
587  int length=strat->sl;
588  int i;
589  int an = 0;
590  int en= length;
591
592  if ((len>setL[length])
593      || ((len==setL[length]) && (pLmCmp(set[length],p)== -1)))
594    return length+1;
595
596  loop
597    {
598      if (an >= en-1)
599        {
600          if ((len<setL[an])
601              || ((len==setL[length]) && (pLmCmp(set[an],p) == 1))) return an;
602          return en;
603        }
604      i=(an+en) / 2;
605      if ((len<setL[i])
606          || ((len==setL[i]) && (pLmCmp(set[i],p) == 1))) en=i;
607      //else if ((len>setL[i])
608      //|| ((len==setL[i]) && (pLmCmp(set[i],p) == -1))) an=i;
609      else an=i;
610    }
611}
612/*2
613 *if the leading term of p
614 *divides the leading term of some S[i] it will be canceled
615 */
616static inline void clearS (poly p, unsigned long p_sev,int l, int* at, int* k,
617                           kStrategy strat)
618{
619  assume(p_sev == pGetShortExpVector(p));
620  if (!pLmShortDivisibleBy(p,p_sev, strat->S[*at], ~ strat->sevS[*at])) return;
621  if (l>=strat->lenS[*at]) return;
622  PrintS("!");mflush();
623  //pDelete(&strat->S[*at]);
624  deleteInS((*at),strat);
625  (*at)--;
626  (*k)--;
627}
628void add_to_basis(poly h, int i_pos, int j_pos,calc_dat* c)
629{
630  BOOLEAN R_found=FALSE;
631  void* hp;
632  PrintS("s");
633  ++(c->n);
634  ++(c->S->ncols);
635  int i,j;
636  i=c->n-1;
637  c->T_deg=(int*) omrealloc(c->T_deg,c->n*sizeof(int));
638  c->T_deg[i]=pFDeg(h);
639  hp=omrealloc(c->rep, c->n *sizeof(int));
640  if (hp!=NULL){
641    c->rep=(int*) hp;
642  } else {
643    exit(1);
644  }
645  c->short_Exps=(long *) omrealloc(c->short_Exps ,c->n*sizeof(long));
646
647  hp=omrealloc(c->lengths, c->n *sizeof(int));
648  if (hp!=NULL){
649    c->lengths=(int*) hp;
650  } else {
651    exit(1);
652  }
653  c->misses=(int*) omrealloc(c->misses,c->n*sizeof(int));
654  c->misses[i]=0;
655  c->lengths[i]=pLength(h);
656  hp=omrealloc(c->states, c->n * sizeof(char*));
657  if (hp!=NULL){
658    c->states=(char**) hp;
659  } else {
660    exit(1);
661  }
662  c->deg=(int**) omrealloc(c->deg, c->n * sizeof(int*));
663  c->rep[i]=i;
664  hp=omalloc(i*sizeof(char));
665  if (hp!=NULL){
666    c->states[i]=(char*) hp;
667  } else {
668    exit(1);
669  }
670  c->deg[i]=(int*) omalloc(i*sizeof(int));
671  hp=omrealloc(c->S->m,c->n*sizeof(poly));
672  if (hp!=NULL){
673    c->S->m=(poly*) hp;
674  } else {
675    exit(1);
676  }
677  c->S->m[i]=h;
678  c->short_Exps[i]=p_GetShortExpVector(h,c->r);
679  for (j=0;j<i;j++){
680    c->deg[i][j]=pLcmDeg(c->S->m[i],c->S->m[j]);
681    if (c->rep[j]==j){
682      //check product criterion
683
684      c->states[i][j]=UNCALCULATED;
685
686
687      //lies I[i] under I[j] ?
688      if(p_LmShortDivisibleBy(c->S->m[i],c->short_Exps[i],c->S->m[j],~(c->short_Exps[j]),c->r)){
689        c->rep[j]=i;
690        PrintS("R"); R_found=TRUE;
691
692        c->misses[i_pos]--;
693        c->misses[j_pos]--;
694        for(int z=0;z<j;z++){
695          if (c->states[j][z]==UNCALCULATED){
696            c->states[j][z]=UNIMPORTANT;
697          }
698        }
699        for(int z=j+1;z<i;z++){
700          if (c->states[z][j]==UNCALCULATED){
701            c->states[z][j]=UNIMPORTANT;
702          }
703        }
704      }
705    }
706    else {
707      c->states[i][j]=UNIMPORTANT;
708    }
709    if ((c->lengths[i]==1) && (c->lengths[j]==1))
710      c->states[i][j]=HASTREP;
711    else if (pHasNotCF(c->S->m[i],c->S->m[j]))
712      c->states[i][j]=HASTREP;
713
714  }
715  if (c->skipped_i>0){
716    c->continue_i=c->skipped_i;
717    c->continue_j=0;
718    c->skipped_i=-1;
719  }
720
721  //i=posInS(c->strat,c->strat->sl,h,0 /*ecart*/);
722    i=simple_posInS(c->strat,h,c->lengths[c->n-1]);
723
724    LObject P; memset(&P,0,sizeof(P));
725    P.tailRing=c->r;
726    P.p=h; /*p_Copy(h,c->r);*/
727    P.FDeg=pFDeg(P.p,c->r);
728    if (!rField_is_Zp(c->r)) pCleardenom(P.p);
729    //enterT(P,c->strat,-1);
730    c->strat->enterS(P,i,c->strat);
731    c->strat->lenS[i]=/*pLength(c->strat->S[i]);*/ c->lengths[c->n-1];
732    pNorm(c->strat->S[i]);
733    if (0 /*R_found && (i<c->strat->sl)*/)
734      {
735        i++;
736        unsigned long h_sev = pGetShortExpVector(h);
737        loop
738          {
739            if (i>c->strat->sl) break;
740            clearS(h,h_sev,c->lengths[c->n-1], &i,&(c->strat->sl),c->strat);
741            i++;
742          }
743      }
744    if (c->lengths[c->n-1]==1)
745      shorten_tails(c,c->S->m[c->n-1]);
746    //you should really update c->lengths, c->strat->lenS, and the oder of polys in strat if you sort after lengths
747
748    //for(i=c->strat->sl; i>0;i--)
749    //  if(c->strat->lenS[i]<c->strat->lenS[i-1]) printf("fehler bei %d\n",i);
750}
751#if 0
752static poly redNF (poly h,kStrategy strat)
753{
754  int j = 0;
755  int z = 3;
756  unsigned long not_sev;
757
758  if (0 > strat->sl)
759    {
760      return h;
761    }
762  not_sev = ~ pGetShortExpVector(h);
763  loop
764    {
765      if (pLmShortDivisibleBy(strat->S[j], strat->sevS[j], h, not_sev))
766        {
767          //if (strat->interpt) test_int_std(strat->kIdeal);
768          /*- compute the s-polynomial -*/
769#ifdef KDEBUG
770          if (TEST_OPT_DEBUG)
771            {
772              PrintS("red:");
773              wrp(h);
774              PrintS(" with ");
775              wrp(strat->S[j]);
776            }
777#endif
778          h = ksOldSpolyRed(strat->S[j],h,strat->kNoether);
779#ifdef KDEBUG
780          if (TEST_OPT_DEBUG)
781            {
782              PrintS("\nto:");
783              wrp(h);
784              PrintLn();
785            }
786#endif
787          if (h == NULL) return NULL;
788          z++;
789          if (z>=10)
790            {
791              z=0;
792              pNormalize(h);
793            }
794          /*- try to reduce the s-polynomial -*/
795          j = 0;
796          not_sev = ~ pGetShortExpVector(h);
797        }
798      else
799        {
800          if (j >= strat->sl) return h;
801          j++;
802        }
803    }
804}
805#else
806static poly redNF2 (poly h,calc_dat* c , int &len)
807{
808  len=0;
809  if (h==NULL) return NULL;
810  int j;
811
812  kStrategy strat=c->strat;
813  len=pLength(h);
814  int len_upper_bound=len;
815  if (0 > strat->sl)
816    {
817      return h;
818    }
819  LObject P(h);
820  P.SetShortExpVector();
821  P.bucket = kBucketCreate(currRing);
822  kBucketInit(P.bucket,P.p,len /*pLength(P.p)*/);
823  //int max_pos=simple_posInS(strat,P.p);
824  loop
825    {
826      j=kFindDivisibleByInS(strat->S,strat->sevS,strat->sl,&P);
827      if (j>=0)
828        {
829          poly sec_copy=NULL;
830          //pseudo code
831          bool must_replace_in_basis=(len_upper_bound<=strat->lenS[j]);//first test
832          if (must_replace_in_basis)
833            {
834              //second test
835              if (pLmEqual(P.p,strat->S[j]))
836                {
837                  PrintS("b");
838                  sec_copy=kBucketClear(P.bucket);
839                  kBucketInit(P.bucket,pCopy(sec_copy),pLength(sec_copy));
840                }
841              else
842                {
843                  must_replace_in_basis=false;
844                }
845            }
846          nNormalize(pGetCoeff(P.p));
847#ifdef KDEBUG
848          if (TEST_OPT_DEBUG)
849            {
850              PrintS("red:");
851              wrp(h);
852              PrintS(" with ");
853              wrp(strat->S[j]);
854            }
855#endif
856          len_upper_bound=len_upper_bound+strat->lenS[j]-2;
857          number coef=kBucketPolyRed(P.bucket,strat->S[j],
858                                     strat->lenS[j]/*pLength(strat->S[j])*/,
859  strat->kNoether);
860 nDelete(&coef);
861 h = kBucketGetLm(P.bucket);
862 //pseudo code
863 if (must_replace_in_basis){
864   int old_pos_in_c=-1;
865   poly p=strat->S[j];
866   int z;
867   
868   int new_length=pLength(sec_copy);
869   Print("%i",strat->lenS[j]-new_length);
870   len_upper_bound=new_length +strat->lenS[j]-2;//old entries length
871   int new_pos=simple_posInS(c->strat,strat->S[j],new_length);
872   strat->S[j]=sec_copy;
873   c->strat->lenS[j]=new_length;
874   
875   for (z=c->n;z;z--)
876     {
877       if(p==c->S->m[z-1])
878         {
879           c->S->m[z-1]=sec_copy;
880           old_pos_in_c=z-1;
881           c->lengths[z-1]=new_length;
882           if (new_length==1)
883             {
884               int i;
885               for ( i=0;i<z-1;i++)
886                 {
887                   if (c->lengths[i]==1)
888                     c->states[z-1][i]=HASTREP;
889                 }
890               for ( i=z;i<c->n;i++){
891                 if (c->lengths[i]==1)
892                   c->states[i][z-1]=HASTREP;
893               }
894               shorten_tails(c,sec_copy);
895             }
896           break;
897         }
898     }
899   pDelete(&p);
900
901   //   replace_quietly(c,j,sec_copy);
902   // have to do many additional things for consistency
903   {
904     
905     
906     
907
908     int old_pos=j;
909     assume(new_pos<=old_pos);
910 
911
912     c->strat->lenS[old_pos]=new_length;
913     int i=0;
914     for(i=new_pos;i<old_pos;i++){
915       if (strat->lenS[i]<=new_length)
916         new_pos++;
917       else
918         break;
919     }
920     if (new_pos<old_pos)
921       move_forward_in_S(old_pos,new_pos,c->strat);
922   
923
924   }
925 }
926 if (h==NULL) return NULL;
927 P.p=h;
928 P.t_p=NULL;
929 P.SetShortExpVector();
930#ifdef KDEBUG
931 if (TEST_OPT_DEBUG)
932   {
933     PrintS("\nto:");
934     wrp(h);
935     PrintLn();
936   }
937#endif
938        }
939    else
940      {
941        kBucketClear(P.bucket,&(P.p),&len);
942        kBucketDestroy(&P.bucket);
943        pNormalize(P.p);
944        return P.p;
945      }
946    }
947}
948static poly redNF (poly h,kStrategy strat, int &len)
949{
950  len=0;
951  if (h==NULL) return NULL;
952  int j;
953
954  len=pLength(h);
955  if (0 > strat->sl)
956    {
957      return h;
958    }
959  LObject P(h);
960  P.SetShortExpVector();
961  P.bucket = kBucketCreate(currRing);
962  kBucketInit(P.bucket,P.p,len /*pLength(P.p)*/);
963  //int max_pos=simple_posInS(strat,P.p);
964  loop
965    {
966      j=kFindDivisibleByInS(strat->S,strat->sevS,strat->sl,&P);
967      if (j>=0)
968        {
969          nNormalize(pGetCoeff(P.p));
970#ifdef KDEBUG
971          if (TEST_OPT_DEBUG)
972            {
973              PrintS("red:");
974              wrp(h);
975              PrintS(" with ");
976              wrp(strat->S[j]);
977            }
978#endif
979          number coef=kBucketPolyRed(P.bucket,strat->S[j],
980                                     strat->lenS[j]/*pLength(strat->S[j])*/,
981  strat->kNoether);
982 nDelete(&coef);
983 h = kBucketGetLm(P.bucket);
984 if (h==NULL) return NULL;
985 P.p=h;
986 P.t_p=NULL;
987 P.SetShortExpVector();
988#ifdef KDEBUG
989 if (TEST_OPT_DEBUG)
990   {
991     PrintS("\nto:");
992     wrp(h);
993     PrintLn();
994   }
995#endif
996        }
997    else
998      {
999        kBucketClear(P.bucket,&(P.p),&len);
1000        kBucketDestroy(&P.bucket);
1001        pNormalize(P.p);
1002        return P.p;
1003      }
1004    }
1005}
1006#endif
1007#ifdef REDTAIL_S
1008poly redNFTail (poly h,const int sl,kStrategy strat, int len)
1009{
1010  if (h==NULL) return NULL;
1011  if (pNext(h)==NULL) return h;
1012  pTest(h);
1013  if (0 > sl)
1014    return h;
1015
1016  int j;
1017  poly res=h;
1018  poly act=res;
1019  LObject P(pNext(h));
1020  pNext(res)=NULL;
1021  P.bucket = kBucketCreate(currRing);
1022  len--;
1023  h=P.p;
1024  if (len <=0) len=pLength(h);
1025  kBucketInit(P.bucket,h /*P.p*/,len /*pLength(P.p)*/);
1026  pTest(h);
1027  loop
1028    {
1029      P.p=h;
1030      P.t_p=NULL;
1031      P.SetShortExpVector();
1032      loop
1033        {
1034          j=kFindDivisibleByInS(strat->S,strat->sevS,sl,&P);
1035          if (j>=0)
1036            {
1037#ifdef REDTAIL_PROT
1038              PrintS("r");
1039#endif
1040              nNormalize(pGetCoeff(P.p));
1041#ifdef KDEBUG
1042              if (TEST_OPT_DEBUG)
1043                {
1044                  PrintS("red tail:");
1045                  wrp(h);
1046                  PrintS(" with ");
1047                  wrp(strat->S[j]);
1048                }
1049#endif
1050              number coef;
1051              pTest(strat->S[j]);
1052              coef=kBucketPolyRed(P.bucket,strat->S[j],
1053                                  strat->lenS[j]/*pLength(strat->S[j])*/,strat->kNoether);
1054              pMult_nn(res,coef);
1055              nDelete(&coef);
1056              h = kBucketGetLm(P.bucket);
1057              pTest(h);
1058              if (h==NULL)
1059                {
1060#ifdef REDTAIL_PROT
1061                  PrintS(" ");
1062#endif
1063                  return res;
1064                }
1065              pTest(h);
1066              P.p=h;
1067              P.t_p=NULL;
1068              P.SetShortExpVector();
1069#ifdef KDEBUG
1070              if (TEST_OPT_DEBUG)
1071                {
1072                  PrintS("\nto tail:");
1073                  wrp(h);
1074                  PrintLn();
1075                }
1076#endif
1077            }
1078          else
1079            {
1080#ifdef REDTAIL_PROT
1081              PrintS("n");
1082#endif
1083              break;
1084            }
1085        } /* end loop current mon */
1086      poly tmp=pHead(h /*kBucketGetLm(P.bucket)*/);
1087      act->next=tmp;pIter(act);
1088      poly tmp2=pHead(h);
1089      pNeg(tmp2);
1090      int ltmp2=1;
1091      pTest(tmp2);
1092      kBucket_Add_q(P.bucket, tmp2, &ltmp2);
1093
1094      h = kBucketGetLm(P.bucket);
1095      if (h==NULL)
1096        {
1097#ifdef REDTAIL_PROT
1098          PrintS(" ");
1099#endif
1100          return res;
1101        }
1102      pTest(h);
1103    }
1104}
1105#endif
1106
1107void do_this_spoly_stuff(int i,int j,calc_dat* c){
1108  poly f=c->S->m[i];
1109  poly g=c->S->m[j];
1110  poly h=ksOldCreateSpoly(f, g, NULL, c->r);
1111  poly hr=NULL;
1112#ifdef FULLREDUCTIONS
1113  if (h!=NULL)
1114    {
1115      int len;
1116
1117      hr=redNF2(h,c,len);
1118#ifdef HALFREDUCTIONS
1119      int real_sl=c->strat->sl;
1120      int l;
1121      for(l=0;l<c->n;l++){
1122        if (c->strat->lenS[l]>4)
1123          break;
1124      }
1125      c->strat->sl=l;
1126#endif
1127      if (hr!=NULL)
1128#ifdef REDTAIL_S
1129        hr = redNFTail(hr,c->strat->sl,c->strat,len);
1130#else
1131      hr = redtailBba(hr,c->strat->sl,c->strat);
1132#endif
1133#ifdef HALFREDUCTIONS
1134      c->strat->sl=real_sl;
1135#endif
1136    }
1137#else
1138  if (h!=NULL)
1139    {
1140      int len;
1141      hr=redNF2(h,c,len);
1142    }
1143#endif
1144  c->normal_forms++;
1145  if (hr==NULL)
1146    {
1147      PrintS("-");
1148      c->misses_counter++;
1149      c->misses[i]++;
1150      c->misses[j]++;
1151      c->misses_series++;
1152    }
1153  else
1154    {
1155      c->misses_series=0;
1156#ifdef HEAD_BIN
1157      hr=p_MoveHead(hr,c->HeadBin);
1158#endif
1159      add_to_basis(hr, i,j,c);
1160    }
1161}
1162
1163ideal t_rep_gb(ring r,ideal arg_I){
1164  ideal I_temp=idCopy(arg_I); //kInterRed(arg_I);
1165  ideal I=idCompactify(I_temp);
1166  idDelete(&I_temp);
1167  calc_dat* c=(calc_dat*) omalloc(sizeof(calc_dat));
1168  c->r=currRing;
1169  c->S=I;
1170  initial_data(c);
1171  while (find_next_pair(c)){
1172    int i,j;
1173    int z;
1174    i=c->found_i;
1175    j=c->found_j;
1176    replace_pair(i,j,c);
1177    if (!(c->states[max(i,j)][min(i,j)]==HASTREP)){
1178      do_this_spoly_stuff(i,j,c);
1179    }
1180    //else
1181    //  PrintS("f");
1182
1183    now_t_rep(i,j,c);
1184    now_t_rep(c->found_i,c->found_j,c);
1185#ifdef RANDOM_WALK
1186    c->continue_i=my_rand(c->n-1)+1;
1187    c->continue_j=my_rand(c->continue_i);
1188#endif
1189    int_pair_node* h=c->soon_free;
1190    while(h!=NULL)
1191      {
1192        int_pair_node* s=h;
1193        now_t_rep(h->a,h->b,c);
1194
1195        h=h->next;
1196        omfree(s);
1197      }
1198
1199
1200  }
1201  omfree(c->rep);
1202  for(int z=0;z<c->n;z++){
1203    omfree(c->states[z]);
1204  }
1205  omfree(c->states);
1206  omfree(c->lengths);
1207  printf("calculated %d NFs\n",c->normal_forms);
1208  omfree(c);
1209  return(I);
1210}
1211void now_t_rep(int arg_i, int arg_j, calc_dat* c){
1212  int i,j;
1213  if (arg_i==arg_j){
1214    return;
1215  }
1216  if (arg_i>arg_j){
1217    i=arg_j;
1218    j=arg_i;
1219  } else {
1220    i=arg_i;
1221    j=arg_j;
1222  }
1223  c->states[j][i]=HASTREP;
1224}
1225bool has_t_rep(const int & arg_i, const  int & arg_j, calc_dat* state){
1226
1227  if (arg_i==arg_j)
1228    {
1229      return (true);
1230    }
1231  if (arg_i>arg_j)
1232    {
1233      return (state->states[arg_i][arg_j]==HASTREP);
1234    } else
1235      {
1236        return (state->states[arg_j][arg_i]==HASTREP);
1237      }
1238}
1239int pLcmDeg(poly a, poly b)
1240{
1241  int i;
1242  int n=0;
1243  for (i=pVariables; i; i--)
1244    {
1245      n+=max( pGetExp(a,i), pGetExp(b,i));
1246    }
1247  return n;
1248
1249}
1250int pMinDeg3(poly f){
1251  if (f==NULL){
1252    return(-1);
1253
1254  }
1255
1256  poly h=f->next;
1257  int n=pFDeg(h);
1258  int i=0;
1259  while((i<2)){
1260    if (h==NULL)
1261      return(n);
1262    h=h->next;
1263    if (h!=NULL){
1264      n=min(n,pFDeg(h));
1265    }
1266    i++;
1267  }
1268
1269  return n;
1270}
1271
1272
1273void shorten_tails(calc_dat* c, poly monom)
1274{
1275
1276  for(int i=0;i<c->n;i++)
1277    {
1278      //enter tail
1279      if (c->rep[i]!=i) continue;
1280      if (c->S->m[i]==NULL) continue;
1281      poly tail=c->S->m[i]->next;
1282      poly prev=c->S->m[i];
1283      bool did_something=false;
1284      while((tail!=NULL)&& (pLmCmp(tail, monom)>=0))
1285        {
1286          if (p_LmDivisibleBy(monom,tail,c->r))
1287            {
1288              did_something=true;
1289              prev->next=tail->next;
1290              tail->next=NULL;
1291              p_Delete(& tail,c->r);
1292              tail=prev;
1293              //PrintS("Shortened");
1294              c->lengths[i]--;
1295            }
1296          prev=tail;
1297          tail=tail->next;
1298        }
1299      if (did_something)
1300        {
1301          int new_pos=simple_posInS(c->strat,c->S->m[i],c->lengths[i]);
1302          int old_pos=-1;
1303          //assume new_pos<old_pos
1304          for (int z=new_pos;z<=c->strat->sl;z++)
1305            {
1306              if (c->strat->S[z]==c->S->m[i])
1307                {
1308                  old_pos=z;
1309                  break;
1310                }
1311            }
1312          if (old_pos== -1)
1313            for (int z=new_pos-1;z>=0;z--)
1314              {
1315                if (c->strat->S[z]==c->S->m[i])
1316                  {
1317                    old_pos=z;
1318                    break;
1319                  }
1320              }
1321          assume(old_pos>=0);
1322          assume(pLength(c->strat->S[old_pos])==c->lengths[i]);
1323          c->strat->lenS[old_pos]=c->lengths[i];
1324          if (new_pos<old_pos)
1325            move_forward_in_S(old_pos,new_pos,c->strat);
1326          if (c->lengths[i]==1)
1327            {
1328             
1329              int j;
1330              for ( j=0;j<i;j++)
1331                {
1332                  if (c->lengths[j]==1)
1333                    c->states[i][j]=HASTREP;
1334                }
1335              for ( j=i+1;j<c->n;j++){
1336                if (c->lengths[j]==1)
1337                  c->states[j][i]=HASTREP;
1338              }
1339              shorten_tails(c,c->S->m[i]);
1340            }
1341        }
1342    }
1343}
Note: See TracBrowser for help on using the repository browser.