source: git/Singular/tgb.cc @ cae0b6

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