1 /* -*- mode:C++ ; compile-command: "g++ -I.. -g -c gausspol.cc" -*- */
2 /*
3  *  Copyright (C) 2000,2014 B. Parisse, Institut Fourier, 38402 St Martin d'Heres
4  *
5  *  This program is free software; you can redistribute it and/or modify
6  *  it under the terms of the GNU General Public License as published by
7  *  the Free Software Foundation; either version 3 of the License, or
8  *  (at your option) any later version.
9  *
10  *  This program is distributed in the hope that it will be useful,
11  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
12  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
13  *  GNU General Public License for more details.
14  *
15  *  You should have received a copy of the GNU General Public License
16  *  along with this program. If not, see <http://www.gnu.org/licenses/>.
17  */
18 #ifndef _GIAC_GAUSSPOL_H_
19 #define _GIAC_GAUSSPOL_H_
20 #include "first.h"
21 #include "poly.h"
22 #include "gen.h"
23 #ifdef HAVE_PTHREAD_H
24 #include <pthread.h>
25 #endif
26 
27 #ifndef NO_NAMESPACE_GIAC
28 namespace giac {
29 #endif // ndef NO_NAMESPACE_GIAC
30 
31 #ifdef USE_GMP_REPLACEMENTS
32 #undef HAVE_GMPXX_H
33 #undef HAVE_LIBMPFR
34 #endif
35 
36   extern const int primes[100] ;
37 
38   class nfactor {
39   public:
40     gen fact;
41     int mult;
nfactor()42     nfactor():fact(1),mult(0) {}
nfactor(const nfactor & n)43     nfactor(const nfactor &n) : fact(n.fact),mult(n.mult) {}
nfactor(const gen & n,int m)44     nfactor(const gen & n, int m) : fact(n),mult(m) {}
45   };
46   std::vector<nfactor> trivial_n_factor(gen &n);
47   vecteur cyclotomic(int n);
48 
49   // gen pow(const gen & n,int k);
50   typedef std::vector< monomial<gen> > monomial_v;
51   typedef tensor<gen> polynome;
52 
53   // function on polynomials
54   polynome gen2polynome(const gen & e,int dim);
55   // check type of coefficients
56   int coefftype(const polynome & p,gen & coefft);
57   // remove MODulo coefficients
58   void unmodularize(const polynome & p,polynome & res);
59   polynome unmodularize(const polynome & p);
60   void modularize(polynome & d,const gen & m);
61   // remove EXT, also checks that pmin is the min poly
62   bool unext(const polynome & p,const gen & pmin,polynome & res);
63   bool ext(polynome & res,const gen & pmin);
64   void ext(const polynome & p,const gen & pmin,polynome & res);
65   // arithmetic
66   bool is_one(const polynome & p);
67   bool operator < (const polynome & f,const polynome & g);
68   bool operator < (const facteur<polynome> & f,const facteur<polynome> & g);
69   polynome firstcoeff(const polynome & p);
70   void Add_gen ( std::vector< monomial<gen> >::const_iterator & a,
71 		 std::vector< monomial<gen> >::const_iterator & a_end,
72 		 std::vector< monomial<gen> >::const_iterator & b,
73 		 std::vector< monomial<gen> >::const_iterator & b_end,
74 		 std::vector< monomial<gen> > & new_coord,
75 		 bool (* is_strictly_greater)( const index_m &, const index_m &)) ;
76   polynome operator + (const polynome & th,const polynome & other);
77   void Sub_gen ( std::vector< monomial<gen> >::const_iterator & a,
78 		 std::vector< monomial<gen> >::const_iterator & a_end,
79 		 std::vector< monomial<gen> >::const_iterator & b,
80 		 std::vector< monomial<gen> >::const_iterator & b_end,
81 		 std::vector< monomial<gen> > & new_coord,
82 		 bool (* is_strictly_greater)( const index_m &, const index_m &)) ;
83   polynome operator - (const polynome & th,const polynome & other);
84   // Fast multiplication using hash maps, might also use an int for reduction
85   // but there is no garantee that res is smod-ed modulo reduce
86   // Use reduce=0 for non modular
87   void mulpoly (const polynome & th, const polynome & other,polynome & res,const gen & reduce);
88   polynome operator * (const polynome & th, const polynome & other) ;
89   polynome & operator *= (polynome & th, const polynome & other) ;
90   void Mul_gen ( std::vector< monomial<gen> >::const_iterator & ita,
91 		 std::vector< monomial<gen> >::const_iterator & ita_end,
92 		 std::vector< monomial<gen> >::const_iterator & itb,
93 		 std::vector< monomial<gen> >::const_iterator & itb_end,
94 		 std::vector< monomial<gen> > & new_coord,
95 		 bool (* is_strictly_greater)( const index_t &, const index_t &),
96 		 const std::pointer_to_binary_function < const monomial<gen> &, const monomial<gen> &, bool> m_is_greater
97 		 ) ;
98   void mulpoly(const polynome & th,const gen & fact,polynome & res);
99   polynome operator * (const polynome & th, const gen & fact) ;
100   inline polynome operator * (const gen & fact, const polynome & th){ return th*fact; }
101   // a*b+c*d
102   gen foisplus(const polynome & a,const polynome & b,const polynome & c,const polynome & d);
103   gen foisplus(const gen & a,const gen & b,const gen & c,const gen & d);
104   polynome operator - (const polynome & th) ;
105   polynome operator / (const polynome & th,const polynome & other);
106   polynome operator / (const polynome & th,const gen & fact );
107   polynome operator % (const polynome & th,const polynome & other);
108   polynome operator % (const polynome & th, const gen & modulo);
109   polynome re(const polynome & th);
110   polynome im(const polynome & th);
111   polynome conj(const polynome & th);
112   polynome poly1_2_polynome(const vecteur & v, int dimension);
113   void polynome2poly1(const polynome & p,int var,vecteur & v);
114   vecteur polynome12poly1(const polynome & p);
115   int inner_POLYdim(const vecteur & v);
116   vecteur polynome2poly1(const polynome & p,int var);
117   vecteur polynome2poly1(const polynome & p); // for algebraic ext.
118   gen polynome2poly1(const gen & e,int var);
119   void poly12polynome(const vecteur & v, int var,polynome & p,int dimension=0);
120   polynome poly12polynome(const vecteur & v,int var,int dimension=0);
121   polynome poly12polynome(const vecteur & v);
122   gen untrunc(const gen & e,int degree,int dimension);
123   gen vecteur2polynome(const vecteur & v,int dimension);
124   bool divrem1(const polynome & a,const polynome & b,polynome & quo,polynome & r,int exactquo=0,bool allowrational=false) ;
125   bool divrem (const polynome & th, const polynome & other, polynome & quo, polynome & rem, bool allowrational = false );
126   bool divremmod (const polynome & th,const polynome & other, const gen & modulo,polynome & quo, polynome & rem);
127   bool exactquotient(const polynome & a,const polynome & b,polynome & quo,bool allowrational=true);
128   bool powpoly (const polynome & th, int u,polynome & res);
129   polynome pow(const polynome & th,int n);
130   bool is_positive(const polynome & p);
131   polynome pow(const polynome & p,const gen & n);
132   polynome powmod(const polynome &p,int n,const gen & modulo);
133   polynome gcd(const polynome & p,const polynome & q);
134   void gcd(const polynome & p,const polynome & q,polynome & d);
135   void lcmdeno(const polynome & p, gen & res);
136   gen lcoeffn(const polynome & p);
137   gen lcoeff1(const polynome & p);
138   polynome ichinrem(const polynome &p,const polynome & q,const gen & pmod,const gen & qmod);
139   // set i to i+(j-i)*B mod A, inplace operation
140   void ichrem_smod_inplace(mpz_t * Az,mpz_t * Bz,mpz_t * iz,mpz_t * tmpz,gen & i,const gen & j);
141   polynome resultant(const polynome & p,const polynome & q);
142   bool resultant_sylvester(const polynome &p,const polynome &q,matrice & S,polynome & res);
143   bool resultant_sylvester(const polynome &p,const polynome &q,vecteur &pv,vecteur &qv,matrice & S,gen & determinant);
144   polynome lgcd(const polynome & p);
145   gen ppz(polynome & p,bool divide=true);
146   void lgcdmod(const polynome & p,const gen & modulo,polynome & pgcd);
147   polynome gcdmod(const polynome &p,const polynome & q,const gen & modulo);
148   polynome content1mod(const polynome & p,const gen & modulo,bool setdim=true);
149   void contentgcdmod(const polynome &p, const polynome & q, const gen & modulo, polynome & cont,polynome & prim);
150   polynome pp1mod(const polynome & p,const gen & modulo);
151   // modular gcd via PSR, used when not enough eval points available
152   // a and b must be primitive and will be scratched
153   void psrgcdmod(polynome & a,polynome & b,const gen & modulo,polynome & prim);
154   // Find non zeros coeffs of p, res contains the positions of non-0 coeffs
155   int find_nonzero(const polynome & p,index_t & res);
156   polynome pzadic(const polynome &p,const gen & n);
157   bool gcd_modular_algo(polynome &p,polynome &q, polynome &d,bool compute_cof);
158   bool listmax(const polynome &p,gen & n );
159   bool gcdheu(const polynome &p,const polynome &q, polynome & p_simp, gen & np_simp, polynome & q_simp, gen & nq_simp, polynome & d, gen & d_content ,bool skip_test=false,bool compute_cofactors=true);
160   polynome gcdpsr(const polynome &p,const polynome &q,int gcddeg=0);
161   void simplify(polynome & p,polynome & q,polynome & p_gcd);
162   polynome simplify(polynome &p,polynome &q);
163   void egcdlgcd(const polynome &p1, const polynome & p2, polynome & u,polynome & v,polynome & d);
164   void egcdpsr(const polynome &p1, const polynome & p2, polynome & u,polynome & v,polynome & d);
165   void egcd(const polynome &p1, const polynome & p2, polynome & u,polynome & v,polynome & d);
166   // Input a,b,c,u,v,d such that a*u+b*v=d,
167   // Output u,v,C such that a*u+b*v=c*C
168   void egcdtoabcuv(const tensor<gen> & a,const tensor<gen> &b, const tensor<gen> &c, tensor<gen> &u,tensor<gen> &v, tensor<gen> & d, tensor<gen> & C);
169 
170   bool findabcdelta(const polynome & p,polynome & a,polynome &b,polynome & c,polynome & delta);
171   bool findde(const polynome & p,polynome & d,polynome &e);
172   factorization sqff(const polynome &p );
173   // factorize a square-free univariate polynomial
174   bool sqfffactor(const polynome &p, vectpoly & v,bool with_sqrt,bool test_composite,bool complexmode);
175   bool sqff_evident(const polynome & p,factorization & f,bool withsqrt,bool complexmode);
176   // factorization over Z[i]
177   bool cfactor(const polynome & p, gen & an,factorization & f,bool withsqrt,gen & extra_div);
178   // add a dimension in front of pcur for algebraic extension variable
179   bool algext_convert(const polynome & pcur,const gen & e,polynome & p_y);
180   // convert minimal polynomial of algebraic extension
181   void algext_vmin2pmin(const vecteur & v_mini,polynome & p_mini);
182 
183   // factorization over an algebraic extension
184   // the main variable of G is the algebraic extension variable
185   // the minimal polynomial of this variable is p_mini
186   // G is assumed to be square-free
187   // See algorithm 3.6.4 in Henri Cohen book starting at step 3
188   bool algfactor(const polynome & G,const polynome & p_mini,int & k,factorization & f,bool complexmode,gen & extra_div,polynome & Gtry);
189   // sqff factorization over a finite field
190   factorization squarefree_fp(const polynome & p,unsigned n,unsigned exposant);
191   // univariate factorization over a finite field, once sqff
192   bool sqff_ffield_factor(const factorization & sqff_f,int n,environment * env,factorization & f);
193   // p is primitive wrt the main var
194   bool mod_factor(const polynome & p_orig,polynome & p_content,int n,factorization & f);
195 
196   // factorization over Z[e] where e is an algebraic extension
197   bool ext_factor(const polynome &p,const gen & e,gen & an,polynome & p_content,factorization & f,bool complexmode,gen &extra_div);
198   // factorization over Z[coeff_of_p]
199   bool factor(const polynome &p,polynome & p_content,factorization & f,bool isprimitive,bool withsqrt,bool complexmode,const gen & divide_by_an,gen & extra_div);
200   void unitarize(const polynome &pcur, polynome &unitaryp, polynome & an);
201   polynome ununitarize(const polynome & unitaryp, const polynome & an);
202   void partfrac(const polynome & num, const polynome & den, const std::vector< facteur< polynome > > & v , std::vector < pf <gen> > & pfde_VECT, polynome & ipnum, polynome & ipden, bool rational=true );
203   pf<gen> intreduce_pf(const pf<gen> & p_cst, std::vector< pf<gen> > & intde_VECT ,bool residue=false);
204   // Sturm sequences
205   vecteur vector_of_polynome2vecteur(const vectpoly & v);
206   vecteur sturm_seq(const polynome & p,polynome & cont);
207 
208   // prototype of factorization of univariate sqff unitary polynomial
209   // provided e.g. by smodular
210   bool factorunivsqff(const polynome & q,environment * env,vectpoly & v,int & ithprime,int debug,int modfactor_primes);
211   // find linear factor only
212   int linearfind(const polynome & q,environment * env,polynome & qrem,vectpoly & v,int & ithprime);
213   // prototype of modular 1-d gcd algorithm
214   bool gcd_modular_algo1(polynome &p,polynome &q,polynome &d,bool compute_cof);
215   polynome smod(const polynome & th, const gen & modulo);
216   void smod(const polynome & th, const gen & modulo,polynome & res);
217   bool gcdmod_dim1(const polynome &p,const polynome & q,const gen & modulo,polynome & d,polynome & pcof,polynome & qcof,bool compute_cof,bool & real);
218 
219   // evaluate p at v by replacing the last variables of p with values of v
220   gen peval(const polynome & p,const vecteur & v,const gen &m,bool simplify_at_end=false,std::vector<int_unsigned> * pptr=0);
221   int total_degree(const polynome & p);
222 
223   // build a multivariate poly
224   // with normal coeff from a multivariate poly with multivariate poly coeffs
225   polynome unsplitmultivarpoly(const polynome & p,int inner_dim);
226   polynome splitmultivarpoly(const polynome & p,int inner_dim);
227   polynome split(const polynome & p,int inner_dim);
228 
229 
230   template <class U>
convert_myint(const polynome & p,const index_t & deg,std::vector<T_unsigned<my_mpz,U>> & v)231   bool convert_myint(const polynome & p,const index_t & deg,std::vector< T_unsigned<my_mpz,U> >  & v){
232     typename std::vector< monomial<gen> >::const_iterator it=p.coord.begin(),itend=p.coord.end();
233     v.clear();
234     v.reserve(itend-it);
235     U u;
236     my_mpz tmp;
237     index_t::const_iterator itit,ditbeg=deg.begin(),ditend=deg.end(),dit;
238     T_unsigned<my_mpz,U> gu;
239     for (;it!=itend;++it){
240       u=0;
241       itit=it->index.begin();
242       for (dit=ditbeg;dit!=ditend;++itit,++dit)
243 	u=u*unsigned(*dit)+unsigned(*itit);
244       gu.u=u;
245       if (it->value.type==_ZINT)
246 	mpz_set(gu.g.ptr,*it->value._ZINTptr);
247       else {
248 	if (it->value.type!=_INT_)
249 	  return false;
250 	mpz_set_si(gu.g.ptr,it->value.val);
251       }
252       v.push_back(gu);
253     }
254     return true;
255   }
256 
257 
258 #ifdef HAVE_GMPXX_H
259   template <class U>
convert_myint(const polynome & p,const index_t & deg,std::vector<T_unsigned<mpz_class,U>> & v)260   bool convert_myint(const polynome & p,const index_t & deg,std::vector< T_unsigned<mpz_class,U> >  & v){
261     typename std::vector< monomial<gen> >::const_iterator it=p.coord.begin(),itend=p.coord.end();
262     v.clear();
263     v.reserve(itend-it);
264     U u;
265     index_t::const_iterator itit,ditbeg=deg.begin(),ditend=deg.end(),dit;
266     for (;it!=itend;++it){
267       u=0;
268       itit=it->index.begin();
269       for (dit=ditbeg;dit!=ditend;++itit,++dit)
270 	u=u*unsigned(*dit)+unsigned(*itit);
271       T_unsigned<mpz_class,U> gu;
272       gu.u=u;
273       if (it->value.type==_ZINT){
274 	mpz_set(gu.g.get_mpz_t(),*it->value._ZINTptr);
275       }
276       else {
277 	if (it->value.type!=_INT_)
278 	  return false;
279 	gu.g=it->value.val;
280       }
281       v.push_back(gu);
282     }
283     return true;
284   }
285 #endif
286 
coeff_type(const std::vector<T_unsigned<gen,U>> & p,unsigned & maxint)287   template<class U> int coeff_type(const std::vector< T_unsigned<gen,U> > & p,unsigned & maxint){
288     maxint=0;
289     typename std::vector< T_unsigned<gen,U> >::const_iterator it=p.begin(),itend=p.end();
290     if (it==itend)
291       return -1;
292     int t=it->g.type,tt;
293     register int tmp;
294     for (++it;it!=itend;++it){
295       tt=it->g.type;
296       if (tt!=t)
297 	return -1;
298       if (!tt){
299 	if (it->g.val>0)
300 	  tmp=it->g.val;
301 	else
302 	  tmp=-it->g.val;
303 	if (maxint<tmp)
304 	  maxint=tmp;
305       }
306     }
307     return t;
308   }
309 
310   bool is_integer_poly(const polynome & p,bool intonly);
311 
312   template <class U>
convert_double(const polynome & p,const index_t & deg,std::vector<T_unsigned<double,U>> & v)313   bool convert_double(const polynome & p,const index_t & deg,std::vector< T_unsigned<double,U> >  & v){
314     typename std::vector< monomial<gen> >::const_iterator it=p.coord.begin(),itend=p.coord.end();
315     v.clear();
316     v.reserve(itend-it);
317     T_unsigned<double,U> gu;
318     U u;
319     index_t::const_iterator itit,ditbeg=deg.begin(),ditend=deg.end(),dit;
320     for (;it!=itend;++it){
321       u=0;
322       itit=it->index.begin();
323       for (dit=ditbeg;dit!=ditend;++itit,++dit)
324 	u=u*unsigned(*dit)+unsigned(*itit);
325       gu.u=u;
326       if (it->value.type!=_DOUBLE_)
327 	return false;
328       gu.g=it->value._DOUBLE_val;
329       v.push_back(gu);
330     }
331     return true;
332   }
333 
334   template <class U>
335   bool convert_int32(const polynome & p,const index_t & deg,std::vector< T_unsigned<int,U> > & v,int modulo=0){
336     typename std::vector< monomial<gen> >::const_iterator it=p.coord.begin(),itend=p.coord.end();
337     v.clear();
338     v.reserve(itend-it);
339     U u;
340     index_t::const_iterator itit,oldit,ititend,ditbeg=deg.begin(),ditend=deg.end(),dit;
341     for (;it!=itend;++it){
342       u=0;
343       oldit=itit=it->index.begin();
344       for (dit=ditbeg;dit!=ditend;++itit,++dit)
345 	u=u*unsigned(*dit)+unsigned(*itit);
346       if (it->value.type==_INT_){
347 	if (modulo)
348 	  v.push_back(T_unsigned<int,U>(it->value.val % modulo,u));
349 	else
350 	  v.push_back(T_unsigned<int,U>(it->value.val,u));
351       }
352       else {
353 	if (modulo && it->value.type==_ZINT)
354 	  v.push_back(T_unsigned<int,U>(smod(it->value,modulo).val,u));
355 	else
356 	  return false;
357       }
358       int nterms=*(itit-1);
359       if (nterms<=1 || nterms>=itend-it)
360 	continue;
361       itit = (it+nterms)->index.begin();
362       ititend = itit + p.dim-1;
363       if (*(ititend))
364 	continue;
365       for (;itit!=ititend;++oldit,++itit){
366 	if (*itit!=*oldit)
367 	  break;
368       }
369       if (itit!=ititend)
370 	continue;
371       // for dense poly, make all terms with the same x1..xn-1 powers
372       for (;nterms;--nterms){
373 	++it;
374 	--u;
375 	if (it->value.type==_INT_){
376 	  if (modulo)
377 	    v.push_back(T_unsigned<int,U>(it->value.val % modulo,u));
378 	  else
379 	    v.push_back(T_unsigned<int,U>(it->value.val,u));
380 	}
381 	else {
382 	  if (modulo && it->value.type==_ZINT)
383 	    v.push_back(T_unsigned<int,U>(smod(it->value,modulo).val,u));
384 	  else
385 	    return false;
386 	}
387       }
388     }
389     return true;
390   }
391 
392   template <class U>
convert_int(const polynome & p,const index_t & deg,std::vector<T_unsigned<longlong,U>> & v,longlong & maxp)393   bool convert_int(const polynome & p,const index_t & deg,std::vector< T_unsigned<longlong,U> >  & v,longlong & maxp){
394     typename std::vector< monomial<gen> >::const_iterator it=p.coord.begin(),itend=p.coord.end();
395     v.clear();
396     v.reserve(itend-it);
397     T_unsigned<longlong,U> gu;
398     U u;
399     maxp=0;
400     longlong tmp;
401     mpz_t tmpz;
402     mpz_init(tmpz);
403     index_t::const_iterator itit,ditbeg=deg.begin(),ditend=deg.end(),dit;
404     for (;it!=itend;++it){
405       u=0;
406       itit=it->index.begin();
407       for (dit=ditbeg;dit!=ditend;++itit,++dit)
408 	u=u*unsigned(*dit)+unsigned(*itit);
409       gu.u=u;
410       if (it->value.type==_INT_)
411 	gu.g=it->value.val;
412       else {
413 	if (it->value.type!=_ZINT || mpz_sizeinbase(*it->value._ZINTptr,2)>62){
414 	  mpz_clear(tmpz);
415 	  return false;
416 	}
417 	mpz2longlong(it->value._ZINTptr,&tmpz,gu.g);
418       }
419       tmp=gu.g>0?gu.g:-gu.g;
420       if (tmp>maxp)
421 	maxp=tmp;
422       v.push_back(gu);
423     }
424     mpz_clear(tmpz);
425     return true;
426   }
427 
428 #ifdef INT128
429   template <class U>
convert_int(const polynome & p,const index_t & deg,std::vector<T_unsigned<int128_t,U>> & v,int128_t & maxp)430   bool convert_int(const polynome & p,const index_t & deg,std::vector< T_unsigned<int128_t,U> >  & v,int128_t & maxp){
431     typename std::vector< monomial<gen> >::const_iterator it=p.coord.begin(),itend=p.coord.end();
432     v.clear();
433     v.reserve(itend-it);
434     T_unsigned<int128_t,U> gu;
435     U u;
436     maxp=0;
437     int128_t tmp;
438     mpz_t tmpz;
439     mpz_init(tmpz);
440     index_t::const_iterator itit,ditbeg=deg.begin(),ditend=deg.end(),dit;
441     for (;it!=itend;++it){
442       u=0;
443       itit=it->index.begin();
444       for (dit=ditbeg;dit!=ditend;++itit,++dit)
445 	u=u*unsigned(*dit)+unsigned(*itit);
446       gu.u=u;
447       if (it->value.type==_INT_)
448 	gu.g=it->value.val;
449       else {
450 	if (it->value.type!=_ZINT || mpz_sizeinbase(*it->value._ZINTptr,2)>124){
451 	  mpz_clear(tmpz);
452 	  return false;
453 	}
454 	mpz2int128(it->value._ZINTptr,&tmpz,gu.g);
455       }
456       tmp=gu.g>0?gu.g:-gu.g;
457       if (tmp>maxp)
458 	maxp=tmp;
459       v.push_back(gu);
460     }
461     mpz_clear(tmpz);
462     return true;
463   }
464 #endif
465 
convert_longlong(const std::vector<T_unsigned<gen,U>> & p,std::vector<T_unsigned<longlong,U>> & pd)466   template<class U> void convert_longlong(const std::vector< T_unsigned<gen,U> > & p,std::vector< T_unsigned<longlong,U> > & pd){
467     typename std::vector< T_unsigned<gen,U> >::const_iterator it=p.begin(),itend=p.end();
468     pd.reserve(itend-it);
469     for (;it!=itend;++it)
470       pd.push_back(T_unsigned<longlong,U>(it->g.val,it->u));
471   }
472 
convert_from(const std::vector<T_unsigned<T,U>> & p,std::vector<T_unsigned<gen,U>> & pd)473   template<class T,class U> void convert_from(const std::vector< T_unsigned<T,U> > & p,std::vector< T_unsigned<gen,U> > & pd){
474     typename std::vector< T_unsigned<T,U> >::const_iterator it=p.begin(),itend=p.end();
475     pd.reserve(itend-it);
476     for (;it!=itend;++it){
477       if (it->g)
478 	pd.push_back(T_unsigned<gen,U>(gen(it->g),it->u));
479     }
480   }
481 
482   // mode=0: fill both, =1 fill the gen part, =2 fill the index_m part
483   template<class T,class U>
484   void convert_from(typename std::vector< T_unsigned<T,U> >::const_iterator it,typename std::vector< T_unsigned<T,U> >::const_iterator itend,const index_t & deg,typename std::vector< monomial<gen> >::iterator jt,int mode=0){
485     if (mode==1){
486       for (;it!=itend;++jt,++it){
487 	jt->value=gen(it->g);
488       }
489       return;
490     }
491     index_t::const_reverse_iterator ditbeg=deg.rbegin(),ditend=deg.rend(),dit;
492     int pdim=int(deg.size());
493     U u,prevu=0;
494     int k;
495     int count=0;
496 #if defined(GIAC_NO_OPTIMIZATIONS) || ((defined(VISUALC) || defined(__APPLE__)) && !defined(GIAC_VECTOR)) || defined __clang__ // || defined NSPIRE || defined(FXCG)
497     if (0){ count=0; }
498 #else
499     if (pdim<=POLY_VARS){
500       deg_t i[POLY_VARS+1];
501       i[0]=2*pdim+1;
502       deg_t * iitbeg=i+1,*iit,*iitback=i+pdim,*iitbackm1=iitback-1;
503       for (iit=iitbeg;iit!=iitback;++iit)
504 	*iit=0;
505       *iitback=0;
506       for (--prevu;it!=itend;++it,++jt){
507 	u=it->u;
508 	if (prevu<=u+*iitback){
509 	  *iitback -= deg_t(prevu-u);
510 	  prevu=u;
511 	}
512 	else {
513 	  if (pdim>1 && (*iitbackm1)>0 && prevu<=u+*ditbeg+*iitback){
514 	    --(*iitbackm1);
515 	    *iitback += deg_t((u+(*ditbeg))-prevu);
516 	    prevu=u;
517 	  }
518 	  else
519 	  {
520 	    prevu=u;
521 	    for (k=pdim,dit=ditbeg;dit!=ditend;++dit,--k){
522 	      // qr=div(u,*dit);
523 	      i[k]=u % (*dit); // qr.rem;
524 	      u= u / (*dit); // qr.quot;
525 	      count += pdim;
526 	    }
527 	  }
528 	}
529 	jt->index=i;
530 	if (mode)
531 	  continue;
532 	jt->value=gen(it->g);
533 	// p.coord.push_back(monomial<gen>(gen(it->g),i));
534       }
535     }
536 #endif
537     else {
538       index_t i(pdim);
539       index_t::iterator /*iitbeg=i.begin(),*/ iitback=i.end()-1,iitbackm1=iitback-1;
540       for (--prevu;it!=itend;++it,++jt){
541 	u=it->u;
542 	if (prevu<=u+*iitback){
543 	  *iitback -= short(prevu-u);
544 	  prevu=u;
545 	}
546 	else {
547 	  if (pdim>1 && (*iitbackm1)>0 && prevu<=u+*ditbeg+*iitback){
548 	    --(*iitbackm1);
549 	    *iitback += short((u+(*ditbeg))-prevu);
550 	    prevu=u;
551 	    // cerr << "/" << u << ":" << i << '\n';
552 	  }
553 	  else
554           {
555 	    prevu=u;
556 	    for (k=pdim-1,dit=ditbeg;dit!=ditend;++dit,--k){
557 	      // qr=div(u,*dit);
558 	      i[k]=u % (*dit); // qr.rem;
559 	      u= u / (*dit); // qr.quot;
560 	      count += pdim;
561 	      // i[k]=u % unsigned(*dit);
562 	      // u = u/unsigned(*dit);
563 	    }
564 	  }
565 	}
566 	jt->index=i;
567 	if (mode)
568 	  continue;
569 	jt->value=gen(it->g);
570 	// p.coord.push_back(monomial<gen>(gen(it->g),i));
571       }
572     }
573     if (debug_infolevel>5)
574       CERR << "Divisions: " << count << '\n';
575   }
576 
577 
578   template<class T,class U>
579   struct convert_t {
580     typename std::vector< T_unsigned<T,U> >::const_iterator it,itend;
581     const index_t * degptr;
582     typename std::vector< monomial<gen> >::iterator jt;
583     int mode;
584   };
585 
586   template<class T,class U>
do_convert_from(void * ptr)587   void * do_convert_from(void * ptr){
588     convert_t<T,U> * argptr = (convert_t<T,U> *) ptr;
589     convert_from<T,U>(argptr->it,argptr->itend,*argptr->degptr,argptr->jt,argptr->mode);
590     return 0;
591   }
592 
593   extern int threads;
594 
595   // conversion in parallel will not be much faster if T!=gen because
596   // we must allocate memory for each gen
597   template<class T,class U>
598   void convert_from(const std::vector< T_unsigned<T,U> > & v,const index_t & deg,polynome & p,bool threaded=false,bool coeff_apart=false){
599     typename std::vector< T_unsigned<T,U> >::const_iterator it=v.begin(),itend=v.end();
600     p.dim=int(deg.size());
601     // p.coord.clear(); p.coord.reserve(itend-it);
602     p.coord=std::vector< monomial<gen> >(itend-it);
603     std::vector< monomial<gen> >::iterator jt=p.coord.begin();
604     int nthreads=threads;
605     if (nthreads==1 || !threaded || p.dim>POLY_VARS){
606       convert_from<T,U>(it,itend,deg,jt,0);
607       return;
608     }
609 #if defined(HAVE_PTHREAD_H) && !defined(EMCC) // && !defined(__clang__)
610     unsigned taille=itend-it;
611     if (nthreads>1
612 	&& int(taille)>nthreads*1000
613 	){
614       pthread_t tab[nthreads];
615       std::vector< convert_t<T,U> > arg(nthreads);
616       if (coeff_apart){
617 	convert_from<T,U>(it,itend,deg,jt,1); // convert first coefficients
618 	if (debug_infolevel>5)
619 	  CERR << CLOCK()*1e-6 << " end coefficients conversion" << '\n';
620       }
621       for (int i=0;i<nthreads;i++){
622 	convert_t<T,U> tmp={it+i*(taille/nthreads),it+(i+1)*taille/nthreads,&deg,jt+i*(taille/nthreads),coeff_apart?2:0};
623 	if (i==nthreads-1){
624 	  tmp.itend=itend;
625 	  convert_from<T,U>(tmp.it,tmp.itend,deg,tmp.jt,tmp.mode);
626 	}
627 	else {
628 	  arg[i]=tmp;
629 	  int res=pthread_create(&tab[i],(pthread_attr_t *) NULL,do_convert_from<T,U>,(void *) &arg[i]);
630 	  if (res) // thread not created
631 	    convert_from<T,U>(tmp.it,tmp.itend,deg,tmp.jt,tmp.mode);
632 	}
633       }
634       for (int i=0;i<nthreads-1;++i){
635 	void * ptr;
636 	pthread_join(tab[i],&ptr);
637       }
638       return;
639     } // end if (nthreads>1)
640 #endif
641     convert_from<T,U>(it,itend,deg,jt,0);
642   }
643 
644 #ifndef NO_NAMESPACE_GIAC
645 } // namespace giac
646 #endif // ndef NO_NAMESPACE_GIAC
647 
648 #endif // _GIAC_GAUSSPOL_H_
649