Borderbasix

corealgo3sdpa.hpp
Go to the documentation of this file.
1 
2 #ifndef ALREADY_corealgo3
3 #define ALREADY_corealgo3
4 #include<algorithm>
5 #include<iostream>
6 #include<string>
7 #include<fstream>
8 #include<iomanip>
9 #include<iterator>
10 #include <borderbasix/mpoly.hpp>
11 
12 
13 #include "linalg/clapack.h"
14 
15 namespace coinor {
16 // extern "C" {
17 #define NOSHORTS
18 using namespace std;
19 
20 #include "sdpinclude.hpp"
21 //#include "solvermosek.hpp"
22 //#include "declarations.hpp"
23 //}
24 }
25 
26 template<typename typP,typename typPk,typename Base>
27 void Crochetchoix(const typP & P,typPk & Pk, typPk & redspol,int k,
28  Base &b);
29 
30 template<typename typPk>
31 void stab(typPk & Pk)
32 {
33  typedef typename typPk::value_type pol;
34  for(typename typPk::iterator iter=Pk.begin();iter!=Pk.end();iter++)
35  {
36  pol tmpres=mon(0);
37  for(typename pol::iterator iterpol=iter->begin();iterpol!=iter->end();
38  iterpol++)
39  {
40  if(!Iszero(iterpol->GetCoeff()))
41  tmpres+=*iterpol;
42  }
43  *iter=tmpres;
44  }
45 }
46 
47 template<typename T>
48 void my_merge(T &l1,const T &l2)
49 {
50  typename T::const_iterator iter=l2.begin();
51  for(;iter!=l2.end();iter++)
52  l1.push_back(*iter);
53 }
54 template<typename mon>
55 mon lcm(const mon &m1,const mon &m2)
56 {
57  mon m(1);
58  int maxi=max(lvar(m1.rep),lvar(m2.rep));
59  for(int i=0;i<=maxi;i++)
60  m*=mon(i,max(m1.GetDegree(i),m2.GetDegree(i)));
61  return m;
62 }
63 
64 static int max_degppcm;
65 
66 template<typename typppcm,typename Base>
67 int compute_ppcm(typppcm &ppcm,const Base & b)
68 {
69  typedef typename Base::mon mon;
70  //typppcm tmpppcm;
71  max_degppcm=0;
72  int maxdeg=-1;
73  list<mon> list1,list2,tmplist2;
74 
75 
76  for(typename Base::const_reverse_iterator iter=b.rbegin()
77  ;iter==b.rbegin();iter++)
78  {
79 #ifndef oldppcm
80  for(int i=0;i<iter->taille2;i++)
81  for(int j=i+1;j<iter->taille2;j++)
82  {
83  mon tmpmon=lcm(iter->refuse[i],iter->refuse[j]);
84 
85  if((!ppcm.count(tmpmon))&&
86  (tmpmon.GetDegree()!=
87  (iter->refuse[i].GetDegree()+iter->refuse[j].GetDegree())))
88  {
89  ppcm[tmpmon]=1;
90  if(tmpmon.GetDegree()>max_degppcm)
91  max_degppcm=tmpmon.GetDegree();
92  }
93 
94  // mdebug()<<"ppcm rajoute "<<tmpmon<<" "<<iter->refuse[i]<<" "<<iter->refuse[j];
95  }
96 #else
97  for(int i=0;i<iter->taille2;i++)
98  for(int j=i+1;j<iter->taille2;j++)
99  {
100  mon tmpmon=lcm(iter->refuse[i],iter->refuse[j]);
101 
102  if((tmpmon.GetDegree()!=
103  (iter->refuse[i].GetDegree()+iter->refuse[j].GetDegree())))
104  {
105  int ok=1;
106  for(int k=0;k<iter->taille2;k++)
107  {
108  if((k!=i) && (k!=j)
109  && (isdivisible(iter->refuse[k],tmpmon,b)))
110  {
111  mon lc_i=lcm(iter->refuse[k],iter->refuse[i]);
112  mon lc_j=lcm(iter->refuse[k],iter->refuse[j]);
113  if((lc_i.rep!=tmpmon.rep)
114  && (isdivisible(lc_i,tmpmon,b)) &&
115  (lc_j.rep!=tmpmon.rep)
116  &&(isdivisible(lc_j,tmpmon,b)))
117  ok=0;
118  }
119  }
120  if(ok)
121  list2.push_back(tmpmon);
122  }
123  }
124 #endif
125  }
126 #ifdef oldppcm
127  maxdeg=0;
128  // mdebug()<<"liste2.size "<<list2.size();
129  for(typename list<mon>::iterator iter2=list2.begin();
130  iter2!=list2.end();iter2++)
131  {
132  // mdebug()<<"je passe ici "<<*iter2;
133  ppcm[*iter2]=1;
134  if(iter2->GetDegree()>maxdeg)
135  maxdeg=iter2->GetDegree();
136  }
137  max_degppcm=maxdeg;
138 #endif
139  return maxdeg;
140 }
141 
142 int GB_stop(int k)
143 {
144  int res=0;
145  if(k>max_degppcm) res=1;
146  return res;
147 }
148 
149 template<typename mon,typename Base,typename typppcm,typename typMk>
150 int goodspol(const mon &m1,const mon &m2, typppcm & ppcm,const Base &b
151  ,const typMk &Spoldejavu)
152 {
153  //list<mon> ppcm;
154  int res=0,comp=0;//,nobplus=1;
155  mon m=m1;
156  for(int i=0;i<b.nbvar();i++)
157  if(m1.GetDegree(i)!=m2.GetDegree(i))
158  {
159  comp++;
160  if (m1.GetDegree(i)-m2.GetDegree(i)==1)
161  res+=1;
162  if (m1.GetDegree(i)-m2.GetDegree(i)==-1)
163  {
164  res+=1;
165  m*=mon(i,1);
166  }
167  };
168  //compute_ppcm(ppcm,b);
169  //if((res==2)&&(comp==2))
170  // for(int i=0;i<b.nbvar();i++)
171  // if(m.GetDegree(i)>0)
172  // nobplus=nobplus && !IsinB(my_divmon(m,i,b),b);
173  //mdebug()<<"m1 "<<m1<<" m2 "<<m2<<" resultat "<<((res==2)&&(comp==2));
174  //return ((res==2)&&(comp==2)&&nobplus);
175  //mdebug()<<m1<<" < "<<m2<<" res "<<(member(ppcm,m)&&(res==2)&&(comp==2));
176 
177  // mdebug()<<"m "<<m;
178  //il y a un pbm a ne pas refaire les spol deja vu je ne comprends
179  //pas pourquoi voir example caprass
180 #ifdef NO_TWICE
181  return (res==2)&&(comp==2)&&ppcm[m]&&!member(Spoldejavu,m);
182 #else
183  return (res==2)&&(comp==2)&&ppcm.count(m);
184 #endif
185 }
186 
187 template<typename typPk,typename Base,typename typMk>
188 typPk Spolordsup(const typPk &Pk,const Base & b, const typMk & Spoldejavu)
189 {
190  typedef typename typPk::value_type::monom_t mon;
191  //typedef typename typ::order_t O;
192  map<mon,int> ppcm;
193  typPk res;
194  int maxdeg;
195  maxdeg=compute_ppcm(ppcm,b);
196  // mdebug()<<"maxdeg"<<maxdeg;
197  //mdebug()<<"test goodspol";
198  for(typename typPk::const_iterator iter1=Pk.begin();iter1!=Pk.end();iter1++)
199 
200  for(typename typPk::const_iterator iter2=iter1;iter2!=Pk.end();iter2++)
201  { /*mdebug()<<"iter1->ind"<<iter1->ind;
202  mdebug()<<"iter2->ind"<<iter2->ind;
203  mdebug()<<"goodspol"<<goodspol(iter1->ind,iter2->ind,ppcm,b,Spoldejavu);*/
204  if (goodspol(iter1->ind,iter2->ind,ppcm,b,Spoldejavu))
205  {
206  // mdebug()<<"iter1"<<*iter1;
207  // mdebug()<<"iter2"<<*iter2;
208  res.push_back(*iter1);
209  res.push_back(*iter2);
210  }
211  }
212  return res;
213 }
214 template<typename typP>
215 int minmaxdegP(typP P,int &maxdeg)
216 {
217  int res=Degree(P.front());
218  maxdeg=res;
219  for(typename typP::iterator iter=P.begin();iter!=P.end();iter++)
220  {
221  res=(res>Degree(*iter))?Degree(*iter):res;
222  maxdeg=(maxdeg<Degree(*iter))?Degree(*iter):maxdeg;
223  }
224  return res;
225 }
226 template<typename typP>
227 typP Crochet(const typP &P, int k)
228 {
229  typP res;
230  for(typename typP::const_iterator iter=P.begin();iter!=P.end();iter++)
231  {
232  if(Degree(*iter)==k) res.push_back(*iter);
233  }
234  return res;
235 }
236 
237 template<typename mon, typename T, typename monomial_server>
238 std::list<pol<mon,T> > Crochet(const std::list<pol<mon,T> > &P, int k,
239  monomial_server &serv)
240 {
241  std::list<pol<mon,T> > res;
242  for(typename std::list<pol<mon,T> >::const_iterator iter=P.begin()
243  ;iter!=P.end();iter++)
244  {
245  if(Degree(*iter,serv)==k) res.push_back(*iter);
246  }
247  return res;
248 }
249 
250 template<typename typP, typename Base>
251 void Setnbvars(const typP &tmpPk,Base& b)
252 {
253  int tmpvar=0;
254  for(typename typP::const_iterator iter=tmpPk.begin();iter!=tmpPk.end();iter++)
255  for(typename typP::value_type::const_iterator iterpol=iter->begin();
256  iterpol!=iter->end();iterpol++)
257  {
258  // mdebug()<<"iterpol->rep"<<iterpol->rep;
259  // mdebug()<<"lvar(iterpol->rep)"<<lvar(iterpol->rep);
260  if(tmpvar<lvar(iterpol->rep)) tmpvar=lvar(iterpol->rep);
261  }
262  b.nbvars=tmpvar+1;
263 }
264 
265 template<typename typP>
266 int Setnvars(const typP &tmpPk)
267 {
268  int tmpvar=0;
269  for(typename typP::const_iterator iter=tmpPk.begin();iter!=tmpPk.end();iter++)
270  for(typename typP::value_type::const_iterator iterpol=iter->begin();
271  iterpol!=iter->end();iterpol++)
272  {
273  // mdebug()<<"iterpol->rep"<<iterpol->rep;
274  //mdebug()<<"lvar(iterpol->rep)"<<lvar(iterpol->rep);
275  if(tmpvar<lvar(iterpol->rep)) tmpvar=lvar(iterpol->rep);
276  }
277  return tmpvar+1;
278 
279 }
280 
281 template<typename typP, typename typMk, typename Base,
282  typename monomial_server>
283 void DoBaseMk(const typP &tmpPk, typMk &Mk, Base &b, monomial_server &serv)
284 {
285  typedef typename typMk::value_type mon;
286  typename Base::predicat_t tmp;
287  int i=0;
288  //int tmpvar=0;
289 
290  for(typename typP::const_iterator iter=tmpPk.begin();iter!=tmpPk.end();iter++)
291  {
292  mon tmp(choice(*iter,b,serv));
293  int yes=1;
294  // mdebug()<<"pol dans DoBaseMk "<<endl<<*iter<<endl<<"g choisi "<<tmp;
295  tmp.SetCoeff(1);//*=mon(1/tmp.GetCoeff());
296  for(typename typMk::iterator iter=Mk.begin();iter!=Mk.end();iter++)
297  if(*iter==tmp) {yes=0;break;};
298  if(yes) Mk.push_back(tmp);
299  }
300  //mdebug()<<"Mk.size()"<<Mk.size();
301  tmp.taille2=Mk.size();
302  tmp.refuse=(mon*)MAC_REV_MALLOC<mon>(tmp.taille2*sizeof(mon));
303  tmp.taille1=1;
304  tmp.accept=(mon*)MAC_REV_MALLOC<mon>(sizeof(mon));
305  // mdebug()<<"mon(1) "<<mon(1);
306  tmp.accept[0]=mon(1);
307  //mdebug()<<"tmp.accept "<<tmp.accept[0];
308  i=0;//pour en etre bien sur
309  for(typename typMk::iterator iter=Mk.begin();iter!=Mk.end();iter++,i++)
310  {
311  tmp.refuse[i]=*iter;
312  };
313  if(!b.def.empty())
314  {
315  typename Base::predicat_t verytmp;
316  for(typename Base::iterator iter=b.begin();iter!=b.end();iter++)
317  {
318  verytmp=*iter;
319  MAC_REV_FREE<mon>(verytmp.refuse,verytmp.taille2*sizeof(mon));
320  MAC_REV_FREE<mon>(verytmp.accept,verytmp.taille1*sizeof(mon));
321  }
322  b.erase(b.begin(),b.end());
323  }
324  b.def.push_back(tmp);
325  return;
326 }
327 
328 template <typename typPk,typename typMk,typename typmat >
329 typmat*
330 MatrixOf(const typPk & L,
331  const typMk & lm,const typmat &mat)
332 {
333  typedef typename typPk::value_type POL;
334  typedef typename typMk::value_type monom_t;
335  typmat *res=NULL;
336  map<monom_t,int> index;
337  int l=1;
338  if (res!=NULL) {delete res;};
339  for(typename typMk::const_iterator m = lm.begin(); m!=lm.end(); m++){
340  // mdebug()<<"*m "<<*m;
341  index[*m]=l;
342  l++;
343  }
344  int c;
345  res=new typmat(L.size(),lm.size());
346  l=0;
347 
348  //mdebug()<<"L.size() "<<L.size()<<" lm.size() "<<lm.size();
349  for(typename typPk::const_iterator pi = L.begin(); pi != L.end(); ++pi){
350  for(typename POL::const_iterator m = pi->begin(); m!=pi->end(); m++)
351  {
352 
353  if((index.count(*m)))
354  {
355  c=index[*m];
356  (*res)(l,c-1,m->GetCoeff());
357  }
358  }
359  l++;
360  }
361  //mdebug()<<"fin";
362  return res;
363 }
364 
365 template<typename typMat,typename typP,typename typPk,typename typMk, typename Base,typename typdump, class monomial_server, typename workspace>
366 void init_algo(typP &P,int &nouveaumon, int &maxdeg, int & k,
367  typPk &Pk, typMk & Mk, Base & b, typdump & dump, workspace &w,
368  monomial_server &serv, Base & b2)
369 {
370  typedef typename typPk::value_type pol;
371  typedef typename pol::coeff_t coeff;
372  typedef typename pol::monom_t mon;
373  const typename typP::value_type::order_t (tmpord)(const mon,const mon);
374  typP tmpPk=P;
375  pol debugpol;
376  typMat *Mat;
377  int *pr,*pc,flag=1;
378  serv.setorder(my_ord);
379  nouveaumon=1;
380  Setnbvars(P,b);
381  Setnbvars(P,b2);
382 
383 
384  //mdebug()<<"Init Space";
385  init_space(&w);
386  //mdebug()<<"Init Space fini";
387 #ifdef unimul
388  serv.init_multiple(b.nbvars);
389 #endif
390 #ifdef uninf
391  serv.init_nf(b.nbvars);
392 #endif
393  flag=1;
394  while(flag)
395  {
396  mdebug()<<"Dans While";
397  Mk.erase(Mk.begin(),Mk.end());
398  mdebug()<<"Dans While 1";
399  k=minmaxdegP(tmpPk,maxdeg);
400  mdebug()<<"Dans While 2"<<"k"<<k<<"hola";
401  tmpPk=Crochet(tmpPk,k);
402  // mdebug()<<"Dans While 3";
403  DoBaseMk(tmpPk,Mk,b,serv);
404  DoBaseMk(tmpPk,Mk,b2,serv);
405  // mdebug()<<"Dans While 4";
406  Mat=MatrixOf(tmpPk,Mk,typMat());
407  // mdebug()<<"Dans While 5";
408  tmpPk=Solve(*Mat,tmpPk,&pr,&pc);
409  // mdebug()<<"Dans While 6";
410  selectnozero(tmpPk);
411  stab(tmpPk);
412  // mdebug()<<"my_rank(*Mat)"<<my_rank(*Mat);
413  // mdebug()<<"tmpPk.size()"<<tmpPk.size();
414  if(my_rank(*Mat)==(int)tmpPk.size()) flag=0;
415  Mat->destroystore();
416  delete(Mat);
417  free(pr);
418  free(pc);
419  }
420  // mdebug()<<"While fini";
421  for(typename typP::iterator iter=tmpPk.begin();iter!=tmpPk.end();iter++)
422  {
423  Pk.push_back(convert<pol>(*iter,b,serv));
424  mdebug()<<"tmpPk"<<*iter;
425  }
426 
427  Dump(Pk,dump);
428 
429 
430  return ;
431 }
432 
433 template<typename typMat,typename typP,typename typPk,typename typMk, typename Base,typename typdump, class monomial_server, typename workspace>
434 void init_algo1(typP &P,int &nouveaumon, int &maxdeg, int & k,
435  typPk &Pk, typMk & Mk, Base & b, typdump & dump, workspace &w,
436  monomial_server &serv)
437 {
438  typedef typename typPk::value_type pol;
439  typedef typename pol::coeff_t coeff;
440  typedef typename pol::monom_t mon;
441  const typename typP::value_type::order_t (tmpord)(const mon,const mon);
442  typP tmpPk=P;
443  pol debugpol;
444  typMat *Mat;
445  int *pr,*pc,flag=1;
446  serv.setorder(my_ord);
447  nouveaumon=1;
448  Setnbvars(P,b);
449 
450 
451  mdebug()<<"Init Space";
452  init_space(&w);
453  mdebug()<<"Init Space fini";
454 #ifdef unimul
455  serv.init_multiple(b.nbvars);
456 #endif
457 #ifdef uninf
458  serv.init_nf(b.nbvars);
459 #endif
460  flag=1;
461  while(flag)
462  {
463  mdebug()<<"Dans While";
464  Mk.erase(Mk.begin(),Mk.end());
465  mdebug()<<"Dans While 1";
466  k=minmaxdegP(tmpPk,maxdeg);
467  mdebug()<<"Dans While 2";
468  tmpPk=Crochet(tmpPk,k);
469  mdebug()<<"Dans While 3";
470  DoBaseMk(tmpPk,Mk,b,serv);
471  mdebug()<<"Dans While 4";
472  Mat=MatrixOf(tmpPk,Mk,typMat());
473  mdebug()<<"Dans While 5";
474  tmpPk=Solve(*Mat,tmpPk,&pr,&pc);
475  mdebug()<<"Dans While 6";
476  selectnozero(tmpPk);
477  stab(tmpPk);
478  mdebug()<<"my_rank(*Mat)"<<my_rank(*Mat);
479  mdebug()<<"tmpPk.size()"<<tmpPk.size();
480  if(my_rank(*Mat)==(int)tmpPk.size()) flag=0;
481  Mat->destroystore();
482  delete(Mat);
483  free(pr);
484  free(pc);
485  }
486  mdebug()<<"While fini";
487  for(typename typP::iterator iter=tmpPk.begin();iter!=tmpPk.end();iter++)
488  {
489  Pk.push_back(convert<pol>(*iter,b,serv));
490  mdebug()<<"tmpPk"<<*iter;
491  }
492 
493  Dump(Pk,dump);
494 
495 
496  return ;
497 }
498 
499 template<typename pol,typename coeff, typename workspace>
500 void inline my_plus(pol &res,const pol & nf,const coeff toto,workspace *w)
501 {
502  typedef typename pol::monom_t mon;
503  typedef typename mon::coeff_t mcoeff_t;
504  //symbolic size
505  int i=0,j=0,k=0,p=0,q=0;
506  int tmpsize=0,tmp;
507  mcoeff_t *tmpnf,*tmpspace,zero=0;
508  unsigned char *tmpnfind;
509 #ifdef DEB
510  mdebug()<<"dans myplus je calcul"<<endl<<invconv<Poly>(res)<<endl<<
511  "toto "<<toto<<endl<<"nf"<<endl<<invconv<Poly>(nf);
512 #endif
513  if(res.sizenf<nf.sizenf)
514  {
515  // mdebug()<<"res.sizenf<nf.sizenf ";
516  if(w->sizenf<nf.sizenf)
517  increase_space(w);
518 
519  //#pragma omp parallel for shared(w,toto)
520  for(i=0;i<res.sizenf;i++)
521  {
522  if(res.nfind[i]|nf.nfind[i])
523  {
524  w->tabnf[i]=res.nfind[i]|nf.nfind[i];
525  for(int j=0;j<8;j++)
526  {
527  if ((res.nfind[i]>>j)&1)
528  w->tabcoeff[8*i+j]=res.nf[p++];
529  if((nf.nfind[i]>>j)&1)
530  w->tabcoeff[8*i+j]-=toto*nf.nf[q++];
531  if(Iszero(w->tabcoeff[8*i+j]))
532  w->tabnf[i]&=~(1<<j);
533  }
534  }
535  }
536  //#pragma omp parallel for shared(w,toto)
537  for(i=res.sizenf;i<nf.sizenf;i++)
538  {
539  if(nf.nfind[i])
540  for(int j=0;j<8;j++)
541  if ((nf.nfind[i]>>j)&1)
542  if (!Iszero(nf.nf[q++]))
543  w->tabcoeff[8*i+j]=((mcoeff_t)-1)*toto*(nf.nf[q-1]);
544  else
545  w->tabnf[i]&=~(1<<j);
546  w->tabnf[i]=nf.nfind[i];
547  }
548  for(int i=0;i<nf.sizenf;i++)
549  tmpsize+=nbbits[w->tabnf[i]];
550  res.nf=(COEFF*)MAC_REV_REALLOC<COEFF>(res.nf,res.size*sizeof(COEFF)
551  ,tmpsize*sizeof(COEFF));
552  res.size=tmpsize;
553  res.nfind=(unsigned char*)realloc(res.nfind,nf.sizenf);
554  res.sizenf=nf.sizenf;
555 
556  for(int i=0;i<nf.sizenf;i++)
557  {
558  if(w->tabnf[i])
559  for(int j=0;j<8;j++)
560  if ((w->tabnf[i]>>j)&1)
561  //if (!Iszero(w->tabcoeff[8*i+j]))
562  {
563  res.nf[k++]=w->tabcoeff[8*i+j];
564  w->tabcoeff[8*i+j]=zero;
565  }
566  //else
567  // w->tabnf[i]&=~(1<<j);
568  res.nfind[i]=w->tabnf[i];
569  w->tabnf[i]=0;
570  }
571  }
572  else
573  {
574  if(w->sizenf<res.sizenf)
575  increase_space(w);
576  //#pragma omp parallel for shared(w,toto)
577  for(i=0;i<nf.sizenf;i++)
578 #if 1
579  {
580  if(res.nfind[i])
581  {
582  unsigned char mask=1;
583  w->tabnf[i]=res.nfind[i];
584  if(res.nfind[i]==255)
585  for(int j=0;j<8;j++)
586  w->tabcoeff[8*i+j]=res.nf[p++];
587  else
588  for(int j=0;j<8;j++,mask<<=1)
589  {
590  //if ((res.nfind[i]>>j)&1)
591  if (res.nfind[i]&mask)
592  w->tabcoeff[8*i+j]=res.nf[p++];
593  }
594  }
595  if(nf.nfind[i])
596  {
597  unsigned char mask=1;
598  w->tabnf[i]|=nf.nfind[i];
599  if(nf.nfind[i]==255)
600  for(int j=0;j<8;j++,mask<<=1)
601  {
602  w->tabcoeff[8*i+j]-=toto*nf.nf[q++];
603  if(Iszero(w->tabcoeff[8*i+j]))
604  w->tabnf[i]&=~mask;
605  }
606  else
607  for(int j=0;j<8;j++,mask<<=1)
608  {
609  if(nf.nfind[i]&mask)
610  w->tabcoeff[8*i+j]-=toto*nf.nf[q++];
611  if(Iszero(w->tabcoeff[8*i+j]))
612  w->tabnf[i]&=~mask;
613  }
614  }
615  }
616 #endif
617 #if 0
618  if(res.nfind[i]|nf.nfind[i])
619  {
620  w->tabnf[i]=res.nfind[i]|nf.nfind[i];
621  for(int j=0;j<8;j++)
622  {
623  if ((res.nfind[i]>>j)&1)
624  {
625  w->tabcoeff[8*i+j]=res.nf[p++];
626  }
627  if((nf.nfind[i]>>j)&1)
628  w->tabcoeff[8*i+j]-=toto*nf.nf[q++];
629  if(Iszero(w->tabcoeff[8*i+j]))
630  w->tabnf[i]&=~(1<<j);
631  }
632  }
633 #endif
634  //#pragma omp parallel for shared(w,toto)
635  for(i=nf.sizenf;i<res.sizenf;i++)
636  {
637  if(res.nfind[i])
638  for(int j=0;j<8;j++)
639  if ((res.nfind[i]>>j)&1)
640  if (!Iszero(res.nf[p++]))
641  w->tabcoeff[8*i+j]=res.nf[p-1];
642  else
643  w->tabnf[i]&=~(1<<j);
644  w->tabnf[i]=res.nfind[i];
645  }
646  for(int i=0;i<res.sizenf;i++)
647  tmpsize+=nbbits[w->tabnf[i]];
648  res.nf=(COEFF*)MAC_REV_REALLOC<COEFF>(res.nf,res.size*sizeof(COEFF)
649  ,tmpsize*sizeof(COEFF));
650  res.size=tmpsize;
651  // res.nfind=(unsigned char*)realloc(res.nfind,nf.sizenf);
652  //res.sizenf=nf.sizenf;
653  for(int i=0;i<res.sizenf;i++)
654  {
655  if(w->tabnf[i])
656  for(int j=0;j<8;j++)
657  if ((w->tabnf[i]>>j)&1)
658  //if (!Iszero(w->tabcoeff[8*i+j]))
659  {
660  res.nf[k++]=w->tabcoeff[8*i+j];
661  w->tabcoeff[8*i+j]=zero;
662  }
663  //else
664  // w->tabnf[i]&=~(1<<j);
665 
666  res.nfind[i]=w->tabnf[i];
667  w->tabnf[i]=0;
668  }
669  }
670 
671 #ifdef DEB
672  mdebug()<<"k "<<k<<" tmpsize "<<tmpsize;
673  int cpt=0;
674  for(int i=0;i<res.sizenf;i++)
675  if(res.nfind[i])
676  for(int j=0;j<8;j++)
677  if ((res.nfind[i]>>j)&1)
678  if (Iszero(res.nf[cpt++]))
679  mdebug()<<"j'ai fait n'importe quoi ici";
680  mdebug()<<"les resultat "<<endl<<invconv<Poly>(res);
681 #endif
682  return;
683 }
684 template<typename pol,typename coeff>
685 void inline my_plus(pol &res,const pol & nf,const coeff toto,int removeme)
686 {
687  typedef typename pol::monom_t mon;
688  typedef typename mon::coeff_t mcoeff_t;
689  //symbolic size
690  int i=0,j=0,k=0,p=0,q=0;
691  int tmpsize=0,tmp;
692  mcoeff_t *tmpnf,*tmpspace;
693  unsigned char *tmpnfind;
694 #ifdef DEB
695  mdebug()<<"dans myplus je calcul"<<endl<<invconv<Poly>(res)<<endl<<
696  "toto "<<toto<<endl<<"nf"<<endl<<invconv<Poly>(nf);
697 #endif
698  if(res.sizenf<nf.sizenf)
699  {
700  tmpspace=(mcoeff_t*)
701  MAC_REV_MALLOC<mcoeff_t>(sizeof(mcoeff_t)*nf.sizenf*8);
702  for(int i=0;i<nf.sizenf*8;i++)
703  tmpspace[i]=0;
704  for(int i=0;i<res.sizenf;i++)
705  tmpsize+=nbbits[res.nfind[i]|nf.nfind[i]];
706  for(int i=0;i<nf.sizenf;i++)
707  tmpsize+=nbbits[nf.nfind[i]];
708  tmpnf=(mcoeff_t*)
709  MAC_REV_MALLOC<mcoeff_t>(sizeof(mcoeff_t)*tmpsize);
710  tmpnfind=(unsigned char*)
711  MAC_REV_MALLOC<unsigned char>(nf.sizenf);
712  // memset(tmpnfind,0,nf.sizenf);
713 
714  for(int i=0;i<res.sizenf;i++)
715  {
716  if(res.nfind[i]|nf.nfind[i])
717  for(int j=0;j<8;j++)
718  {
719  if ((res.nfind[i]>>j)&1)
720  tmpspace[8*i+j]=res.nf[p++];
721  if((nf.nfind[i]>>j)&1)
722  tmpspace[8*i+j]-=toto*nf.nf[q++];
723  }
724  tmpnfind[i]=res.nfind[i]|nf.nfind[i];
725  }
726 
727  for(int i=res.sizenf;i<nf.sizenf;i++)
728  {
729  if(nf.nfind[i])
730  for(int j=0;j<8;j++)
731  if ((nf.nfind[i]>>j)&1)
732  tmpspace[8*i+j]=((mcoeff_t)-1)*toto*(nf.nf[q++]);
733  tmpnfind[i]=nf.nfind[i];
734  }
735  for(int i=0;i<nf.sizenf;i++)
736  if(tmpnfind[i])
737  for(int j=0;j<8;j++)
738  if ((tmpnfind[i]>>j)&1)
739  if (!Iszero(tmpspace[8*i+j]))
740  tmpnf[k++]=tmpspace[8*i+j];
741  else
742  tmpnfind[i]&=~(1<<j);
743  MAC_REV_FREE<mcoeff_t>(tmpspace,nf.sizenf*8*sizeof(mcoeff_t));
744  res.sizenf=nf.sizenf;
745  }
746  else
747  {
748  tmpspace=(mcoeff_t*)
749  MAC_REV_MALLOC<mcoeff_t>(sizeof(mcoeff_t)*res.sizenf*8);
750  for(int i=0;i<res.sizenf*8;i++)
751  tmpspace[i]=0;
752  for(int i=0;i<nf.sizenf;i++)
753  tmpsize+=nbbits[res.nfind[i]|nf.nfind[i]];
754  for(int i=0;i<res.sizenf;i++)
755  tmpsize+=nbbits[res.nfind[i]];
756  tmpnf=(mcoeff_t*)
757  MAC_REV_MALLOC<mcoeff_t>(sizeof(mcoeff_t)*tmpsize);
758  tmpnfind=(unsigned char*)
759  MAC_REV_MALLOC<unsigned char>(res.sizenf);
760 
761  for(int i=0;i<nf.sizenf;i++)
762  {
763  if(res.nfind[i]|nf.nfind[i])
764  for(int j=0;j<8;j++)
765  {
766  if((res.nfind[i]>>j)&1)
767  tmpspace[8*i+j]=res.nf[p++];
768  if((nf.nfind[i]>>j)&1)
769  tmpspace[8*i+j]-=toto*nf.nf[q++];
770  }
771  tmpnfind[i]=res.nfind[i]|nf.nfind[i];
772  }
773 
774  for(int i=nf.sizenf;i<res.sizenf;i++)
775  {
776  if(res.nfind[i])
777  for(int j=0;j<8;j++)
778  if((res.nfind[i]>>j)&1)
779  tmpspace[8*i+j]=res.nf[p++];
780  tmpnfind[i]=res.nfind[i];
781  }
782 
783  for(int i=0;i<res.sizenf;i++)
784  if(tmpnfind[i])
785  for(int j=0;j<8;j++)
786  if ((tmpnfind[i]>>j)&1)
787  if (!Iszero(tmpspace[8*i+j]))
788  tmpnf[k++]=tmpspace[8*i+j];
789  else
790  tmpnfind[i]&=~(1<<j);
791 
792  MAC_REV_FREE<mcoeff_t>(tmpspace,res.sizenf*8*sizeof(mcoeff_t));
793  }
794 #ifdef DEB
795  mdebug()<<"k "<<k<<" tmpsize "<<tmpsize;
796 #endif
797  MAC_REV_FREE<mcoeff_t>(res.nf,res.size*sizeof(mcoeff_t));
798  //res.nf=tmpnf;
799  MAC_REV_FREE<unsigned char>(res.nfind,res.sizenf);
800  //res.nfind=tmpnfind;
801 
802  res.nf=tmpnf;
803  res.nf=(mcoeff_t*)MAC_REV_REALLOC<mcoeff_t>(res.nf,res.size*sizeof(mcoeff_t)
804  ,k*sizeof(mcoeff_t));
805 
806  //mdebug()<<"tmpnfind "<<tmpnfind;
807  res.nfind=tmpnfind;
808  // mdebug()<<"res.nfind "<<res.nfind;
809  res.size=k;
810 #ifdef DEB
811  int cpt=0;
812  for(int i=0;i<res.sizenf;i++)
813  if(tmpnfind[i])
814  for(int j=0;j<8;j++)
815  if ((tmpnfind[i]>>j)&1)
816  if (Iszero(res.nf[cpt++]))
817  mdebug()<<"j'ai fait n'importe quoi ici";
818  mdebug()<<"les resultat "<<endl<<invconv<Poly>(res);
819 #endif
820  return;
821 }
822 
823 //#undef DEB
824 //renvoit l'indice du monome x_i*m_j
825 template<typename mon,typename Base>
826 int indice2(int i,int j,const Base & b)
827 {
828  int res=-1;
829  mon tmp;
830  int2mon(j,tmp);
831  tmp*=mon(i,1);//on suppose qd meme que j correspond a qqchose
832  //mdebug()<<"indice2 IsinB("<<mon(i,1)*tmp<<",b) "<<IsinB(mon(i,1)*tmp,b);
833  if(IsinB(tmp,b))
834  {
835  res=mon2int(tmp);
836  // mdebug()<<"le monome "<<mon(i,1)*tmp<<" l'indice "<<res;
837  }
838  return res;
839 }
840 
841 //on lui donne l'ensemble des formes normales, dump, l'indice de b(cf ci dessus)
842 // le numero de la variable par laquelle on multiplie
843 //et elle renvoie un pol(cf def du type) sur la forme normale de var*m_{indmon}
844 #ifndef uninf
845 template<typename typpol,typename typdump,typename Base>
846 typpol nf(int var ,int indmon,const typdump &dump,const Base & b)
847 {
848  //on cherche dans dump...
849  //a optimiser plus tard
850  typedef typename typdump::value_type dumpstruct;
851  typedef typpol pol;
852  typedef typename pol::monom_t mon;
853  pol res;
854  res.nf=NULL;
855  res.size=-1;
856  mon tmpmon;//=mon(1,var)*int2mon<mon>(indmon);
857  int deg=0;//=tmpmon.GetDegree()+1;
858  typename typdump::const_iterator iter;
859  int2mon(indmon,tmpmon);
860  deg=tmpmon.GetDegree()+1;
861  //mdebug()<<"mon("<<var<<",1) "<<mon(var,1)<<" et deg "<<deg;
862  tmpmon*=mon(var,1);
863  for(iter=dump.begin();
864  (iter!=dump.end()) && (iter->k!=deg);iter++);
865  // mdebug()<<"iter-k "<<iter->k;
866  // if(iter!=dump.end()) mdebug()<<"iter->k en sortie "<<iter->k<<" et deg "<<deg;
867  //mdebug()<<"valuer a la sortie "<<(--iter)++->k;
868  // mdebug()<<"iter!=dump.end() "<<(iter!=dump.end())
869  // <<" dump.empty() "<<dump.empty();
870  if(iter!=dump.end())//on cherche pas n'importe quoi qd meme
871  {
872  int i;
873  for(i=0;i<iter->size;i++)
874  {
875  // mdebug()<<"dans nf iter->nf[i].ind "<<iter->nf[i].ind;
876  if(iter->nf[i].ind.rep==tmpmon.rep)
877  break;
878  }
879  if(i!=iter->size)
880  //Tue Dec 4 16:39:25 MET 2001
881  // densify(iter->nf[i],res);
882 
883  res=(iter->nf[i]);
884  }
885  //mdebug()<<"la gueule de res "<<res;
886  //if (res!=NULL) mdebug()<<"suite de res "<<res->ind<<" "<<res->nf;
887  return res;
888 }
889 #endif
890 
891 void getperm(int *perm, int sizeperm, int *indextab)
892 {
893  int median=indextab[sizeperm/2],i,j=sizeperm-1;
894  if(1)
895  for(i=0;i<=j;i++)
896  if(indextab[i]>median)
897  {
898  int tmp=indextab[j];
899  indextab[j]=indextab[i];
900  indextab[i]=tmp;
901  j--;i--;
902  }
903  perm[sizeperm/2]=sizeperm-i;
904  getperm(perm,sizeperm-i-1,indextab);
905  getperm(perm+sizeperm-i,i,indextab+sizeperm-i);
906 }
907 
908 //effectue la multiplication du pol T par la variable x_i
909 template<typename typpol,typename typdump,typename Base,typename mon,
910  typename workspace, typename monomial_server>
911 typpol inline mult(int i,const typpol &pol,const typdump &dump,
912  const Base &b,mon **exce,int &exsize,workspace &w,
913  monomial_server &serv)
914 {
915  typedef typename typpol::coeff_t coeff;
916  static const mon monzero=mon(0);
917  static COEFF zero=0;
918  typpol ret;
919  (*exce)=(mon*)MAC_REV_MALLOC<mon>((pol.size+1)*sizeof(mon));
920  int *stockindicetab=(int*)malloc(pol.size*sizeof(int));
921  int noquo=0,maxind=0,p=0,q=0;
922  coeff *tmpspace;
923  // for(int i=0;i<w.sizenf;i++)
924  // p+=nbbits[w.tabnf[i]];
925  for(int j=0;j<pol.sizenf;++j)
926  if (pol.nfind[j])
927  {
928  for(int k=0;k<8;k++)
929  if((pol.nfind[j]>>k)&1)
930  {
931  int stockindice2=mulind(8*j+k,i,b);
932  if(stockindice2<0)
933  stockindice2=-1*(8*j+k)-1;
934  stockindicetab[p++]=stockindice2;
935 
936  // mdebug()<<"stockindice2 "<<stockindice2<<" pol.size "<<
937  //pol.size<<" J "<<j;
938  if (stockindice2>=maxind)
939  maxind=stockindice2;
940  }
941  }
942  //mdebug()<<"p "<<p<<" pol.size "<<pol.size;
943  for(;maxind/8+1>=w.sizenf;increase_space(&w));
944  //tmpspace=(coeff*)MAC_REV_MALLOC<coeff>((maxind/8+1)*8*sizeof(coeff));
945  //ret.sizenf=(maxind/8+1);
946  //ret.nfind=(unsigned char*)MAC_REV_MALLOC<unsigned char>(ret.sizenf);
947  //memset(ret.nfind,0,ret.sizenf);
948 
949 
950 
951  //for(int j=0;j<ret.sizenf*8;j++)
952  // tmpspace[j]=0;
953 #ifdef DEB
954  q=0;
955  for(int j=0;j<p;j++)
956  q+=nbbits[pol.nfind[j]];
957  mdebug()<<" a '"<<q;
958 #endif
959  int j;
960  //#pragma omp parallel for shared(stockindicetab,w,serv)
961  for(j=0;j<p;j++)
962  if(stockindicetab[j]>=0)
963  {
964 
965 #ifdef DEB
966  mon tmpmm;
967  int2mon(stockindicetab[j],tmpmm);
968  mdebug()<<"le monome la "<<tmpmm<<" "<<pol.nf[j];
969 #endif
970  w.tabcoeff[stockindicetab[j]]+=pol.nf[j];
971  if(!Iszero(w.tabcoeff[stockindicetab[j]]))
972  w.tabnf[stockindicetab[j]/8]|=1<<(stockindicetab[j]%8);
973  }
974  else
975  {
976  // mdebug()<<"ici aussi";
977  typpol tmp;
978  tmp=serv.nf(i,-1*(stockindicetab[j]+1),dump,b);
979  if (tmp.size!=-1)
980  {
981  int compteurtmp=0;
982  //On fait un plus dense dans tmpspace
983  for(;tmp.sizenf>=w.sizenf;increase_space(&w));
984  // if(tmp.sizenf>res.sizenf)
985  // {
986  // ret.nfind=(unsigned char*)realloc(ret.nfind,tmp.sizenf);
987  // for(int k=ret.sizenf;k<tmp.sizenf;k++)
988  // ret.nfind[k]=0;
989  // tmpspace=(coeff*)MAC_REV_REALLOC<coeff>(tmpspace,
990  // ret.sizenf*8*sizeof(coeff),
991  // tmp.sizenf*8*sizeof(coeff));
992  // for(int k=ret.sizenf*8;k<tmp.sizenf*8;k++)
993  // tmpspace[k]=0;
994  // ret.sizenf=tmp.sizenf;
995 
996  // //normalement on ne devrait JAMAIS passer par la
997  // }
998 #if 1
999  int k,l;
1000  unsigned *tmpptr=(unsigned*)tmp.nfind;
1001  unsigned mask=1;
1002  for(k=0;k<tmp.sizenf/sizeof(unsigned);k++)
1003  if(tmpptr[k])
1004  {
1005  ((unsigned*)w.tabnf)[k]|=tmpptr[k];
1006  mask=1;
1007  for(l=0;l<8*sizeof(unsigned);l++,mask<<=1)
1008  {
1009  if(tmpptr[k]&mask)
1010  w.tabcoeff[8*sizeof(unsigned)*k+l]
1011  -=pol.nf[j]*tmp.nf[compteurtmp++];
1012  // if(Iszero(w.tabcoeff[8*k*sizeof(unsigned)+l]))
1013  // ((unsigned*)w.tabnf)[k]&=~mask;
1014  }
1015  }
1016  for(k=tmp.sizenf-(tmp.sizenf%sizeof(unsigned))
1017  ;k<tmp.sizenf;k++)
1018  if(tmp.nfind[k])
1019  {
1020  w.tabnf[k]|=tmp.nfind[k];
1021  mask=1;
1022  for(l=0;l<8;l++,mask<<=1)
1023  {
1024  if(tmp.nfind[k]&mask)
1025  w.tabcoeff[8*k+l]-=pol.nf[j]*tmp.nf[compteurtmp++];
1026  //if(Iszero(w.tabcoeff[8*k+l]))
1027  // w.tabnf[k]&=~mask;
1028  }
1029  }
1030 #endif
1031 #if 0
1032  for(int k=0;k<tmp.sizenf;k++)
1033  {
1034  if(tmp.nfind[k])
1035  {
1036  if(tmp.nfind[k]==255)
1037  {
1038  unsigned char mask=1;
1039  w.tabnf[k]|=tmp.nfind[k];
1040  for(int l=0;l<8;l++,mask<<=1)
1041  {
1042  w.tabcoeff[8*k+l]-=
1043  pol.nf[j]*tmp.nf[compteurtmp++];
1044  if(Iszero(w.tabcoeff[8*k+l]))
1045  w.tabnf[k]&=~mask;
1046  }
1047  }
1048  else
1049  {
1050  unsigned char mask=1;
1051  w.tabnf[k]|=tmp.nfind[k];
1052  for(int l=0;l<8;l++,mask<<=1)
1053  if(tmp.nfind[k]&mask)
1054  {
1055  w.tabcoeff[8*k+l]-=
1056  pol.nf[j]*tmp.nf[compteurtmp++];
1057  if(Iszero(w.tabcoeff[8*k+l]))
1058  w.tabnf[k]&=~mask;
1059  }
1060  }
1061  }
1062 
1063  }
1064 #endif
1065  }
1066  else
1067  {
1068  mon tmpmon;
1069  (*exce)[noquo]=mon(i,1);
1070  int2mon(-1*(stockindicetab[j]+1),tmpmon);
1071  (*exce)[noquo]*=tmpmon;//int2mon<mon>(j);
1072  (*exce)[noquo]*=pol.nf[j];
1073  ++noquo;
1074  };
1075  }
1076  ret.ind=pol.ind;
1077  ret.ind*=mon(i,1);
1078  (*exce)[noquo++]=ret.ind;
1079  exsize=noquo;
1080  (*exce)=(mon*)MAC_REV_REALLOC<mon>((*exce),(pol.size+1)*sizeof(mon)
1081  ,(noquo)*sizeof(mon));
1082 #ifdef DEB
1083  for(int j=0;j<exsize;j++)
1084  mdebug()<<" "<<(*exce)[j];
1085  mdebug();
1086 #endif
1087  ret.ind.SetCoeff(noquo+1);
1088  ret.size=0;
1089  p=0;
1090  for(int j=0;j<w.sizenf;j++)
1091  if(w.tabnf[j])
1092  {
1093  ret.size+=nbbits[w.tabnf[j]];
1094  p=j;
1095  }
1096 #ifdef DEB
1097  mdebug()<<" val p "<<p<<" et pol.size "<<pol.size;
1098 #endif
1099  ret.nf=(coeff*)MAC_REV_MALLOC<coeff>(ret.size*sizeof(coeff));
1100  ret.sizenf=p+1;
1101  ret.nfind=(unsigned char *)malloc(ret.sizenf);
1102  p=0;
1103  for(int j=0;j<ret.sizenf;j++)
1104  {
1105  ret.nfind[j]=w.tabnf[j];
1106  if(w.tabnf[j])
1107  for(int k=0;k<8;k++)
1108  if((w.tabnf[j]>>k)&1)
1109  {
1110  if(!Iszero(w.tabcoeff[8*j+k]))
1111  {
1112  ret.nf[p++]=w.tabcoeff[8*j+k];
1113  w.tabcoeff[8*j+k]=zero;
1114  }
1115  else
1116  {
1117  ret.nfind[j]&=~(1<<k);
1118  }
1119  }
1120  w.tabnf[j]=0;
1121  }
1122  ret.nf=(coeff*)MAC_REV_REALLOC<coeff>(ret.nf,ret.size*sizeof(coeff),
1123  p*sizeof(coeff));
1124  ret.size=p;
1125  //MAC_REV_FREE<coeff>(tmpspace,ret.sizenf*8*sizeof(coeff));
1126  free(stockindicetab);
1127 #ifdef DEB
1128  mdebug()<<"ret.size "<<ret.size;
1129  mdebug()<<invconv<Poly>(ret);
1130  mdebug()<<"-----------------------------------";
1131  int cpt=0;
1132  for(int i=0;i<ret.sizenf;i++)
1133  if(ret.nfind[i])
1134  for(int j=0;j<8;j++)
1135  if ((ret.nfind[i]>>j)&1)
1136  if (Iszero(ret.nf[cpt++]))
1137  mdebug()<<"nimpici";
1138 #endif
1139  return ret;
1140 }
1141 
1142 // renvoit mon/x_i .... mais suppose qu'on ne fait pas n'importe quoi
1143 template <typename Mon,typename Base>
1144 Mon my_divmon(const Mon &mon,int i,const Base& b)
1145 {
1146  Mon res(1);
1147  for(int j=0;j<b.nbvar();j++)
1148  if(i==j)
1149  {
1150  res*=Mon(j,mon.GetDegree(j)-1);
1151  }
1152  else
1153  {
1154  res*=Mon(j,mon.GetDegree(j));
1155  }
1156  return res;
1157 }
1158 
1159 template<typename typmon,typename Base>
1160 int GoodDir(const typmon &m,int i, const Base & b)
1161 {
1162  int j,res=0;
1163  mon tmpm=m*mon(i,1);
1164  //mdebug()<<"dans GoodDir les monomes "<<m<<"********"<<typmon(i,1);
1165  //mdebug()<<"b.nbvar() "<<b.nbvar();
1166  //#pragma omp parallel for shared(tmpm,res)
1167  for(j=0;j<b.nbvar();j++)
1168  {
1169  if(tmpm.GetDegree(j)>=1)
1170  if(IsinB(my_divmon(tmpm,j,b),b)) {res|=1;};
1171  }
1172  //mdebug()<<" le verdict de GoodDir "<<res;
1173  return res;
1174 }
1175 
1176 //fonction adhoc qui merge une liste de polyalp avec une liste de pol
1177 //en cours de calcul cad liste de pol + liste excedentaire
1178 template<typename typPk,typename typP,typename typexceed,
1179  typename typsize ,typename Base,typename monomial_server>
1180 void conv_merge(typPk &Pk,typP &P,typexceed &ex,typsize &sizeex,
1181  const Base &b, monomial_server &serv)
1182 {
1183  typedef typename typP::value_type polalp;
1184  typedef typename polalp::order_t ord;
1185  typedef typename polalp::monom_t mon;
1186  typedef typename typPk::value_type pol;
1187  typedef typename pol::coeff_t coeff;
1188  coeff * tabcoeff=NULL;
1189  unsigned char *tabindices=NULL;
1190  int tabsize=0,tabnfsize=0;
1191 
1192  for(typename typP::iterator iter=P.begin();iter!=P.end();iter++)
1193  {
1194  //mdebug()<<"dans for";
1195  int comp=0;int i=0;
1196  pol tmp;
1197  int maxind=0;
1198  tmp.ind=mon(1);
1199  for(typename polalp::iterator iterpol=iter->begin();
1200  iterpol!=iter->end();iterpol++)
1201  {
1202  // mdebug()<<"iterpol"<<*iterpol;
1203  // mdebug()<<"IsinB(*iterpol,b)"<<IsinB(*iterpol,b);
1204  // mdebug()<<"serv.mon2int(*iterpol)"<<serv.mon2int(*iterpol);
1205  if (IsinB(*iterpol,b) && (serv.mon2int(*iterpol)>maxind))
1206  maxind=serv.mon2int(*iterpol);
1207  }
1208  // mdebug()<<"maxind"<<maxind;
1209  tabcoeff=(coeff*)MAC_REV_MALLOC<coeff>((maxind+1)*sizeof(coeff));
1210  tabindices=(unsigned char*)malloc(maxind/8+1);
1211 
1212  tabsize=maxind+1;
1213  tabnfsize=maxind/8+1;
1214  memset((void*)tabindices,0,tabnfsize);
1215  maxind=0;
1216  //ici il faudrait un type template
1217  for(i=0;i<tabsize;i++)
1218  {
1219  tabcoeff[i]=0;
1220  }
1221  //memset((void*)tabcoeff,0,tabsize*sizeof(coeff));
1222  //ici aussi
1223  i=0;
1224  mon * tmpex=NULL;
1225  tmpex=(mon*)MAC_REV_REALLOC<mon>((void*)tmpex,0,iter->size()
1226  *sizeof(mon));
1227 
1228  for(typename polalp::iterator iterpol=iter->begin();
1229  iterpol!=iter->end();iterpol++)
1230  {
1231  //#if 1
1232  // mdebug()<<"IsinB("<<*iterpol<<","<<"b)"<<IsinB(*iterpol,b);
1233  //#endif
1234  if (IsinB(*iterpol,b))
1235  {
1236 #ifdef BONRAJOUT
1237  mdebug()<<"dedans";
1238 #endif
1239  int tmp=serv.mon2int(*iterpol);
1240  tabcoeff[tmp]=iterpol->GetCoeff();
1241  tabindices[tmp/8]|=1<<(tmp%8);
1242  i++;
1243  maxind++;
1244  }
1245  else
1246  {
1247 #if 1
1248  mdebug()<<"y en a un qu'est pas dans B "<<*iterpol;
1249 #endif
1250  if(!Iszero(iterpol->GetCoeff()))
1251  tmpex[comp++]=*iterpol;
1252  }
1253  }
1254  tmpex=(mon*)MAC_REV_REALLOC<mon>((void*)tmpex,iter->size()*sizeof(mon)
1255  ,(comp+1)*sizeof(mon));
1256  tmp.size=maxind;
1257  tmp.nf=(coeff*)MAC_REV_MALLOC<coeff>(tmp.size*sizeof(coeff));
1258  tmp.nfind=tabindices;
1259  tmp.sizenf=tabnfsize;
1260  tmp.ind=mon(comp+1);
1261  maxind=0;
1262  for(int i=0;i<tabnfsize;i++)
1263  if(tabindices[i])
1264  {
1265  for(int j=0;j<8;j++)
1266  if (((tabindices[i]>>j)&1) && !Iszero(tabcoeff[8*i+j]))
1267  tmp.nf[maxind++]=tabcoeff[8*i+j];
1268  else
1269  tmp.nfind[i]&=~(1<<j);
1270  }
1271  tmp.nf=(coeff*)MAC_REV_REALLOC<coeff>(tmp.nf,tmp.size*sizeof(coeff)
1272  ,maxind*sizeof(coeff));
1273  tmp.size=maxind;
1274  MAC_REV_FREE<coeff>(tabcoeff,tabsize*sizeof(coeff));
1275  // mdebug()<<"tmp"<<tmp.nf;
1276  // mdebug()<<"tmpex"<<*tmpex;
1277  // mdebug()<<"comp"<<comp;
1278  Pk.push_back(tmp);
1279  ex.push_back(tmpex);
1280  sizeex.push_back(comp);
1281 #ifdef BONRAJOUT
1282  mdebug()<<"ici dans conv merge";
1283  for(int i=0;i<comp;i++)
1284  mdebug()<<tmpex[i]<<" ";
1285  mdebug();
1286 #endif
1287  }
1288 }
1289 //member ne test que les exposants pas les coefficients
1290 template<typename typMk>
1291 int member(const typMk & tmpMk,const typename typMk::value_type mon)
1292 {
1293  int res=0;
1294  for(typename typMk::const_iterator iter=tmpMk.begin();
1295  iter!=tmpMk.end();iter++)
1296  if(iter->rep==mon.rep) {res=1;break;}
1297  return res;
1298 }
1299 
1300 template<typename typMat,typename typPk,typename typP,
1301  typename typMk,typename Base,typename typdump,typename typrecall,
1302  typename workspace, typename monomial_server>
1303 int newPkMkMatrixof(typMat &mat,typPk &Pk, typP &P,typMk &Mk,
1304  int &k,Base &b, typdump &dump,typrecall &recall
1305  ,workspace &w, monomial_server &serv)
1306 {
1307  typedef typename typMk::value_type mon;
1308  typedef typename typP::value_type::order_t ord;
1309  typedef typename typPk::value_type pol;
1310  typedef typename pol::coeff_t coeff;
1311  typedef typename typP::value_type polyalp;
1312  typPk tmpPk,workvect;
1313  typMk tmpMk;
1314  list<mon *> exceed;
1315  list<int> sizeexceed;
1316  int nbproj=0,flagdoplus=1,*ncol;
1317  //prise en compte des pol pas encore introduits
1318  mdebug()<<"k on entering newpkmkmatrixof "<<k;
1319  typP ttmpP=Crochet(P,k+1);//
1320  typrecall verytmp=Crochet(recall,k+1,serv);
1321  recall.erase(recall.begin(),recall.end());
1322  if (!ttmpP.empty())
1323  {
1324  int tmpdeg;
1325  proj(ttmpP,dump,b,serv);
1326  selectnozero(ttmpP);
1327 #if 1
1328  mdebug()<<"---------";
1329  for(typename typP::iterator iter=ttmpP.begin();iter!=ttmpP.end();iter++)
1330  mdebug()<<*iter;
1331 #endif
1332 #ifdef NO_TWICE_FALL
1333  for(typename typP::iterator iter=P.begin();iter!=P.end();iter++)
1334  if(Degree(*iter)==k+1) P.erase(iter--);
1335  for(typename typP::iterator iter=ttmpP.begin()
1336  ;iter!=ttmpP.end();iter++)
1337  P.push_back(*iter);
1338 #endif
1339  // mdebug()<<"pol mal red "<<*iter;
1340  //mdebug()<<"passage par ici";
1341  nbproj=ttmpP.size();
1342  if (!ttmpP.empty())
1343  {
1344  // mdebug()<<"apres select non zero "<<ttmpP.size();
1345  tmpdeg=minmaxdegP(ttmpP,tmpdeg);
1346  // mdebug()<<"apres minmax deg "<<tmpdeg<<" "<<k;
1347  if(tmpdeg<k+1)
1348  {
1349  flagdoplus=0;
1350  mdebug()<<"un des pol de pretendu degre le bon se reduit mal";
1351 
1352  //en fait a un polynome de degre inferieur
1353  typPk tmpconv;
1354  //mdebug()<<"recover";
1355  recover(tmpPk,Mk,dump,tmpdeg);
1356  //Listedejafait(Mk,Spoldejavu);
1357  for (typename typP::iterator iter=ttmpP.begin();iter!=ttmpP.end();iter++)
1358  mdebug()<<"pol ICI "<<*iter;
1359  // mdebug()<<"pol ICI "<<invconv<Poly>(*iter,serv);
1360  // mdebug()<<"forget "<<tmpdeg;
1361  recall=verytmp;
1362  forget(b,dump,tmpdeg,recall);
1363  serv.compress(tmpPk,dump,b,tmpdeg);
1364  verytmp.erase(verytmp.begin(),verytmp.end());
1365 
1366 
1367 
1368  //forget(b,dump,tmpdeg);
1369 
1370  ttmpP=Crochet(ttmpP,tmpdeg);
1371 #ifdef NO_TWICE
1372  for(typename typP::iterator iter=P.begin()
1373  ;iter!=P.end();iter++)
1374  if(Degree(*iter)==tmpdeg)
1375  P.erase(iter--);
1376 #endif
1377  nbproj=ttmpP.size();
1378  //for(typename typP::iterator iter=ttmpP.begin();iter!=ttmpP.end();iter++)
1379  // tmpconv.push_back(convert<pol>(*iter,b));
1380  //ttmpP.erase(ttmpP.begin(),ttmpP.end());//samedi 24 aout
1381  //nbproj=1;// samedi 24 aout
1382  //Crochetchoix(P,Pk,tmpconv,tmpdeg,b);
1383  //ttmpP.push_back(invconv<typename typP::value_type>(Pk.back()));
1384  //Pk.pop_back();//on elimine de Pk le pol de proj choisi
1385  //pour se ramner au cas ou proj != 0 et P ok
1386  Mk.erase(Mk.begin(),Mk.end());
1387  for (typename typPk::iterator iterpk=tmpPk.begin();
1388  iterpk!=tmpPk.end();iterpk++)
1389  {
1390  // mdebug()<<"iterpk->ind "<<iterpk->ind
1391  // <<"et Isinb(iterpk-ind,b) "<<IsinB(iterpk->ind,b);
1392  if (IsinB(iterpk->ind,b))
1393  {
1394  for(typename Base::iterator iterb=b.begin()
1395  ;iterb!=b.end();iterb++)
1396  {
1397  iterb->refuse=(mon*)MAC_REV_REALLOC<mon>
1398  (iterb->refuse,(iterb->taille2)*sizeof(mon),
1399  (iterb->taille2+1)*sizeof(mon));
1400  iterb->refuse[iterb->taille2]=iterpk->ind;
1401  iterb->refuse[iterb->taille2].SetCoeff(1);
1402  iterb->taille2++;
1403  }
1404  }
1405  }
1406  for(typename typPk::iterator iter=tmpPk.begin()
1407  ;iter!=tmpPk.end();iter++)
1408  {
1409  mon *tmp=(mon*)MAC_REV_MALLOC<mon>(sizeof(mon));
1410  *tmp=iter->ind;
1411  tmpMk.push_back(*tmp);
1412  exceed.push_back(tmp);
1413  sizeexceed.push_back(1);
1414  // mdebug()<<"pour Mk iter->ind "<<iter->ind;
1415  }
1416  //on vient de chuter de
1417  //degre et on a reconstitue les bons mk et pk
1418  k=tmpdeg-1;
1419  //Dump(Pk,dump);
1420  }
1421  }
1422  }
1423  if(flagdoplus)
1424  {
1425 #ifdef unimul
1426  serv.reinit_multiple();
1427 #endif
1428 #ifdef uninf
1429  serv.reinit_nf();
1430 #endif
1431 
1432  //list<mon> dejavu;
1433  //construction des pol ord sup
1434  //
1435  // Problewme de gestion des accumulateurs ici
1436  int i;
1437  //#pragma omp parallel for shared(tmpPk,tmpMk,exceed,sizeexceed) private(i)
1438  for(i=0;i<b.nbvar();i++)
1439  {
1440 
1441  for(typename typPk::iterator iter=Pk.begin();iter!=Pk.end();iter++)
1442  {
1443  //mdebug()<<"avant gooddir l'ind "<<iter->ind<<" et i "<<i;
1444  if(GoodDir((*iter).ind,i,b))
1445  {
1446  mon tmpmon=iter->ind*mon(i,1);
1447  tmpmon.SetCoeff(1);
1448  //mdebug()<<"tmpmon ici "<<tmpmon;
1449  if(!member(tmpMk,tmpmon))
1450  {
1451  //
1452  //ici aussi rendre les op atomic ne suffit pas
1453  //il faut aussi faire l'ajout.
1454  //#pragma omp atomic
1455  typename typrecall::iterator iterrec;
1456 
1457  mon * tmp=NULL;
1458  int tmpsize;
1459  //mdebug()<<"avant mult";
1460  pol stockmult=mult(i,*iter,dump,b,&tmp,tmpsize,w
1461  ,serv);
1462 
1463 
1464 
1465  //attention le coefficient de ind est la taille de
1466  // l'exceed correspondant
1467 #ifdef DEB
1468  mdebug()<<"apres mult "<<tmpsize;
1469 #endif
1470  //#pragma omp critical
1471  {
1472  tmpMk.push_back(tmpmon);
1473  tmpPk.push_back(stockmult);
1474  exceed.push_back(tmp);
1475  sizeexceed.push_back(tmpsize);
1476  }
1477  }
1478  }
1479  }
1480  }
1481  // recall=verytmp;
1482  }
1483  // mdebug()<<"tmpPK.size "<<tmpPk.size();
1484  //Constructyion de pol terminee
1485  //mdebug()<<"ttmpP.size "<<ttmpP.size();
1486  //mdebug()<<"avant conv merge";
1487  conv_merge(tmpPk,ttmpP,exceed,sizeexceed,b,serv);
1488  //mdebug()<<"apres conv merge";
1489  //mdebug()<<"tmpMK==empty?? "<<tmpMk.empty();
1490 #ifdef DEB
1491  for(typename typPk::iterator iter=tmpPk.begin();iter!=tmpPk.end();iter++)
1492  mdebug()<<"invocn"<<invconv<Poly>(*iter,serv);
1493 #endif
1494  // {
1495  // mdebug()<<"iter->ind.GetCoeff() "<<iter->ind.GetCoeff()-1
1496  // <<" et ind "<<iter->ind;
1497  // }
1498  //mdebug()<<"tmpPK.size "<<tmpPk.size();
1499 
1500  //Tue Nov 27 15:57:25 MET 2001
1501  /*
1502  for(typename typPk::iterator iter=tmpPk.begin();iter!=tmpPk.end();iter++)
1503 
1504  {
1505  mon tmp=iter->ind;
1506  tmp.SetCoeff(1);
1507  if((iter->ind.GetDegree()!=0)&&(!member(tmpMk,tmp)))
1508  {
1509  //mdebug()<<"iter->ind.GetCoeff()"<<iter->ind.GetCoeff()-(coeff)1<<" "
1510  // <<tmp;
1511  tmpMk.push_back(tmp);
1512  //tmpMk.back().SetCoeff(1);
1513  }
1514  }
1515  */
1516  //construction de la matrice maintenant
1517  //positionnement des pointeurs de matrices pour gagner en tps
1518  {
1519  int nnz=0;
1520  //-----------------
1521  //for(typename typMk::iterator iter=tmpMk.begin();iter!=tmpMk.end();iter++)
1522  // mdebug()<<"mon "<<*iter;
1523  //-----------------
1524  // mdebug()<<"nb proj "<<nbproj;
1525  //for(list<int>::iterator iter=sizeexceed.begin();iter!=sizeexceed.end()
1526  // ;iter++)
1527  // {
1528  // nnz+=*iter;
1529  // mdebug()<<"on a donc de monme pour celui la "<<
1530  // (int)((iter->ind).GetCoeff()-1)<<" d'ou nnz "<<nnz;
1531  // }
1532  //for(typename typPk::iterator iter=tmpPk.begin();iter!=tmpPk.end();iter++)
1533  // {
1534  // int tmp=0;
1535  // for(;(iter->ind).GetCoeff()>(coeff)tmp;tmp++);
1536  // nnz+=tmp-1;
1537  // mdebug()<<"on a donc de monme pour celui la "<<
1538  // (int)((iter->ind).GetCoeff()-1)<<" d'ou nnz "<<nnz;
1539  // }
1540  // nnz+=tmpPk.size()-nbproj;//le 10/05/2001 avant: nnz+=tmpMk.size()-nbproj;
1541  //le -nbproj rajoute le 23/01/2001!!
1542  // #define DEBUG
1543 #ifdef DEBUG
1544  typename typPk::iterator tmpiter=tmpPk.begin();
1545  list<int>::iterator itersize=sizeexceed.begin();
1546  for(typename list<mon*>::iterator iter=exceed.begin();iter!=exceed.end()
1547  ;tmpiter++,iter++,itersize++)
1548  {
1549  mdebug()<<"pour un exceeed la taille attendue est "
1550  <<*itersize;
1551  for(int i=0;*itersize>i;i++)//;tmpiter->ind.GetCoeff()>i+1;i++)
1552  {
1553  mdebug()<<"le gus que je regarde "<<(*iter)[i];
1554  if(!member(tmpMk,(*iter)[i]))
1555  {
1556  mdebug()<<(*iter)[i];
1557  mdebug()<<"dans ce pol "<<invconv<Poly>(*tmpiter);
1558  mdebug()<<"l'ind "<<tmpiter->ind;
1559  exit(-10);
1560  }
1561  }
1562  }
1563 #endif
1564  // #undef DEBUG
1565  mat.nrow=tmpPk.size();
1566  mat.ncol=tmpMk.size();
1567  ncol=(int*)calloc(mat.ncol,sizeof(int));
1568  //memset((void*)ncol,0,(*mat)->ncol*sizeof(int));
1569  map<mon,int,ord> mapint;
1570  typename typMk::iterator itermk=tmpMk.begin();
1571  typename typPk::iterator iterpk=tmpPk.begin();
1572  typename list<mon*>::iterator iterex=exceed.begin();
1573  list<int>::iterator itersize=sizeexceed.begin();
1574  int line_num=0;
1575  for(int i=0;itermk!=tmpMk.end();i++,itermk++)
1576  {
1577  mapint[*itermk]=i;
1578  //mdebug()<<"mapint "<<*itermk<<" "<<i;
1579  }
1580  iterpk=tmpPk.begin();
1581  nnz=0;
1582  for(;iterpk!=tmpPk.end();iterpk++,iterex++,itersize++,line_num++)
1583  {
1584  nnz+=*itersize;
1585 #ifdef DEB
1586  mdebug()<<"itesize "<<*itersize;
1587 #endif
1588  // if(iterpk->ind.GetDegree()>0)
1589  // {
1590  // mon tmpmon=iterpk->ind;
1591  // tmpmon.SetCoeff(1);
1592  // ncol[mapint[tmpmon]]++;
1593  // mdebug()<<"on est a la ligne "<<line_num<<" ai augmente "<<mapint[(iterpk->ind*(1/iterpk->ind.GetCoeff()))]<<" a cause de son monome dom ";
1594  // mdebug()<<"le voila "<<(iterpk->ind*(1/iterpk->ind.GetCoeff()));
1595  //}
1596 
1597  {
1598  // coeff tmpcoeff=(iterpk->ind).GetCoeff()-(coeff)1;
1599  for(int i=0;*itersize>i;i++)
1600  {
1601  //mdebug()<<"iterpk icici ->GetCoeff()"<<iterpk->ind.GetCoeff();
1602  if (mapint.count((*iterex)[i]))
1603  {
1604 #ifdef DEB
1605  mdebug()<<"iterex ici "<<(*iterex)[i];
1606 #endif
1607  ++ncol[mapint[(*iterex)[i]]];
1608  // mdebug()<<" ai augmente "<<mapint[(*iterex)[i]]<<"a cause de iterex ";
1609  }
1610 #ifdef FAST
1611  else
1612  {
1613  int j,tmpsize;
1614  mdebug()<<"OUILLE OUILLE OUILLE "<<(*iterex)[i];
1615  mat.ncol++;
1616  mat.nrow++;
1617  mapint[(*iterex)[i]]=mat.ncol-1;
1618  ncol=(int*)realloc(ncol,mat.ncol*sizeof(int));
1619  ncol[mapint[(*iterex)[i]]]=1;
1620  //reste a calculer un proj de ce monome et le tour sera joue
1621  Poly tmppolyalp;
1622  mon tmpmon=(*iterex)[i];
1623  pol tmppol;
1624  mon *tmp;
1625  tmpMk.push_back(tmpmon);
1626  for(j=0;j<b.nbvar()&&(!tmpmon.GetDegree(j))
1627  &&IsinB(my_divmon(tmpmon,j,b),b);j++);
1628  tmppolyalp=Poly(my_divmon(tmpmon,j,b));
1629  //proj(tmppolyalp,dump,b);
1630  tmppol=convert<pol>(tmppolyalp,b);
1631  tmpPk.push_back(mult(j,tmppol,dump,b,&tmp,tmpsize,w));
1632  exceed.push_back(tmp);
1633  sizeexceed.push_back(tmpsize);
1634  mat.nnz+=tmpsize;
1635  }
1636 #endif
1637  }
1638  }
1639  }
1640  //rajoute le 23/01/2001
1641  //enleve le 17.05/2001
1642  // for(int i=mat.nrow-nbproj;i<mat.nrow;i++)
1643  // ncol[i]--;
1644  //fin rajout c pour prendre en compte que lors de
1645  //de la proj le monome de tete se trouve dans exceed
1646  //fin enlevement
1647  // mdebug()<<"nnz of matrix "<<nnz;
1648  //en general la matrice est rectangulaire
1649  void * toto=MAC_REV_MALLOC<coeff>(nnz*sizeof(typename typMat::coeff_t));
1650  //+nnz*sizeof(int)+((*mat)->ncol+1)*sizeof(int));
1651  if (mat.nzval!=NULL) MAC_REV_FREE<typename typMat::coeff_t>
1652  (mat.nzval,mat.nnz
1653  *sizeof(typename typMat::coeff_t));
1654  mat.nnz=nnz;
1655  mat.nzval=toto;
1656 
1657 
1658  //for(int i=0;i<nnz;i++)
1659  // ((typename typMat::coeff_t*)mat.nzval)[i]=(coeff)1000000000000.0;
1660 
1661 
1662  // mstore->rowind=(int*)((double*)toto+nnz);
1663  if(mat.rowind!=NULL) free(mat.rowind);
1664  mat.rowind=(int*)malloc(nnz*sizeof(int));
1665  //mstore->colptr=(int*)(mstore->rowind+nnz);//
1666  if(mat.colptr!=NULL) free(mat.colptr);
1667  mat.colptr=(int*)malloc((mat.ncol+1)*sizeof(int));
1668  // mdebug()<<"matrix of size "<<mat.nrow<<"x"<<mat.ncol;
1669  //matrice cree rest plus qu'a ecrire dedans ca
1670  typename typMat::coeff_t* nzval=(typename typMat::coeff_t*)mat.nzval;
1671  int *rowind=mat.rowind;
1672  int *colptr=mat.colptr;
1673  //mise en place des colptr
1674  *colptr=0;
1675  for(int i=0;i<mat.ncol-1;i++)
1676  {
1677  colptr[i+1]=colptr[i]+ncol[i];
1678  //mdebug()<<"nb element colonne "<<i<<" "<<ncol[i];
1679  //mdebug()<<"les valeurs de debut de colonne "<<8*colptr[i+1];
1680  }
1681  colptr[mat.ncol]=mat.nnz;
1682 
1683  //mdebug()<<"mat.nnz "<<mat.nnz<<" "<<mat.nnz*8;
1684  //fin mise en place
1685  //la on fait des manip symboliques pour savoir ou ecrire
1686  //les coeffs de chaques colonnes ca evite d'avoir a les
1687  //recopier 50000 fois pour construire la matrice
1688  //c bcp plus rapide
1689  itermk=tmpMk.begin();
1690  iterpk=tmpPk.begin();
1691  iterex=exceed.begin();
1692  itersize=sizeexceed.begin();
1693  line_num=0;
1694  //mdebug()<<"nbproj "<<nbproj;
1695  for(;0&&(line_num < (int)(tmpPk.size()-nbproj));
1696  iterpk++,iterex++,itersize++,line_num++)
1697  //les polynomes projetes sont foutut au bout de la liste!!
1698  //d'ou le tmpPk.size()-nbproj
1699  {
1700  mon tmpmon=iterpk->ind;
1701  tmpmon.SetCoeff(1);
1702  //const unsigned int tmpplus=line_num+1,tmp=line_num;
1703  const unsigned int tmpplus
1704  =mapint[tmpmon]+1,tmp=tmpplus-1;
1705  // mdebug()<<"ici line_num"<<line_num<<" le mon dom est col "<<tmp;
1706  //ici couille car le ind n'est pas forcement le mon diag
1707  nzval[colptr[tmpplus]-ncol[tmp]]=1;
1708  //fin couille changer un truc
1709  // mdebug()<<"mat.nnz "<<mat.nnz;
1710  //mdebug()<<"colptr[...]"<<colptr[tmpplus];
1711  //mdebug()<<"ncol[....]"<<ncol[tmp];
1712  //mdebug()<<"1 a la pos "<<colptr[tmpplus]-ncol[tmp];
1713  rowind[colptr[tmpplus]-ncol[tmp]]=line_num;
1714  //mdebug()<<"1 a la ligne "<<rowind[colptr[tmpplus]-ncol[tmp]];
1715  ncol[tmp]--;
1716  {
1717  //coeff tmpcoeff=iterpk->ind.GetCoeff()-(coeff)1;
1718  for(int i=0;*itersize>i;i++)
1719  {
1720  const coeff stock=(*iterex)[i].GetCoeff();
1721  (*iterex)[i].SetCoeff(1);
1722  const unsigned int tmpplusex=mapint[(*iterex)[i]]+1;
1723  const unsigned int tmpex=mapint[(*iterex)[i]];
1724  nzval[colptr[tmpplusex]-ncol[tmpex]]
1725  =stock;
1726  rowind[colptr[tmpplusex]-ncol[tmpex]]=line_num;
1727  ncol[tmpex]--;
1728  mdebug()<<(*iterex)[i]<<" valeur de ncol["<<tmpex<<"] "<<ncol[tmpex];
1729  }
1730  }
1731  MAC_REV_FREE<mon>(*iterex,*itersize*sizeof(mon));
1732  }
1733  for(;(iterpk!=tmpPk.end());iterpk++,iterex++,itersize++,line_num++)
1734  {
1735  //La c le traitement des projetes
1736  //const unsigned int tmpplus=line_num+1,tmp=line_num;
1737  {
1738  //mdebug()<<"==================================================";
1739  //coeff tmpcoeff=iterpk->ind.GetCoeff()-(coeff)1;
1740  for(int i=0;*itersize>i;i++)
1741  {
1742  // mdebug()<<" le monome "<<(*iterex)[i];
1743  const coeff stock=(*iterex)[i].GetCoeff();
1744  (*iterex)[i].SetCoeff(1);
1745  const unsigned int tmpplusex=mapint[(*iterex)[i]]+1;
1746  const unsigned int tmpex=mapint[(*iterex)[i]];
1747  nzval[colptr[tmpplusex]-ncol[tmpex]]
1748  =stock;
1749  rowind[colptr[tmpplusex]-ncol[tmpex]]=line_num;
1750  ncol[tmpex]--;
1751  // mdebug()<<"tmpex "<<tmpex;
1752  // mdebug()<<"tmpplusex "<<tmpplusex;
1753  // mdebug()<<"colptr[tmpplusex] "<<colptr[tmpplusex];
1754  // mdebug()<<"ncol[tmpex] "<<ncol[tmpex];
1755  // mdebug()<<"linenum "<<line_num;
1756  // mdebug()<<"nzval "<<stock;
1757  // mdebug()<<"rowind "<<line_num;
1758  // mdebug()<<"colptr[tmpplusex]-ncol[tmpex] "<<colptr[tmpplusex]-ncol[tmpex]+1;
1759  }
1760  }
1761  MAC_REV_FREE<mon>(*iterex,*itersize*sizeof(mon));
1762  //free(*iterex);
1763  }
1764  free(ncol);
1766  for(typename typPk::iterator iterpk=tmpPk.begin()
1767  ;iterpk!=tmpPk.end();iterpk++)
1768  iterpk->ind.SetCoeff(1);
1769  //iterpk->ind*=1/iterpk->ind.GetCoeff();
1770  // mdebug()<<"la matrice dans newpkmkmatrixof";
1771  //mdebug()<<" nombre de non zero "<<mat.nnz;
1772 #ifdef DEB
1773  affdebug(mat);
1774 #endif
1775  //mdebug()<<"le colptr ";
1776  //for(int i=0;i<mat.ncol;i++)
1777  // mdebug()<<mat.colptr[i]<<" ";
1778  // mdebug()<<"les indices de lignes et le nnz";
1779  // for(int i=0;i<mat.nnz;i++)
1780  // mdebug()<<mat.rowind[i]<<" "<<((typename typMat::coeff_t*)mat.nzval)[i];
1781  // mdebug()<<"fin matrice";
1782  Pk.swap(tmpPk);
1783  Mk.swap(tmpMk);
1784  }
1785  // for(int i=0;i<mat.nnz;i++)
1786  // mdebug()<<"coeff mat "<<((typename typMat::coeff_t*)mat.nzval)[i]<<" et i "<<i;
1787  return flagdoplus;
1788 }
1789 //newpkmkmatrixof termine reste plus qu'a tester .... ;-(
1790 
1791 template<typename typdump,typename mon,typename pol>
1792 void findcorresppol(const mon & m,const typdump &dump,pol &res)
1793 {
1794  int i,d=m.GetDegree();
1795  typename typdump::const_iterator iter=dump.begin();
1796  for(;(iter!=dump.end())&&(iter->k!=d);iter++);
1797  {
1798  // mdebug()<<"dump"<<iter->k;
1799  //mdebug()<<"nf"<<iter->nf[i].ind.rep;
1800  //mdebug()<<"m.rep"<<m.rep;
1801  //localise la dump struct de bon degre
1802  if(iter!=dump.end())
1803  {
1804  for(i=0;(i<iter->size)&&(iter->nf[i].ind.rep!=m.rep);i++);
1805  // mdebug()<<"monomes passe en revue pour chercher "<<m<<" : "<<iter->nf[i].ind;;
1806  //trouve dans la dump struct le pol qui marche
1807  if (i<iter->size)//estce que je connais?
1808  //oui super
1809  {
1810  //mdebug()<<"g trouve un gus qui convient";
1811  //for(int j=0;j<iter->nf[i].size;j++)
1812  // mdebug()<<"et ses coeff sont "<<iter->nf[i].nf[j];
1813  //Tue Dec 4 16:12:10 MET 2001
1814  //densify(iter->nf[i],res);}
1815  //Tue Dec 4 16:12:10 MET 2001
1816  res=iter->nf[i];}
1817  else {
1818  //non ben alors il se reecrit en lui meme
1819  res.ind=mon(0);//m*(1/m.GetCoeff());
1820  res.size=-1;//pour proj avant ct -1
1821  res.nf=NULL;
1822  res.nfind=NULL;
1823  res.sizenf=0;
1824  }
1825  }
1826  }
1827  return;
1828 }
1829 template<typename iterator, typename mon,typename Base>
1830 void varmult(iterator iter1,iterator iter2,mon &var1, mon &var2,const Base &b)
1831 {
1832  for(int i=0;i<b.nbvar();i++)
1833  {
1834  if(iter1->ind.GetDegree(i)>iter2->ind.GetDegree(i))
1835  var2=mon(i,1);
1836  if(iter2->ind.GetDegree(i)>iter1->ind.GetDegree(i))
1837  var1=mon(i,1);
1838  }
1839 }
1840 
1841 template <typename mon,typename pol, typename typdump,typename Base,
1842  typename monomial_server>
1843 pol halfred(const mon &targetind, const pol & p, const typdump & dump
1844  , const Base & b, monomial_server & serv)
1845 {
1846  typedef typename pol::coeff_t coeff;
1847  static const mon monzero=mon(0);
1848  static const typename pol::coeff_t zero(0);
1849  int *place=(int*)MAC_REV_MALLOC<int>(p.size*sizeof(int)),comp=0;
1850  pol res;
1851 #ifdef uninf
1852  int varref;
1853 #else
1854  mon varref;
1855 #endif
1856  int i,maxref=0,ppp=0,nbbitsun=0;
1857  for(int j=0;j<p.size;j++)
1858  place[j]=0;
1859  // coeff *tmpspace;
1860  for(i=0;i<b.nbvar();i++)
1861  {
1862  // mdebug()<<"targetind.GetDegree(i)"<<targetind.GetDegree(i);
1863  // mdebug()<<"p.ind.GetDegree(i)"<<p.ind.GetDegree(i);
1864  if(targetind.GetDegree(i)>p.ind.GetDegree(i))
1865  break;
1866  }
1867 #ifdef uninf
1868  varref=i;
1869 #else
1870  varref=mon(i,1);
1871 #endif
1872 
1873  // for(i=0;i<p.size;i++)
1874  // mdebug()<<p.nf[i]<<" ";
1875  //mdebug();
1876  //mdebug()<<"ppp "<<ppp<<" p.size "<<p.size;
1877  //mdebug()<<"p en entree"<<invconv<Poly>(p);
1878  //mdebug()<<"mulind(0,1,b) "<<mulind(0,1,b);
1879  //mdebug()<<"target "<<targetind<<" et ind "<<p.ind<<" et var "<<varref;
1880  // mdebug()<<"p.size";
1881  for(i=0;i<p.sizenf;i++)
1882  {
1883  if(p.nfind[i])
1884  {
1885 #ifdef DEB
1886  mdebug()<<nbbitsun<<" p.nfind[i] "<<(int)p.nfind[i];
1887 #endif
1888  nbbitsun+=nbbits[p.nfind[i]];
1889  for(int j=0;j<8;j++)
1890  {
1891  if((p.nfind[i]>>j)&1)
1892  {
1893 #ifdef unimul
1894  // mdebug()<<"8*i+j"<<8*i+j;
1895  // mdebug()<<"varref"<<varref;
1896  // mdebug()<<"comp"<<comp;
1897  // mdebug()<<"place[comp]"<<place[comp];
1898  place[comp]=serv.mulind(8*i+j,varref,b);
1899  // mdebug()<<"place[comp]"<<place[comp];
1900  if(place[comp]<0)
1901  place[comp]=-1*(8*i+j+1);
1902  if (place[comp]>maxref)
1903  maxref=place[comp];
1904  // mdebug()<<"maxref"<<maxref;
1905 
1906 #else
1907  //ATTENTION UNIQUEMENT UNIMUL SINON PAS DEBUGGE
1908  int indtmpmon;
1909  mon tmpmon;
1910  int2mon(8*i+j,tmpmon);
1911  tmpmon*=varref;
1912  if(IsinB(tmpmon,b))
1913  {
1914  //int j;
1915  indtmpmon=mon2int(tmpmon);
1916  //for(j=0;j<p.size&&p.nfind[j]!=indtmpmon;j++);
1917  //if(j<p.size)
1918  maxref++;
1919  place[comp]=indtmpmon;
1920  }
1921  else
1922  {
1923  place[comp]=-1*(8*i+j);
1924  }
1925 #endif
1926  comp++;
1927  }
1928  }
1929  }
1930  }
1931  // mdebug()<<"despues del for";
1932  coeff *tmpspace=(coeff*)MAC_REV_MALLOC<coeff>((maxref/8+1)*8*sizeof(coeff));
1933 #ifdef DEB
1934  mdebug()<<"comp "<<comp<<" et p.size "<<p.size;
1935 #endif
1936  res.sizenf=(maxref/8+1);
1937  for(int j=0;j<8*res.sizenf;j++)
1938  tmpspace[j]=0;
1939  res.nfind=(unsigned char*)MAC_REV_MALLOC<unsigned char>(res.sizenf);
1940  memset(res.nfind,0,res.sizenf);
1941  // mdebug()<<"por aqui";
1942  for(int j=0;j<nbbitsun;j++)
1943  if(place[j]>=0)
1944  {
1945 #ifdef DEB
1946  mon tmpmm;
1947  int2mon(place[j],tmpmm);
1948  mdebug()<<"le monome la "<<tmpmm<<" "<<p.nf[j];
1949 #endif
1950  tmpspace[place[j]]+=p.nf[j];
1951  res.nfind[place[j]/8]|=1<<(place[j]%8);
1952  }
1953  else
1954  {
1955  pol tmp;
1956 #ifdef unimul
1957  tmp=serv.nf(varref,-1*place[j]-1,dump,b);
1958 #else
1959  mon tmpmon;
1960  int2mon(-1*place[j],tmpmon);
1961  tmpmon*=mon(i,1);
1962  tmp=nf<typpol>(tmpmon,dump,b);
1963 #endif
1964 
1965  if (tmp.size>=0)
1966  {
1967  int compteurtmp=0;
1968  //mdebug()<<" je reduit prr lui"<<endl<<invconv<Poly>(tmp);
1969  //On fait un plus dense dans tmpspace
1970  if(tmp.sizenf>res.sizenf)
1971  {
1972  res.nfind=(unsigned char*)realloc(res.nfind,tmp.sizenf);
1973  for(int k=res.sizenf;k<tmp.sizenf;k++)
1974  res.nfind[k]=0;
1975  tmpspace=(coeff*)MAC_REV_REALLOC<coeff>(tmpspace,
1976  res.sizenf*8*sizeof(coeff),
1977  tmp.sizenf*8*sizeof(coeff));
1978  for(int k=res.sizenf*8;k<tmp.sizenf*8;k++)
1979  tmpspace[k]=0;
1980  res.sizenf=tmp.sizenf;
1981 
1982  }
1983  for(int k=0;k<tmp.sizenf;k++)
1984  {
1985  if(tmp.nfind[k])
1986  {
1987  for(int l=0;l<8;l++)
1988  if((tmp.nfind[k]>>l)&1)
1989  tmpspace[8*k+l]-=p.nf[j]*tmp.nf[compteurtmp++];
1990  }
1991  res.nfind[k]|=tmp.nfind[k];
1992  }
1993  }
1994  else
1995  {
1996  mon tmpmon;
1997  serv.int2mon(-1*place[j]-1,tmpmon);
1998  tmpmon*=mon(varref,1);
1999  for(int i=0;i<serv.sizeplacemon();i++)
2000  mdebug()<<i<<" "<<serv.int2mon_[i];
2001  mdebug()<<"monom "<<tmpmon<<" "<<-1*place[j]
2002  <<" "<<p.nf[j];
2003  cerr<<"Je n'aurrai pas du passer par la ";
2004  exit(10);
2005  };
2006  }
2007  // mdebug()<<"despues del segundo for";
2008  res.ind=mon(0);
2009  comp=0;
2010  for(int j=0;j<res.sizenf;j++)
2011  if(res.nfind[j])
2012  comp+=nbbits[res.nfind[j]];
2013  res.nf=(coeff*)MAC_REV_MALLOC<coeff>(comp*sizeof(coeff));
2014  res.size=comp;
2015  comp=0;
2016  // mdebug()<<"antes del ultimo for";
2017  for(int j=0;j<res.sizenf;j++)
2018  if(res.nfind[j])
2019  for(int k=0;k<8;k++)
2020  if((res.nfind[j]>>k)&1)
2021  if(tmpspace[8*j+k]!=0)
2022  {res.nf[comp++]=tmpspace[8*j+k];
2023 
2024  }
2025  else
2026  res.nfind[j]&=~(1<<k);
2027  // mdebug()<<"despues del ultimo for";
2028  res.nf=(coeff*)MAC_REV_REALLOC<coeff>(res.nf,res.size*sizeof(coeff)
2029  ,comp*sizeof(coeff));
2030  // mdebug()<<"res.nf";
2031  res.size=comp;
2032  // mdebug()<<"res.size";
2033  // for(int j=0;j<8*res.sizenf;j++)
2034  // tmpspace[j]=0;
2035  // mdebug()<<"tmpspace.size"<<res.sizenf*8*sizeof(coeff)<<res.sizenf*8;
2036  MAC_REV_FREE<coeff>(tmpspace,res.sizenf*8*sizeof(coeff));
2037  // mdebug()<<"antes de free";
2038  free(place);
2039  //mdebug()<<" apres reduction "<<invconv<Poly>(res);
2040  return res;
2041 }
2042 
2043 template <typename typPk, typename typdump,typename Base,
2044  typename workspace ,typename monomial_server>//,typename typMk>
2045 void ReduceSpol(typPk &S,const typdump & dump,typPk &redspol,const Base &b
2046  ,workspace &w, monomial_server &serv)
2047 //, typMk & Spoldejavu)
2048 {
2049  typedef typename typPk::value_type pol;
2050  typedef typename pol::monom_t mon;
2051  typedef typename pol::coeff_t coeff;
2052  list<mon> leadmon;
2053  list<typPk> membermon;
2054  mdebug()<<"en entree de RedSpol S.size "<<S.size()/2;
2055  mdebug()<<"reduce";
2056  for(typename typPk::iterator iterS=S.begin();iterS!=S.end();iterS++,iterS++)
2057  {
2058  typename typPk::iterator iterS2=iterS;
2059  typename list<mon>::iterator itermon=leadmon.begin();
2060  typename list<typPk>::iterator itermemb=membermon.begin();
2061  mon var1,var2,tmpleadmon;
2062  iterS2++;
2063  varmult(iterS,iterS2,var1,var2,b);
2064  tmpleadmon=iterS->ind*var1;
2065  //mdebug()<<"tmpleadmon"<<tmpleadmon;
2066  // Spoldejavu.push_back(tmpleadmon);
2067  for(;itermon!=leadmon.end()&&itermon->rep!=tmpleadmon.rep;
2068  itermon++,itermemb++);
2069  if(itermon!=leadmon.end())
2070  {
2071  int yes1=1,yes2=1;
2072  for(typename typPk::iterator tmpiter=itermemb->begin();
2073  tmpiter!=itermemb->end();tmpiter++)
2074  {
2075  if(tmpiter->ind.rep==iterS->ind.rep) yes1=0;
2076  if(tmpiter->ind.rep==iterS2->ind.rep) yes2=0;
2077  }
2078  if(yes1) itermemb->push_back(*iterS);
2079  if(yes2) itermemb->push_back(*iterS2);
2080  }
2081  else
2082  {
2083  // mdebug()<<"iterS"<<iterS->ind;
2084  //mdebug()<<"iterS2"<<iterS2->ind;
2085  typPk verytmp;
2086  leadmon.push_back(tmpleadmon);
2087  //mdebug()<<"le tmpleadmon "<<tmpleadmon;
2088  membermon.push_back(verytmp);
2089  membermon.back().push_back(*iterS);
2090  membermon.back().push_back(*iterS2);
2091  }
2092  }
2093  mdebug()<<"sale del primer for";
2094  //maintenant on a mis en place la structure de Redspol
2095  //reste plus qu'a calculer les bons et a reduire tout ca!!
2096  {
2097  typename list<mon>::iterator itermon=leadmon.begin();
2098  typename list<typPk>::iterator itermemb=membermon.begin();
2099  // mdebug()<<"antes del otro for";
2100  for(;itermon!=leadmon.end();itermon++,itermemb++)
2101  {
2102  typename typPk::value_type ref;
2103  typename typPk::iterator tmpiter=itermemb->begin();
2104  mdebug()<<"antes halfred";
2105  mdebug()<<"itermon"<<*itermon;
2106  // mdebug()<<"tmpiter"<<tmpiter->ind;
2107  ref=halfred(*itermon,*tmpiter,dump,b,serv);
2108  // mdebug()<<"ref.size "<<ref.size<<" "<<ref.ind;
2109  //mdebug()<<"ref en entier"<<endl<<invconv<Poly>(ref);
2110  tmpiter++;
2111  for(;tmpiter!=itermemb->end();tmpiter++)
2112  {
2113  typename typPk::value_type redcur;
2114  redcur=halfred(*itermon,*tmpiter,dump,b,serv);
2115  my_plus(redcur,ref,(coeff)1,&w);
2116  if(!Ispolzero(redcur))
2117  {
2118  redspol.push_back(redcur);
2119 
2120  // mdebug()<<"le leadmon qui donne un non zero "<<*itermon;
2121  }
2122  else
2123  {
2124  MAC_REV_FREE<coeff>(redcur.nf,redcur.size*sizeof(coeff));
2125  MAC_REV_FREE<unsigned char>(redcur.nfind,redcur.sizenf);
2126  }
2127  }
2128  MAC_REV_FREE<coeff>(ref.nf,ref.size*sizeof(coeff));
2129  MAC_REV_FREE<unsigned char>(ref.nfind,ref.sizenf);
2130  }
2131  }
2132 }
2133 
2134 template<typename typPk,typename monomial_server>
2135 int mindeg(const typPk &l,monomial_server &serv)
2136 {
2137  typedef typename typPk::value_type::monom_t mon;
2138  typedef typename mon::coeff_t coeff;
2139  int res=-1;
2140  // mdebug()<<"dans mindeg res au debut"<<res;
2141  for(typename typPk::const_iterator iter=l.begin() ;iter!=l.end();iter++)
2142  {
2143  int mmax=-1;
2144  // mdebug()<<"tour de boucle de mindeg ";
2145  for(int i=0;i<iter->sizenf;i++)
2146  if(iter->nfind[i])
2147  for(int j=0;j<8;j++)
2148  if((iter->nfind[i]>>j)&1)
2149  {
2150  mon tmpmon;
2151  serv.int2mon(8*i+j,tmpmon);
2152  mmax=(mmax<tmpmon.GetDegree())?tmpmon.GetDegree():mmax;
2153  }
2154  // mdebug()<<"mmax ici "<<mmax;
2155  if (iter->ind.GetCoeff()==(coeff)1)
2156  mmax=(mmax<iter->ind.GetDegree())?iter->ind.GetDegree():mmax;
2157  //mdebug()<<"les gens que je regarde "<<iter->ind;
2158  if(mmax!=-1)
2159  {
2160  if(res==-1) res=mmax;
2161  res=(mmax>res)?res:mmax;
2162  }
2163  }
2164  //mdebug()<<"resultat de mindeg "<<res;
2165  return res;
2166 }
2167 
2168 template<typename coeff,typename typiterpk, typename workspace>
2169 void inline update_crochetchoix(typiterpk &iterPk, typiterpk &iter,int tmpind
2170  ,workspace &w)
2171 {
2172 
2173  if((iterPk->sizenf>tmpind/8)&&(((iterPk->nfind[tmpind/8])>>(tmpind%8))&1))
2174  {
2175  // mdebug()<<"g trouve un gus qui a un coeff du mon choisi ";
2176  //mdebug()<<"le voila "<<invconv<Poly>(*iterPk);
2177  int i=0;
2178  coeff tmp;
2179  for(int j=0;8*j<=tmpind;j++)
2180  i+=nbbits[iterPk->nfind[j]];
2181  for(int j=7;j>=(tmpind % 8);j--)
2182  if((iterPk->nfind[tmpind/8]>>j)&1)
2183  i--;
2184 
2185 #ifdef DEB
2186  mdebug()<<"i ici "<<i<<" et iterPk->size "<<iterPk->size;
2187 #endif
2188  tmp=iterPk->nf[i];
2189  // mdebug()<<"le coeff "<<tmp<<" tmpind "<<tmpind;
2190  //mdebug()<<"les voisins "<<iterPk->nf[i-1]<<" "<<iterPk->nf[i+1];
2191 
2192  //mdebug()<<"les voisins plus loin"<<iterPk->nf[i-2]<<" "<<iterPk->nf[i+2];
2193  for(int j=i;j<iterPk->size-1;j++){
2194  iterPk->nf[j]=iterPk->nf[j+1];
2195  }
2196  iterPk->nfind[tmpind/8]&=~(1<< (tmpind%8));
2197  iterPk->nf=(coeff*)MAC_REV_REALLOC<coeff>(iterPk->nf,
2198  iterPk->size*sizeof(coeff),(iterPk->size-1)*sizeof(coeff));
2199  iterPk->size--;
2200  //iterPk->nf[i]=0;
2201  my_plus(*iterPk,*iter,tmp,&w);
2202 
2203  //mdebug()<<"apres plus le coeff vaut "
2204  // <<iterPk->nf[mon2int(iter->ind)];
2205  }
2206  //mdebug()<<"apres mise a jour"<<endl<<invconv<Poly>(*iterPk);
2207 }
2208 template<typename typP,typename typPk,typename typMk
2209  ,typename Base,typename typdump, typename typrecall,
2210  typename workspace ,typename monomial_server>
2211 void NewCrochetchoix(typP & P,typPk & Pk,typMk & Mk, typPk & redspol
2212  ,int &k,typdump &dump, Base &b, typrecall & recall
2213  ,workspace &w, monomial_server &serv)
2214 //,typMk &Spoldejavu)
2215 //template<typename typP,typename typPk,typename typMk,typename Base,typename typdump>
2216 //void NewCrochetchoix(const typP & P,typPk & Pk,typMk & Mk, typPk & redspol
2217 // ,int &k,typdump &dump, Base &b)
2218 {
2219  typedef typename typP::value_type polyalp;
2220  typedef typename polyalp::monom_t mon;
2221  typedef typename polyalp::coeff_t coeff;
2222  //for(typename typPk::iterator iter=redspol.begin();iter!=redspol.end();iter++)
2223  // P.push_back(invconv<Poly>(*iter));
2224  //mdebug()<<"dans newcrochetchoix";
2225  /*
2226 
2227  //mdebug()<<"les sizes "<<iter->size;
2228  {
2229  int cpt=0;
2230  mdebug()<<invconv<typename typP::value_type>(*iter);
2231  for(int i=0;i<iter->sizenf;i++)
2232  for(int j=0;j<8;j++)
2233  if (((iter->nfind[i]>>j)&1) && (iter->nf[cpt++]==0))
2234  {
2235  mdebug()<<"ici ca couille a la place "<<8*i+j;
2236  exit(1);
2237  }
2238  }
2239  */
2240  typPk res=Crochet(redspol,k,serv),tmpres;
2241  //mdebug()<<"redspol.size "<<redspol.size();
2242  //mdebug()<<"res.size "<<res.size();
2243  for(typename typPk::iterator iter=res.begin();iter!=res.end();iter++)
2244  {
2245  for(typename typPk::iterator iterredspol=redspol.begin();
2246  iterredspol!=redspol.end();iterredspol++)
2247  if(iter->nf==iterredspol->nf)
2248  {
2249  redspol.erase(iterredspol);
2250  break;
2251  }
2252  }
2253  NewDestroyredspol(redspol,Pk);
2254  redspol.erase(redspol.begin(),redspol.end());
2255  //for(typename typPk::iterator iter=res.begin();iter!=res.end();iter++)
2256  // mdebug()<<"res->sizenf<< "<<iter->sizenf;
2257  for(typename typPk::iterator iter=res.begin();iter!=res.end();iter++)
2258  {
2259  mdebug()<<"un nouveau tour de boucle ";
2260  if (!Ispolzero(*iter))
2261  {
2262  int tmpind,i;
2263  iter->ind=choice(*iter,b,serv);
2264  //mdebug()<<"le pol pas zero "<<(invconv<polyalp>(*iter));
2265  //mdebug()<<"iter -> size "<<iter->size<<" "<<iter->ind;
2266  //for (int i=0;i<iter->size;i++)
2267  // mdebug()<<iter->nf[i]<<" "<<iter->nfind[i];
2268  if(iter->ind.GetDegree()<k)
2269  {
2270  //mdebug()<<"le pol qui fait chuter"<<(invconv<polyalp>(*iter));
2271  k=iter->ind.GetDegree();
2272  iter->ind=mon(0);
2273  //mdebug()<<"passage par bout de code en plus ";
2274  my_merge(res,tmpres);
2275  my_merge(redspol,res);
2276  tmpres.erase(tmpres.begin(),tmpres.end());
2277  res.erase(res.begin(),res.end());
2278  res=Crochet(redspol,k,serv);
2279  //mdebug()<<"newres.size "<<res.size();
2280  for(typename typPk::iterator iterres=res.begin()
2281  ;iterres!=res.end();iterres++)
2282  {
2283  for(typename typPk::iterator iterredspol=redspol.begin();
2284  iterredspol!=redspol.end();iterredspol++)
2285  if(iterres->nf==iterredspol->nf)
2286  {
2287  redspol.erase(iterredspol);
2288  break;
2289  }
2290  }
2291  for(typename typPk::iterator iterredspol=redspol.begin();
2292  iterredspol!=redspol.end();iterredspol++)
2293  if(Degree(*iterredspol,serv)<=k)
2294  {
2295  tmpres.push_back(*iterredspol);
2296  redspol.erase(iterredspol--);
2297  }
2298  NewDestroyredspol(redspol,Pk);
2299  redspol.erase(redspol.begin(),redspol.end());
2300  my_merge(redspol,tmpres);
2301  tmpres.erase(tmpres.begin(),tmpres.end());
2302  iter=res.begin();
2303  iter->ind=choice(*iter,b,serv);
2304 
2305  Dump(Pk,dump);//c pour detruire Pk dans forget
2306  //ca ne respecte pas la structure du dump!!!
2307  recover(Pk,Mk,dump,k);
2308  // Listedejafait(Mk,Spoldejavu);
2309  forget(b,dump,k,recall);
2310  //forget(b,dump,k);
2311  // mdebug()<<"un coup de compress";
2312  //ce compresss la vire le 06 10 2006
2313  serv.compress(Pk,redspol,dump,b,k);
2314  //pose un probleme de oherence pour le int2mon
2315  //si on le laisse
2316  //swap(*iter,*res.rbegin());//on veut ne plus
2317 
2318 
2319  //iter=--res.end();//refaire d'autre tours
2320  }
2321  //mdebug()<<"g choisi dans Crochetchoix "<<iter->ind;
2322  //mdebug()<<"a t'il un mon2int ?"<<IsinB(iter->ind,b);
2323  //mon tmp10,tmp11;
2324  // int2mon(10,tmp10);int2mon(11,tmp11);
2325  //mdebug()<<"qu'en est il de int2mon(10) et 11 "<<tmp10<<" "<<tmp11;
2326  //
2327  //
2328  //
2329  //UGLY HACK
2330  remiseenforme(*iter);
2331  //SALE HACK
2332  //
2333  //
2334  //
2335  //
2336  //
2337  *iter/=iter->ind.GetCoeff();//pour la "nf"
2338  //mdebug()<<"le coeff "<<iter->ind;
2339  iter->ind*=(coeff)1/iter->ind.GetCoeff();//pour l'ind
2340  tmpind=serv.mon2int(iter->ind);
2341  i=0;
2342  for(int j=0;8*j<=tmpind;j++)
2343  i+=nbbits[iter->nfind[j]];
2344  for(int j=7;j>=(tmpind % 8);j--)
2345  if((iter->nfind[tmpind/8]>>j)&1)
2346  i--;
2347  for(int j=i;j<iter->size-1;j++)
2348  iter->nf[j]=iter->nf[j+1];
2349  iter->nf=(coeff*)MAC_REV_REALLOC<coeff>(iter->nf,
2350  iter->size*sizeof(coeff),
2351  (iter->size-1)*sizeof(coeff));
2352  iter->nfind[tmpind/8]&=~(1<<(tmpind%8));
2353  iter->size--;
2354 #ifdef DEB
2355  mdebug()<<"apres tambouille"<<endl<<invconv<Poly>(*iter);
2356 #endif
2357 
2358  //
2359  //
2360  //
2361  //UGLY HACK
2362  remiseenforme(*iter);
2363  //SALE HACK
2364  //
2365  //
2366  //
2367  //
2368  //
2369 
2370  tmpres.push_back(*iter);
2371  //on met a 0 le coeff du mon passe en ind
2372  //maintenant on a une regle en forme comme il faut reste a:
2373  //-reduire les regles deja existantes par rapport a celle ci
2374  //-reduire les autres redspol par rapport a celle ci
2375  for(typename typPk::iterator iterPk=Pk.begin();iterPk!=Pk.end();iterPk++)
2376  update_crochetchoix<coeff>(iterPk,iter,tmpind,w);
2377  typename typPk::iterator iterredspol=iter;
2378  iterredspol++;
2379  for(;iterredspol!=res.end();iterredspol++)
2380  {
2381  //mdebug()<<"redspol.size() "<<redspol.size();
2382  //mdebug()<<"iterredspol-ind "<<iterredspol->ind<<" "<<iterredspol->nf;
2383  update_crochetchoix<coeff>(iterredspol,iter,tmpind,w);
2384  }
2385  //au tour de tmpres...
2386  typename typPk::iterator iterend=tmpres.end();
2387  iterend--;
2388  for(iterredspol=tmpres.begin();iterredspol!=iterend;iterredspol++)
2389  {
2390  //mdebug()<<"redspol.size() "<<redspol.size();
2391  //mdebug()<<"iterredspol-ind "<<iterredspol->ind<<" "<<iterredspol->nf;
2392  update_crochetchoix<coeff>(iterredspol,iter,tmpind,w);
2393  }
2394 
2395  //voila c fait
2396  //maintenant on met a jour la Base b.
2397  //comment on met a jour la base b
2398  mdebug()<<"avant met a jour";
2399  for(typename Base::iterator iterb=b.begin();iterb!=b.end();iterb++)
2400 
2401  {
2402  mdebug()<<"met a jour";
2403  mdebug()<<"iterb->taille1"<<iterb->taille1;
2404  mdebug()<<"iterb->taille2"<<iterb->taille2;
2405  int tmptaille1=0;
2406  mon *tmp=(mon*)MAC_REV_MALLOC<mon>(iterb->taille1*sizeof(mon));
2407  for(int i=0;i<iterb->taille1;i++)
2408  if(!isdivisible(iter->ind,iterb->accept[i],b))
2409  tmp[tmptaille1++]=iterb->accept[i];
2410  tmp=(mon*)MAC_REV_REALLOC<mon>(tmp,iterb->taille1*sizeof(mon)
2411  ,tmptaille1*sizeof(mon));
2412  MAC_REV_FREE<mon>(iterb->accept,iterb->taille1*sizeof(mon));
2413  iterb->accept=tmp;
2414  iterb->taille1=tmptaille1;
2415  if(iterb->taille1==0)
2416  {
2417  mdebug()<<"passe par ici";
2418  MAC_REV_FREE<mon>(iterb->refuse,iterb->taille2*sizeof(mon));
2419  b.erase(iterb--);
2420  continue;
2421  }
2422  iterb->refuse=(mon*)MAC_REV_REALLOC<mon>
2423  (iterb->refuse,(iterb->taille2)*sizeof(mon),
2424  (iterb->taille2+1)*sizeof(mon));
2425  iterb->refuse[iterb->taille2]=iter->ind;
2426  mdebug()<<"reffuse rajoute "<<iter->ind;
2427  iterb->refuse[iterb->taille2].SetCoeff(1);
2428  iterb->taille2++;
2429  }
2430  res.erase(iter--);
2431  }
2432  //break;//comme ca on fait pas 50 tours.
2433  }
2434  for (typename typPk::iterator iterpk=Pk.begin();iterpk!=Pk.end();iterpk++)
2435  {
2436  //mdebug()<<"iterpk->ind "<<iterpk->ind<<"et Isinb(iterpk-ind,b) "<<IsinB(iterpk->ind,b);
2437  //mdebug()<<"le pol "<<invconv<Poly>(*iterpk);
2438  if (IsinB(iterpk->ind,b))
2439  {
2440  for(typename Base::iterator iterb=b.begin();iterb!=b.end();iterb++)
2441  {
2442  iterb->refuse=(mon*)MAC_REV_REALLOC<mon>
2443  (iterb->refuse,(iterb->taille2)*sizeof(mon),
2444  (iterb->taille2+1)*sizeof(mon));
2445  iterb->refuse[iterb->taille2]=iterpk->ind;
2446  iterb->refuse[iterb->taille2].SetCoeff(1);
2447  iterb->taille2++;
2448  }
2449  }
2450  }
2451 #ifdef NO_TWICE_FALL
2452  for(typename typPk::iterator iter=tmpres.begin();iter!=tmpres.end();iter++)
2453  P.push_back(invconv<polyalp>(*iter));
2454 #endif
2455 #ifdef BONRAJOUT
2456  for(typename typPk::iterator iter=tmpres.begin();iter!=tmpres.end();iter++)
2457  P.push_back(invconv<polyalp>(*iter));
2458  mdebug()<<"a la fin de \A7New";
2459 #endif
2460 #ifdef DEB
2461  for(typename typPk::iterator iter=tmpres.begin();iter!=tmpres.end();iter++)
2462  mdebug()<<"tmpres"<<endl<<invconv<polyalp>(*iter);
2463  mdebug()<<"tmpres.size "<<tmpres.size();
2464 #endif
2465  my_merge(Pk,tmpres);//les gens dans res\tmpres ne sont pas detruits ici!!
2466  my_merge(redspol,res);
2467 
2468 
2469  /* for(typename typPk::iterator iter=Pk.begin();iter!=Pk.end();iter++)
2470  //mdebug()<<"les sizes "<<iter->size;
2471  {
2472  int cpt=0;
2473  mdebug()<<invconv<typename typP::value_type>(*iter);
2474  for(int i=0;i<iter->sizenf;i++)
2475  for(int j=0;j<8;j++)
2476  if (((iter->nfind[i]>>j)&1) && (iter->nf[cpt++]==0))
2477  {
2478  mdebug()<<"ici ca couille a la place "<<8*i+j;
2479  exit(1);
2480  }
2481  }
2482  */
2483  return;
2484 }
2485 //#undef DEB
2486 template<typename typPk>
2487 void NewDestroyredspol(typPk &redspol,const typPk &Pk)
2488 {
2489  // typename typPk::const_iterator iterpk=--Pk.end();
2490  for(typename typPk::iterator iter=redspol.begin();
2491  iter!=redspol.end();iter++)
2492  {
2493  // if(iter->nf!=iterpk->nf)
2494  MAC_REV_FREE<typename typPk::value_type::coeff_t>(iter->nf,iter->size*sizeof(typename typPk::value_type::coeff_t));
2495  MAC_REV_FREE<unsigned char>(iter->nfind,iter->sizenf);
2496  }
2497 }
2498 
2499 template<typename typMk>
2500 void Listedejafait(const typMk &Mk,typMk &Spoldejavu)
2501 {
2502  Spoldejavu.erase(Spoldejavu.begin(),Spoldejavu.end());
2503  for(typename typMk::const_iterator iter1=Mk.begin();iter1!=Mk.end();
2504  iter1++)
2505  for(typename typMk::const_iterator iter2=iter1;iter2!=Mk.end();
2506  iter2++)
2507  {
2508  typename typMk::value_type tmp=lcm(*iter1,*iter2);
2509  if (tmp.GetDegree()==iter1->GetDegree()+1)
2510  Spoldejavu.push_back(tmp);
2511  }
2512 }
2513 
2514 //on suppose stockmon vide
2515 
2516 
2517 template<typename Base, typename typstock, typename typdump>
2518 void enumeratemoninv(typstock &stockmon, int k,
2519  const Base &b, const typdump &dump)
2520 {
2521  typedef typename typstock::value_type mon;
2522  int maxdeg=0;
2523  typstock current;
2524  current.push_back(mon(1));
2525  stockmon.push_back(mon(1));
2526  mdebug()<<"k"<<k;
2527  //
2528  //
2529  //on a le degre max : k
2530  // on part de 1 et on decrit le reste en multipliant par les variables
2531  //
2532  //
2533  for(int j=0;j<k;j++)
2534  {
2535  typstock tmpcurrent=current;
2536  current.erase(current.begin(),current.end());
2537  for(int i=b.nbvar()-1;i>=0;i--)
2538 
2539  {
2540  for(typename list<mon>::iterator iter=tmpcurrent.begin();
2541  iter!=tmpcurrent.end();iter++)
2542  {
2543  mon tmpmon=(*iter)*mon(i,1);
2544  tmpmon.SetCoeff(1);
2545  current.push_back(tmpmon);
2546  }
2547  }
2548  //maintenant on tri;
2549  // mdebug()<<"tri";
2550  // mdebug()<<current.size();
2551  for(typename list<mon>::iterator iter=current.begin();iter!=current.end();iter++)
2552  if(!member(stockmon,*iter))
2553  {
2554  stockmon.push_back(*iter);
2555  }
2556  else
2557  {
2558  current.erase(iter--);
2559  }
2560  }
2561  return;
2562 }
2563 
2564 template<typename Base, typename typstock, typename typdump>
2565 void enumeratemon(typstock &stockmon, int k,
2566  const Base &b, const typdump &dump)
2567 {
2568  typedef typename typstock::value_type mon;
2569  int maxdeg=0;
2570  typstock current;
2571  current.push_back(mon(1));
2572  stockmon.push_back(mon(1));
2573  mdebug()<<"k"<<k;
2574  //
2575  //
2576  //on a le degre max : k
2577  // on part de 1 et on decrit le reste en multipliant par les variables
2578  //
2579  //
2580  for(int j=0;j<k;j++)
2581  {
2582  typstock tmpcurrent=current;
2583  current.erase(current.begin(),current.end());
2584  for(int i=0;i<b.nbvar();i++)
2585 
2586  {
2587  for(typename list<mon>::iterator iter=tmpcurrent.begin();
2588  iter!=tmpcurrent.end();iter++)
2589  {
2590  mon tmpmon=(*iter)*mon(i,1);
2591  tmpmon.SetCoeff(1);
2592  current.push_back(tmpmon);
2593  }
2594  }
2595  //maintenant on tri;
2596  // mdebug()<<"tri";
2597  // mdebug()<<current.size();
2598  for(typename list<mon>::iterator iter=current.begin();iter!=current.end();iter++)
2599  if(!member(stockmon,*iter))
2600  {
2601  stockmon.push_back(*iter);
2602  }
2603  else
2604  {
2605  current.erase(iter--);
2606  }
2607  }
2608  //mdebug()<<"ttttttttt "<<stockmon.size();
2609  return;
2610 }
2611 
2612 
2613 
2614 template<typename Base, typename typstock, typename typdump>
2615 void enumeratemonvar(typstock &stockmon, int k,
2616  const Base &b, const typdump &dump,int n)
2617 {
2618  typedef typename typstock::value_type mon;
2619  int maxdeg=0;
2620  typstock current;
2621  current.push_back(mon(1));
2622  stockmon.push_back(mon(1));
2623  mdebug()<<"k"<<k;
2624  //
2625  //
2626  //on a le degre max : k
2627  // on part de 1 et on decrit le reste en multipliant par les variables
2628  //
2629  //
2630  for(int j=0;j<k;j++)
2631  {
2632  typstock tmpcurrent=current;
2633  current.erase(current.begin(),current.end());
2634  for(int i=0;i<n;i++)
2635 
2636  {
2637  for(typename list<mon>::iterator iter=tmpcurrent.begin();
2638  iter!=tmpcurrent.end();iter++)
2639  {
2640  mon tmpmon=(*iter)*mon(i,1);
2641  tmpmon.SetCoeff(1);
2642  current.push_back(tmpmon);
2643  }
2644  }
2645  //maintenant on tri;
2646  // mdebug()<<"tri";
2647  // mdebug()<<current.size();
2648  for(typename list<mon>::iterator iter=current.begin();iter!=current.end();iter++)
2649  if(!member(stockmon,*iter))
2650  {
2651  stockmon.push_back(*iter);
2652  }
2653  else
2654  {
2655  current.erase(iter--);
2656  }
2657  }
2658  //mdebug()<<"ttttttttt "<<stockmon.size();
2659  return;
2660 }
2661 
2662 template <class M> struct Dlexpolprod
2663 {
2664  static bool less (const M &, const M &);
2665  bool operator() (const M &, const M &) const;
2666 };
2667 template <class M>
2668 bool Dlexpolprod<M>::operator()(const M & m1, const M & m2) const {
2669  return less(m1,m2);
2670 }
2671 
2672 template <class M>
2673 inline bool Dlexpolprod<M>::less(const M & p1, const M & p2){
2674  typedef typename M::monom_t mon;
2675  mon m1=*(p1.begin());
2676  mon m2=*(p2.begin());
2677  int dm1=0,dm2=0,firstdiff=0,tmpdeg1,tmpdeg2;
2678  //int dm1=Degree(m1),dm2=Degree(m2);
2679  int n1 = lvar(m1.rep), n2 =lvar(m2.rep);
2680  int n = max(n1,n2);
2681  if(n1 <0 && n2>0) return false;
2682  if(n2 <0 && n1>0) return true;
2683  for(int i=0;i<=n;i++)
2684  {
2685  tmpdeg1=m1.GetDegree(i);
2686  tmpdeg2=m2.GetDegree(i);
2687  if(tmpdeg1>tmpdeg2 && firstdiff==0) firstdiff=1;
2688  if(tmpdeg2>tmpdeg1 && firstdiff==0) firstdiff=-1;
2689  dm1+=tmpdeg1;
2690  dm2+=tmpdeg2;
2691  }
2692  if (dm1 < dm2) return false;
2693  else if (dm1 > dm2) return true;
2694  return (firstdiff==1);
2695 }
2696 
2697 template<typename typprod,typename Base, typename typstock, typename typdump>
2698 void computeprod(typprod &prod,typstock &stockmon, Base & b,
2699  typdump & dump)
2700 {
2701  //
2702  for(typename typstock::iterator i=stockmon.begin();
2703  i!=stockmon.end();i++)
2704  {
2705  for(typename typstock::iterator j=i;j!=stockmon.end();j++)
2706  prod.push_back(*i * *j);
2707  }
2708  //Maintenant on tri et on ne selectionne que ceux qui sont en dehors de B
2709  //
2710  //Aie ici il faut define DLEX sur des polynomes!
2711  prod.sort(Dlexpolprod<typename typprod::value_type>());
2712  prod.unique();
2713  // for(typename typprod::iterator i=prod.begin();i!=prod.end();i++)
2714  // if (IsinB(*(i->begin()),b)) prod.erase(i--);
2715 }
2716 
2717 
2718 template<typename T>
2719 int dgesvd_(char *jobu, char *jobvt, long *m, long *n,
2720  T *a, long *lda, T *s, T *u, long *
2721  ldu, T *vt, long *ldvt, T *work, long *lwork,
2722  long *info);
2723 
2724 template<typename typP, typename typdump,typename Base, typename typserv>
2725 void generickernel(typP &newpol, int k,
2726  const typdump &dump, Base &b, typserv &serv)
2727 {
2728  using namespace coinor;
2729  typedef typename typP::value_type::order_t ord;
2730  //
2731  //Allocation du problem et des
2732  //contraintes
2733  struct blockmatrix<COEFF> C;
2734  struct sparseblock<COEFF> *blockptr;
2735  COEFF *bb;
2736  int i;
2737  struct constraintmatrix<COEFF> *constraints;
2738 
2739  //
2740  //On introduit un variable par monome de degre d en dehors de B
2741  //On va contruire l'operateur de de Hankel associ\E9
2742  //
2743 
2744  //
2745  //On commence par enumerer les monomes de degre plus petit que k
2746  //
2747  //
2748  typP stockprod,stockproj;//dans stockprod on a des polynomes
2749  //
2750  //
2751  //
2752  int nbconstraint=0;
2753  list<typename typserv::monom_t> stockmon,stockksur2;
2754  map<mon,int,ord> mapintligne,mapintcolonne,maphankel;
2755  enumeratemon(stockmon,k,b,dump);
2756  //
2757  //ICI on cree un index pour le monomes de degre k
2758  //il n'a a priori rien a voir avec int2mon/mon2int du reste du code
2759  //sic.....
2760  enumeratemon(stockksur2,k/2,b,dump);
2761  //
2762  //pour le suite il nous faut les monomes de B de degre <=k/2
2763  //
2764 
2765  for(typename list<typename typserv::monom_t>::iterator
2766  iter=stockksur2.begin(); iter!=stockksur2.end();iter++)
2767  if(!IsinB(*iter,b)) stockksur2.erase(iter--);
2768 #ifdef DEBUGSDP
2769  mdebug()<<"monome en degre k/2 ";
2770  for(typename list<typename typserv::monom_t>::iterator
2771  iter=stockksur2.begin(); iter!=stockksur2.end();iter++)
2772  mdebug()<<*iter<<" ";
2773  mdebug();
2774 #endif
2775  i=1;
2776  for(typename list<typename typserv::monom_t>::iterator
2777  iter=stockksur2.begin();iter!=stockksur2.end();iter++,i++)
2778  {
2779  int j=1;
2780  for(typename list<typename typserv::monom_t>::iterator
2781  iter2=stockksur2.begin();iter2!=stockksur2.end();iter2++,j++)
2782  {
2783  typename typserv::monom_t mono=*iter * (*iter2);
2784  // mdebug()<<"monome "<<mono<<" "<<i<<" "<<j;
2785  if (!mapintligne.count(mono))
2786  {
2787  mapintligne[mono]=i;
2788  mapintcolonne[mono]=j;
2789  }
2790  }
2791  }
2792  //mdebug()<<"pass\E9";
2793 
2794  //
2795  // On genere les formes normales des produits de monomes de
2796  // de tous les monomes de degre <=k/2 dez B.
2797  //
2798 
2799  computeprod(stockprod,stockksur2,b,dump);
2800  stockproj=stockprod;
2801  proj(stockproj,dump,b,serv);
2802 
2803 #ifdef DEBUGSDP
2804  mdebug()<<"produits ";
2805  for(typename typP::iterator iter=stockprod.begin();iter!=stockprod.end();
2806  iter++)
2807  mdebug()<<*iter<<" ";
2808  mdebug();
2809  mdebug()<<"proj ";
2810  for(typename typP::iterator iter=stockproj.begin();iter!=stockproj.end();
2811  iter++)
2812  mdebug()<<*iter<<" ";
2813  mdebug();
2814 #endif
2815  //
2816  //
2817  //computerelations();//
2818  //
2819  //
2820  //
2821  //
2822  // \A7Ici on va avoir besoin d'une structure contenant un monome de tete
2823  // hors de B et sa forme normal
2824 
2825  for(typename typP::iterator iter=stockprod.begin();
2826  iter!=stockprod.end();iter++)
2827  if(!IsinB(*(iter->begin()),b))
2828  nbconstraint++;
2829 #ifndef NONEW
2830  constraints=new constraintmatrix<COEFF>[nbconstraint+1];
2831 #else
2832  constraints=(struct constraintmatrix<COEFF> *)
2833  malloc((nbconstraint+1)*sizeof(struct constraintmatrix<COEFF>));
2834 #endif
2835  /*
2836  *
2837  * Les indices partent de 1 pas de 0!
2838  *
2839  */
2840  //
2841  //On construit la liste des contraintes de la fin vers le debut
2842  //
2843  //Ca c la fin
2844  constraints[0].blocks=NULL;//not used dixit the manual
2845  //
2846  i=1;
2847  typename typP::iterator iterstockprod=stockprod.begin();
2848  for(typename typP::iterator iter=stockproj.begin();
2849  iter!=stockproj.end();iter++,iterstockprod++)
2850  {
2851  if(!IsinB(*(iterstockprod->begin()),b))
2852  {
2853  int k;
2854  int decal;
2855  //
2856  //On construit la liste des contraintes de la fin vers le debut
2857  //
2858  //Ca c la fin
2859  constraints[i].blocks=NULL;
2860  //
2861 #ifndef NONEW
2862  blockptr=new sparseblock<COEFF>[1];
2863 #else
2864  blockptr=(struct sparseblock<COEFF> *)malloc(sizeof(struct sparseblock<COEFF>));
2865 #endif
2866  blockptr->blocknum=1;//
2867  // ?????? relations.size()-i;
2868 #if 0
2869  blockptr->blocksize=stockksur2.size();
2870  blockptr->numentries=2*iter->size()+2;//pour l'ind
2871  decal=iter->size()+1;//pour l'ind
2872 #else
2873  blockptr->blocksize=stockksur2.size();
2874  blockptr->numentries=iter->size()+1;//pour l'ind
2875  decal=0;
2876 #endif
2877  blockptr->constraintnum=i;
2878  blockptr->next=NULL;
2879  blockptr->nextbyblock=NULL;
2880 #ifndef NONEW
2881  blockptr->entries=new COEFF[(blockptr->numentries+1)];
2882 #else
2883  blockptr->entries=(COEFF *)
2884  malloc((blockptr->numentries+1)*sizeof(COEFF));
2885 #endif
2886 #ifndef NOSHORTS
2887  blockptr->iindices=(short unsigned int *)
2888  malloc((blockptr->numentries+1)*sizeof(short unsigned int));
2889  blockptr->jindices=(short unsigned int *) malloc((blockptr->numentries+1)*sizeof(short unsigned int));
2890 #else
2891  blockptr->iindices=( int *)
2892  malloc((blockptr->numentries+1)*sizeof( int));
2893  blockptr->jindices=(int *) malloc((blockptr->numentries+1)*sizeof( int));
2894 #endif
2895 
2896  k=1;
2897  // mdebug()<<" matrix contrainte "<<endl<<*iter;
2898  for(typename typP::value_type::iterator iterpol=iter->begin();
2899  iterpol!=iter->end();k++,iterpol++)
2900  {
2901  if(mapintligne[*iterpol]<mapintcolonne[*iterpol])
2902  {
2903  blockptr->iindices[k]=mapintligne[*iterpol];
2904  blockptr->jindices[k]=mapintcolonne[*iterpol];
2905  blockptr->entries[k]=((COEFF)-1)*iterpol->GetCoeff();
2906  // mdebug()<<blockptr->iindices[k]<<" "<<blockptr->jindices[k]
2907  // <<" "<<blockptr->entries[k];
2908  //
2909  }
2910  else
2911  {
2912  //blockptr->numentries--;
2913  //decal--;
2914  if (mapintligne[*iterpol]!=mapintcolonne[*iterpol] )
2915  {
2916  blockptr->jindices[k]=mapintligne[*iterpol];
2917  blockptr->iindices[k]=mapintcolonne[*iterpol];
2918  blockptr->entries[k]=((COEFF)-1)*iterpol->GetCoeff();
2919  }
2920  else
2921  {
2922  blockptr->jindices[k]=mapintligne[*iterpol];
2923  blockptr->iindices[k]=mapintcolonne[*iterpol];
2924  blockptr->entries[k]=((COEFF)-2)*iterpol->GetCoeff();
2925  }
2926  }
2927 
2928  }
2929 
2930  // mdebug()<<"index en plus";
2931  if(mapintligne[*(iterstockprod->begin())]<
2932  mapintcolonne[*(iterstockprod->begin())] )
2933  {
2934  blockptr->iindices[k]=mapintligne[*(iterstockprod->begin())];
2935  blockptr->jindices[k]=mapintcolonne[*(iterstockprod->begin())];
2936  blockptr->entries[k]=1;
2937  // mdebug()<<blockptr->iindices[k]<<" "<<blockptr->jindices[k]
2938  // <<" "<<blockptr->entries[k];
2939  }
2940  else
2941  {
2942  if (mapintligne[*(iterstockprod->begin())]!=
2943  mapintcolonne[*(iterstockprod->begin())])
2944  {
2945  blockptr->iindices[k]
2946  =mapintcolonne[*(iterstockprod->begin())];
2947 
2948  blockptr->jindices[k]
2949  =mapintligne[*(iterstockprod->begin())];
2950  blockptr->entries[k]=1;
2951  // mdebug()<<"indices ici "<<mapintligne[*(iterstockprod->begin())]<<" "<<
2952  // mapintcolonne[*(iterstockprod->begin())];
2953  }
2954  else
2955  {
2956  blockptr->iindices[k]=mapintligne[*(iterstockprod->begin())];
2957  blockptr->jindices[k]=mapintcolonne[*(iterstockprod->begin())];
2958  blockptr->entries[k]=2;
2959  }
2960  }
2961  /*
2962  * Insert block into the linked list of A1 blocks.
2963  */
2964  // blockptr->blocksize--;
2965  blockptr->next=NULL;
2966  //=constraints[i].blocks;
2967  constraints[i].blocks=blockptr;
2968  // mdebug()<<"i "<<i<<" "<<stockprod.size();
2969  i++;
2970 
2971  }
2972  }
2973  //
2974  // On donne les relations Hankel
2975  //
2976  //
2977  //
2978  //
2979  //
2980  //
2981  int nbhankelrel=0,i_counter=1,j_counter;
2982  for(typename list<typename typserv::monom_t>::iterator i=stockksur2.begin();
2983  i!=stockksur2.end();i++,i_counter++)
2984  {
2985  j_counter=1;
2986  for(typename list<typename typserv::monom_t>::iterator j=stockksur2.begin();
2987  j!=stockksur2.end();j++,j_counter++)
2988  if(!maphankel.count(*i * *j))
2989  {
2990  maphankel[*i * *j]=j_counter*stockksur2.size()+i_counter;
2991  }
2992  else
2993  if((maphankel[*i * *j]%stockksur2.size()!=j_counter)
2994  ||(maphankel[*i * *j]/stockksur2.size()!=i_counter))
2995  {
2996  int decal=0;
2997  nbhankelrel++;
2998 #ifndef NONEW
2999  constraintmatrix<COEFF> *tmpconstraints=new constraintmatrix<COEFF>[nbhankelrel+nbconstraint+1];
3000  for(int ii=0;ii<nbhankelrel+nbconstraint;ii++)
3001  tmpconstraints[ii]=constraints[ii];
3002  delete[] constraints;
3003  constraints=tmpconstraints;
3004  blockptr=new sparseblock<COEFF>[1];
3005 #else
3006  constraints=(struct constraintmatrix<COEFF> *)
3007  realloc(constraints,
3008  (nbhankelrel+nbconstraint+1)
3009  *sizeof(struct constraintmatrix<COEFF>));
3010 
3011  blockptr=(struct sparseblock<COEFF> *)malloc(sizeof(struct sparseblock<COEFF>));
3012 #endif
3013  blockptr->blocknum=1;//
3014  blockptr->blocksize=stockksur2.size();
3015 #if 0
3016  blockptr->numentries=4;
3017 #else
3018  blockptr->numentries=2;
3019 #endif
3020  blockptr->constraintnum=nbconstraint+nbhankelrel;
3021  blockptr->next=NULL;
3022  blockptr->nextbyblock=NULL;
3023 #ifndef NONEW
3024  blockptr->entries=new COEFF[(blockptr->numentries+1)];
3025 #else
3026  blockptr->entries=(COEFF *)
3027  malloc((blockptr->numentries+1)*sizeof(COEFF));
3028 #endif
3029 #ifndef NOSHORTS
3030  blockptr->iindices=(short unsigned int *) malloc((blockptr->numentries+1)*sizeof(short unsigned int));
3031  blockptr->jindices=(short unsigned int *) malloc((blockptr->numentries+1)*sizeof(short unsigned int));
3032 #else
3033  blockptr->iindices=(int *) malloc((blockptr->numentries+1)*sizeof(int));
3034  blockptr->jindices=(int *) malloc((blockptr->numentries+1)*sizeof(int));
3035 #endif
3036  //blockptr->numentries=blockptr->blocksize-1; //On a rajoute un avant
3037  if((maphankel[*i * *j]%stockksur2.size())<
3038  (maphankel[*i * *j]/stockksur2.size()))
3039  {
3040  blockptr->iindices[1]=maphankel[*i * *j]%stockksur2.size();
3041  blockptr->jindices[1]=maphankel[*i * *j]/stockksur2.size();
3042  blockptr->entries[1]=1;
3043  }
3044  else
3045  {
3046  if((maphankel[*i * *j]/stockksur2.size())!=
3047  (maphankel[*i * *j]%stockksur2.size()))
3048  {
3049  blockptr->jindices[1]=maphankel[*i * *j]%stockksur2.size();
3050  blockptr->iindices[1]=maphankel[*i * *j]/stockksur2.size();
3051  blockptr->entries[1]=1;
3052  }
3053  else
3054  {
3055  blockptr->jindices[1]=maphankel[*i * *j]%stockksur2.size();
3056  blockptr->iindices[1]=maphankel[*i * *j]/stockksur2.size();
3057  blockptr->entries[1]=2;
3058  }
3059 
3060  }
3061  if(i_counter<j_counter)
3062  {
3063  blockptr->iindices[2]=i_counter;
3064  blockptr->jindices[2]=j_counter;
3065  blockptr->entries[2]=-1;
3066  }
3067  else
3068  {
3069  if(i_counter!=j_counter)
3070  {
3071  blockptr->iindices[2]=j_counter;
3072  blockptr->jindices[2]=i_counter;
3073  blockptr->entries[2]=-1;
3074  }
3075  else
3076  {
3077  blockptr->iindices[2]=j_counter;
3078  blockptr->jindices[2]=i_counter;
3079  blockptr->entries[2]=-2;
3080  }
3081  }
3082 #if 0
3083  if((maphankel[*i * *j]/stockksur2.size())!=
3084  (maphankel[*i * *j]%stockksur2.size()))
3085  {
3086  blockptr->iindices[3]=maphankel[*i * *j]/stockksur2.size();
3087  blockptr->jindices[3]=maphankel[*i * *j]%stockksur2.size();;
3088  blockptr->entries[3]=1;
3089  }
3090  else
3091  {
3092  decal++;
3093  blockptr->numentries--;
3094  blockptr->entries[1]=2;
3095  }
3096 
3097  if(i_counter!=j_counter)
3098  {
3099  blockptr->iindices[4-decal]=j_counter;
3100  blockptr->jindices[4-decal]=i_counter;
3101  blockptr->entries[4-decal]=-1;
3102  }
3103  else
3104  {
3105  blockptr->numentries--;
3106  blockptr->entries[2]=-2;
3107 
3108  }
3109 #endif
3110  // mdebug()<<"indice "<<maphankel[*i * *j]%stockksur2.size()
3111  // <<" "<<maphankel[*i * *j]/stockksur2.size();
3112  // mdebug()<<"indice "<<i_counter<<" "<<j_counter;
3113  // mdebug();
3114  // mdebug()<<"indice "<<l<<" "<<1;
3115  //mdebug()<<"indice "<<l-k<<" "<<k+1;
3116  /*
3117  * Insert block into the linked list of A1 blocks.
3118  */
3119 
3120  blockptr->next=NULL;//constraints[i].blocks;
3121  constraints[nbconstraint+nbhankelrel].blocks=blockptr;
3122 
3123 
3124  }
3125  }
3126  //
3127  //
3128  //
3129  // les relations de symetrie sont implicites X est recherche SDP
3130  // ONrajoute la contrainte sur le coeff en haut a gauche
3131  //
3132  //
3133  nbhankelrel++;
3134 #ifndef NONEW
3135  constraintmatrix<COEFF> *tmpconstraints=
3136  new constraintmatrix<COEFF>[(nbhankelrel+nbconstraint+1)];
3137  for(int ii=0;ii<nbhankelrel+nbconstraint;ii++)
3138  tmpconstraints[ii]=constraints[ii];
3139  delete [] constraints;
3140  constraints=tmpconstraints;
3141 #else
3142  constraints=(struct constraintmatrix<COEFF> *)
3143  realloc(constraints,
3144  (nbhankelrel+nbconstraint+1)
3145  *sizeof(struct constraintmatrix<COEFF>));
3146 #endif
3147 
3148 #ifndef NONEW
3149  blockptr=new sparseblock<COEFF>[1];
3150 #else
3151  blockptr=(struct sparseblock<COEFF> *)malloc(sizeof(struct sparseblock<COEFF>));
3152 #endif
3153  blockptr->blocknum=1;//
3154  blockptr->blocksize=stockksur2.size();
3155  blockptr->numentries=1;
3156  blockptr->constraintnum=nbconstraint+nbhankelrel;
3157  blockptr->next=NULL;
3158  blockptr->nextbyblock=NULL;
3159 #ifndef NONEW
3160  blockptr->entries=new COEFF[blockptr->numentries+1];
3161 #else
3162  blockptr->entries=(COEFF *)
3163  malloc((blockptr->numentries+1)*sizeof(COEFF));
3164 #endif
3165 #ifndef NOSHORTS
3166  blockptr->iindices=(short unsigned int *) malloc((blockptr->numentries+1)*sizeof(short unsigned int));
3167  blockptr->jindices=(short unsigned int *) malloc((blockptr->numentries+1)*sizeof(short unsigned int));
3168 #else
3169  blockptr->iindices=(int *) malloc((blockptr->numentries+1)*sizeof(int));
3170  blockptr->jindices=(int *) malloc((blockptr->numentries+1)*sizeof(int));
3171 #endif
3172  //blockptr->numentries=blockptr->blocksize-1; //On a rajoute un avant
3173 
3174  blockptr->iindices[1]=1;
3175  blockptr->jindices[1]=1;
3176  blockptr->entries[1]=1;
3177  blockptr->next=NULL;//constraints[i].blocks;
3178  constraints[nbconstraint+nbhankelrel].blocks=blockptr;
3179  //
3180  //On optimise la fonction "0" sous l es contraintes lineaires precedents
3181  //
3182  struct blockmatrix<COEFF> X,Z;
3183  COEFF *y;
3184  COEFF pobj,dobj;
3185  int tmpsize=(int)stockksur2.size();
3186 #ifndef NONEW
3187  bb=new COEFF[nbhankelrel+nbconstraint+1];
3188 #else
3189  bb=(COEFF *)malloc((nbhankelrel+nbconstraint+1)*sizeof(COEFF));
3190 #endif //oon a un block de 0
3191  C.nblocks=1;
3192 #ifndef NONEW
3193  C.blocks=new blockrec<COEFF>[2];
3194 #else
3195  C.blocks=(struct blockrec<COEFF> *)malloc(2*sizeof(struct blockrec<COEFF>));
3196 #endif
3197 #ifndef CCCCC
3198  C.blocks[1].blockcategory=MATRIX;
3199  C.blocks[1].blocksize=stockksur2.size();
3200 #ifndef NONEW
3201  C.blocks[1].data.mat=new COEFF[stockksur2.size()*stockksur2.size()+1];
3202 #else
3203  C.blocks[1].data.mat=(COEFF *)malloc(
3204  (stockksur2.size()*stockksur2.size()+1)*sizeof(COEFF));
3205 #endif
3206  for(int k=0;k<=(stockksur2.size()*stockksur2.size());k++)
3207  C.blocks[1].data.mat[k]=0;
3208  //
3209  C.blocks[1].data.mat[1]=0;
3210 #else
3211  C.blocks[1].blockcategory=DIAG;
3212  C.blocks[1].blocksize=stockksur2.size();
3213  C.blocks[1].data.vec=(COEFF *)malloc(
3214  (stockksur2.size()+1)*sizeof(COEFF));
3215  for(int k=0;k<=(stockksur2.size());k++)
3216  C.blocks[1].data.vec[k]=0;
3217 
3218 #endif
3219  for(i=0;i<=nbhankelrel+nbconstraint;i++)
3220  bb[i]=0;
3221  bb[nbconstraint+nbhankelrel]=(COEFF)1.;//on double qd on est sur la diagonale
3222  //mdebug()<<"stockprod.size() "<<stockprod.size();
3223  // mdebug()<<"nbconstraint "<<nbconstraint;
3224  // mdebug()<<"tmpsize "<<tmpsize;
3225  // mdebug()<<"Affichage contraintes";
3226  for(i=1;i<nbconstraint+nbhankelrel+1;i++)
3227  {
3228  // mdebug()<<"le b associe "<<bb[i];
3229  int tsize=constraints[i].blocks->blocksize;
3230  COEFF* tmp = new COEFF[tsize*tsize];
3231  // mdebug()<<"taille "<<tsize;
3232  for(int j;j<tsize*tsize;j++)
3233  tmp[j]=0;
3234  for(int j=1;j<=constraints[i].blocks->numentries;j++)
3235  tmp[constraints[i].blocks->iindices[j]-tsize-1+
3236  constraints[i].blocks->jindices[j]*tsize]=
3237  constraints[i].blocks->entries[j];
3238 #ifdef DEBUGSDP
3239  for(int j=0;j<tsize;j++)
3240  {
3241  for(int jj=0;jj<tsize;jj++)
3242  mdebug()<<tmp[tsize*(jj)+j]<<" ";
3243  mdebug();
3244  }
3245  mdebug();
3246 #endif
3247  }
3248  char nom[50];
3249  static int d=0;
3250  sprintf(nom,"ficprob_%d",d++);
3251  write_prob(nom,tmpsize,
3252  nbhankelrel+nbconstraint,C,bb,constraints);
3253  //Precision(53);
3254  mdebug()<<"SDP of size : "<<stockksur2.size()<<"x"<<nbhankelrel+
3255  nbconstraint;
3256  initsoln(tmpsize,nbhankelrel+nbconstraint,
3257  C,bb,constraints,&X,&y,&Z);
3258  struct timeval tvini, tvfin;
3259  gettimeofday(&tvini,NULL);
3260  int return_code=easy_sdp((int)stockksur2.size(),nbhankelrel+
3261  nbconstraint,
3262  C,bb,constraints,(COEFF)(0.0),&X,&y,&Z,&pobj,&dobj);
3263  //Precision(350);
3264  gettimeofday(&tvfin,NULL);
3265  tempsCDSP+=tvfin.tv_sec-tvini.tv_sec+1e-6*(tvfin.tv_usec-tvini.tv_usec);
3266  // write_sol("prob.sol",tmpsize,nbhankelrel+nbconstraint,X,y,Z);
3267  //
3268  //
3269  //
3270  // A ce niveau la il faut verifier le rang
3271  // Le cas echeant recuperer le kernel
3272  //
3273  mdebug()<<"pobj "<<pobj<<" dobj "<<dobj<<" diff "<<pobj-dobj;
3274  // mdebug()<<"return code "<<return_code;
3275  if (1||return_code==8)//c le code de retour indiquant une singularite
3276  {
3277  long n=stockksur2.size();
3278  COEFF *matdense=new COEFF[n*n];
3279  COEFF *V=new COEFF[n*n];
3280  COEFF *S=new COEFF[n];
3281  COEFF *work=new COEFF[16*n];
3282  long lwork=16*n,info;
3283  long k=0;//servira pour le rang0
3284  for(int j=0;j<n*n;j++)
3285  matdense[j]=0;
3286  /*
3287  Warning !!!
3288  je ne traite pas les differetn cas des blocks de X
3289  il peuvent etre soit DIAG soit MATRIX
3290  soit PACKEDMATRIX (on peut s'en passez il y a invalid dasn writesol
3291  */
3292 
3293  //
3294  //on est pas de rang plein
3295  //
3296  //il faut calculer le noyau
3297  //
3298  //
3299 #ifdef DEBUGSDP
3300  for(int i=1;i<=nbconstraint+nbhankelrel;i++)
3301  mdebug()<<"y "<<i<<" "<<y[i];
3302  mdebug()<<"nombre de blocks de X "<<X.nblocks;
3303 #endif
3304  for(int i=1;i<=X.nblocks;i++)
3305  {
3306  //Pour chque block il faut regarder son noyau
3307  struct blockrec<COEFF> tmp=X.blocks[i];
3308  //
3309  //On recopie le bloc en dense
3310  //
3311 #ifdef DEBUGSDP
3312  mdebug()<<"n _________________________________ "<<n;
3313  mdebug()<<"blockcat "<<tmp.blockcategory<<" "<<k;
3314 #endif
3315  switch(tmp.blockcategory)
3316  {
3317  case DIAG:
3318  for(int j=1;j<=tmp.blocksize;j++)
3319  matdense[n*(k+j-1)+(k+j-1)]=tmp.data.vec[j];
3320 #ifdef DEBUGSDP
3321  for(int j=1;j<=tmp.blocksize;j++)
3322  {
3323  for(int jj=1;jj<=tmp.blocksize;jj++)
3324  mdebug()<<matdense[n*(k+jj-1)+j-1]<<" ";
3325  mdebug();
3326  }
3327  mdebug();
3328 #endif
3329  break;
3330  case MATRIX:
3331  for(int j=1;j<=tmp.blocksize;j++)
3332  for(int jj=1;jj<=tmp.blocksize;jj++)
3333  matdense[n*(k+jj-1)+j-1]=
3334  tmp.data.mat[tmp.blocksize*(jj-1)+j-1];
3335 #ifdef DEBUGSDP
3336  for(int j=1;j<=tmp.blocksize;j++)
3337  {
3338  for(int jj=1;jj<=tmp.blocksize;jj++)
3339  mdebug()<<matdense[n*(k+jj-1)+j-1]<<" ";
3340  mdebug();
3341  }
3342 #endif
3343  }
3344  k+=tmp.blocksize;
3345 
3346  }
3347  //
3348  //On calcul le noyau
3349  //
3350  dgesvd_((char*)"N",(char*)"A",&n,&n,matdense,
3351  &n,S,(COEFF*)NULL,&n,V,&n,work,&lwork,&info);
3352 #if 1
3353  mdebug()<<"info svd "<<info;
3354  mdebug()<<"val singulieres ";
3355  for(k=0;k<n;k++)
3356  mdebug()<<S[k]<<" ";
3357  mdebug();
3358 #endif
3359  for(k=1;k<n;k++)
3360  //if(S[0]>S[k]*1e6) break;//MAGIC
3361  if((abs(S[k-1])==0) || (abs(S[k-1])>abs(S[k])*1e6)) break;//MAGIC
3362  //#ifdef DEBUGSDP
3363  mdebug()<<"----k ici "<<k<<" et n "<<n;
3364  //#endif
3365  for(int j=1;j<=n;j++)
3366  {
3367  for(int jj=1;jj<=n;jj++)
3368  mdebug()<<V[n*(jj-1)+j-1]<<" ";
3369  mdebug();
3370  }
3371  //
3372  //On construit les nouvelles equations
3373  //
3374  for(;k<n;k++)
3375  {
3376 #ifdef DEBUGSDP
3377  mdebug()<<"equation dans le noyau";
3378 #endif
3379  typedef typename typP::value_type polyalp;
3380  polyalp tmp(0);
3381  for(typename list<typename typserv::monom_t>::iterator
3382  iter=stockksur2.begin();
3383  iter!=stockksur2.end();iter++)
3384  {
3385 
3386  //mdebug()<<"dans construct eq "<<k<<" map "<<mapintcolonne[*iter];
3387  //mdebug()<<V[k+n*(mapintcolonne[*iter]-1)]<<" "<<*iter<<" ";
3388  if(!Iszero(abs(V[k+n*(mapintcolonne[*iter]-1)])))
3389 
3390  {
3391  // mdebug()<<"tmp avant"<<endl<<tmp;
3392  //mdebug()<<"ajout "<<endl<<polyalp(*iter*V[k+n*(mapintcolonne[*iter]-1)]);
3393  tmp+=polyalp(*iter*(COEFF)V[k+n*(mapintcolonne[*iter]-1)]);
3394  //mdebug()<<"tmp aapre"<<endl<<tmp;
3395  }
3396  //iter->first;// est le monome
3397  //iter->second;//est l'entier associ\E9
3398  }
3399  // mdebug();
3400  // mdebug()<<"g trouve ca comme nvlle eq "<<endl<<tmp;
3401  newpol.push_back(tmp);
3402 
3403  }
3404  delete[] S;
3405  delete[] V;
3406  delete[] work;
3407  delete[] matdense;
3408  }
3409  for(int i=1;i<nbconstraint;i++)
3410  {
3411  free(constraints[i].blocks->iindices);
3412  free(constraints[i].blocks->jindices);
3413  delete [] constraints[i].blocks->entries;
3414  delete [] constraints[i].blocks;
3415  }
3416  delete [] constraints;
3417  delete [] bb;
3418  delete [] C.blocks[1].data.mat;
3419  delete [] C.blocks;
3420  // il faut desallouer les contraintes mieux (les matrices pointees
3421  //il faut desalluer C
3422  //Precision(53);
3423  //mdebug()<<"en sortie "<<newpol.empty();
3424 }
3425 
3426 //#include"genekersdpa.hpp"
3427 
3428 template<typename typmatmap, typename typserv>
3429 void write_to_sdpa_fic(typmatmap &matrixmap, int &dim, typserv &serv,
3430  int d)
3431 {
3432  ofstream tmpfic;
3433  typename typserv::monom_t m(1);
3434  COEFF *tmpval=matrixmap[m];
3435  char nomdufic[50];
3436  sprintf(nomdufic,"tmpfile_%d",d);
3437  tmpfic.open(nomdufic);
3438  tmpfic<<setprecision(80);
3439  tmpfic<<"\" Ceci est un message \"";
3440  tmpfic<<matrixmap.size()-1<<" = mDIM";
3441  tmpfic<<1<<" = nBLOCK";
3442  tmpfic<<dim<<" = bLOCKsTRUCT";
3443  tmpfic<<"{";
3444  for(int i=0;i<((long)matrixmap.size())-2;i++)
3445  tmpfic<<"0,";
3446  tmpfic<<"0}";
3447  tmpfic<<"{";
3448  for(int i=0;i<dim-1;i++)
3449  {
3450  tmpfic<<"{";
3451  for(int j=0;j<dim-1;j++)
3452  {
3453  //mdebug()<<(-1*matrixmap[m][i*dim+j])<<" ";
3454  tmpfic<<(COEFF)(-1*matrixmap[m][i*dim+j])<<",";
3455  }
3456  mdebug();
3457  tmpfic<<(COEFF)(-1*matrixmap[m][i*dim+dim-1])<<"},";
3458  }
3459  tmpfic<<"{";
3460 
3461  for(int j=0;j<dim-1;j++)
3462  tmpfic<<-1*matrixmap[m][(dim-1)*dim+j]<<",";
3463  tmpfic<<-1*matrixmap[m][dim*dim-1]<<"}}";
3464 
3465  matrixmap.erase(matrixmap.find(m));
3466  for(typename typmatmap::iterator itermat=matrixmap.begin();
3467  itermat!=matrixmap.end();itermat++)
3468  {
3469  //mdebug()<<"matrice de "<<itermat->first;
3470  COEFF *tt=itermat->second;
3471  tmpfic<<"{";
3472  for(int i=0;i<dim-1;i++)
3473  {
3474  tmpfic<<"{";
3475  for(int j=0;j<dim-1;j++)
3476  tmpfic<<tt[i*dim+j]<<",";
3477  tmpfic<<tt[i*dim+dim-1]<<"},";
3478  }
3479  tmpfic<<"{";
3480 
3481  for(int j=0;j<dim-1;j++)
3482  tmpfic<<tt[(dim-1)*dim+j]<<",";
3483  tmpfic<<tt[dim*dim-1]<<"}}";
3484 
3485  }
3486  matrixmap[m]=tmpval;
3487  tmpfic.close();
3488 }
3489 
3490 template<typename typmatmap, typename coeff, typename typserv>
3491 void read_sol_from_sdpa(coeff *&sol, typmatmap &matmap,int dim, typserv &serv)
3492 {
3493  typedef typename typserv::monom_t mon;
3494  ifstream tmpfic;
3495  string ligne,ligne2;
3496  typename string::iterator debut;
3497  int i;
3498  coeff *xVec,*xMat;
3499  stringstream tmpstr;
3500  system(
3501  "perl -i.bak -0777 -ne 'if (/^xMat *= *(.*)(?=yMat)/ms) {$toto=$1; $toto=~s/[[:blank:]{}]\\+?//g; $toto=~s/,\\+?/ /g; print $toto}' lenom" // "perl -0777 -i.bak -ne '/^xMat *= *(.*)(?=yMat)/ms and $toto=$1 and $toto=~s/[[:blank:]{}]\\+?//g and $toto=~s/,\\+?/ /g and print $toto' lenom"
3502  );
3503  tmpfic.open("lenom");
3504 #if 0
3505  tmpfic>>ligne;//chaine Xvec
3506  mdebug()<<"ligne xVec :"<<ligne;
3507  tmpfic>>ligne;// les coordonn\E9es interessantes
3508  mdebug()<<"ligne coordq :"<<ligne;
3509  dim=count(ligne.begin(),ligne.end(),'i')+1;
3510  xVec=new coeff[dim];
3511  sol=new coeff[dim*dim];
3512  i=0;
3513  int tmpindex;
3514  for(string::iterator iter=ligne.begin()+1;iter!=ligne.end();
3515  iter=ligne.begin()+tmpindex+1)
3516  {
3517  stringstream tmpstr;
3518  mdebug()<<"iter "<<iter-ligne.begin();
3519  tmpindex=ligne.find(",",iter-ligne.begin());
3520  if (tmpindex==string::npos) tmpindex=ligne.size()-1;
3521  mdebug()<<"tmpindex "<<tmpindex;
3522  copy(iter,ligne.begin()+tmpindex,ostream_iterator<char>(tmpstr,""));
3523  tmpstr>>ligne2;
3524  mdebug()<<"eeeee "<<ligne2<<" "<<ligne.find(",",iter-ligne.begin());
3525  tmpstr<<ligne2;
3526  tmpstr>>xVec[i++];
3527  mdebug()<<"i "<<i;
3528  }
3529 #endif
3530  //en theorie maintenant on a xVec
3531 #if 0
3532  tmpfic>>ligne;//chaine xMat
3533  mdebug()<<"ligne xmat :"<<ligne;
3534  tmpfic>>ligne;//chaine '{'
3535  for(i=0;i<dim;i++)
3536  {
3537  tmpfic>>ligne;// les coordonn\E9es interessantes
3538  int j=0;
3539  for(string::iterator iter=ligne.begin()+ligne.find_last_of('{');
3540  iter!=ligne.end();
3541  iter=ligne.begin()+ligne.find(",",iter-ligne.begin()))
3542  {
3543  stringstream tmpstr;
3544  // copy(iter,ligne.begin()+ligne.find(',',iter-ligne.begin())
3545  // ,ostream_iterator<char>(tmpstr,""));
3546  int toto=ligne.find(",",iter-ligne.begin());
3547  mdebug()<<"ligne "<<ligne;
3548  mdebug()<<"toto "<<toto<<" "<<string::npos;
3549  if(toto==string::npos) toto=ligne.size()-2;
3550  mdebug()<<"toto "<<toto<<" "<<string::npos;
3551  tmpstr<<ligne.substr(iter-ligne.begin(),toto);
3552  tmpstr>>xMat[dim*j++ + i];
3553  }
3554  }
3555 #endif
3556  //maintenant on a relu le xVec et le Xmat
3557  //Reconstruction de la solution
3558  sol=new coeff[dim*dim];
3559  for(int i=0;i<dim*dim;i++)
3560  {
3561  tmpfic>>sol[i];
3562  }
3563  tmpfic.close();
3564  mdebug()<<"Solution ";
3565  for(int i=0;i<dim*dim;i++)
3566  mdebug()<<sol[i]<<" ";
3567  mdebug();
3568  // copy(sol,sol+dim*dim,ostream_iterator<string>(cout," "));
3569 #if 0
3570  for(int i=0;i<dim*dim;i++)
3571  sol[i]+=matmap[mon(1)][i];
3572  matmap.erase(matmap.find(mon(1)));
3573  int j=0;
3574  for(typename typmatmap::iterator itermat=matmap.begin();
3575  itermat!=matmap.end();itermat++)
3576  {
3577  for(int i=0;i<dim*dim;i++)
3578  sol[i]-=xVec[j]*itermat->second[i];
3579  j++;
3580  }
3581 #endif
3582 }
3583 
3584 template<typename typP, typename typdump,typename Base, typename typserv>
3585 void generickernelsdpa(typP &newpol, int k,
3586  const typdump &dump, Base &b, typserv &serv)
3587 {
3588  typedef typename typP::value_type::order_t ord;
3589  typedef typename typserv::monom_t Monomial;
3590  //
3591  //Allocation du problem et des
3592  //contraintes
3593  COEFF * C;
3594  map<Monomial,COEFF*,ord> matrixmap;
3595  int i;
3596 
3597  //
3598  //On introduit un variable par monome de degre d en dehors de B
3599  //On va contruire l'operateur de de Hankel associ\E9
3600  //
3601 
3602  //
3603  //On commence par enumerer les monomes de degre plus petit que k
3604  //
3605  //
3606  typP stockprod,stockproj;//dans stockprod on a des polynomes
3607  //
3608  //
3609  //
3610  int nbconstraint=0;
3611  typedef list<typename typserv::monom_t> typstock;
3612  list<typename typserv::monom_t> stockmon,stockksur2;
3613  map<mon,int,ord> mapintligne,mapintcolonne,maphankel;
3614  enumeratemon(stockmon,k,b,dump);
3615  //
3616  //ICI on cree un index pour le monomes de degre k
3617  //il n'a a priori rien a voir avec int2mon/mon2int du reste du code
3618  //sic.....
3619  enumeratemon(stockksur2,k/2,b,dump);
3620  //
3621  //pour le suite il nous faut les monomes de B de degre <=k/2
3622  //
3623 
3624  for(typename list<typename typserv::monom_t>::iterator
3625  iter=stockksur2.begin(); iter!=stockksur2.end();iter++)
3626  if(!IsinB(*iter,b)) stockksur2.erase(iter--);
3627  i=1;
3628  for(typename list<typename typserv::monom_t>::iterator
3629  iter=stockksur2.begin();iter!=stockksur2.end();iter++,i++)
3630  {
3631  int j=1;
3632  for(typename list<typename typserv::monom_t>::iterator
3633  iter2=stockksur2.begin();iter2!=stockksur2.end();iter2++,j++)
3634  {
3635  typename typserv::monom_t mono=*iter * (*iter2);
3636  // mdebug()<<"monome "<<mono<<" "<<i<<" "<<j;
3637  if (!mapintligne.count(mono))
3638  {
3639  mapintligne[mono]=i;
3640  mapintcolonne[mono]=j;
3641  }
3642  }
3643  }
3644  //mdebug()<<"pass\E9";
3645 
3646  //
3647  // On genere les formes normales des produits de monomes de
3648  // de tous les monomes de degre <=k/2 dez B.
3649  //
3650  int dim=stockksur2.size();
3651  for(typename typstock::iterator i=stockmon.begin();
3652  i!=stockmon.end();i++)
3653  if(IsinB(*i,b))
3654  {
3655  mdebug()<<"ici "<<*i<<" dim "<<dim;
3656  matrixmap[*i]=new COEFF [dim*dim];
3657  for(int ii=0;ii<dim*dim;ii++)
3658  matrixmap[*i][ii]=0;
3659  }
3660 
3661 
3662  int ii=0;
3663  mdebug()<<"coiucou";
3664  for(typename typstock::iterator i=stockksur2.begin();
3665  i!=stockksur2.end();i++,ii++)
3666  {
3667  int jj=ii;
3668  mdebug()<<"coiucou";
3669  for(typename typstock::iterator j=i;j!=stockksur2.end();jj++,j++)
3670  {
3671 
3672  typename typserv::monom_t tmpmon=(*j)*(*i);
3673  //mdebug()<<"coiucoufffffff "<<tmpmon;
3674  if (IsinB(tmpmon,b))
3675  {
3676  matrixmap[tmpmon][ii*dim+jj]+=1;
3677  if(ii!=jj)
3678  matrixmap[tmpmon][jj*dim+ii]+=1;
3679  }
3680  else
3681  {
3682  typP tmplist;
3683  tmplist.push_back(tmpmon);
3684  proj(tmplist,dump,b,serv);
3685  // mdebug()<<"ici tmpmon "<<tmpmon;
3686  //mdebug()<<"tmplist.front "<<tmplist.front();
3687  for(typename typP::value_type::iterator iter=tmplist.front().begin();
3688  iter!=tmplist.front().end();iter++)
3689  {
3690  Monomial verytmpmon=*iter;
3691  verytmpmon.SetCoeff(1);
3692  //mdebug()<<"tmpmon "<<verytmpmon<<" "<<*iter;
3693  matrixmap[verytmpmon][ii*dim+jj]+=iter->GetCoeff();
3694  if(ii!=jj)
3695  matrixmap[verytmpmon][jj*dim+ii]+=iter->GetCoeff();
3696 
3697  }
3698  }
3699  }
3700  }
3701 
3702  if(dim>1)
3703  {
3704  struct timeval tvini, tvfin;
3705  COEFF *sol;
3706  static int d=0;
3707  static char nom[1024];
3708  gettimeofday(&tvini,NULL);
3709  mdebug()<<"coucou dim "<<dim;
3710  write_to_sdpa_fic(matrixmap,dim,serv,d);
3711  sprintf(nom,"sdpa_gmp -dd tmpfile_%d -o lenom",d++);
3712  system(nom);
3713  read_sol_from_sdpa(sol,matrixmap,dim,serv);
3714 
3715 
3716  gettimeofday(&tvfin,NULL);
3717  tempsCDSP+=tvfin.tv_sec-tvini.tv_sec+1e-6*(tvfin.tv_usec-tvini.tv_usec);
3718 
3719  //
3720  //mdebug()<<"coiucou";
3721  // A ce niveau la il faut verifier le rang
3722  // Le cas echeant recuperer le kernel
3723  //
3724  // mdebug()<<"return code "<<return_code;
3725  long n=stockksur2.size();
3726  COEFF *V=new COEFF[n*n];
3727  COEFF *S=new COEFF[n];
3728  COEFF *work=new COEFF[16*n];
3729  long lwork=16*n,info;
3730  long kk=0;//servira pour le rang 0
3731  dgesvd_((char*)"N",(char*)"A",&n,&n,sol,
3732  &n,S,(COEFF*)NULL,&n,V,&n,work,&lwork,&info);
3733 #if 1
3734  mdebug()<<"info svd "<<info;
3735  mdebug()<<"val sing";
3736  for(kk=0;kk<n;kk++)
3737  mdebug()<<S[kk]<<" ";
3738  mdebug();
3739 #endif
3740  mdebug()<<"coiucou";
3741  for(kk=1;kk<n;kk++)
3742  //if(S[0]>S[k]*1e6) break;//MAGIC
3743  if((abs(S[kk-1])==0) ||(abs(S[kk-1])>abs(S[kk])*1e6)) break;//MAGIC
3744 #if 1
3745  mdebug()<<"kk ici "<<kk<<" et n "<<n;
3746 #endif
3747  for(int j=0;j<n;j++)
3748  {
3749  for(int jj=0;jj<n;jj++)
3750  mdebug()<<V[n*(jj)+j]<<" ";
3751  mdebug();
3752  }
3753  //
3754  //On construit les nouvelles equations
3755  //
3756  for(;kk<(int)n;kk++)
3757  {
3758 #if 1
3759  mdebug()<<"equation dans le noyau";
3760 #endif
3761  typedef typename typP::value_type polyalp;
3762  polyalp tmp(0);
3763  for(typename list<typename typserv::monom_t>::iterator
3764  iter=stockksur2.begin();
3765  iter!=stockksur2.end();iter++)
3766  {
3767 
3768  //mdebug()<<"dans construct eq "<<k<<" map "<<mapintcolonne[*iter];
3769  //mdebug()<<V[k+n*(mapintcolonne[*iter]-1)]<<" "<<*iter<<" ";
3770  if(!Iszero(abs(V[kk+n*(mapintcolonne[*iter]-1)])))
3771 
3772  {
3773  // mdebug()<<"tmp avant"<<endl<<tmp;
3774  //mdebug()<<"ajout "<<endl<<polyalp(*iter*V[k+n*(mapintcolonne[*iter]-1)]);
3775  tmp+=polyalp(*iter*(COEFF)V[kk+n*(mapintcolonne[*iter]-1)]);
3776  //mdebug()<<"tmp aapre"<<endl<<tmp;
3777  }
3778  //iter->first;// est le monome
3779  //iter->second;//est l'entier associ\E9
3780  }
3781  // mdebug();
3782  mdebug()<<"g trouve ca comme nvlle eq \n "<<tmp;
3783  newpol.push_back(tmp);
3784 
3785  }
3786 
3787  delete[] S;
3788  delete[] V;
3789  delete[] work;
3790  delete[] sol;
3791  }
3792 }
3793 
3794 template<typename typdump,typename Base, typename typserv>
3795 int testsdp(int k, const typdump &dump,const Base &b,const typserv &serv)
3796 {
3797  return !(k%2);
3798 }
3799 template<typename typP,typename typPk,typename typdump,typename Base
3800  , template<typename T> class tmptypmat ,typename monomial_server>
3801 void algo(typP P,typdump &dump,Base &b, monomial_server &serv)
3802 {
3803 #ifdef PEQUAN
3804  int toto=0;
3805  FILE *fic;
3806  char name[8];
3807 #endif
3808  typedef typename typP::value_type polyalp;
3809  typedef list<typename typPk::value_type::monom_t> typMk;
3810  typedef typename typMk::value_type mon;
3811  typedef typename typPk::value_type pol;
3812 
3813  typMk Mk;
3814  typMk Spoldejavu;
3815  typPk Pk,S;
3816 
3817  //typP recall
3818  typPk recall;
3819  //Base b;
3820  workspace<typename typP::value_type::coeff_t> w;
3821  typedef tmptypmat<typename typP::value_type::coeff_t> typMat;
3822  typMat Mat,L;
3823  int onemoretime=1,nouveaumon,allpolused=0,maxdeg,k,*pr,*pc;
3824  mdebug()<<"avant init";
3825 
3826  init_algo1<typMat>(P,nouveaumon,maxdeg,k,Pk,Mk,b,dump,w,serv);
3827  mdebug()<<"apres init";
3828  mdebug()<<"onemoretime "<<onemoretime;
3829  mdebug()<<"nouveaumon "<<nouveaumon;
3830  mdebug()<<"k "<<k;
3831  mdebug()<<"maxdeg "<<maxdeg;
3832 #ifndef GROEBNER
3833 
3834 
3835  // while((critere(b,k)==-1)||(k<=maxdeg))
3836 #ifdef ONEMORETIME
3837  mdebug()<<"one";
3838  while(nouveaumon||(k<=maxdeg)||onemoretime)
3839 #else
3840  while(nouveaumon||(k<=maxdeg))
3841 #endif
3842 #else
3843  while(!GB_stop(k)||(k<=maxdeg))
3844 #endif
3845 
3846  {
3847 #ifdef ONEMORETIME
3848  mdebug()<<"two";
3849  mdebug()<<"nouveaumon";
3850  if(!nouveaumon)
3851  mdebug()<<"si";
3852  onemoretime=0;
3853 #endif
3854  //if (toto>=5) {exit(0);}
3855  //toto++;
3856 
3857  int killspol;
3858  //test si on est dans les conditions d'application
3859  gettimeofday(&initclock,NULL);
3860  mdebug()<<"onemoretime "<<onemoretime;
3861  mdebug()<<"nouveaumon"<<nouveaumon<<"\n";
3862  mdebug()<<"k "<<k;
3863  mdebug()<<"maxdeg "<<maxdeg;
3864  mdebug()<<"testsdp(k,dump,b,serv) "<<testsdp(k,dump,b,serv);
3865  if(testsdp(k,dump,b,serv)||(!onemoretime))
3866  //envoie la routine de calcul de sdp
3867  {
3868  int mindeg_newpol;
3869  typP newpol;
3870  mdebug()<<"un nouveau tour de boucle et k vaut"<<k;
3871  if(!onemoretime)//sdpa!!!!!!!!
3872  // generickernel(newpol,2*k,dump,b,serv);
3873  generickernel_testT(newpol,2*k,dump,b,serv);
3874  // generickernelsdpa(newpol,2*k,dump,b,serv);
3875  else
3876  // generickernel(newpol,k,dump,b,serv);
3877  generickernel_testT(newpol,k,dump,b,serv);
3878  //generickernelsdpa(newpol,k,dump,b,serv);
3879  gettimeofday(&tmpclock,NULL);
3880  tempsdssdp+=tmpclock.tv_sec-initclock.tv_sec+1e-6*(tmpclock.tv_usec-initclock.tv_usec);
3881 
3882  //sdprealrad(newpol,dump,b,serv);
3883  //penser a decroitre k le cas echeant avec un coup de recover
3884  mdebug()<<"newpol "<<newpol.empty();
3885  if(!newpol.empty())
3886  {
3887  mindeg_newpol=Degree(newpol.front());
3888  // mdebug()<<"mindeg "<<mindeg_newpol<<" k "<<k;
3889  for(typename typP::iterator iter=newpol.begin();
3890  iter!=newpol.end();iter++)
3891  mindeg_newpol=(mindeg_newpol<Degree(*iter)
3892  ?mindeg_newpol:Degree(*iter));
3893  for(typename typP::iterator iter=newpol.begin();
3894  iter!=newpol.end();iter++)
3895  if(Degree(*iter)==mindeg_newpol)
3896  P.push_back(*iter);
3897 
3898  if(k>=mindeg_newpol)
3899  {
3900  list<mon *> exceed;
3901  list<int> sizeexceed;
3902  typPk tmpPk;
3903  conv_merge(tmpPk,newpol,exceed,sizeexceed,b,serv);
3904  // mdebug()<<"merde merde k "<<k<<" mindeg_newpol "<<mindeg_newpol;
3905  k=mindeg_newpol;
3906  // forget(b,dump,k,recall);
3907  // recover(Pk,Mk,dump,k);
3908  recover(Pk,Mk,dump,k);
3909 
3910  //Listedejafait(Mk,Spoldejavu);
3911  // mdebug()<<"Pk.size() apres recover "<<Pk.size();
3912  if((k>0)||(k==-1))
3913  forget(b,dump,k,recall);
3914  else
3915  forget0(b,dump,k,recall);
3916  // mdebug()<<"forget fait";
3917  //forget(b,dump,k);
3918  NewCrochetchoix(P,Pk,Mk,tmpPk,k,dump,b,recall,w,serv);
3919  serv.compress(Pk,dump,b,k);
3920  mdebug()<<"compress";
3921  Dump(Pk,dump);
3922  continue;
3923  }
3924  }
3925  }
3926 #ifdef BONRAJOUT
3927  for(typename Base::const_iterator p=b.begin();p!=b.end();p++)
3928  {
3929  for(int i=0;i<p->taille1;i++)
3930  {
3931  mdebug()<<"p.accept[i] "<<p->accept[i];
3932  }
3933  for(int i=0;i<p->taille2;i++)
3934  {
3935  mdebug()<<"p.refuse[i] "<<p->refuse[i];
3936  }
3937  }
3938 #endif
3939 #if 1//def DEB
3940 
3941  //mdebug()<<"mon a la fin de la boulcle de l'algo "<<Mk.size();
3942  //for(typename typMk::iterator iter=Mk.begin();iter!=Mk.end();iter++)
3943  // mdebug()<<*iter;
3944  mdebug()<<"poly a la fin de la boulcle de l'algo "<<Pk.size();
3945  for(typename typPk::iterator iter=Pk.begin();iter!=Pk.end();iter++)
3946  //mdebug()<<"les sizes "<<iter->size;
3947  {
3948  int cpt=0;
3949  mdebug()<<invconv<typename typP::value_type>(*iter,serv);
3950  for(int i=0;i<iter->sizenf;i++)
3951  for(int j=0;j<8;j++)
3952  if (((iter->nfind[i]>>j)&1) && (iter->nf[cpt++]==0))
3953  {
3954  mdebug()<<"ici ca couille a la place"<<8*i+j;
3955  mdebug()<<iter->sizenf<<" "<<iter->size;
3956  exit(1);
3957  }
3958  }
3959 #endif
3960 
3961  // mdebug()<<"Spolordsup";
3962  S=Spolordsup(Pk,b,Spoldejavu);
3963  //S.erase(S.begin(),S.end());
3964  //S=Spolordsup(Pk,b);
3965  //mdebug();
3966  //mdebug()<<endl<<"newpkmkmatrixof";
3967  //mdebug()<<"Pk.size() "<<Pk.size();
3968  //killspol=!newPkMkMatrixof(Mat,Pk,P,Mk,k,b,dump,recall,Spoldejavu);
3969  killspol=!newPkMkMatrixof(Mat,Pk,P,Mk,k,b,dump,recall,w,serv);
3970 #ifdef PEQUAN
3971  sprintf(name,"mat%d",toto);
3972  fic=fopen(name,"w");
3973  printmatfile_pequan(Mat,fic);
3974  fclose(fic);
3975  toto++;
3976 #endif
3977  // if (toto>=5) {mdebug()<<"killredspol "<<killspol; exit(0);}
3978  // killspol=!newPkMkMatrixof(Mat,Pk,P,Mk,k,b,dump);
3979  if(killspol) S.erase(S.begin(),S.end());
3980  if (k+1>=maxdeg) allpolused=1;
3981  if(Mk.size()==0) nouveaumon=0;
3982  //mdebug()<<"j'en suis a k="<<k;
3983  //mdebug()<<"Ici dans algo Mat";
3984  //#ifdef BONRAJOUT
3985  //affdebug(Mat);
3986  //#endif
3987  //Solve affecte a Mat la valeur U
3988  //mdebug()<<"solve";
3989  //for(typename typPk::iterator iter=Pk.begin();iter!=Pk.end();iter++)
3990  // mdebug()<<"avant solve les ind "<<iter->ind;
3991  //for(typename typPk::iterator iter=Pk.begin();iter!=Pk.end();iter++)
3992  // mdebug()<<"avant solve les pol "<<invconv<Poly>(*iter);
3993  //if(k==3) exit(-10);
3994  //for(typename typMk::iterator iter=Mk.begin();iter!=Mk.end();iter++)
3995  // mdebug()<<"monomes dans Mk "<<*iter;
3996 
3997  if (Mat.ncol!=0)
3998  {
3999  //mdebug()<<"debut du solve";
4000  Solve(Mat,L,Pk,&pr,&pc,w);
4001  // mdebug()<<"fin du solve";
4002  //for(int kk=0;kk<Mk.size();kk++)
4003  // mdebug()<<"permc "<<kk<<" "<<pc[kk];
4004  inverse(pc,Mk.size());
4005  Dopermu(Mk,pc);
4006  {
4007  //ici on remet les ind a la bonne valuer
4008  //a inclure dans une petite fonction inline apres
4009  int j=0;
4010  typename typMk::iterator itermon=Mk.begin();
4011  typename typPk::iterator iter;
4012  for(iter=Pk.begin();
4013  itermon!=Mk.end();iter++,itermon++,j++)
4014  {
4015  if(!Iszero(Mat(j,j)))
4016  {
4017  // mdebug()<<"ligne "<<j<<" je mets le monome "<<*itermon;
4018  iter->ind=*itermon*(iter->ind.GetCoeff());
4019  }
4020  else
4021  {
4022  iter->ind=mon(0);
4023  }
4024  }
4025  //for(;itermon!=Mk.end();iter++,itermon++)
4026  // iter->ind=mon(0);
4027  }
4028 
4029  free(pr);
4030  free(pc);
4031  }
4032 
4033  //appliquer les permutation a MK!!!!!!!!
4034  // for(typename typPk::iterator iter=Pk.begin();iter!=Pk.end();iter++)
4035  // mdebug()<<"apres solve les ind "<<iter->ind;
4036  //for(typename typPk::iterator iter=Pk.begin();iter!=Pk.end();iter++)
4037  // mdebug()<<"apres solve les pol "<<invconv<Poly>(*iter);
4038  //
4039  //
4040  //APPLIQUER LES PERMUTYATIONS A MK!!
4041  //
4042  //
4043  //
4044  //
4045  if ((unsigned)my_rank(Mat)==Mk.size())
4046  //matrice de rang plein
4047  {
4048  typMk tmpSpoldejavu;
4049  typPk redspol,secd=secondmembre(Mat,Pk);;
4050  //secondmembre2(Mat,Pk,secd);
4051  //typename typPk::iterator iter=Pk.begin();
4052  mdebug()<<"matrice de rang plein ";
4053  //for(unsigned int i=0;i<Mk.size();i++,iter++);
4054  //Pk.erase(iter,Pk.end());
4055  //mdebug()<<"Dump ";
4056  Dump(Pk,dump);
4057  //mdebug()<<"Reduce Spol ";
4058 
4059  ReduceSpol(S,dump,redspol,b,w,serv);//tmpSpoldejavu);
4060  //mdebug()<<"secd.size() "<<secd.size();
4061  my_merge(redspol,secd);
4062  if(redspol.size()==0)
4063  {
4064  //a voir.... Dump(Pk,dump);
4065  k++;
4066  my_merge(Spoldejavu,tmpSpoldejavu);
4067  L.destroystore();
4068  continue;
4069  //la matrice a rang plein
4070  //ET les Spol\secondmembre de ce niveau se ,reduisent a 0
4071  }
4072  else
4073  {
4074  //on a des spol\secondmembre qui se reduisent pas a 0
4075  //cad un defaut de commuatation
4076  int tmpdeg=mindeg(redspol,serv);
4077  max_degppcm=k+2;
4078  nouveaumon=1;
4079  mdebug()<<"redspol.size() "<<redspol.size();
4080  //mdebug()<<invconv<typename typP::value_type>(redspol.front());
4081  if (tmpdeg==(k+1))
4082  my_merge(Spoldejavu,tmpSpoldejavu);
4083  k=tmpdeg;
4084  //mdebug()<<"le mindeg ici vaut "<<k;
4085  //mdebug()<<"Pk.size() avant recover "<<Pk.size();
4086  recover(Pk,Mk,dump,k);
4087 
4088  //Listedejafait(Mk,Spoldejavu);
4089  //mdebug()<<"Pk.size() apres recover "<<Pk.size();
4090  forget(b,dump,k,recall);
4091  //forget(b,dump,k);
4092 #ifdef DEB
4093  mdebug()<<"Pk.size() avant Crochetchoix "<<Pk.size();
4094 #endif
4095  //NewCrochetchoix(P,Pk,Mk,redspol,k,dump,b,recall,Spoldejavu);
4096  NewCrochetchoix(P,Pk,Mk,redspol,k,dump,b,recall,w,serv);
4097  //NewCrochetchoix(P,Pk,Mk,redspol,k,dump,b);
4098 #ifdef DEB
4099  mdebug()<<"Pk.size apres Crochetchoix "<<Pk.size()<<" et k "<<k;
4100 #endif
4101  NewDestroyredspol(redspol,Pk);
4102  serv.compress(Pk,dump,b,k);
4103  Dump(Pk,dump);
4104  L.destroystore();
4105  //misa jour de B dans crochet choix!!
4106  continue;
4107  }
4108  }
4109  else
4110  //la matrice n'est pas de rang plein
4111  {
4112  mdebug()<<"matrice de rang pas plein ";
4113  typPk secd;
4114 #ifdef BONRAJOUT
4115  AddB_dim(Mat,L,P,Pk,Mk,dump,b);
4116 #else
4117  AddB(Mat,Pk,Mk,dump,b,serv);
4118 #endif
4119  //
4120  //secondmembre supprime des trucs dans Pk!!
4121  //
4122  secd=secondmembre(Mat,Pk);
4123  rajoute_voisins(Pk,Mk,dump,b,w,serv);
4124  //ne selectionne que les non nuls
4125  if(secd.size()==0)
4126  {
4127  typPk redspol;
4128  //for(int i=Pk.size();i>rank(Mat);i--)
4129  //{
4130  // if(Pk.back().size!=0)
4131  // MAC_REV_FREE<typename typMat::coeff_t>(Pk.rbegin()->nf
4132  // ,Pk.rbegin()->size*sizeof(typename typMat::coeff_t));
4133  // Pk.pop_back();
4134  //}//on vire les 0 avant de dumper
4135  // AddB(Mat,Pk,Mk,dump,b);
4136  Dump(Pk,dump);
4137  ReduceSpol(S,dump,redspol,b,w,serv);//,Spoldejavu);
4138 #ifdef NO_TWICE
4139  for(typename typPk::iterator iter=redspol.begin();iter!=redspol.end();iter++)
4140  P.push_back(invconv<polyalp>(*iter));
4141 #endif
4142  if(redspol.size()==0)
4143  {
4144  // forget(b,dump,k);
4145  //AddB(Mat,Pk,Mk,b);
4146  //Dump(Pk,dump);
4147  //ici on modifie Pk dans AddB
4148  //<=> forget(b,dump,k) AddB(..) Dump(Pk,dump);
4149  k++;
4150  L.destroystore();
4151  continue;
4152  //la matrice a pas rang plein second membre =0
4153  //ET les Spol de ce niveau se ,reduisent a 0
4154  }
4155  else
4156  {
4157  //on a des spol qui se reduisent pas a 0
4158  //cad un defaut de commutation
4159  k=mindeg(redspol,serv);
4160  recover(Pk,Mk,dump,k);
4161  //Listedejafait(Mk,Spoldejavu);
4162  forget(b,dump,k,recall);
4163  //forget(b,dump,k);
4164  mdebug()<<"il ya un passage par mat rang pas"
4165  <<"plein et pol non zero";
4166  //NewCrochetchoix(P,Pk,Mk,redspol,k,dump,b,recall,Spoldejavu);
4167  NewCrochetchoix(P,Pk,Mk,redspol,k,dump,b,recall,w,serv);
4168  //NewCrochetchoix(P,Pk,Mk,redspol,k,dump,b);
4169  NewDestroyredspol(redspol,Pk);
4170  serv.compress(Pk,dump,b,k);
4171  Dump(Pk,dump);
4172  L.destroystore();
4173  continue;
4174  }
4175  }
4176  else
4177  {
4178  //on a des pol qui se reduisent pas a 0
4179  //cad un defaut de commuatation dans les coins
4180  typPk redspol;
4181  //for(int i=Pk.size();i>rank(Mat);i--)
4182  // Pk.pop_back();
4183  //mdebug()<<"PK. size "<<Pk.size()<<" Mat.nbcol "<<Mat.ncol;
4184  //mdebug()<<"rank(Mat) "<<rank(Mat);
4185  //int j=0;
4186  //mdebug()<<"avant recover et meme avant le dump ";
4187  //for(typename typPk::iterator iter = secd.begin();iter!=secd.end()
4188  // ;iter++,j++)
4189  // mdebug()<<"pol de secd num "<<j<<endl<<invconv<Poly>(*iter);
4190  // AddB(Mat,Pk,Mk,b);
4191  //mdebug()<<"apres ADDB le dump ";
4192  //for(typename typPk::iterator iter = Pk.begin();iter!=Pk.end()
4193  // ;iter++)
4194  // mdebug()<<invconv<Poly>(*iter);
4195  Dump(Pk,dump);
4196  //mdebug()<<"secd.size()" <<secd.size();
4197  ReduceSpol(S,dump,redspol,b,w,serv);//,Spoldejavu);
4198  my_merge(secd,redspol);
4199 #ifdef NO_TWICE
4200  for(typename typPk::iterator iter=redspol.begin();iter!=redspol.end();iter++)
4201  P.push_back(invconv<polyalp>(*iter));
4202 #endif
4203  k=mindeg(secd,serv);
4204  recover(Pk,Mk,dump,k);
4205  //mdebug()<<"k "<<k;
4206  //Listedejafait(Mk,Spoldejavu);
4207  //mdebug()<<"apres recover";
4208  // for(typename typPk::iterator iter = Pk.begin();iter!=Pk.end()
4209  // ;iter++)
4210  // mdebug()<<invconv<Poly>(*iter);
4211  //mdebug()<<"les Spol non reduits";
4212  //for(typename typPk::iterator iter = secd.begin();iter!=secd.end()
4213  // ;iter++)
4214  // mdebug()<<invconv<Poly>(*iter);
4215  forget(b,dump,k,recall);
4216 
4217  //forget(b,dump,k);
4218  //mdebug()<<"il ya un passage par mat pas rang plein et Spol non zero";
4219  //for(typename typPk::iterator iter = secd.begin();iter!=secd.end()
4220  // ;iter++)
4221  // mdebug()<<invconv<Poly>(*iter);
4222  //NewCrochetchoix(P,Pk,Mk,redspol,k,dump,b,recall,Spoldejavu);
4223  NewCrochetchoix(P,Pk,Mk,secd,k,dump,b,recall,w,serv);
4224  //NewCrochetchoix(P,Pk,Mk,secd,k,dump,b);
4225  NewDestroyredspol(secd,Pk);
4226 #ifdef BONRAJOUT
4227  mdebug()<<"Apres NewCrochetchoix";
4228 
4229  for(typename Base::const_iterator p=b.begin();p!=b.end();p++)
4230  {
4231  for(int i=0;i<p->taille1;i++)
4232  {
4233  mdebug()<<"p.accept[i] "<<p->accept[i];
4234  }
4235  for(int i=0;i<p->taille2;i++)
4236  {
4237  mdebug()<<"p.refuse[i] "<<p->refuse[i];
4238  }
4239  }
4240 #endif
4241  //for(typename typPk::iterator iter = Pk.begin();iter!=Pk.end()
4242  // ;iter++)
4243  // mdebug()<<invconv<Poly>(*iter);
4244  serv.compress(Pk,dump,b,k);
4245  Dump(Pk,dump);
4246  L.destroystore();
4247  continue;
4248  }
4249 
4250  }
4251  }
4252 
4253  Mat.destroystore();
4254  //while nouveau mon etc...
4255  // Pk=recomposePk(dump);
4256  //return Pk;
4257  // return;
4258 #ifdef DEBUG
4259  for(typename Base::const_iterator p=b.begin();p!=b.end();p++)
4260  {
4261  for(int i=0;i<p->taille1;i++)
4262  {
4263  mdebug()<<"p.accept[i] "<<p->accept[i];
4264  }
4265  for(int i=0;i<p->taille2;i++)
4266  {
4267  mdebug()<<"p.refuse[i] "<<p->refuse[i];
4268  }
4269  }
4270 #endif
4271  destroy_space(&w);
4272  //mdebug()<<"sortie d'algo flash et dimension ";//<<critere(b,k);
4273  return ;//dump;
4274  //comme ca on conserve le tri par degre des formes normales et
4275  //donc la construction des operateurs de multiplication se
4276  //fait plus vite.. :-)
4277 
4278 }
4279 
4280 template<typename typP,typename typPk,typename typdump,typename Base
4281  , template<typename T> class tmptypmat ,typename monomial_server,typename coeff>
4282 void algo(typP P,typdump &dump,Base &b, monomial_server &serv,typP F, coeff *&fmin, Base &b2, typP restr, coeff *&solfin,int solver, int &method, string param, double &tresh, int borderbasis)
4283 {
4284 #ifdef PEQUAN
4285  int toto=0;
4286  FILE *fic;
4287  char name[8];
4288 #endif
4289  typedef typename typP::value_type polyalp;
4290  typedef list<typename typPk::value_type::monom_t> typMk;
4291  typedef typename typMk::value_type mon;
4292  typedef typename typPk::value_type pol;
4293 
4294  typMk Mk;
4295  typMk Spoldejavu;
4296  typPk Pk,S;
4297  typPk recall;
4298  typP Pini,restr2;
4299  workspace<typename typP::value_type::coeff_t> w;
4300  typedef tmptypmat<typename typP::value_type::coeff_t> typMat;
4301  typMat Mat,L;
4302  int onemoretime=1,nouveaumon,allpolused=0,maxdeg,k,maxdegb,kb,*pr,*pc,optimo=0,k1,rest,m=0,kk,drest,dmax=0,sup,flat,noeq=0,numeq=0;
4303  int nogap=0,deg=0,nok=1, noreduced=0, kbon,kini,basesdifferent,fois=0,premierfois=0,knew=0,bonbases=0,bonborder=0,elemnoyau=0,nvar,solsize;
4304  fmin =new coeff[1];
4305 
4306  typP newpol,newpol1;
4307  COEFF *sol2,*sol3,*solfinal;
4308 
4309 
4310  for(typename typP::iterator iter2=P.begin();
4311  iter2!=P.end();iter2++)
4312  {
4313  numeq++;
4314  Pini.push_back(*iter2);
4315  }
4316  init_algo<typMat>(P,nouveaumon,maxdeg,k,Pk,Mk,b,dump,w,serv,b2);
4317 
4318  //1006
4319  if (k!=100)
4320  {
4321  if (!restr.empty())
4322  {
4323  mdebug()<<"K"<<k;
4324  mdebug()<<"maxdeg"<<maxdeg;
4325  if (maxdeg==Degree(F.front())) maxdeg*=2;
4326  else k=max(maxdeg,Degree(F.front()));//+1;
4327  if (maxdeg%2==1){
4328  k++;
4329  maxdeg++;
4330  }
4331  }
4332  else maxdeg=max(maxdeg,Degree(F.front()));
4333  sup=maxdeg;
4334  // if (!restr.empty())
4335  // //drest=minmaxdegP(restr,maxdeg);
4336  // if (2*maxdeg>Degree(F.front()))
4337  // {
4338 
4339  // maxdeg*=2;
4340  // sup=maxdeg;
4341  // mdebug()<<"maxdeg"<<maxdeg;
4342  // }
4343 
4344  }
4345  else
4346  {
4347 
4348  if (!restr.empty())
4349  drest=minmaxdegP(restr,maxdeg);
4350  else maxdeg=1;//cas de la base complete
4351 
4352  noeq=1;
4353  if (2*maxdeg>Degree(F.front()))
4354  {
4355 // k=maxdeg+2;
4356 // maxdeg=2*k;
4357 // sup=2*k;
4358  k=maxdeg;
4359  maxdeg=k;
4360  sup=k;
4361  mdebug()<<"se mete por aqui maxdeg"<<maxdeg;
4362  }
4363  else{
4364  mdebug()<<"K"<<k;
4365  if (method==1 && (2*maxdeg!=Degree(F.front()))) k=Degree(F.front());
4366  else k=Degree(F.front())+1;
4367  if (k%2==1) k++;
4368  maxdeg=k;
4369  sup=k;
4370  // k=maxdeg+2;
4371  // maxdeg=2*k;
4372  // sup=2*k;
4373  }
4374  }
4375 
4376 
4377  nvar=b.nbvars;
4378 
4379  mdebug()<<"no hay constraint"<<restr.empty()<<"k"<<k;
4380 
4381 #ifndef GROEBNER
4382 
4383 #ifdef ONEMORETIME
4384  mdebug()<<"one";
4385  while(nouveaumon||(k<=maxdeg)||onemoretime)
4386 #else
4387  while(nouveaumon||(k<=maxdeg)|| (!optimo) || nok )
4388 
4389 #endif
4390 #else
4391  while(!GB_stop(k)||(k<=maxdeg))
4392 #endif
4393 
4394  {
4395 #ifdef ONEMORETIME
4396  mdebug()<<"two";
4397  mdebug()<<"nouveaumon";
4398  if(!nouveaumon)
4399  mdebug()<<"si";
4400  onemoretime=0;
4401 #endif
4402 
4403 
4404  int killspol;
4405  //test si on est dans les conditions d'application
4406  gettimeofday(&initclock,NULL);
4407  mdebug()<<"onemoretime "<<onemoretime;
4408  mdebug()<<"nouveaumon"<<nouveaumon;
4409  mdebug()<<"k "<<k;
4410  mdebug()<<"maxdeg "<<maxdeg;
4411  mdebug()<<"testsdp(k,dump,b,serv) "<<testsdp(k,dump,b,serv);
4412  for(typename Base::iterator p=b.begin();p!=b.end();p++)
4413  {
4414  mdebug()<<"p->taille1"<<p->taille1;
4415  mdebug()<<"p->taille2"<<p->taille2;
4416  for(int i=0;i<p->taille1;i++)
4417  {
4418  mdebug()<<"p.accept[i] "<<p->accept[i];
4419  }
4420  for(int i=0;i<p->taille2;i++)
4421  {
4422  mdebug()<<"p.refuse[i] "<<p->refuse[i];
4423  }
4424  }
4425  mdebug()<<"base apres";
4426 
4427 
4428  mdebug()<<"maxdeg"<<maxdeg;
4429  mdebug()<<"optimo"<<optimo;
4430  mdebug()<<"nok"<<nok;
4431  mdebug()<<"deg"<<deg;
4432  mdebug()<<"sup"<<sup;
4433  basesdifferent=0;
4434 
4435  if((testsdp(k,dump,b,serv) && (!optimo)) && (k!=1) || (k==deg) )
4436 
4437  {
4438  int mindeg_newpol;
4439 
4440 
4441  mdebug()<<"k cuando emepzamos"<<k<<"sup"<<sup;
4442 
4443 
4444  //la variable premierfois veux dire si on a fait ou pas SDPA
4445  //on commence avec le bases egales dans cet if on copy notre base b
4446  //en la base b2
4447  if (premierfois==0 && (k<=sup) && (!basesdifferent))
4448  {
4449 
4450  // if (k>sup)
4451  // {
4452  // bonbases=1;
4453  // if (k-1!=sup) k-=2;
4454  // }
4455 
4456  mdebug()<<"un nouveau tour de boucle et k vaut avec compare bases"<<k;
4457  //on compare les bases
4458  for(typename Base::iterator p=b.begin();p!=b.end();p++)
4459  for(typename Base::iterator p2=b2.begin();p2!=b2.end();p2++)
4460  {
4461 
4462  if (p2->taille2!=p->taille2)
4463  basesdifferent=1;
4464  else
4465  for(int i=0;i<p2->taille2;i++)
4466  {
4467  if (p2->refuse[i]!=p->refuse[i])
4468  basesdifferent=1;
4469  }
4470  }
4471  mdebug()<<"basesdifferent"<<basesdifferent;
4472  if (basesdifferent)
4473  {
4474  fois++;
4475  for(typename Base::iterator p=b.begin();p!=b.end();p++)
4476  for(typename Base::iterator p2=b2.begin();p2!=b2.end();p2++)
4477  {
4478 
4479 
4480  p2->taille2 =p->taille2 ;
4481  p2->refuse=(mon*)MAC_REV_MALLOC<mon>(p2->taille2*sizeof(mon));
4482 
4483  for(int i=0;i<p2->taille1;i++)
4484  {
4485  p2->accept[i]=p->accept[i];
4486  }
4487 
4488  for(int i=0;i<p2->taille2;i++)
4489  {
4490 
4491  p2->refuse[i]=p->refuse[i];
4492 
4493  }
4494 
4495  }
4496  }
4497 
4498  }
4499 
4500 
4501 
4502  if ((borderbasis==0 && !basesdifferent && k>=maxdeg ) || (borderbasis==1 && k>=maxdeg) )
4503  // if ((borderbasis==0 && k>=maxdeg ) || (borderbasis==1 && k>=maxdeg) )
4504 
4505  {
4506 
4507 
4508  if (nogap)
4509  {
4510 
4511  if (k==deg)
4512  {
4513  newpol1.erase(newpol1.begin(),newpol1.end());
4514 
4515 
4516  for(typename typP::iterator iter=newpol.begin();iter!=newpol.end();
4517  iter++)
4518  {
4519  mdebug()<<*iter<<" ";
4520  mdebug();
4521  }
4522  }
4523 
4524 
4525  newpol1.erase(newpol1.begin(),newpol1.end());
4526 
4527 
4528 
4529 
4530  //generickernel_testF(newpol1,k,dump,b,serv,F,*&sol3,restr,*&solfinal,solver,method,flat,param,tresh);
4531  mdebug()<<"tresh"<<tresh;
4532  generickernelNEW(newpol1,k,dump,b,serv,F,*&sol3,restr,*&solfinal,solver,method,flat,param,tresh);
4533 
4534  if (solfinal[0]!=0)
4535  {
4536  mdebug()<<"solfinal[0]"<<solfinal[0];
4537  mdebug()<<"nvar"<<nvar;
4538  //solsize=solfinal[0]*nvar+1;
4539  solsize=solfinal[0]*nvar+2;
4540 
4541 
4542  mdebug()<<"solsize"<<solsize;
4543  solfin =new coeff[solsize];
4544  for (int l=0; l<solsize;l++)
4545  {
4546  solfin[l]=solfinal[l];
4547  mdebug()<<"solution final"<<solfin[l]<<"l"<<l;
4548  }
4549 
4550  /*if (solver==4) fmin[0]=sol2[0];
4551  else */ //fmin[0]=sol2[2];
4552  fmin[0]=solfin[1];
4553 
4554  return;
4555  }
4556 
4557 
4558 
4559  // mdebug()<<"solucionprimal k"<< sol2[2] << endl;
4560  // mdebug()<<"solucionprimal k+1"<< sol3[2] << endl;
4561  // mdebug()<<"diferenciasprimales"<<abs(sol2[2]-sol3[2])*1e+3;
4562  mdebug()<<"k1"<<k1;
4563  mdebug()<<"k"<<k;
4564  // mdebug()<<"gap"<<abs(sol3[2]-sol3[3])*1e+4;
4565  mdebug()<<"deg"<<deg;
4566  mdebug()<<"!newpol.empty()"<<!newpol.empty();
4567  mdebug()<<"newpol1.empty()"<<newpol1.empty();
4568 
4569 
4570  if (((abs(sol2[2]-sol3[2])*1e+4<1) && (!newpol.empty()) && (k1!=k))
4571  || ((k==deg)))
4572  {
4573  if (k==deg && newpol1.empty())
4574  {
4575  nok=0;
4576  continue;
4577  }
4578  else
4579  {
4580  elemnoyau++;
4581  if (k!=deg)
4582  {
4583  optimo=1;
4584  fmin[0]=sol3[2];
4585  deg=k+2;
4586 
4587 
4588  int numkerf=0;
4589  mdebug()<<"noyau bon";
4590  for(typename typP::iterator iter=newpol.begin();iter!=newpol.end();
4591  iter++,numkerf++)
4592  {
4593  mdebug()<<*iter<<" ";
4594  mdebug();
4595  }
4596 
4597  }
4598 
4599  if (k==deg && (!newpol.empty()) )
4600  {
4601  newpol.erase(newpol.begin(),newpol.end());
4602  newpol=newpol1;
4603 
4604  }
4605 
4606  mindeg_newpol=Degree(newpol.front());
4607  for(typename typP::iterator iter=newpol.begin();
4608  iter!=newpol.end();iter++)
4609  mindeg_newpol=(mindeg_newpol<Degree(*iter)
4610  ?mindeg_newpol:Degree(*iter));
4611 
4612  P.erase(P.begin(),P.end());
4613  for(typename typP::iterator iter=newpol.begin();
4614  iter!=newpol.end();iter++)
4615  if(Degree(*iter)==mindeg_newpol)
4616  P.push_back(*iter);
4617  for(typename typP::iterator iter2=Pini.begin();
4618  iter2!=Pini.end();iter2++)
4619  {
4620  //ajoutons le relations initials ( ideal gradient, langrange, kkt)
4621  //il y a des cases ou les ajouter produce error numerique
4622 
4623  P.push_back(*iter2);
4624  }
4625  for(typename typP::iterator iter3=P.begin();
4626  iter3!=P.end();iter3++)
4627  mindeg_newpol=(mindeg_newpol<Degree(*iter3)
4628  ?mindeg_newpol:Degree(*iter3));
4629  if(k-2>=mindeg_newpol)
4630  {
4631  list<mon *> exceed;
4632  list<int> sizeexceed;
4633  typPk tmpPk;
4634  conv_merge(tmpPk,newpol,exceed,sizeexceed,b,serv);
4635  k=mindeg_newpol;
4636 
4637  recover(Pk,Mk,dump,k);
4638  if (k!=deg)
4639  {
4640  if((k>0)||(k==-1))
4641  forget(b,dump,k,recall);
4642  else
4643  forget0(b,dump,k,recall);
4644  }
4645 
4646  NewCrochetchoix(P,Pk,Mk,tmpPk,k,dump,b,recall,w,serv);
4647 
4648 
4649  mdebug();
4650  serv.compress(Pk,dump,b,k);
4651  mdebug()<<"compress";
4652  Dump(Pk,dump);
4653  continue;
4654  }
4655  }
4656  }
4657  else
4658  {
4659 
4660 
4661  sol2=sol3;
4662  newpol=newpol1;
4663 
4664 
4665  }
4666 
4667  }
4668  else
4669  {
4670 
4671  k1=k;
4672  newpol.erase(newpol.begin(),newpol.end());
4673 
4674 
4675  mdebug()<<"premier k"<<k;
4676 
4677  //generickernel_testF(newpol,k,dump,b,serv,F,*&sol2,restr,*&solfinal,solver,method,flat,param,tresh);
4678  generickernelNEW(newpol,k,dump,b,serv,F,*&sol2,restr,*&solfinal,solver,method,flat,param,tresh);
4679  mdebug()<<"Number of solutions"<<solfinal[0];
4680  if (solfinal[0]!=0 && solfinal[0]!=-1 && method==0)
4681  {
4682 
4683  mdebug()<<"nvar"<<nvar;
4684  //solsize=solfinal[0]*nvar+1;
4685  solsize=solfinal[0]*nvar+2;
4686 
4687  mdebug()<<"solsize"<<solsize;
4688  solfin =new coeff[solsize];
4689  for (int k=0; k<solsize;k++)
4690  {
4691  solfin[k]=solfinal[k];
4692  mdebug()<<"solution final"<<solfin[k];
4693  }
4694  /*if (solver==4) fmin[0]=sol2[0];
4695  else /*///fmin[0]=sol2[2];
4696  fmin[0]=solfinal[1];
4697 
4698 
4699  return;
4700  }
4701 
4702 
4703 
4704  if (solfinal[0]==-1)
4705  {
4706  method=1;
4707  return;
4708  }
4709 
4710 
4711  premierfois=1;
4712  int numker=0;
4713 
4714  mdebug()<<"apres decomposer";
4715  for(typename typP::iterator iter=newpol.begin();iter!=newpol.end();
4716  iter++,numker++)
4717  mdebug()<<*iter<<" ";
4718  mdebug();
4719 
4720  if (method==1 && (!newpol.empty()) && flat==1)
4721  {
4722  mdebug()<<"pasamos";
4723  mindeg_newpol=Degree(newpol.front());
4724  for(typename typP::iterator iter=newpol.begin();
4725  iter!=newpol.end();iter++)
4726  mindeg_newpol=(mindeg_newpol<Degree(*iter)
4727  ?mindeg_newpol:Degree(*iter));
4728 
4729  P.erase(P.begin(),P.end());
4730  for(typename typP::iterator iter=newpol.begin();
4731  iter!=newpol.end();iter++)
4732  if(Degree(*iter)==mindeg_newpol)
4733  P.push_back(*iter);
4734 
4735  for(typename typP::iterator iter3=P.begin();
4736  iter3!=P.end();iter3++)
4737  mindeg_newpol=(mindeg_newpol<Degree(*iter3)
4738  ?mindeg_newpol:Degree(*iter3));
4739  if(k-2>=mindeg_newpol)
4740  {
4741  mdebug()<<"pasamos1";
4742  list<mon *> exceed;
4743  list<int> sizeexceed;
4744  typPk tmpPk;
4745  conv_merge(tmpPk,newpol,exceed,sizeexceed,b,serv);
4746  k=mindeg_newpol;
4747 
4748  recover(Pk,Mk,dump,k);
4749  if (k!=deg)
4750  {
4751  if((k>0)||(k==-1))
4752  forget(b,dump,k,recall);
4753  else
4754  forget0(b,dump,k,recall);
4755  }
4756 
4757  NewCrochetchoix(P,Pk,Mk,tmpPk,k,dump,b,recall,w,serv);
4758 
4759 
4760  serv.compress(Pk,dump,b,k);
4761 
4762  Dump(Pk,dump);
4763  //continue;
4764  }
4765  fmin[0]=solfinal[1];
4766  break;
4767  }
4768 
4769 
4770  }
4771 
4772  }
4773  }
4774 
4775 
4776 
4777 
4778 #ifdef BONRAJOUT
4779  for(typename Base::const_iterator p=b.begin();p!=b.end();p++)
4780  {
4781  for(int i=0;i<p->taille1;i++)
4782  {
4783  mdebug()<<"p.accept[i] "<<p->accept[i];
4784  }
4785  for(int i=0;i<p->taille2;i++)
4786  {
4787  mdebug()<<"p.refuse[i] "<<p->refuse[i];
4788  }
4789  }
4790 #endif
4791  //#ifdef DEB
4792 
4793  //mdebug()<<"mon a la fin de la boulcle de l'algo "<<Mk.size();
4794  //for(typename typMk::iterator iter=Mk.begin();iter!=Mk.end();iter++)
4795  // mdebug()<<*iter;
4796  mdebug()<<"poly a la fin de la boulcle de l'algo "<<Pk.size();
4797  for(typename typPk::iterator iter=Pk.begin();iter!=Pk.end();iter++)
4798  //mdebug()<<"les sizes "<<iter->size;
4799  {
4800  int cpt=0;
4801  // mdebug()<<invconv<typename typP::value_type>(*iter,serv);
4802  for(int i=0;i<iter->sizenf;i++)
4803  for(int j=0;j<8;j++)
4804  if (((iter->nfind[i]>>j)&1) && (iter->nf[cpt++]==0))
4805  {
4806  mdebug()<<"ici ca couille a la place"<<8*i+j;
4807  mdebug()<<iter->sizenf<<" "<<iter->size;
4808  exit(1);
4809  }
4810  }
4811  //#endif
4812 
4813  S=Spolordsup(Pk,b,Spoldejavu);
4814 
4815 
4816  mdebug()<<"newpkmkmatrixof";
4817  // mdebug()<<"k"<<k;
4818 
4819  killspol=!newPkMkMatrixof(Mat,Pk,P,Mk,k,b,dump,recall,w,serv);
4820 #ifdef PEQUAN
4821  sprintf(name,"mat%d",toto);
4822  fic=fopen(name,"w");
4823  printmatfile_pequan(Mat,fic);
4824  fclose(fic);
4825  toto++;
4826 #endif
4827  mdebug()<<"apres kill";
4828  if(killspol) S.erase(S.begin(),S.end());
4829  if (k+1>=maxdeg) allpolused=1;
4830  // mdebug()<<"Mk.size()"<<Mk.size();
4831  if(Mk.size()==0) nouveaumon=0;
4832  //Solve affecte a Mat la valeur U
4833  //mdebug()<<"solve";
4834  // for(typename typPk::iterator iter=Pk.begin();iter!=Pk.end();iter++)
4835  // mdebug()<<"avant solve les ind "<<iter->ind;
4836  // for(typename typMk::iterator iter=Mk.begin();iter!=Mk.end();iter++)
4837  // mdebug()<<"monomes dans Mk "<<*iter;
4838  // mdebug()<<"Mat.ncol"<<Mat.ncol;
4839  if (Mat.ncol!=0)
4840  {
4841  mdebug()<<"debut du solve";
4842  Solve(Mat,L,Pk,&pr,&pc,w);
4843  mdebug()<<"fin du solve";
4844  //for(int kk=0;kk<Mk.size();kk++)
4845  // mdebug()<<"permc "<<kk<<" "<<pc[kk];
4846  inverse(pc,Mk.size());
4847  Dopermu(Mk,pc);
4848  {
4849  //ici on remet les ind a la bonne valuer
4850  //a inclure dans une petite fonction inline apres
4851  int j=0;
4852  typename typMk::iterator itermon=Mk.begin();
4853  typename typPk::iterator iter;
4854  for(iter=Pk.begin();
4855  itermon!=Mk.end();iter++,itermon++,j++)
4856  {
4857  if(!Iszero(Mat(j,j)))
4858  {
4859  // mdebug()<<"ligne "<<j<<" je mets le monome "<<*itermon;
4860  iter->ind=*itermon*(iter->ind.GetCoeff());
4861  }
4862  else
4863  {
4864  iter->ind=mon(0);
4865  }
4866  }
4867 
4868  }
4869 
4870  free(pr);
4871  free(pc);
4872  }
4873 
4874  //appliquer les permutation a MK!!!!!!!!
4875 
4876  //
4877  //APPLIQUER LES PERMUTYATIONS A MK!!
4878  //
4879  //
4880  //
4881  //
4882  if ((unsigned)my_rank(Mat)==Mk.size())
4883  //matrice de rang plein
4884  {
4885  typMk tmpSpoldejavu;
4886  typPk redspol,secd=secondmembre(Mat,Pk);;
4887 
4888  Dump(Pk,dump);
4889 
4890 
4891  ReduceSpol(S,dump,redspol,b,w,serv);//tmpSpoldejavu);
4892 
4893 
4894  my_merge(redspol,secd);
4895  if(redspol.size()==0)
4896  {
4897 
4898  k++;
4899  my_merge(Spoldejavu,tmpSpoldejavu);
4900  L.destroystore();
4901  continue;
4902  //la matrice a rang plein
4903  //ET les Spol\secondmembre de ce niveau se ,reduisent a 0
4904  }
4905  else
4906  {
4907  //on a des spol\secondmembre qui se reduisent pas a 0
4908  //cad un defaut de commuatation
4909  int tmpdeg=mindeg(redspol,serv);
4910  max_degppcm=k+2;
4911  nouveaumon=1;
4912 
4913  if (tmpdeg==(k+1))
4914  my_merge(Spoldejavu,tmpSpoldejavu);
4915  k=tmpdeg;
4916 
4917  recover(Pk,Mk,dump,k);
4918 
4919 
4920 
4921  forget(b,dump,k,recall);
4922 
4923 #ifdef DEB
4924  mdebug()<<"Pk.size() avant Crochetchoix "<<Pk.size();
4925 #endif
4926 
4927  NewCrochetchoix(P,Pk,Mk,redspol,k,dump,b,recall,w,serv);
4928 
4929 #ifdef DEB
4930  mdebug()<<"Pk.size apres Crochetchoix "<<Pk.size()<<" et k "<<k;
4931 #endif
4932  NewDestroyredspol(redspol,Pk);
4933  serv.compress(Pk,dump,b,k);
4934  Dump(Pk,dump);
4935 
4936  L.destroystore();
4937 
4938  //misa jour de B dans crochet choix!!
4939  continue;
4940  }
4941  }
4942  else
4943  //la matrice n'est pas de rang plein
4944  {
4945  mdebug()<<"matrice de rang pas plein ";
4946  for(typename Base::const_iterator p=b.begin();p!=b.end();p++)
4947  {
4948  for(int i=0;i<p->taille1;i++)
4949  {
4950  mdebug()<<"p.accept[i] "<<p->accept[i];
4951  }
4952  for(int i=0;i<p->taille2;i++)
4953  {
4954  mdebug()<<"p.refuse[i] "<<p->refuse[i];
4955  }
4956  }
4957  typPk secd;
4958 #ifdef BONRAJOUT
4959  AddB_dim(Mat,L,P,Pk,Mk,dump,b);
4960 #else
4961  AddB(Mat,Pk,Mk,dump,b,serv);
4962 
4963 #endif
4964  mdebug()<<"apres addb";
4965  for(typename Base::const_iterator p=b.begin();p!=b.end();p++)
4966  {
4967  for(int i=0;i<p->taille1;i++)
4968  {
4969  mdebug()<<"p.accept[i] "<<p->accept[i];
4970  }
4971  for(int i=0;i<p->taille2;i++)
4972  {
4973  mdebug()<<"p.refuse[i] "<<p->refuse[i];
4974  }
4975  }
4976  //
4977  //secondmembre supprime des trucs dans Pk!!
4978  //
4979  secd=secondmembre(Mat,Pk);
4980  rajoute_voisins(Pk,Mk,dump,b,w,serv);
4981  //ne selectionne que les non nuls
4982  if(secd.size()==0)
4983  {
4984  typPk redspol;
4985  Dump(Pk,dump);
4986 
4987  ReduceSpol(S,dump,redspol,b,w,serv);
4988 #ifdef NO_TWICE
4989  for(typename typPk::iterator iter=redspol.begin();iter!=redspol.end();iter++)
4990  P.push_back(invconv<polyalp>(*iter));
4991 #endif
4992  if(redspol.size()==0)
4993  {
4994  //ici on modifie Pk dans AddB
4995  //<=> forget(b,dump,k) AddB(..) Dump(Pk,dump);
4996  k++;
4997  L.destroystore();
4998  continue;
4999  //la matrice a pas rang plein second membre =0
5000  //ET les Spol de ce niveau se ,reduisent a 0
5001  }
5002  else
5003  {
5004  //on a des spol qui se reduisent pas a 0
5005  //cad un defaut de commutation
5006  k=mindeg(redspol,serv);
5007  recover(Pk,Mk,dump,k);
5008 
5009  forget(b,dump,k,recall);
5010 
5011  mdebug()<<"il ya un passage par mat rang pas"
5012  <<"plein et pol non zero";
5013 
5014  NewCrochetchoix(P,Pk,Mk,redspol,k,dump,b,recall,w,serv);
5015 
5016  NewDestroyredspol(redspol,Pk);
5017  serv.compress(Pk,dump,b,k);
5018  Dump(Pk,dump);
5019 
5020  L.destroystore();
5021  continue;
5022  }
5023  }
5024  else
5025  {
5026  mdebug()<<"entre en bucle une de los pol ne reduise pas a 0";
5027  for(typename Base::const_iterator p=b.begin();p!=b.end();p++)
5028  {
5029  for(int i=0;i<p->taille1;i++)
5030  {
5031  mdebug()<<"p.accept[i] "<<p->accept[i];
5032  }
5033  for(int i=0;i<p->taille2;i++)
5034  {
5035  mdebug()<<"p.refuse[i] "<<p->refuse[i];
5036  }
5037  }
5038 
5039  //on a des pol qui se reduisent pas a 0
5040  //cad un defaut de commuatation dans les coins
5041  typPk redspol;
5042 
5043  Dump(Pk,dump);
5044 
5045 
5046  ReduceSpol(S,dump,redspol,b,w,serv);//,Spoldejavu);
5047  my_merge(secd,redspol);
5048 #ifdef NO_TWICE
5049  for(typename typPk::iterator iter=redspol.begin();iter!=redspol.end();iter++)
5050  P.push_back(invconv<polyalp>(*iter));
5051 #endif
5052  k=mindeg(secd,serv);
5053  recover(Pk,Mk,dump,k);
5054 
5055  forget(b,dump,k,recall);
5056 
5057 
5058  for(typename Base::const_iterator p=b.begin();p!=b.end();p++)
5059  {
5060  for(int i=0;i<p->taille1;i++)
5061  {
5062  mdebug()<<"p.accept[i] "<<p->accept[i];
5063  }
5064  for(int i=0;i<p->taille2;i++)
5065  {
5066  mdebug()<<"p.refuse[i] "<<p->refuse[i];
5067  }
5068  }
5069  NewCrochetchoix(P,Pk,Mk,secd,k,dump,b,recall,w,serv);
5070 
5071  NewDestroyredspol(secd,Pk);
5072 
5073  serv.compress(Pk,dump,b,k);
5074  Dump(Pk,dump);
5075 
5076  L.destroystore();
5077  continue;
5078  }
5079  }
5080  }
5081 
5082  Mat.destroystore();
5083 
5084 #ifdef DEBUG
5085  for(typename Base::const_iterator p=b.begin();p!=b.end();p++)
5086  {
5087  for(int i=0;i<p->taille1;i++)
5088  {
5089  mdebug()<<"p.accept[i] "<<p->accept[i];
5090  }
5091  for(int i=0;i<p->taille2;i++)
5092  {
5093  mdebug()<<"p.refuse[i] "<<p->refuse[i];
5094  }
5095  }
5096 #endif
5097 
5098  destroy_space(&w);
5099  // delete[] &sol2;
5100  // delete[] &sol3;
5101 
5102  mdebug()<<"sortie d'algo flash et dimension ";//<<critere(b,k);
5103  mdebug()<<"TPS SDP "<<tempsCDSP;
5104  return ;//dump;
5105  //comme ca on conserve le tri par degre des formes normales et
5106  //donc la construction des operateurs de multiplication se
5107  //fait plus vite.. :-)
5108 
5109 }
5110 
5111 
5112 
5113 #endif //ALREADY_corealgo3
monomial_server::int2mon
void int2mon(const int &i, mon &res)
Definition: placemon3.hpp:664
tmpclock
struct timeval initclock tmpclock
Definition: solver_bb_floating.cpp:67
workspace
Definition: workspace.hpp:11
remiseenforme
void remiseenforme(typpol &p)
Definition: ugly.hpp:2
lcm
mon lcm(const mon &m1, const mon &m2)
Definition: corealgo3sdpa.hpp:55
my_plus
void my_plus(pol &res, const pol &nf, const coeff toto, workspace *w)
Definition: corealgo3sdpa.hpp:500
dumpstruct
Definition: dump.hpp:2
init_algo1
void init_algo1(typP &P, int &nouveaumon, int &maxdeg, int &k, typPk &Pk, typMk &Mk, Base &b, typdump &dump, workspace &w, monomial_server &serv)
Definition: corealgo3sdpa.hpp:434
MPoly
Multivariate polynomials.
Definition: MPoly.hpp:28
AddB
void AddB(const typmat &mat, typPk &Pk, const typMk &Mk, const typexceed &exceed, const typsizeexceed &sizeexceed, const typdump &dump, typB &b, monomial_server &serv)
Definition: pol2terspecial.hpp:2470
affdebug
void affdebug(const typmat &m)
Definition: io.hpp:30
conv_merge
void conv_merge(typPk &Pk, typP &P, typexceed &ex, typsize &sizeex, const Base &b, monomial_server &serv)
Definition: corealgo3sdpa.hpp:1180
max
#define max(a, b)
Definition: alp_f2c.H:167
member
int member(const typMk &tmpMk, const typename typMk::value_type mon)
Definition: corealgo3sdpa.hpp:1291
Solve
list< polbis< mono, T > > & Solve(typmat &S, std::list< polbis< mono, T > > &TildePk, int **pr, int **pc, workspace &w)
Definition: pol2bisspecial.hpp:609
Dlexpolprod::operator()
bool operator()(const M &, const M &) const
Definition: corealgo3sdp.hpp:2484
monomial_server::mon2int
int mon2int(const mon &mm)
Definition: placemon3.hpp:679
isdivisible
int isdivisible(const mon &m1, const mon &m2, int nbvar)
Definition: isdivisible.hpp:5
proj
void proj(T &Crochet, const typdump &dump, const Base &b)
Definition: proj.hpp:51
lvar
int lvar(const dynamicexp< X, E > &A)
Definition: dynamicexp.hpp:178
free
void free(void *)
monomial_server::setorder
void setorder(ptrfun tmpord)
Definition: placemon3.hpp:301
read_sol_from_sdpa
void read_sol_from_sdpa(coeff *&sol, typmatmap &matmap, int dim, typserv &serv)
Definition: corealgo3sdpa.hpp:3491
monomial_server::int2mon_
mon * int2mon_
Definition: placemon3.hpp:21
computeprod
void computeprod(typprod &prod, typstock &stockmon, Base &b, typdump &dump)
Definition: corealgo3sdpa.hpp:2698
coinor
Definition: corealgo3sdp.hpp:14
forget0
void forget0(Base &b, typdump &dump, int k, typrecall &recall)
Definition: dump3sdp.hpp:197
findcorresppol
void findcorresppol(const mon &m, const typdump &dump, pol &res)
Definition: corealgo3sdpa.hpp:1792
pol
Definition: pol.hpp:6
generickernel
void generickernel(typP &newpol, int k, const typdump &dump, Base &b, typserv &serv)
Definition: corealgo3sdpa.hpp:2725
generickernelNEW
void generickernelNEW(typP &newpol, int k, const typdump &dump, Base &b, typserv &serv, typP F, coeff *&sol2, typP rest, coeff *&solfinal, int solver, int &method, int &flat, string param, double &tresh)
Definition: generkeropt.hpp:18
Base::iterator
list< predicat >::iterator iterator
Definition: types.hpp:20
GoodDir
int GoodDir(const typmon &m, int i, const Base &b)
Definition: corealgo3sdpa.hpp:1160
monomial_server::init_multiple
void init_multiple(int n)
Definition: placemon3.hpp:59
std
Definition: Scl_mpfr.hpp:827
Ispolzero
int Ispolzero(const pol &p)
Definition: Iszero.hpp:30
MatrixOf
typmat * MatrixOf(const typPk &L, const typMk &lm, const typmat &mat)
Definition: corealgo3sdpa.hpp:330
predicat::accept
mon * accept
Definition: types.hpp:6
my_ord
int my_ord(const mon &m1, const mon &m2)
Definition: ord.hpp:2
pobj
MSKsoltypee MSKrealt * pobj
Definition: mosek.h:3621
tempsdssdp
double tempsdssdp
Definition: solver_bb_floating.cpp:68
Crochet
typP Crochet(const typP &P, int k)
Definition: corealgo3sdpa.hpp:227
nf
typpol nf(int var, int indmon, const typdump &dump, const Base &b)
Definition: corealgo3sdpa.hpp:846
flag
long flag
Definition: alp_f2c.H:52
inverse
void inverse(int *&, int)
Base::const_iterator
list< predicat >::const_iterator const_iterator
Definition: types.hpp:21
update_crochetchoix
void update_crochetchoix(typiterpk &iterPk, typiterpk &iter, int tmpind, workspace &w)
Definition: corealgo3sdpa.hpp:2169
y
MSKcallbackcodee MSKsoltypee MSKprostae MSKsolstae MSKstakeye MSKstakeye MSKstakeye MSKrealt MSKrealt MSKrealt * y
Definition: mosek.h:2689
malloc
void * malloc(YYSIZE_T)
recover
void recover(typPk &Pk, typMk &Mk, const typdump &dump, int k)
Definition: dump.hpp:192
goodspol
int goodspol(const mon &m1, const mon &m2, typppcm &ppcm, const Base &b, const typMk &Spoldejavu)
Definition: corealgo3sdpa.hpp:150
workspace::tabcoeff
coeff * tabcoeff
Definition: workspace.hpp:12
testsdp
int testsdp(int k, const typdump &dump, const Base &b, const typserv &serv)
Definition: corealgo3sdpa.hpp:3795
predicat
Definition: types.hpp:3
Dump
void Dump(const typPk &Pk, typdump &dump)
Definition: dump.hpp:48
pol::coeff_t
T coeff_t
Definition: pol.hpp:9
dim
MSKint32t MSKCONST MSKint32t * dim
Definition: mosek.h:2407
predicat::taille1
int taille1
Definition: types.hpp:5
selectnozero
void selectnozero(typP &ttmp)
Definition: Iszero.hpp:40
Base::nbvar
int nbvar()
Definition: types.hpp:37
init_algo
void init_algo(typP &P, int &nouveaumon, int &maxdeg, int &k, typPk &Pk, typMk &Mk, Base &b, typdump &dump, workspace &w, monomial_server &serv, Base &b2)
Definition: corealgo3sdpa.hpp:366
mon2int
int mon2int(const mon &mm)
Definition: placemon.hpp:294
monomial_server::reinit_multiple
void reinit_multiple()
Definition: placemon3.hpp:85
write_to_sdpa_fic
void write_to_sdpa_fic(typmatmap &matrixmap, int &dim, typserv &serv, int d)
Definition: corealgo3sdpa.hpp:3429
j
MSKconetypee MSKrealt MSKint32t MSKint32t j
Definition: mosek.h:2421
workspace::tabnf
unsigned char * tabnf
Definition: workspace.hpp:13
stab
void stab(typPk &Pk)
Definition: corealgo3sdpa.hpp:31
destroy_space
void destroy_space(workspace< coeff > *w)
Definition: workspace.hpp:34
Dlexpolprod::less
static bool less(const M &, const M &)
Definition: corealgo3sdp.hpp:2489
int2mon
void int2mon(const int &i, mon &res)
Definition: placemon.hpp:288
mdebug
Definition: mdebug.hpp:11
k
MSKint32t k
Definition: mosek.h:2713
Monom::rep
R rep
Definition: Monom.hpp:30
rajoute_voisins
void rajoute_voisins(typPk &Pk, const typMk &Mk, const typdump &dump, Base &b, workspace &w, monomial_server &serv)
Definition: IsinB3.hpp:111
Monom::coeff_t
C coeff_t
Definition: Monom.hpp:26
init_space
void init_space(workspace< coeff > *w)
Definition: workspace.hpp:43
halfred
pol halfred(const mon &targetind, const pol &p, const typdump &dump, const Base &b, monomial_server &serv)
Definition: corealgo3sdpa.hpp:1843
minmaxdegP
int minmaxdegP(typP P, int &maxdeg)
Definition: corealgo3sdpa.hpp:215
nbbits
int nbbits[256]
Definition: pol2ter.hpp:4
param
MSKdparame param
Definition: mosek.h:2743
mult
typpol mult(int i, const typpol &pol, const typdump &dump, const Base &b, mon **exce, int &exsize, workspace &w, monomial_server &serv)
Definition: corealgo3sdpa.hpp:911
algo
void algo(typP P, typdump &dump, Base &b, monomial_server &serv)
Definition: corealgo3sdpa.hpp:3801
my_divmon
Mon my_divmon(const Mon &mon, int i, const Base &b)
Definition: corealgo3sdpa.hpp:1144
monomial_server
Definition: placemon3.hpp:14
indice2
int indice2(int i, int j, const Base &b)
Definition: corealgo3sdpa.hpp:826
varmult
void varmult(iterator iter1, iterator iter2, mon &var1, mon &var2, const Base &b)
Definition: corealgo3sdpa.hpp:1830
index
MSKCONST char MSKint32t MSKint32t * index
Definition: mosek.h:2860
Degree
dynamicexp< X, E >::degree_t Degree(const dynamicexp< X, E > &t)
Definition: dynamicexp.hpp:91
printmatfile_pequan
void printmatfile_pequan(const typmat< Scl< MPF > > &m, FILE *fic)
Definition: io.hpp:381
Base::end
iterator end()
Definition: types.hpp:30
IsinB
int IsinB(const mon &m, const Base &b)
Definition: IsinB3.hpp:68
Base::def
list< predicat > def
Definition: types.hpp:16
Crochetchoix
void Crochetchoix(const typP &P, typPk &Pk, typPk &redspol, int k, Base &b)
enumeratemon
void enumeratemon(typstock &stockmon, int k, const Base &b, const typdump &dump)
Definition: corealgo3sdpa.hpp:2565
Setnvars
int Setnvars(const typP &tmpPk)
Definition: corealgo3sdpa.hpp:266
mon
Mon mon
Definition: solver_bb_floating.cpp:136
Base::nbvars
int nbvars
Definition: types.hpp:17
monomial_server::sizeplacemon
int sizeplacemon()
Definition: placemon3.hpp:313
increase_space
void increase_space(workspace< coeff > *w)
Definition: workspace.hpp:19
enumeratemonvar
void enumeratemonvar(typstock &stockmon, int k, const Base &b, const typdump &dump, int n)
Definition: corealgo3sdpa.hpp:2615
predicat::taille2
int taille2
Definition: types.hpp:7
Base
Definition: types.hpp:14
Listedejafait
void Listedejafait(const typMk &Mk, typMk &Spoldejavu)
Definition: corealgo3sdpa.hpp:2500
Monom::GetDegree
exponent_t GetDegree() const
Definition: Monom.hpp:70
dmax
#define dmax(a, b)
Definition: alp_f2c.H:169
ReduceSpol
void ReduceSpol(typPk &S, const typdump &dump, typPk &redspol, const Base &b, workspace &w, monomial_server &serv)
Definition: corealgo3sdpa.hpp:2045
enumeratemoninv
void enumeratemoninv(typstock &stockmon, int k, const Base &b, const typdump &dump)
Definition: corealgo3sdpa.hpp:2518
c
MSKrealt * c
Definition: mosek.h:2678
i
MSKint32t MSKint32t MSKint32t i
Definition: mosek.h:2278
Iszero
int Iszero(const Scl< MPQ > &c)
Definition: Iszero.hpp:14
dobj
MSKsoltypee MSKrealt MSKrealt MSKrealt MSKrealt MSKrealt MSKrealt MSKrealt * dobj
Definition: mosek.h:3621
Poly
MPoly< vector< Mon >, my_Dlex< Mon > > Poly
Definition: solver_bb_floating.cpp:138
workspace::sizenf
int sizenf
Definition: workspace.hpp:15
abs
#define abs(x)
Definition: f2c.H:161
Mon
Monom< COEFF, dynamicexp<'x'> > Mon
Definition: solver_bb_floating.cpp:134
dgesvd_
int dgesvd_(char *jobu, char *jobvt, long *m, long *n, T *a, long *lda, T *s, T *u, long *ldu, T *vt, long *ldvt, T *work, long *lwork, long *info)
mpoly.hpp
pol::size
int size
Definition: pol.hpp:10
a
MSKstreamtypee MSKint32t MSKint32t MSKint32t MSKint32t MSKint32t MSKint32t MSKint32t MSKint32t MSKint32t a
Definition: mosek.h:3833
Dlexpolprod
Definition: corealgo3sdp.hpp:2478
compute_ppcm
int compute_ppcm(typppcm &ppcm, const Base &b)
Definition: corealgo3sdpa.hpp:67
COEFF
#define COEFF
Definition: symbo2.hpp:8
forget
void forget(Base &b, typdump &dump, int k, typrecall &recall)
Definition: dump.hpp:75
GB_stop
int GB_stop(int k)
Definition: corealgo3sdpa.hpp:142
predicat::refuse
mon * refuse
Definition: types.hpp:8
generickernelsdpa
void generickernelsdpa(typP &newpol, int k, const typdump &dump, Base &b, typserv &serv)
Definition: corealgo3sdpa.hpp:3585
monomial_server::mulind
int mulind(int indexbase, int var, const typBase &b)
Definition: placemon3.hpp:126
AddB_dim
void AddB_dim(const typmat &mat, const typmat &L, typP &P, typPk &Pk, const typMk &Mk, const typdump &dump, typB &b)
Definition: IsinB3.hpp:545
NewCrochetchoix
void NewCrochetchoix(typP &P, typPk &Pk, typMk &Mk, typPk &redspol, int &k, typdump &dump, Base &b, typrecall &recall, workspace &w, monomial_server &serv)
Definition: corealgo3sdpa.hpp:2211
newPkMkMatrixof
int newPkMkMatrixof(typMat &mat, typPk &Pk, typP &P, typMk &Mk, int &k, Base &b, typdump &dump, typrecall &recall, workspace &w, monomial_server &serv)
Definition: corealgo3sdpa.hpp:1303
Setnbvars
void Setnbvars(const typP &tmpPk, Base &b)
Definition: corealgo3sdpa.hpp:251
Base::begin
iterator begin()
Definition: types.hpp:25
my_merge
void my_merge(T &l1, const T &l2)
Definition: corealgo3sdpa.hpp:48
Base::const_reverse_iterator
list< predicat >::const_reverse_iterator const_reverse_iterator
Definition: types.hpp:23
Base::mon
predicat::monom_t mon
Definition: types.hpp:19
pol
#define pol
Definition: pol2ter.hpp:3
tempsCDSP
double tempsCDSP
Definition: solver_bb_floating.cpp:69
Spolordsup
typPk Spolordsup(const typPk &Pk, const Base &b, const typMk &Spoldejavu)
Definition: corealgo3sdpa.hpp:188
mindeg
int mindeg(const typPk &l, monomial_server &serv)
Definition: corealgo3sdpa.hpp:2135
Monom::SetCoeff
void SetCoeff(const C &c)
Definition: Monom.hpp:68
DoBaseMk
void DoBaseMk(const typP &tmpPk, typMk &Mk, Base &b, monomial_server &serv)
Definition: corealgo3sdpa.hpp:283
monomial_server::compress
void compress(typPk &Pk, typdump &dump, const Base &b, int k)
Definition: placemon3.hpp:1076
pol::ind
mon ind
Definition: pol.hpp:11
Monom
Multivariate monomials.
Definition: Monom.hpp:21
Base::erase
iterator erase(iterator first, iterator last)
Definition: types.hpp:35
pol::nf
T * nf
Definition: pol.hpp:12
Base::rbegin
reverse_iterator rbegin()
Definition: types.hpp:27
getperm
void getperm(int *perm, int sizeperm, int *indextab)
Definition: corealgo3sdpa.hpp:891
NewDestroyredspol
void NewDestroyredspol(typPk &redspol, const typPk &Pk)
Definition: corealgo3sdpa.hpp:2487
name
MSKint32t MSKint32t char * name
Definition: mosek.h:2853
Home  |  Download & Install |  Contributions
Generated on Tue Sep 29 2015 11:43:20 for Borderbasix by 1.8.9.1