source: git/kernel/F4.cc @ 210e07

spielwiese
Last change on this file since 210e07 was 599326, checked in by Kai Krüger <krueger@…>, 14 years ago
Anne, Kai, Frank: - changes to #include "..." statements to allow cleaner build structure - affected directories: omalloc, kernel, Singular - not yet done: IntergerProgramming git-svn-id: file:///usr/local/Singular/svn/trunk@13032 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 18.5 KB
Line 
1#include <kernel/mod2.h>
2#include <kernel/F4.h>
3#include <kernel/tgb_internal.h>
4#include <kernel/tgbgauss.h>
5/****************************************
6*  Computer Algebra System SINGULAR     *
7****************************************/
8/* $Id$ */
9/*
10* ABSTRACT: F4 implementation
11*/
12static int posInPolys (poly*  p, int pn, poly qe,slimgb_alg* c)
13{
14  if(pn==0) return 0;
15
16  int length=pn-1;
17  int i;
18  int an = 0;
19  int en= length;
20
21  if (pLmCmp(qe,p[en])==1)
22    return length+1;
23
24  while(1)
25  {
26      //if (an >= en-1)
27    if(en-1<=an)
28    {
29      if (pLmCmp(p[an],qe)==1) return an;
30      return en;
31    }
32    i=(an+en) / 2;
33    if (pLmCmp(p[i],qe)==1)
34      en=i;
35    else an=i;
36  }
37}
38
39static tgb_sparse_matrix* build_sparse_matrix(poly* p,int p_index,poly* done, int done_index, slimgb_alg* c){
40  tgb_sparse_matrix* t=new tgb_sparse_matrix(p_index,done_index,c->r);
41  int i, pos;
42  //  Print("\n 0:%s\n",pString(done[done_index-1]));
43  //Print("\n 1:%s\n",pString(done[done_index-2]));
44  //  assume((!(pLmEqual(done[done_index-1],done[done_index-2]))));
45#ifdef TGB_DEGUG
46  for(i=0;i<done_index;i++)
47{
48  int j;
49  for(j=0;j<i;j++)
50  {
51    assume((!(pLmEqual(done[i],done[j]))));
52  }
53}
54#endif
55  for(i=0;i<p_index;i++)
56{ 
57    // Print("%i ter Eintrag:%s\n",i,pString(p[i]));
58  mac_poly m=NULL;
59  mac_poly* set_this=&m;
60  poly p_i=p[i];
61  while(p_i)
62  {
63
64    int v=-1;
65    pos=posInPolys (done, done_index, p_i,c);
66    if((done_index>pos)&&(pLmEqual(p_i,done[pos])))
67      v=pos;
68    if((pos>0) &&(pLmEqual(p_i,done[pos-1])))
69      v=pos-1;
70    assume(v!=-1);
71      //v is ascending ordered, we need descending order
72    v=done_index-1-v;
73    (*set_this)=new mac_poly_r();
74    (*set_this)->exp=v;
75       
76    (*set_this)->coef=nCopy(p_i->coef);
77    set_this=&(*set_this)->next;
78    p_i=p_i->next;
79
80  }
81  init_with_mac_poly(t,i,m);
82}
83  return t;
84}
85static tgb_matrix* build_matrix(poly* p,int p_index,poly* done, int done_index, slimgb_alg* c){
86  tgb_matrix* t=new tgb_matrix(p_index,done_index);
87  int i, pos;
88  //  Print("\n 0:%s\n",pString(done[done_index-1]));
89  //Print("\n 1:%s\n",pString(done[done_index-2]));
90  assume((!(pLmEqual(done[done_index-1],done[done_index-2]))));
91#ifdef TGB_DEGUG
92  for(i=0;i<done_index;i++)
93{
94  int j;
95  for(j=0;j<i;j++)
96  {
97    assume((!(pLmEqual(done[i],done[j]))));
98  }
99}
100#endif
101  for(i=0;i<p_index;i++)
102{ 
103    // Print("%i ter Eintrag:%s\n",i,pString(p[i]));
104  poly p_i=p[i];
105  while(p_i)
106  {
107
108    int v=-1;
109    pos=posInPolys (done, done_index, p_i,c);
110    if((done_index>pos)&&(pLmEqual(p_i,done[pos])))
111      v=pos;
112    if((pos>0) &&(pLmEqual(p_i,done[pos-1])))
113      v=pos-1;
114    assume(v!=-1);
115      //v is ascending ordered, we need descending order
116    v=done_index-1-v;
117    number nt=t->get(i,v);
118    nDelete(&nt);
119    t->set(i,v,nCopy(p_i->coef));
120    p_i=p_i->next;
121  }
122}
123  return t;
124}
125
126static int retranslate(poly* m,tgb_sparse_matrix* mat,poly* done, slimgb_alg* c){
127  int i;
128  int m_index=0;
129  for(i=0;i<mat->get_rows();i++)
130  {
131    if(mat->zero_row(i))
132    {
133      mat->free_row(i);
134      continue;
135    }
136    m[m_index++]= free_row_to_poly(mat, i, done, mat->get_columns());
137
138  }
139
140  delete mat;
141  return m_index;
142
143}
144//!returns m_index and destroys mat
145 static int retranslate(poly* m,tgb_matrix* mat,poly* done, slimgb_alg* c){
146  int i;
147  int m_index=0;
148  for(i=0;i<mat->get_rows();i++)
149  {
150    if(mat->zero_row(i))
151    {
152      mat->free_row(i);
153      continue;
154    }
155   
156    m[m_index]=pInit();
157    int v=mat->min_col_not_zero_in_row(i);
158    //v=done_index-1-pos; => pos=done_index-1-v=mat->get_columns()-1-v
159    int pos=mat->get_columns()-1-v;
160    int* ev=(int*) omalloc((c->r->N+1)*sizeof(int));
161    pGetExpV(done[pos],ev);
162    pSetExpV(m[m_index],ev);
163    omfree(ev);
164   
165    poly p=m[m_index];
166    pSetCoeff(p,mat->get(i,v));
167    while((v=mat->next_col_not_zero(i,v))!=mat->get_columns())
168    {
169      poly pn=pInit();
170     
171      //v=done_index-1-pos; => pos=done_index-1-v=mat->get_columns()-1-v
172      pos=mat->get_columns()-1-v;
173      ev=(int*) omalloc((c->r->N+1)*sizeof(int));
174      pGetExpV(done[pos],ev);
175      pSetExpV(pn,ev);
176      omfree(ev);
177      pSetCoeff(pn,mat->get(i,v));
178      p->next=pn;
179      p=pn;
180    }
181    p->next=NULL;
182    mat->free_row(i, FALSE);
183    m_index++;
184  }
185
186  delete mat;
187  return m_index;
188
189}
190static int pLmCmp_func(const void* ap1, const void* ap2){
191    poly p1,p2;
192    p1=*((poly*) ap1);
193    p2=*((poly*)ap2);
194
195    return pLmCmp(p1,p2);
196}
197static int monom_poly_crit(const void* ap1, const void* ap2){
198  monom_poly* p1;
199  monom_poly* p2;
200  p1=((monom_poly*) ap1);
201  p2=((monom_poly*)ap2);
202  if(((unsigned long) p1->f)>((unsigned long) p2->f)) return 1;
203  if(((unsigned long) p1->f)<((unsigned long)p2->f)) return -1; 
204
205  return pLmCmp(p1->m,p2->m);
206 
207}
208static int posInMonomPolys (monom_poly*  p, int pn, monom_poly & qe,slimgb_alg* c)
209{
210  if(pn==0) return 0;
211
212  int length=pn-1;
213  int i;
214  int an = 0;
215  int en= length;
216
217  if (monom_poly_crit(&qe,&p[en])==1)
218    return length+1;
219
220  while(1)
221  {
222      //if (an >= en-1)
223    if(en-1<=an)
224    {
225      if (monom_poly_crit(&p[an],&qe)==1) return an;
226      return en;
227    }
228    i=(an+en) / 2;
229    if (monom_poly_crit(&p[i],&qe)==1)
230      en=i;
231    else an=i;
232  }
233}
234static void simplify(monom_poly& h, slimgb_alg* c){
235  mp_array_list* F=c->F;
236  poly_array_list* F_minus=c->F_minus;
237  while(F)
238  {
239    assume(F_minus!=NULL);
240    int i;
241    int posm=posInMonomPolys (F->mp, F->size, h,c);
242#ifdef TGB_DEBUG
243#endif
244     //for(i=0;i<F->size;i++)
245    for(i=si_min(posm,F->size-1);i>=0;i--)
246{
247      //      if((i>=posm)&&(F->mp[i].f!=h.f)) break;
248  if((i<posm)&&(F->mp[i].f!=h.f)) break;
249  if ((h.f==F->mp[i].f) &&(p_LmDivisibleBy(F->mp[i].m,h.m,c->r)))
250  {
251       
252        //      Print("found");
253       
254          //according to the algorithm you should test (!(pIsConstant(F[i].m)))
255          //but I think this is only because of bad formulation
256    int j;
257       
258    poly lm=pLmInit(h.f);
259    pSetCoeff(lm,nInit(1));
260
261    lm=pMult_mm(lm,F->mp[i].m);
262
263    int pos;
264    j=-1;
265    pos=posInPolys (F_minus->p, F_minus->size, lm,c);
266    if((F_minus->size>pos)&&(pLmEqual(lm,F_minus->p[pos])))
267      j=pos;
268    if((pos>0) &&(pLmEqual(lm,F_minus->p[pos-1])))
269      j=pos-1;
270    assume(j!=-1);
271      //        if(j==-1) Print("\n jAltert \n");
272//      for(j=0;j<F_minus->size;j++)
273//      {
274//        if (pLmEqual(F_minus->p[j],lm))
275//          break;
276//      }
277//      assume(j<F_minus->size);
278    pDelete(&lm);
279    if(j>=0)
280    {
281      pExpVectorSub(h.m,F->mp[i].m);
282      h.f=F_minus->p[j];
283    }
284    break;
285  }
286}
287    F=F->next;
288    F_minus=F_minus->next;
289  }
290  assume(F_minus==NULL);
291}
292
293void go_on_F4 (slimgb_alg* c){
294  //set limit of 1000 for multireductions, at the moment for
295  //programming reasons
296  int max_par;
297  if (c->is_homog)
298    max_par=PAR_N_F4;
299  else
300    max_par=200;
301  int done_size=max_par*2;
302  poly* done=(poly*) omalloc(done_size*sizeof(poly));
303  int done_index=0; //done_index must always be smaller than done_size
304  int chosen_size=max_par*2;
305  monom_poly* chosen=(monom_poly*) omalloc(chosen_size*sizeof(monom_poly));
306  int chosen_index=0;
307  //  monom_poly* vgl=(monom_poly*) omalloc(chosen_size*sizeof(monom_poly));
308  int i=0;
309  c->average_length=0;
310  for(i=0;i<c->n;i++){
311    c->average_length+=c->lengths[i];
312  }
313  c->average_length=c->average_length/c->n;
314  i=0;
315  poly* p;
316  int nfs=0;
317  int curr_deg=-1;
318
319  //choose pairs and preprocess symbolically
320  while(chosen_index<max_par)
321  {
322   
323    //    sorted_pair_node* s=c->apairs[c->pair_top];
324    sorted_pair_node* s;
325    if(c->pair_top>=0)
326      s=top_pair(c);//here is actually chain criterium done
327    else
328      s=NULL;
329    if (!s) break;
330    nfs++;
331    if(curr_deg>=0)
332    {
333      if (s->deg >curr_deg) break;
334    }
335
336    else curr_deg=s->deg;
337    quick_pop_pair(c);
338    if(s->i>=0)
339    {
340      //replace_pair(s->i,s->j,c);
341      if(s->i==s->j) 
342      {
343        free_sorted_pair_node(s,c->r);
344        continue;
345      }
346    }
347    assume(s->i>=0);
348    monom_poly h;
349    BOOLEAN only_free=FALSE;
350    if(s->i>=0)
351    {
352     
353      poly lcm=pOne();
354     
355      pLcm(c->S->m[s->i], c->S->m[s->j], lcm);
356      pSetm(lcm);
357      assume(lcm!=NULL);
358      poly factor1=pCopy(lcm);
359      pExpVectorSub(factor1,c->S->m[s->i]);
360      poly factor2=pCopy(lcm);
361      pExpVectorSub(factor2,c->S->m[s->j]);
362
363      if(done_index>=done_size)
364      {
365        done_size+=max_par;
366        done=(poly*) omrealloc(done,done_size*sizeof(poly));
367      }
368      done[done_index++]=lcm;
369      if(chosen_index+1>=chosen_size)
370      {
371        //max_par must be greater equal 2
372        chosen_size+=si_max(max_par,2);
373        chosen=(monom_poly*) omrealloc(chosen,chosen_size*sizeof(monom_poly));
374      }
375      h.f=c->S->m[s->i];
376      h.m=factor1;
377      chosen[chosen_index++]=h;
378      h.f=c->S->m[s->j];
379      h.m=factor2;
380      chosen[chosen_index++]=h;
381    }
382    else
383    {
384      if(chosen_index>=chosen_size)
385      {
386        chosen_size+=max_par;
387        chosen=(monom_poly*) omrealloc(chosen,chosen_size*sizeof(monom_poly));
388      }
389      h.f=s->lcm_of_lm;
390      h.m=pOne();
391      //pSetm(h.m); done by pOne
392      chosen[chosen_index++]=h;
393      //must carefull remember to destroy such a h;
394      poly_list_node* next=c->to_destroy;
395     
396      c->to_destroy=(poly_list_node*) omalloc(sizeof(poly_list_node));
397      c->to_destroy->p=h.f;
398      c->to_destroy->next=next;
399      only_free=TRUE;
400    }
401
402    if(s->i>=0)
403      now_t_rep(s->j,s->i,c);
404
405    if(!(only_free))
406      free_sorted_pair_node(s,c->r);
407    else
408      omfree(s);
409   
410 
411
412  }
413  c->normal_forms+=nfs;
414  if (TEST_OPT_PROT)
415    Print("M[%i, ",nfs);
416  //next Step, simplify all pairs
417  for(i=0;i<chosen_index;i++)
418  {
419    //    PrintS("simplifying ");
420    //Print("from %s",pString(chosen[i].m));
421    //Print(", %s", pString(chosen[i].f));
422    simplify(chosen[i],c);
423    //Print(" to %s",pString(chosen[i].m));
424    //Print(", %s \n", pString(chosen[i].f));
425  }
426 
427  //next Step remove duplicate entries
428  qsort(chosen,chosen_index,sizeof(monom_poly),monom_poly_crit);
429  int pos=0;
430  for(i=1;i<chosen_index;i++)
431  {
432    if(((chosen[i].f!=chosen[pos].f))||(!(pLmEqual(chosen[i].m,chosen[pos].m))))
433      chosen[++pos]=chosen[i];
434    else pDelete(&(chosen[i].m));
435  }
436  if(chosen_index>0)
437    chosen_index=pos+1;
438  //next step process polys
439  int p_size=2*chosen_index;
440  int p_index=0;
441  p=(poly*) omalloc(p_size*sizeof(poly));
442  for(p_index=0;p_index<chosen_index;p_index++)
443  {
444    p[p_index]=ppMult_mm(chosen[p_index].f,chosen[p_index].m);
445  }
446
447  qsort(done, done_index,sizeof(poly),pLmCmp_func);
448  pos=0;
449  //Print("Done_index:%i",done_index);
450  if(done_index>0)
451  {
452    pTest(done[0]);
453  }
454  for(i=1;i<done_index;i++)
455  {
456    pTest(done[i]);
457    if((!(pLmEqual(done[i],done[pos]))))
458      done[++pos]=done[i];
459    else pDelete(&(done[i]));
460  }
461  if(done_index>0)
462    done_index=pos+1;
463#ifdef TGB_DEBUG
464  for(i=0;i<done_index;i++)
465{
466  int j;
467  for(j=0;j<i;j++)
468  {
469    assume((!pLmEqual(done[j],done[i])));
470  }
471}
472#endif
473#ifdef TGB_DEBUG
474  for(i=0;i<done_index;i++)
475{
476  pTest(done[i]);
477}
478#endif
479  poly* m;
480  int m_index=0;
481  int m_size=0;
482  for(i=0;i<p_index;i++)
483{
484  m_size+=pLength(p[i]);
485}
486  m=(poly*) omalloc(m_size*sizeof(poly));
487  //q=(poly*) omalloc(m_size*sizeof(poly));
488
489 
490
491  for(i=0;i<p_index;i++)
492{
493
494  poly p_i=p[i];
495
496  pTest(p[i]);
497
498  while(p_i)
499  {
500
501    m[m_index]=pLmInit(p_i);
502
503    pSetCoeff(m[m_index],nInit(1));
504
505    p_i=p_i->next;
506
507    m_index++;
508
509  }
510}
511
512  int q_size=m_index;
513  poly* q=(poly*) omalloc(q_size*sizeof(poly));
514  int q_index=0;
515  //next Step additional reductors
516
517  while(m_index>0)
518{
519#ifdef TGB_DEBUG
520     
521      for(i=0;i<done_index;i++)
522{
523  pTest(done[i]);
524}
525#endif
526
527    qsort(m, m_index,sizeof(poly),pLmCmp_func);
528
529   
530    pos=0;
531#ifdef TGB_DEBUG
532     
533      for(i=0;i<done_index;i++)
534{
535       
536  pTest(done[i]);
537}
538#endif
539    for(i=1;i<m_index;i++)
540{
541  pTest(m[i]);
542  pTest(m[pos]);
543  if((!(pLmEqual(m[i],m[pos]))))
544    m[++pos]=m[i];
545  else pDelete(&(m[i]));
546}
547    if(m_index>1)
548      m_index=pos+1;
549    if(done_size<done_index+m_index)
550{
551  done_size=done_index+2*m_index;
552  done=(poly*) omrealloc(done,done_size*sizeof(poly));
553}
554    if(chosen_size<chosen_index+m_index)
555{
556  chosen_size=chosen_index+2*m_index;
557  chosen=(monom_poly*) omrealloc(chosen,chosen_size*sizeof(monom_poly));
558}
559    q_index=0;
560    if(q_size<m_index)
561{
562  q_size=2*m_index;
563  omfree(q);
564  q=(poly*) omalloc(q_size*sizeof(poly));
565}
566
567    for(i=0;i<m_index;i++)
568{
569
570#ifdef TGB_DEBUG
571{
572  int my_i;
573  for(my_i=0;my_i<done_index;my_i++)
574  {
575    int my_j;
576    for(my_j=0;my_j<my_i;my_j++)
577    {
578      assume((!pLmEqual(done[my_j],done[my_i])));
579    }
580  }
581}
582#endif
583      BOOLEAN in_done=FALSE;
584      pTest(m[i]);
585
586      pos=posInPolys (done, done_index, m[i],c);
587#ifdef TGB_DEBUG
588{
589  int my_i;
590  for (my_i=0;my_i<done_index;my_i++)
591    pTest(done[my_i]);
592}
593#endif     
594      if(((done_index>pos)&&(pLmEqual(m[i],done[pos]))) ||(pos>0) &&(pLmEqual(m[i],done[pos-1])))
595        in_done=TRUE;
596      if (!(in_done))
597{
598       
599  int S_pos=kFindDivisibleByInS_easy(c->strat,m[i], pGetShortExpVector(m[i]));
600  if(S_pos>=0)
601  {
602    monom_poly h;
603    h.f=c->strat->S[S_pos];
604    h.m=pOne();
605    int* ev=(int*) omalloc((c->r->N+1)*sizeof(int));
606    pGetExpV(m[i],ev);
607    pSetExpV(h.m,ev);
608    omfree(ev);
609    pExpVectorSub(h.m,c->strat->S[S_pos]);
610    simplify(h,c);
611    q[q_index]=ppMult_mm(h.f,h.m);
612    chosen[chosen_index++]=h;
613    q_index++;
614  }
615  pTest(m[i]);
616#ifdef TGB_DEBUG
617{
618  int my_i;
619  for (my_i=0;my_i<done_index;my_i++)
620    pTest(done[my_i]);
621}
622#endif     
623        memmove(&(done[pos+1]),&(done[pos]), (done_index-pos)*sizeof(poly));
624        done[pos]=m[i];
625 done_index++;
626#ifdef TGB_DEBUG
627{
628  int my_i;
629  for (my_i=0;my_i<done_index;my_i++)
630  {
631            //      Print("Position %i pos %i size %i\n",my_i,pos,done_index);
632    pTest(done[my_i]);
633  }
634}
635#endif     
636}
637      else
638        pDelete(&m[i]);
639#ifdef TGB_DEBUG
640{
641  int my_i;
642  for(my_i=0;my_i<done_index;my_i++)
643  {
644    pTest(done[my_i]);
645    int my_j;
646    for(my_j=0;my_j<my_i;my_j++)
647    {
648      assume((!pLmEqual(done[my_j],done[my_i])));
649    }
650  }
651}
652#endif
653}
654    if(p_size<p_index+q_index)
655{
656  p_size=p_index+2*q_index;
657  p=(poly*) omrealloc(p,p_size*sizeof(poly));
658}
659    for (i=0;i<q_index;i++)
660      p[p_index++]=q[i];
661    m_index=0;
662    int sum=0;
663    for(i=0;i<p_index;i++)
664{
665  sum+=pLength(p[i])-1;
666}
667    if (m_size<sum)
668{
669  omfree(m);
670  m=(poly*) omalloc(sum*sizeof(poly));
671}
672    m_size=sum;
673    for(i=0;i<q_index;i++)
674{
675  poly p_i=q[i]->next;
676  while(p_i)
677  {
678    m[m_index]=pLmInit(p_i);
679    pSetCoeff(m[m_index],nInit(1));
680    p_i=p_i->next;
681    m_index++;
682  }
683}
684    //qsort(done, done_index,sizeof(poly),pLmCmp_func);
685#ifdef TGB_DEBUG
686    for(i=0;i<done_index;i++)
687{
688  pTest(done[i]);
689}
690#endif
691}
692#ifdef TGB_DEBUG
693  for(i=0;i<done_index;i++)
694{
695  int j;
696  for(j=0;j<i;j++)
697  {
698    assume((!pLmEqual(done[j],done[i])));
699  }
700}
701#endif
702  omfree(m);
703  omfree(q);
704  if (TEST_OPT_PROT)
705    Print("Mat[%i x %i], ",chosen_index, done_index);
706  //next Step build matrix
707#ifdef TGB_DEBUG
708  for(i=0;i<done_index;i++)
709{
710  int j;
711  for(j=0;j<i;j++)
712  {
713    assume((!pLmEqual(done[j],done[i])));
714  }
715}
716#endif
717  assume(p_index==chosen_index);
718 
719 //  tgb_matrix* mat=build_matrix(p,p_index,done, done_index,c);
720 
721//   simple_gauss2(mat);
722  tgb_sparse_matrix* mat=build_sparse_matrix(p,p_index,done, done_index,c);
723  simple_gauss(mat,c);
724  m_size=mat->get_rows();
725  m=(poly*) omalloc(m_size*sizeof(poly));
726  m_index=retranslate(m,mat,done,c);
727 
728  mat=NULL;
729  for(i=0;i<done_index;i++)
730    pDelete(&done[i]);
731  omfree(done);
732  done=NULL;
733  //next Step addElements to basis
734  int F_plus_size=m_index;
735  poly* F_plus=(poly*)omalloc(F_plus_size*sizeof(poly));
736  int F_plus_index=0;
737  int F_minus_size=m_index;
738  poly* F_minus=(poly*) omalloc(F_minus_size*sizeof(poly));
739  int F_minus_index=0;
740
741  //better algorithm replace p by its monoms, qsort,delete duplicates and binary search for testing if monom is contained in array
742  qsort(p, p_index,sizeof(poly),pLmCmp_func);
743  for(i=0;i<p_index;i++)
744    pDelete(&p[i]->next);
745  for(i=m_index-1;i>=0;--i)
746{
747  int j;
748  int pos=posInPolys (p,p_index, m[i],c);
749  BOOLEAN minus=FALSE;
750  if(((p_index>pos)&&(pLmEqual(m[i],p[pos]))) ||(pos>0) &&(pLmEqual(m[i],p[pos-1])))
751    minus=TRUE;
752   
753  if(minus)
754  {
755    F_minus[F_minus_index++]=m[i];
756    m[i]=NULL;
757  }
758  else
759  {
760    F_plus[F_plus_index++]=m[i];
761    m[i]=NULL;
762  }
763}
764
765//   for(i=m_index-1;i>=0;--i)
766//   {
767//     int j;
768//     BOOLEAN minus=FALSE;
769//     for(j=0;j<p_index;j++)
770//       if (pLmEqual(p[j],m[i]))
771//       {
772//      minus=TRUE;
773//      break;
774//       }
775//     if(minus)
776//     {
777//       F_minus[F_minus_index++]=m[i];
778//       m[i]=NULL;
779//     }
780//     else
781//     {
782//       F_plus[F_plus_index++]=m[i];
783//       m[i]=NULL;
784//     }
785//   }
786  //in this order F_minus will be automatically ascending sorted
787  //to make this sure for foreign gauss
788  //uncomment the following line
789  //  qsort(F_minus, F_minus_index,sizeof(poly),pLmCmp_func);
790  assume((F_plus_index+F_minus_index)==m_index);
791  if (TEST_OPT_PROT)
792    Print("%i]", F_plus_index);
793  for(i=0;i<p_index;i++) 
794    pDelete(&p[i]);
795  omfree(p);
796  p=NULL;
797  omfree(m);
798  m=NULL;
799
800  //the F_minus list must be cleared separately at the end
801  mp_array_list** F_i;
802  poly_array_list** F_m_i;
803  F_i=&(c->F);
804  F_m_i=&(c->F_minus);
805
806  while((*F_i)!=NULL)
807{
808  assume((*F_m_i)!=NULL);
809  F_i=(&((*F_i)->next));
810  F_m_i=(&((*F_m_i)->next));
811}
812  assume((*F_m_i)==NULL);
813  //should resize the array to save memory
814  //F and F_minus
815  qsort(chosen,chosen_index,sizeof(monom_poly),monom_poly_crit);//important for simplify
816  (*F_m_i)=(poly_array_list*) omalloc(sizeof(poly_array_list));
817  (*F_m_i)->size=F_minus_index;
818  (*F_m_i)->p=F_minus;
819  (*F_m_i)->next=NULL;
820  (*F_i)=(mp_array_list*) omalloc(sizeof(poly_array_list));
821  (*F_i)->size=chosen_index;
822  (*F_i)->mp=chosen;
823  (*F_i)->next=NULL;
824 
825  if(F_plus_index>0)
826{
827  int j;
828  int* ibuf=(int*) omalloc(F_plus_index*sizeof(int));
829  sorted_pair_node*** sbuf=(sorted_pair_node***) omalloc(F_plus_index*sizeof(sorted_pair_node**));
830 
831  for(j=0;j<F_plus_index;j++)
832  {
833    int len;
834    poly p=F_plus[j];
835   
836    // delete buf[j];
837    //remember to free res here
838    //    p=redTailShort(p, c->strat);
839    sbuf[j]=add_to_basis_ideal_quotient(p,c,ibuf+j);
840   
841  }
842  int sum=0;
843  for(j=0;j<F_plus_index;j++)
844  {
845    sum+=ibuf[j];
846  }
847  sorted_pair_node** big_sbuf=(sorted_pair_node**) omalloc(sum*sizeof(sorted_pair_node*));
848  int partsum=0;
849  for(j=0;j<F_plus_index;j++)
850  {
851    memmove(big_sbuf+partsum, sbuf[j],ibuf[j]*sizeof(sorted_pair_node*));
852    omfree(sbuf[j]);
853    partsum+=ibuf[j];
854  }
855
856  qsort(big_sbuf,sum,sizeof(sorted_pair_node*),tgb_pair_better_gen2);
857  c->apairs=spn_merge(c->apairs,c->pair_top+1,big_sbuf,sum,c);
858  c->pair_top+=sum;
859  clean_top_of_pair_list(c);
860  omfree(big_sbuf);
861  omfree(sbuf);
862  omfree(ibuf);
863}
864  omfree(F_plus);
865#ifdef TGB_DEBUG
866  int z;
867  for(z=1;z<=c->pair_top;z++)
868{
869  assume(pair_better(c->apairs[z],c->apairs[z-1],c));
870}
871#endif
872
873  return;
874}
Note: See TracBrowser for help on using the repository browser.