source: git/kernel/F4.cc @ fbc7cb

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