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,°,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