1 /* -*- mode:C++ ; compile-command: "g++-3.4 -I. -I.. -I../include -g -c ezgcd.cc -DHAVE_CONFIG_H -DIN_GIAC" -*- */ 2 3 #include "giacPCH.h" 4 /* Multivariate GCD for large data not covered by the heuristic GCD algo 5 * Copyright (C) 2000,7 B. Parisse, Institut Fourier, 38402 St Martin d'Heres 6 * 7 * This program is free software; you can redistribute it and/or modify 8 * it under the terms of the GNU General Public License as published by 9 * the Free Software Foundation; either version 3 of the License, or 10 * (at your option) any later version. 11 * 12 * This program is distributed in the hope that it will be useful, 13 * but WITHOUT ANY WARRANTY; without even the implied warranty of 14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 15 * GNU General Public License for more details. 16 * 17 * You should have received a copy of the GNU General Public License 18 * along with this program. If not, see <http://www.gnu.org/licenses/>. 19 */ 20 using namespace std; 21 #include "threaded.h" 22 #include "ezgcd.h" 23 #include "sym2poly.h" 24 #include "gausspol.h" 25 #include "modpoly.h" 26 #include "monomial.h" 27 #include "derive.h" 28 #include "subst.h" 29 #include "solve.h" 30 #include "giacintl.h" 31 32 #ifndef NO_NAMESPACE_GIAC 33 namespace giac { 34 #endif // ndef NO_NAMESPACE_GIAC 35 add_dim(monomial<gen> & m,int d)36 static void add_dim(monomial<gen> & m,int d){ 37 index_t i(m.index.iref()); 38 for (int j=0;j<d;++j) 39 i.push_back(0); 40 m.index=i; 41 } 42 change_dim(polynome & p,int dim)43 void change_dim(polynome & p,int dim){ 44 vector< monomial<gen> >::iterator it=p.coord.begin(),itend=p.coord.end(); 45 if (p.dim>=dim){ 46 p.dim=dim; 47 for (;it!=itend;++it){ 48 it->index=index_t(it->index.begin(),it->index.begin()+dim); 49 } 50 return; 51 } 52 int delta_dim=dim-p.dim; 53 p.dim=dim; 54 for (;it!=itend;++it) 55 add_dim(*it,delta_dim); 56 } 57 58 // returns q such that p=q [degree] and q has only terms of degree<degree 59 // p=q[N] means that p-q vanishes at v at order N reduce(const polynome & p,const vecteur & v,int degree)60 static polynome reduce(const polynome & p,const vecteur & v,int degree){ 61 int vsize=int(v.size()); 62 if (!vsize) 63 return p; 64 if (v==vecteur(vsize)){ 65 // trivial reduction, remove all terms of total deg >= degree 66 polynome res(p.dim); 67 vector< monomial<gen> >::const_iterator it=p.coord.begin(),itend=p.coord.end(); 68 for (;it!=itend;++it){ 69 if (total_degree(it->index)<degree) 70 res.coord.push_back(*it); 71 } 72 return res; 73 } 74 if (degree<=1){ 75 gen res=peval(p,v,0); 76 if (is_zero(res)) 77 return polynome(p.dim); 78 if (res.type==_POLY){ 79 polynome resp(*res._POLYptr); 80 change_dim(resp,p.dim); 81 return resp; 82 } 83 else 84 return polynome(monomial<gen>(res,0,p.dim)); 85 } 86 polynome pcur(p); 87 polynome y(monomial<gen>(plus_one,1,1,p.dim)); 88 if (!is_zero(v.front())) 89 y.coord.push_back(monomial<gen>(-v.front(),0,1,p.dim)); 90 polynome quo(y.dim),rem(y.dim); 91 pcur.TDivRem1(y,quo,rem); 92 rem=reduce(rem.trunc1(),vecteur(v.begin()+1,v.end()),degree); 93 quo=reduce(quo,v,degree-1); 94 return quo*y+rem.untrunc1(); 95 } 96 reduce_poly(const polynome & p,const vecteur & v,int degree,polynome & res)97 static void reduce_poly(const polynome & p,const vecteur & v,int degree,polynome & res){ 98 res.coord.clear(); 99 res.dim=p.dim; 100 vector< monomial<gen> >::const_iterator it=p.coord.begin(),itend=p.coord.end(); 101 if (is_zero(v)){ 102 index_t::const_iterator jt,jtend; 103 int otherdeg; 104 for (;it!=itend;++it){ 105 jt=it->index.begin()+1; 106 jtend=it->index.end(); 107 for (otherdeg=0;jt!=jtend;++jt){ 108 otherdeg += *jt; 109 } 110 if (otherdeg<degree) 111 res.coord.push_back(*it); 112 } 113 } 114 else { 115 for (;it!=itend;){ 116 int d=it->index.front(); 117 polynome tmp(Tnextcoeff<gen>(it,itend)); 118 res=res+reduce(tmp,v,degree).untrunc1(d); 119 } 120 } 121 } 122 123 // Same as reduce but do it for every coefficient of p with 124 // respect to the main variable reduce_poly(const polynome & p,const vecteur & v,int degree)125 polynome reduce_poly(const polynome & p,const vecteur & v,int degree){ 126 polynome res(p.dim); 127 reduce_poly(p,v,degree,res); 128 return res; 129 } 130 131 // reduce_divrem does a mixed division: euclidean w.r.t. the first var 132 // and ascending power of X-v for the other vars 133 // FIXME: this implementation does not work currently, except if other 134 // depends only on the first var reduce_divrem2(const polynome & a,const polynome & other,const vecteur & v,int n,polynome & quo,polynome & rem,bool allowrational=false)135 static bool reduce_divrem2(const polynome & a,const polynome & other,const vecteur & v,int n,polynome & quo,polynome & rem,bool allowrational=false) { 136 int asize=int(a.coord.size()); 137 if (!asize){ 138 quo=a; 139 rem=a; 140 return true; 141 } 142 int bsize=int(other.coord.size()); 143 if (bsize==0) { 144 #ifdef NO_STDEXCEPT 145 return false; 146 #else 147 setsizeerr(gettext("ezgcd.cc/reduce_divrem2")); 148 #endif 149 } 150 index_m a_max = a.coord.front().index; 151 index_m b_max = other.coord.front().index; 152 quo.coord.clear(); 153 quo.dim=a.dim; 154 rem.dim=a.dim; 155 if ( (bsize==1) && (b_max==b_max*0) ){ 156 rem.coord.clear(); 157 gen b=other.coord.front().value; 158 if (is_one(b)) 159 quo = a ; 160 else { 161 std::vector< monomial<gen> >::const_iterator itend=a.coord.end(); 162 for (std::vector< monomial<gen> >::const_iterator it=a.coord.begin();it!=itend;++it) 163 quo.coord.push_back(monomial<gen>(rdiv(it->value,b,context0),it->index)); 164 } 165 return true; 166 } 167 rem=a; 168 if ( ! (a_max>=b_max) ){ 169 // test that the first power of a_max is < to that of b_max 170 return (a_max.front()<b_max.front()); 171 } 172 gen b(other.coord.front().value); 173 while (a_max >= b_max){ 174 // errors should be trapped here and false returned if error occured 175 gen q(rdiv(rem.coord.front().value,b,context0)); 176 if (!allowrational){ 177 if ( has_denominator(q) || 178 (!is_zero(q*b - rem.coord.front().value)) ) 179 return false; 180 } 181 // end error trapping 182 quo.coord.push_back(monomial<gen>(q,a_max-b_max)); 183 tensor<gen> temp; 184 reduce_poly(other.shift(a_max-b_max,q),v,n,temp); 185 rem = rem-temp; 186 if (rem.coord.size()) 187 a_max=rem.coord.front().index; 188 else 189 break; 190 } 191 return(true); 192 } 193 reduce_divrem(const polynome & a,const polynome & other,const vecteur & v,int n,polynome & quo,polynome & rem)194 bool reduce_divrem(const polynome & a,const polynome & other,const vecteur & v,int n,polynome & quo,polynome & rem) { 195 quo.coord.clear(); 196 quo.dim=a.dim; 197 rem.dim=a.dim; 198 // if ( (a.dim<=1) || (a.coord.empty()) ) 199 return reduce_divrem2(a,other,v,n,quo,rem); 200 #if 0 201 std::vector< monomial<gen> >::const_iterator it=other.coord.begin(); 202 int bdeg=it->index.front(),rdeg; 203 tensor<gen> b0(Tnextcoeff<gen>(it,other.coord.end())); 204 tensor<gen> r(a),q(b0.dim); 205 while ( (rdeg=r.lexsorted_degree()) >=bdeg){ 206 it=r.coord.begin(); 207 tensor<gen> a0(Tnextcoeff<gen>(it,r.coord.end())),tmp(a0.dim); 208 // FIXME: should make ascending power division 209 if (!reduce_divrem(a0,b0,v,n,q,tmp) || !tmp.coord.empty()) 210 return false; 211 q=q.untrunc1(rdeg-bdeg); 212 quo=quo+q; 213 r=r-reduce_poly(q*other,v,n); 214 if (r.coord.empty()) 215 return true; 216 } 217 return true; 218 #endif 219 } 220 221 // increment last index in v up to k, 222 // if last index is k-1 223 // while index[size-pos] is k-pos increment pos 224 // if pos reaches size return false (not possible anymore) 225 // else increment index[size-pos] and set following ones to prev+1 next(vector<int> & v,int dim,int k)226 static bool next(vector<int> & v,int dim,int k){ 227 ++v.back(); 228 if (v.back()!=k) 229 return true; 230 int pos=2; 231 for (;pos<=dim;++pos){ 232 if (v[dim-pos]!=k-pos) 233 break; 234 } 235 if (pos>dim) 236 return false; 237 ++v[dim-pos]; 238 for (--pos;pos>0;--pos){ 239 v[dim-pos]=v[dim-pos-1]+1; 240 } 241 return true; 242 } 243 244 // pcur(x1,...,xk,0,...,0) peval_xk_xn_zero(const polynome & pcur,int k,polynome & pcurx1x2)245 static void peval_xk_xn_zero(const polynome & pcur,int k,polynome & pcurx1x2){ 246 pcurx1x2.coord.clear(); 247 int dim=pcur.dim; 248 pcurx1x2.dim=dim; 249 vector< monomial<gen> >::const_iterator it=pcur.coord.begin(),itend=pcur.coord.end(); 250 for (;it!=itend;++it){ 251 int j=k; 252 index_t::const_iterator i = it->index.begin()+j; 253 for (;j<dim;++j,++i){ 254 if (*i) 255 break; 256 } 257 if (j==dim) 258 pcurx1x2.coord.push_back(*it); 259 } 260 } 261 262 // pcur(x1,...,xk,0,...,0) truncate_xk_xn(polynome & pcur,int k)263 static void truncate_xk_xn(polynome & pcur,int k){ 264 vector< monomial<gen> >::iterator it=pcur.coord.begin(),itend=pcur.coord.end(); 265 for (;it!=itend;++it){ 266 it->index=index_t(it->index.begin(),it->index.begin()+k); 267 } 268 pcur.dim=k; 269 } 270 untruncate_xk_xn(polynome & pcur,int dim)271 static void untruncate_xk_xn(polynome & pcur,int dim){ 272 vector< monomial<gen> >::iterator it=pcur.coord.begin(),itend=pcur.coord.end(); 273 for (;it!=itend;++it){ 274 index_t i (dim); 275 i=it->index.iref(); 276 for (int j=int(i.size());j<dim;++j) 277 i.push_back(0); 278 it->index = i; 279 } 280 pcur.dim=dim; 281 } 282 283 gen _coeff(const gen &,GIAC_CONTEXT); 284 try_sparse_factor(const polynome & pcur,const factorization & v,int mult,factorization & f)285 bool try_sparse_factor(const polynome & pcur,const factorization & v,int mult,factorization & f){ 286 /* Try sparse factorization 287 lcoeff(pcur,x1)^#factors-1 * pcur = product_#factors P_i 288 where P_i has lcoeff(pcur,x1) as leading coeff in x1 289 and same non zeros coeffs pattern as the factors of Fb 290 */ 291 // count number of unknowns 292 factorization::const_iterator vit=v.begin(),vitend=v.end(); 293 int unknowns=0; 294 for (;vit!=vitend;++vit){ 295 if (vit->mult>1) 296 break; // return false might be more appropriate here 297 unknowns += int(vit->fact.coord.size())-1; // lcoeff is known 298 } 299 if (unknowns>=giacmax(5,pcur.lexsorted_degree()/2) || unknowns==0) 300 return false; 301 polynome lcp(Tfirstcoeff(pcur)); 302 int dim=pcur.dim; 303 vecteur lv(dim); 304 for (int i=0;i<dim;++i){ 305 lv[i]=identificateur("x"+print_INT_(i)); 306 } 307 gen mainvar(lv.front()); 308 gen lc=r2sym(lcp,lv,context0); 309 vecteur la(unknowns); 310 for (int i=0;i<unknowns;++i){ 311 la[i]=identificateur("a"+print_INT_(i)); 312 } 313 vecteur la_val(la); 314 int pos=0; 315 // build product(P_i) 316 gen product(1); 317 vecteur Pis; 318 for (vit=v.begin();vit!=vitend;++vit){ 319 const polynome & fact = vit->fact; 320 vector< monomial<gen> >::const_iterator it=fact.coord.begin(),itend=fact.coord.end(); 321 gen Pi=lc*pow(mainvar,it->index.front()); 322 for (++it;it!=itend;++it){ 323 if (pos>=la.size()) 324 return false; 325 Pi += la[pos]*pow(mainvar,it->index.front()); 326 ++pos; 327 } 328 Pis.push_back(Pi); 329 product = product * Pi; 330 } 331 product=product-r2sym(pcur,lv,context0)*pow(lc,int(Pis.size())-1,context0); 332 // solve equation wrt la 333 gen systemeg=_coeff(gen(makevecteur(product,mainvar),_SEQ__VECT),context0); 334 if (systemeg.type!=_VECT) 335 return false; 336 vecteur syst; 337 const_iterateur it=systemeg._VECTptr->begin(),itend=systemeg._VECTptr->end(); 338 for (++it;it!=itend;++it){ 339 if (!is_zero(*it)) 340 syst.push_back(*it); 341 } 342 // to solve syst wrt la, we search all linear equations 343 // if none return false, otherwise solve system, subst 344 while (!syst.empty()){ 345 int N=int(syst.size()); 346 vecteur linear; 347 for (int i=0;i<N;++i){ 348 gen d1=derive(syst[i],la,context0); 349 if (is_zero(derive(d1,la,context0))) 350 linear.push_back(syst[i]); 351 } 352 if (linear.empty()) 353 return false; 354 vecteur indet(lv); 355 lvar(linear,indet); 356 indet=vecteur(indet.begin()+lv.size(),indet.end()); 357 vecteur sols=linsolve(linear,indet,context0); 358 if (sols.size()!=indet.size() || is_undef(sols) || sols.empty()) 359 return false; 360 la_val=subst(la_val,indet,sols,false,context0); 361 gen tmp=recursive_normal(subst(syst,indet,sols,false,context0),context0); 362 if (tmp.type!=_VECT) 363 return false; 364 syst.clear(); 365 const_iterateur it=tmp._VECTptr->begin(),itend=tmp._VECTptr->end(); 366 for (;it!=itend;++it){ 367 if (!is_zero(*it)) 368 syst.push_back(*it); 369 } 370 } 371 // subst la values 372 Pis=subst(Pis,la,la_val,false,context0); 373 for (unsigned int i=0;i<Pis.size();++i){ 374 gen tmp=sym2r(Pis[i],lv,context0),num,den; 375 fxnd(tmp,num,den); 376 if (num.type!=_POLY) 377 return false; 378 const polynome & N=*num._POLYptr; 379 f.push_back(facteur<polynome>(N/lgcd(N),mult)); 380 } 381 return true; 382 } 383 384 // pcur(x,x1,x2,...) with [x1,x2,...]=[t^n1,t^n2,...] eval_tn(const polynome & pcur,const index_t & n,polynome & pt)385 void eval_tn(const polynome & pcur,const index_t & n,polynome & pt){ 386 pt.dim=2; 387 pt.coord.clear(); 388 pt.coord.reserve(pcur.coord.size()); 389 vector< monomial<gen> >::const_iterator it=pcur.coord.begin(),itend=pcur.coord.end(); 390 index_t cur(2); 391 for (;it!=itend;++it){ 392 const index_t & i=it->index.iref(); 393 index_t::const_iterator jt=i.begin(),jtend=i.end(); 394 index_t::const_iterator nt=n.begin(); 395 cur[0]=*jt; 396 int curn=0; 397 for (++jt;jt!=jtend;++jt,++nt) 398 curn += (*jt)*(*nt); 399 cur[1]=curn; 400 pt.coord.push_back(monomial<gen>(it->value,cur)); 401 } 402 pt.tsort(); 403 } 404 405 // return true if none of the coefficients of p with same 1st degree are the same x_degrees(const polynome & p,vector<int> & d)406 bool x_degrees(const polynome & p,vector<int> & d){ 407 d.clear(); 408 vector< monomial<gen> >::const_iterator it=p.coord.begin(),itend=p.coord.end(); 409 int prev=-1; 410 vecteur v; 411 for (;it!=itend;++it){ 412 int cur=it->index.iref().front(); 413 if (cur!=prev){ 414 v=vecteur(1,it->value); 415 d.push_back(cur); 416 prev=cur; 417 } 418 else { 419 if (equalposcomp(v,it->value)) 420 return false; 421 v.push_back(it->value); 422 } 423 } 424 return true; 425 } 426 lex_or_coeff_sort(const monomial<gen> & a,const monomial<gen> & b)427 bool lex_or_coeff_sort(const monomial<gen> & a,const monomial<gen> & b){ 428 if (a.index.front()!=b.index.front()) 429 return a.index.front()>b.index.front(); 430 return is_strictly_greater(a.value,b.value,context0); 431 } 432 try_sparse_factor_bi(polynome & pcur,int mult,factorization & f)433 bool try_sparse_factor_bi(polynome & pcur,int mult,factorization & f){ 434 int dim=pcur.dim; 435 if (dim<=2) 436 return false; 437 /* Try sparse factorization using bivariate images of a factor of 438 pcur(x,x1,x2,...) with [x1,x2,...]=[t^n1,t^n2,...] 439 where n1,n2,...=1,1,... then 2,1,... then 1,2,... 440 */ 441 polynome lcp(Tfirstcoeff(pcur)),lcpt; 442 polynome pt,ptcont; 443 index_t n(dim-1,1); 444 for (;;){ 445 eval_tn(pcur,n,pt); 446 pt=pt/Tlgcd(pt); 447 eval_tn(lcp,n,lcpt); 448 #if POLY_SPARSE_BI 449 factorization ft; 450 gen extra_div_t; 451 factor(pt,ptcont,ft,false,false,false,1,extra_div_t); 452 if (ft.size()==1){ 453 f.push_back(facteur<polynome>(pcur,mult)); 454 return true; 455 } 456 factorization::const_iterator vit=ft.begin(),vitend=ft.end(); 457 #else 458 vecteur lv(makevecteur(vx_var,gen("t",context0))); 459 gen dbg=_poly2symb(makesequence(pt,lv),context0); 460 dbg=_factors(dbg,context0) ; 461 if (dbg.type!=_VECT) return false; 462 vecteur v=*dbg._VECTptr; 463 if (v.size()==2){ 464 f.push_back(facteur<polynome>(pcur,mult)); 465 return true; 466 } 467 iterateur vit=v.begin(),vitend=v.end(); 468 #endif 469 // factor must be distinct from other factors 470 // by one of the degrees in x 471 // select which factor will be reconstructed: 472 // multby=lcpt/lcoeff(factor of ft) must be as simple as possible 473 // Once selected, the factor will be normalized by * by multby 474 vector<int> seldegs; 475 polynome multby,selp; 476 for (;vit!=vitend;++vit){ 477 #if POLY_SPARSE_BI 478 if (vit->mult>1) break; 479 const polynome & p=vit->fact; 480 #else 481 ++vit; 482 if (*vit!=1) break; 483 gen pg=_symb2poly(makesequence(*(vit-1),lv),context0); 484 if (pg.type!=_POLY) break; 485 const polynome & p = *pg._POLYptr; 486 #endif 487 index_t D=p.degree(); 488 double ratio=p.coord.size()/(double(D[0])*D[1]); 489 if (ratio>0.2) 490 return false; 491 vector<int> degs; 492 bool b=x_degrees(p,degs); 493 if (degs==seldegs) break; 494 polynome multbynew=lcpt/Tfirstcoeff(p); 495 if (seldegs.empty() || (b && multbynew.coord.size()<multby.coord.size())){ 496 if (!b){ 497 // some coeffs are the same, dilate randomly 498 // using -1, 1, 2, -2 499 vecteur lv(dim); 500 for (int i=0;i<dim;++i){ 501 lv[i]=identificateur("x"+print_INT_(i)); 502 } 503 gen pcurg=_poly2symb(makesequence(pcur,lv),context0); 504 vecteur lw(lv); 505 vecteur dilate=vranm(dim,4,context0); 506 for (int k=1;k<dim;++k){ 507 int c=dilate[k].val; 508 switch (c){ 509 case 0: 510 dilate[k]=2; 511 break; 512 case 1: case 2: 513 dilate[k]=-1; 514 break; 515 case 3: 516 dilate[k]=2; 517 break; 518 } 519 } 520 for (int k=1;k<dim;++k) 521 lw[k]=dilate[k]*lv[k]; 522 pcurg=subst(pcurg,lv,lw,false,context0); 523 pcurg=_symb2poly(makesequence(pcurg,lv),context0); 524 if (pcurg.type!=_POLY) 525 return false; 526 polynome pcur_dilated=*pcurg._POLYptr; 527 factorization f_dilated; 528 if (!try_sparse_factor_bi(pcur_dilated,mult,f_dilated)) 529 return false; 530 factorization::const_iterator fit=f_dilated.begin(),fitend=f_dilated.end(); 531 for (;fit!=fitend;++fit){ 532 pcurg=_poly2symb(makesequence(fit->fact,lv),context0); 533 for (int k=1;k<dim;++k) 534 lw[k]=lv[k]/dilate[k]; 535 pcurg=subst(pcurg,lv,lw,false,context0); 536 pcurg=_symb2poly(makesequence(pcurg,lv),context0); 537 if (pcurg.type!=_POLY) 538 return false; 539 f.push_back(facteur<polynome>(*pcurg._POLYptr,fit->mult)); 540 } 541 return true; 542 } 543 seldegs=degs; 544 multby=multbynew; 545 selp=multby*p; 546 } 547 } 548 if (vit!=vitend){ 549 ++n[0]; 550 if (n[0]>=4) 551 return false; 552 continue; 553 } 554 // we will deduce x1^ in monomials by comparing with the same factor 555 // of the bivariate factorization with n1=2 instead of n1=1 556 // then x2^ with n1=1 and n2=2 557 // If one bivariate image has less monomials than another one it is an unlucky n, use another one 558 // If one bivariate image has more monomials, then we must throw everything and restart with this bivariate image 559 // Once all monomials are done we should get a factor of pcur 560 // by extracting the primitive part of this factor 561 sort(selp.coord.begin(),selp.coord.end(),lex_or_coeff_sort); 562 polynome curp,recon(selp); recon.dim=pcur.dim; 563 int increment=1,i=0; 564 for (;i<n.size();){ 565 index_t n1(n); 566 n1[i] += increment; 567 int ni=n[i],n1i=n1[i]; 568 eval_tn(pcur,n1,pt); 569 pt=pt/Tlgcd(pt); 570 #if POLY_SPARSE_BI 571 factor(pt,ptcont,ft,false,false,false,1,extra_div_t); 572 vit=ft.begin();vitend=ft.end(); 573 #else 574 dbg=_poly2symb(makesequence(pt,lv),context0); 575 dbg=_factors(dbg,context0) ; 576 if (dbg.type!=_VECT) return false; 577 v=*dbg._VECTptr; 578 iterateur vit=v.begin(),vitend=v.end(); 579 #endif 580 // lcoeff normalization 581 eval_tn(lcp,n1,lcpt); 582 // serch in factorization for seldeg x-degree pattern 583 curp.coord.clear(); 584 for (;vit!=vitend;++vit){ 585 #if POLY_SPARSE_BU 586 if (vit->mult>1){vit=vitend;} break; 587 const polynome & p=vit->fact; 588 #else 589 ++vit; 590 if (*vit!=1) break; 591 gen pg=_symb2poly(makesequence(*(vit-1),lv),context0); 592 if (pg.type!=_POLY) break; 593 const polynome & p = *pg._POLYptr; 594 #endif 595 vector<int> degs; 596 if (!x_degrees(p,degs)) break; 597 if (degs==seldegs){ 598 curp=lcpt/Tfirstcoeff(p)*p; 599 break; 600 } 601 } 602 if (vit==vitend || curp.coord.empty()) break; // not found or not sqrfree 603 // compare with selp 604 if (curp.coord.size()<selp.coord.size()){ // unlucky 605 ++increment; 606 if (increment>3) 607 break; 608 continue; 609 } 610 sort(curp.coord.begin(),curp.coord.end(),lex_or_coeff_sort); 611 if (curp.coord.size()>selp.coord.size()){ 612 // selp was unlucky, restart 613 recon=selp=curp; 614 n=n1; 615 break; 616 } 617 // selp and curp size match, now compare monomial by monomial 618 // and extract x[i] exponent in recon 619 vector< monomial<gen> >::iterator rt=recon.coord.begin(),rtend=recon.coord.end(),st=selp.coord.begin(),ct=curp.coord.begin(); 620 for (;rt!=rtend;++rt,++st,++ct){ 621 if (st->index[0]!=ct->index[0]) 622 break; 623 int idx0=st->index[1]; 624 int idx1=ct->index[1]; 625 index_t I=rt->index.iref(); 626 int delta=(idx1-idx0)/(n1i-ni); 627 if (i==0) 628 I[1]=delta; 629 else 630 I.push_back(delta); 631 if (i==n.size()-2){ 632 for (int j=0;j<=i;++j){ 633 idx1 -= I[j+1]*n1[j]; 634 } 635 I.push_back(idx1/n1[i+1]); 636 } 637 rt->index=I; 638 } 639 if (rt!=rtend) 640 break; 641 increment=1; 642 if (i==n.size()-2) ++i; 643 ++i; 644 } 645 if (i<n.size()){ 646 // restart search 647 ++n[i]; 648 if (n[i]>=4) 649 return false; 650 continue; 651 } 652 recon.tsort(); 653 // divide by reconstructed factor and restart factorization 654 recon=recon/Tlgcd(recon); 655 polynome quo,rem; 656 if (!pcur.TDivRem(recon,quo,rem,false) || !is_zero(rem)) 657 return false; 658 f.push_back(facteur<polynome>(recon,mult)); 659 pcur=quo; 660 return try_sparse_factor_bi(pcur,mult,f); 661 } // end endless for 662 } 663 poly_truncate(const polynome & q,polynome & q1,int j)664 void poly_truncate(const polynome & q,polynome & q1,int j){ 665 q1.coord.clear(); 666 vector< monomial<gen> >::const_iterator jt=q.coord.begin(),jtend=q.coord.end(); 667 for (;jt!=jtend;++jt){ 668 if (jt->index.total_degree()<j) 669 q1.coord.push_back(*jt); 670 } 671 } 672 673 // multiply keep only if total degree < maxdeg mulpoly_truncate(const polynome & p,const polynome & q,polynome & res,int maxdeg)674 void mulpoly_truncate(const polynome & p,const polynome & q,polynome &res,int maxdeg){ 675 res.coord.clear(); 676 int dim=p.dim; 677 polynome p1(dim),q1(dim),tmp(dim); 678 for (int i=0;i<maxdeg;++i){ 679 // p1 total degree i of p, q1 total degree<maxdeg-i of q 680 int j=maxdeg-i; 681 // create p1 and q1 682 p1.coord.clear(); 683 vector< monomial<gen> >::const_iterator it=p.coord.begin(),itend=p.coord.end(); 684 for (;it!=itend;++it){ 685 if (it->index.total_degree()==i) 686 p1.coord.push_back(*it); 687 } 688 poly_truncate(q,q1,j); 689 // multiply, 690 mulpoly(p1,q1,tmp,0); 691 // add to res 692 p1.coord.clear(); 693 tmp.TAdd(res,p1); 694 p1.coord.swap(res.coord); 695 } 696 } 697 698 // keep only monomials of total_degree==j without first degree poly_truncate1(const polynome & q,polynome & q1,int j)699 void poly_truncate1(const polynome & q,polynome & q1,int j){ 700 q1.coord.clear(); 701 vector< monomial<gen> >::const_iterator it=q.coord.begin(),itend=q.coord.end(); 702 index_t::const_iterator jt,jtend; 703 for (;it!=itend;++it){ 704 jt=it->index.begin()+1; 705 jtend=it->index.end(); 706 int otherdeg; 707 for (otherdeg=*jt,++jt;jt!=jtend;++jt){ 708 otherdeg += *jt; 709 } 710 if (otherdeg==j) 711 q1.coord.push_back(*it); 712 } 713 } 714 other_deg(const polynome & p,vector<int> & pdeg)715 void other_deg(const polynome & p,vector<int> & pdeg){ 716 pdeg.reserve(p.coord.size()); pdeg.clear(); 717 vector< monomial<gen> >::const_iterator it=p.coord.begin(),itend=p.coord.end(); 718 for (;it!=itend;++it){ 719 index_t::const_iterator jt,jtend; 720 jt=it->index.begin()+1; 721 //jtend=jt+dim-1; 722 jtend=it->index.end(); 723 int otherdeg; 724 for (otherdeg=*jt,++jt;jt<jtend;++jt){ 725 otherdeg += *jt; 726 } 727 pdeg.push_back(otherdeg); 728 } 729 } 730 731 // multiply keep only if total degree excluding 1st deg == maxdeg mulpoly_truncate1(const polynome & p,const polynome & q,polynome & res,int deg,polynome & p1,polynome & q1,polynome & tmp,vector<int> & pdeg,vector<int> & qdeg)732 void mulpoly_truncate1(const polynome & p,const polynome & q,polynome &res,int deg,polynome & p1,polynome & q1,polynome & tmp,vector<int> & pdeg,vector<int> & qdeg){ 733 bool eq=deg>=0; 734 int maxdeg=eq?deg:-deg; 735 res.coord.clear(); 736 int dim=p.dim; 737 int ps=int(p.coord.size()),qs=int(q.coord.size()); 738 p1.coord.reserve(ps); 739 other_deg(p,pdeg); 740 other_deg(q,qdeg); 741 const vector< monomial<gen> > & pcoord=p.coord; 742 const vector< monomial<gen> > & qcoord=q.coord; 743 for (int i=0;i<=maxdeg;++i){ 744 // p1 total degree <=i of p or ==i if deg>0, 745 // q1 total degree==maxdeg-i of q 746 int j=maxdeg-i; 747 // create p1 and q1 748 p1.coord.clear(); 749 for (int k=0;k<ps;++k){ 750 int otherdeg=pdeg[k]; 751 if (eq?otherdeg==i:otherdeg<=i) 752 p1.coord.push_back(pcoord[k]); 753 } 754 q1.coord.clear(); 755 for (int k=0;k<qs;++k){ 756 if (qdeg[k]==j) 757 q1.coord.push_back(qcoord[k]); 758 } 759 // multiply, 760 mulpoly(p1,q1,tmp,0); 761 // add to res 762 p1.coord.clear(); 763 tmp.TAdd(res,p1); 764 p1.coord.swap(res.coord); 765 } 766 } 767 768 try_hensel_lift_factor(const polynome & pcur,const polynome & F0,const factorization & v0,int mult,factorization & f)769 bool try_hensel_lift_factor(const polynome & pcur,const polynome & F0,const factorization & v0,int mult,factorization & f){ 770 int dim=pcur.dim; 771 int s=int(v0.size()); 772 polynome lcp(Tfirstcoeff(pcur)); 773 if (lcp.coord.back().index.back()!=0) 774 return false; 775 gen lcpb=lcp.coord.back().value; 776 vector<polynome> lcoeffs(s,lcp); 777 bool lcoeff_known=false; 778 factorization::const_iterator F0it=v0.begin(),F0itend=v0.end(); 779 vector<modpoly> F0fact; 780 for (;F0it!=F0itend;++F0it){ 781 if (F0it->mult>1) 782 break; 783 F0fact.push_back(modularize(F0it->fact,0,0)); 784 if (is_undef(F0fact.back())) 785 return false; 786 } 787 if (pcur.dim>2 && lcp.coord.size()>1){ 788 // try bivariate factorization to compute a priori the leadings coefficients 789 // that is factor pcur(x1,x2,0,...,0) = product p_i(x1,x2) 790 // then factor lcoeff(pcur)(x2,x3,..,xn) = product q_j(x2,..,xn) 791 // the lcoeff(pcur) corresponding to lcoeff(p_i)(x2) divides 792 // the product of the q_j such that either q_j(x2,0,...,0) is constant 793 // or gcd(q_j(x2,0,...,0),lcoeff(p_i)(x2)) is not constant 794 // then we know multiples of the lcoeffs, we can therefore replace s, F0it, F0itend, lcoeffs 795 polynome pcurx1x2; 796 peval_xk_xn_zero(pcur,2,pcurx1x2); 797 factorization fx1x2,flcoeff; 798 vector<polynome> flcoeff0; 799 polynome pcurx1x2cont=lgcd(pcurx1x2); 800 gen extra_div=1; 801 // if pcurx1x2cont is not 1, the following code may fail 802 // example f:=3/5*a^2*b^2*c^2+129/20*a^2*b*c^3+18/5*a^2*b*c^2*d-1443/40*a^2*b*c^2+387/10*a^2*c^3*d-387*a^2*c^3-9/20*a^2*c^2*d+9/2*a^2*c^2-273/5*a*b^2*c^2*d+3/5*a*b^2*c^2-11739/20*a*b*c^3*d+129/20*a*b*c^3-18/5*a*b*c^2*d^2+1857/40*a*b*c^2*d-1443/40*a*b*c^2-387/10*a*c^3*d^2+4257/10*a*c^3*d-387*a*c^3+9/20*a*c^2*d^2-99/20*a*c^2*d+9/2*a*c^2+54*b^2*c^2*d^2-54*b^2*c^2*d+1161/2*b*c^3*d^2-1161/2*b*c^3*d-27/4*b*c^2*d^2+27/4*b*c^2*d; factor(f); 803 if ( 804 // is_one(pcurx1x2cont) 805 is_zero(pcurx1x2cont.lexsorted_degree()) 806 ){ 807 truncate_xk_xn(pcurx1x2,2); 808 if (lgcd(pcurx1x2).coord.size()>1) 809 return false; 810 if (!factor(pcurx1x2,pcurx1x2cont,fx1x2,true,false,false,1,extra_div) || extra_div!=1) 811 return false; 812 // fx1x2 contains factorization of pcur(x1,x2,0,...0) 813 // now find factorization of lcoeff(pcur)(x2,...,xn) 814 polynome pcur_lcoeff(Tfirstcoeff(pcur)),pcur_lcoeffcont,pcur_lcoeff_sqfftest; 815 peval_xk_xn_zero(pcur_lcoeff,2,pcur_lcoeff_sqfftest); 816 pcur_lcoeff_sqfftest=pcur_lcoeff_sqfftest.trunc1(); 817 // if (gcd(pcur_lcoeff_sqfftest,pcur_lcoeff_sqfftest.derivative()).lexsorted_degree()) return false; 818 if (!factor(pcur_lcoeff.trunc1(),pcur_lcoeffcont,flcoeff,false,false,false,1,extra_div) || extra_div!=1) 819 return false; 820 factorization::iterator jt=flcoeff.begin(),jtend=flcoeff.end(); 821 polynome constante(pcur_lcoeffcont.untrunc1()*pcurx1x2cont),tmp; 822 for (;jt!=jtend;++jt){ 823 jt->fact=jt->fact.untrunc1(); 824 peval_xk_xn_zero(jt->fact,2,tmp); // should only depend on x2 825 if (Tis_constant(tmp)) 826 constante=constante*pow(jt->fact,jt->mult); 827 //else 828 flcoeff0.push_back(tmp); 829 } // flcoeff0 contains the list of factors of lcoeff(pcur) evaled at 0 830 F0it=fx1x2.begin(); 831 F0itend=fx1x2.end(); 832 s=int(F0itend-F0it); 833 F0fact.clear(); 834 lcoeffs.clear(); 835 modpoly piF(1,1); 836 for (;F0it!=F0itend;++F0it){ 837 if (F0it->mult>1) 838 break; 839 polynome p (F0it->fact); // depends on x1 and x2 840 untruncate_xk_xn(p,dim); 841 peval_xk_xn_zero(p,1,tmp); // make x2=0 842 truncate_xk_xn(tmp,1); 843 modpoly Fi(modularize(tmp,0,0)); 844 if (is_undef(Fi)) 845 return false; 846 if (gcd(piF,Fi,0).size()>1) 847 return false; 848 piF=piF*Fi; 849 F0fact.push_back(Fi); 850 // corresponding lcoeff 851 p=Tfirstcoeff(p); 852 polynome tmp2=constante; 853 for (jt=flcoeff.begin(),jtend=flcoeff.end();jt!=jtend;++jt){ 854 for (int m=jt->mult;m>0;--m){ 855 polynome G(flcoeff0[jt-flcoeff.begin()]); 856 if (Tis_constant(simplify(p,G))) 857 break; 858 else { 859 tmp2 = tmp2 * jt->fact; 860 // mark jt->fact as used 861 --jt->mult; 862 } 863 } 864 } 865 lcoeffs.push_back(tmp2); 866 } 867 lcoeff_known=true; 868 } // if is_one(pcurx1x2cont) 869 } 870 if (F0it!=F0itend) 871 return false; 872 // ok each factor of F0=pcur|0 is square free, they are prime together 873 // if lcp has too much terms it will take too long, because 874 // we must multiply by product(lcoeffs)/lcp 875 // next check was >100 but then heuristic factorization fails 876 // (should also depends on the size of the coeffs and number of variables...) 877 if (!lcoeff_known && pow(lcp,s-1).coord.size()>1000) 878 return false; 879 // we will lift pcur*product(lcoeffs)/lcp = product_i F0fact[i]*lcoeffs[i](b)/lcoeff(F0fact[i]) 880 for (int i=0;i<s;++i){ 881 gen lcoeff=F0fact[i].front(); 882 mulmodpoly(F0fact[i],lcoeffs[i].coord.back().value/lcoeff,F0fact[i]); 883 } 884 polynome pcur_adjusted(pcur); 885 if (!is_one(lcp)){ 886 if (lcoeff_known){ 887 polynome temp(lcoeffs[0]); 888 for (int i=1;i<s;++i){ 889 temp = temp * lcoeffs[i]; 890 } 891 temp = temp / lcp; 892 pcur_adjusted = pcur_adjusted * temp; 893 } 894 else { 895 for (int i=1;i<s;++i){ 896 pcur_adjusted =pcur_adjusted*lcoeffs[i]; 897 } 898 } 899 } 900 // Perhaps the check on lcp should be made here with pcur_adjusted 901 vector<modpoly> u; 902 if (!egcd(F0fact,0,u)) 903 return false;// sum_j U_j * product_{i \neq j} F0fact_i = 1 904 // factor out common deno 905 // sum_j U_j * product_{i \neq j} F0fact_i = D 906 vecteur den(s); 907 gen D(1); 908 for (int i=0;i<s;++i){ 909 lcmdeno(u[i],den[i],context0); 910 D=lcm(D,den[i]); 911 } 912 for (int i=0;i<s;++i) 913 mulmodpoly(u[i],D/den[i],u[i]); 914 vector<polynome> P(s),P0(s),U(s); 915 vecteur b(pcur_adjusted.dim-1); 916 for (int i=0;i<s;++i){ 917 U[i].dim=pcur_adjusted.dim; 918 P[i].dim=pcur_adjusted.dim; 919 P0[i].dim=pcur_adjusted.dim; 920 modpoly::const_iterator it=F0fact[i].begin(),itend=F0fact[i].end(); 921 int deg=int(itend-it)-1; 922 P[i]=lcoeffs[i].trunc1().untrunc1(deg); 923 for (int n=0;it!=itend;++it,++n){ 924 if (!is_zero(*it)){ 925 if (n) 926 P[i].coord.push_back(monomial<gen>(*it,deg-n,1,pcur_adjusted.dim)); 927 P0[i].coord.push_back(monomial<gen>(*it,deg-n,1,pcur_adjusted.dim)); 928 } 929 } 930 U[i].dim=pcur_adjusted.dim; 931 it=u[i].begin(); itend=u[i].end(); 932 deg=int(itend-it)-1; 933 for (int n=0;it!=itend;++it,++n){ 934 if (!is_zero(*it)) 935 U[i].coord.push_back(monomial<gen>(*it,deg-n,1,pcur_adjusted.dim)); 936 } 937 // CERR << Tcontent(U[i]) << '\n'; 938 } 939 polynome quo(dim),rem(dim),tmp(dim); 940 // we have now pcur_adjusted = product P_i + O(total_degree>=1) 941 int Total=pcur_adjusted.total_degree(); 942 // lift to pcur_adjusted = product P_i + O(total_degree>=k+1) 943 // for deg from 1 to total_degree(pcur_adjusted) 944 // P_i += (pcur_adjusted-product P_i) * U_j mod total_degree(k+1) 945 #if 1 // def EZGCD_DEGONLY 946 if (is_zero(b)){ 947 polynome tmp4(dim),tmp5(dim),tmp6(dim),prod(dim); 948 vector<int> tmpi1,tmpi2; 949 for (int deg=1;deg<=Total;++deg){ 950 prod=P[s-2]; 951 for (int i=s-3;i>=0;--i){ 952 // reduce_poly(prod * P[i],b,deg+1,prod); // keep up to deg 953 tmp.coord.clear(); 954 mulpoly_truncate1(prod,P[i],tmp,-deg,tmp4,tmp5,tmp6,tmpi1,tmpi2); 955 prod.coord.swap(tmp.coord); 956 //if (prod!=prod1) CERR << "err " << deg << '\n'; 957 } // end loop on i 958 mulpoly_truncate1(prod,P[s-1],tmp,deg,tmp4,tmp5,tmp6,tmpi1,tmpi2); 959 prod.coord.swap(tmp.coord); 960 poly_truncate1(pcur_adjusted,tmp,deg); 961 prod = tmp - prod; 962 if (prod.coord.empty()){ 963 // check total degrees 964 int tdeg=0; 965 for (int i=0;i<s;++i) 966 tdeg += P[i].total_degree(); 967 if (tdeg==Total){ 968 if (deg!=Total){ 969 prod=P[s-1]; 970 for (int i=s-2;i>=0;--i){ 971 // prod = prod * P[i]; 972 tmp.coord.clear(); 973 mulpoly(prod,P[i],tmp,0); 974 prod.coord.swap(tmp.coord); 975 } 976 // N.B. prod==pcur_adjusted does not always work! 977 if ((prod-pcur_adjusted).coord.empty()) 978 deg=Total; 979 } 980 if (deg==Total){ 981 for (int i=0;i<s;++i){ 982 f.push_back(facteur<polynome>(P[i]/lgcd(P[i]),mult)); 983 } 984 return true; 985 } 986 } 987 continue; 988 } 989 //CERR << Tcontent(prod) << '\n'; 990 for (int i=0;i<s;++i){ 991 // U[i] depends only on 1st var no need to reduce 992 mulpoly(prod,U[i],rem,0); 993 //CERR << "deg " << deg << " " << Tcontent(rem) << '\n'; 994 if (!divrem1(rem,P0[i],quo,tmp,0) && !rem.TDivRem1(P0[i],quo,tmp,true,0)) 995 return false; 996 rem.coord.swap(tmp.coord); // poly_truncate1(tmp,rem,deg); 997 // divide by D 998 vector< monomial<gen> >::const_iterator r1=rem.coord.begin(),r2=rem.coord.end(); 999 Div<gen>(r1,r2,D,rem.coord); 1000 P[i] = P[i] + rem; 1001 } 1002 } 1003 } // end if (is_zero(b)) 1004 else 1005 #endif 1006 for (int deg=1;deg<=Total;++deg){ 1007 polynome prod(P[s-1]); 1008 for (int i=s-2;i>=0;--i){ 1009 // reduce_poly(prod * P[i],b,deg+1,prod); // keep up to deg 1010 tmp.coord.clear(); 1011 mulpoly(prod,P[i],tmp,0); 1012 reduce_poly(tmp,b,deg+1,prod); 1013 } 1014 prod = reduce_poly(pcur_adjusted,b,deg+1) - prod; 1015 if (prod.coord.empty()){ 1016 // check total degrees 1017 int tdeg=0; 1018 for (int i=0;i<s;++i) 1019 tdeg += P[i].total_degree(); 1020 if (tdeg==Total){ 1021 if (deg!=Total){ 1022 prod=P[0]; 1023 for (int i=1;i<s;++i){ 1024 // prod = prod * P[i]; 1025 tmp.coord.clear(); 1026 mulpoly(prod,P[i],tmp,0); 1027 swap(tmp,prod); 1028 } 1029 // N.B. prod==pcur_adjusted does not always work! 1030 if ((prod-pcur_adjusted).coord.empty()) 1031 deg=Total; 1032 } 1033 if (deg==Total){ 1034 for (int i=0;i<s;++i){ 1035 f.push_back(facteur<polynome>(P[i]/lgcd(P[i]),mult)); 1036 } 1037 return true; 1038 } 1039 } 1040 continue; 1041 } 1042 //CERR << Tcontent(prod) << '\n'; 1043 for (int i=0;i<s;++i){ 1044 // U[i] depends only on 1st var no need to reduce 1045 mulpoly(prod,U[i],rem,0); 1046 //CERR << "deg " << deg << " " << Tcontent(rem) << '\n'; 1047 if (!divrem1(rem,P0[i],quo,tmp,0) && !rem.TDivRem1(P0[i],quo,tmp,true,0)) 1048 return false; 1049 reduce_poly(tmp,b,deg+1,rem); 1050 // divide by D 1051 vector< monomial<gen> >::const_iterator r1=rem.coord.begin(),r2=rem.coord.end(); 1052 Div<gen>(r1,r2,D,rem.coord); 1053 P[i] = P[i] + rem; 1054 } 1055 } // end for 1056 // FIXME combine factors 1057 if (s==2){ 1058 f.push_back(facteur<polynome>(pcur,mult)); 1059 return true; 1060 } 1061 int nfact=s; 1062 index_t pcur_deg(pcur_adjusted.degree()); 1063 vector<int> test(1); 1064 for (int k=1;k<=nfact/2;){ 1065 if (debug_infolevel) 1066 COUT << CLOCK() << "Testing combination of " << k << " factors" << '\n'; 1067 // FIXME check on cst coeff 1068 if (1){ 1069 polynome prodP(P[test[0]]); 1070 for (int i=1;i<k;++i){ 1071 mulpoly(prodP,P[test[i]],rem,0); 1072 reduce_poly(rem,b,Total+1,prodP); 1073 } 1074 if (divrem1(pcur_adjusted,prodP,quo,rem,1) && rem.coord.empty()){ 1075 // factor found 1076 pcur_adjusted=quo; 1077 f.push_back(facteur<polynome>(prodP/lgcd(prodP),mult)); 1078 for (int i=k-1;i>=0;--i){ 1079 P.erase(P.begin()+test[i]); 1080 } 1081 nfact -= k; 1082 for (int i=0;i<k;++i) 1083 test[i]=i; 1084 continue; 1085 } 1086 } 1087 if (!next(test,k,nfact)){ 1088 ++k; 1089 test=vector<int>(k); 1090 for (int i=0;i<k;++i) 1091 test[i]=i; 1092 } 1093 } 1094 f.push_back(facteur<polynome>(pcur_adjusted/lgcd(pcur_adjusted),mult)); 1095 return true; 1096 } 1097 1098 // find u,v,d s.t. u*p+v*q=d by Hensel lift try_hensel_egcd(const polynome & p,const polynome & q,polynome & u,polynome & v,polynome & d)1099 bool try_hensel_egcd(const polynome & p,const polynome & q,polynome &u,polynome &v,polynome & d){ 1100 // check # of variables 1101 //if (p.dim<=1 || p.dim!=q.dim) 1102 return false; 1103 // check that 0 is a good evaluation point (same degree, gcd==1) 1104 vecteur b(1,0); 1105 polynome p0(1),q0(1); 1106 find_good_eval(p,q,p0,q0,b,(debug_infolevel>=2)); 1107 if (!is_zero(b)) 1108 return false; 1109 int pdeg=p.lexsorted_degree(),qdeg=q.lexsorted_degree(); 1110 if (p0.lexsorted_degree()!=pdeg || q0.lexsorted_degree()!=qdeg) 1111 return false; 1112 gen g=gcd(pdeg,qdeg); 1113 if (g.type==_POLY && g._POLYptr->lexsorted_degree()) 1114 return false; 1115 // Bezout at other variables==0 1116 polynome u0(1),v0(1),d0(1); 1117 egcd(p0,q0,u0,v0,d0); // d0 must be constant 1118 // now p*u0+q*v0-d0=O(1) where O(k) means of order >= k wrt other variables 1119 // p*uk+q*vk-d0=O(k) -> p*(uk+uk1)+q*(v+vk1)-d0=O(k+1) 1120 // with uk1 and vk1=O(k+1) 1121 // we have p0*uk1+q0*vk1=d0-p*uk-q*vk=yk 1122 // hence uk1=yk*u0/d0 % q0, vk1=yk*v0/d0 % p0 1123 // rational (Pade-like) reconstruction uk=fraction of polynomials 1124 // with max degree wrt other variables <=k/2 1125 // once both fractions corresp. to uk and vk stabilizes, check identity 1126 } 1127 1128 // Hensel linear or quadratic lift 1129 // FIXME Quadratic lift currently works only if lcp is constant 1130 // Lift the equality p(b)=qb*rb [where b is a vecteur like for peval 1131 // assumed to have p.dim-1 coordinates] to p=q*r mod (X-b)^deg 1132 // Assuming that lcoeff(q)=lcp, lcoeff(r)=lcp, lcoeff(p)=lcp^2 1133 // If you want to find factors of a poly P such that P(b)=Qb*Rb, 1134 // if lcp is the leading coeff of P 1135 // then p=P*lcp, qb=Qb*lcp(b)/lcoeff(Qb), rb=Rb*lcp(b)/lcoeff(Rb) hensel_lift(const polynome & p,const polynome & lcp,const polynome & qb,const polynome & rb,const vecteur & b,polynome & q,polynome & r,bool linear_lift,double maxop)1136 bool hensel_lift(const polynome & p, const polynome & lcp, const polynome & qb, const polynome & rb, const vecteur & b,polynome & q, polynome & r,bool linear_lift,double maxop){ 1137 if (maxop) 1138 linear_lift=true; // otherwise please adjust number of operations to do! 1139 double nop=0; 1140 int dim=p.dim; 1141 int deg=total_degree(p); 1142 if ( (qb.dim!=1) || (rb.dim!=1) || (dim==1) ){ 1143 #ifdef NO_STDEXCEPT 1144 return false; 1145 #else 1146 setsizeerr(gettext("Bad dimension for qb or rb or b or degrees")); 1147 #endif 1148 } 1149 polynome qu(1),ru(1),qbd(1); 1150 egcd(qb,rb,qu,ru,qbd); 1151 if (!Tis_constant(qbd)){ 1152 #ifdef NO_STDEXCEPT 1153 return false; 1154 #else 1155 setsizeerr(gettext("qb and rb not prime together!")); 1156 #endif 1157 } 1158 gen qrd(qbd.coord.front().value); 1159 // now we have qu*qb+ru*rb=qrd with 1-d polynomials 1160 change_dim(qu,dim); 1161 change_dim(ru,dim); 1162 // adjust dim & leading coeff of q and r by removing current leading coeff 1163 // and replace by lcp 1164 q=qb; 1165 r=rb; 1166 change_dim(q,dim); 1167 change_dim(r,dim); 1168 polynome q0(q),r0(r); 1169 index_t qshift(q.dim); 1170 qshift[0]=q.lexsorted_degree(); 1171 q=q+(lcp-Tfirstcoeff<gen>(q)).shift(qshift); 1172 qshift[0]=r.lexsorted_degree(); 1173 r=r+(lcp-Tfirstcoeff<gen>(r)).shift(qshift); 1174 polynome p_qr(dim); 1175 for (int n=1;;){ 1176 // qu*q+ru*r=qrd [n] (it's exact at the loop begin) 1177 // p=q*r [n] where [n] means of total valuation >= n 1178 // at the beginning n=1 1179 // enhanced at order 2*n by adding q',r' of valuation >=n 1180 // p-(q+q')*(r+r')=p-q*r - (r'q+q'r)-q'*r' 1181 // hence if we put r', q' such that p-q*r=(r'q+q'r) [2n] 1182 // we are done. Since p-q*r is of order [n], we get the solution 1183 // r'=qu*(p-qr)/qrd and q'=ru*(p-qr)/qrd 1184 if (debug_infolevel) 1185 CERR << "// Hensel " << n << " -> " << deg << '\n'; 1186 if (n>deg) 1187 return false; 1188 if (linear_lift) 1189 ++n; 1190 else 1191 n=2*n; 1192 if (maxop>0){ 1193 nop += double(q.coord.size())*r.coord.size(); 1194 if (debug_infolevel) 1195 CERR << "EZGCD " << nop << ":" << maxop << '\n'; 1196 if (nop>maxop) 1197 return false; 1198 } 1199 p_qr=reduce_poly(p-q*r,b,deg); 1200 if (is_zero(p_qr)) 1201 return true; 1202 if (n>deg) 1203 n=deg; 1204 p_qr=reduce_poly(p_qr,b,n); 1205 polynome qprime(reduce_poly(ru*p_qr,b,n)),qquo(qprime.dim),qrem(qprime.dim); 1206 polynome rprime(reduce_poly(qu*p_qr,b,n)),rquo(rprime.dim),rrem(qprime.dim); 1207 // reduction of qprime and rprime with respect to the main variable 1208 // we know that 1209 // (*) degree(p_qr) < degree(qr) 1210 // where degree is the degree wrt the main variable 1211 // since the leading coeffs of q and r are still adjusted 1212 // Then there is a unique solution to (*) with 1213 // degree(qprime)<degree(q), degree(rprime)<degree(r) 1214 if (linear_lift){ 1215 reduce_divrem(qprime,q0,b,n,qquo,qrem); 1216 reduce_divrem(rprime,r0,b,n,rquo,rrem); 1217 } 1218 else { 1219 reduce_divrem(qprime,q,b,n,qquo,qrem); 1220 reduce_divrem(rprime,r,b,n,rquo,rrem); 1221 } 1222 // reduction of qprime and rprime with respect to the other variables 1223 // maybe we should check that q and r below have integer coeff 1224 q=q+inv(qrd,context0)*qrem; 1225 r=r+inv(qrd,context0)*rrem; 1226 if (!linear_lift && (n<=deg/2)){ 1227 // Now we modify qu and ru so that 1228 // (qu+qu')*q+(ru+ru')*r=qrd [2n] 1229 // therefore qu'*q+ru'*r=qrd-(qu*q+ru*r) [2n] 1230 // hence qu'= qu*[qrd-(qu*q+ru*r)]/qrd, ru'=ru*[qrd-(qu*q+ru*r)]/qrd 1231 p_qr=polynome(monomial<gen>(qrd,0,dim))-reduce_poly(qu*q+ru*r,b,n); 1232 qprime=reduce_poly(qu*p_qr,b,n); 1233 rprime=reduce_poly(ru*p_qr,b,n); 1234 reduce_divrem(qprime,r,b,n,qquo,qrem); 1235 reduce_divrem(rprime,q,b,n,rquo,rrem); 1236 qu=qu+inv(qrd,context0)*qrem; // should check that qu and ru have integer coeff 1237 ru=ru+inv(qrd,context0)*rrem; 1238 } 1239 } 1240 } 1241 1242 // Replace the last coordinates of p with b instead of the first peval_back(const polynome & p,const vecteur & b)1243 gen peval_back(const polynome & p,const vecteur & b){ 1244 int pdim=p.dim,bdim=int(b.size()); 1245 vector<int> cycle(pdim); 1246 int deltad=pdim-bdim; 1247 for (int i=0;i<bdim;++i) 1248 cycle[i]=i+deltad; 1249 for (int i=bdim;i<pdim;++i) 1250 cycle[i]=i-bdim; 1251 polynome pp(p); 1252 pp.reorder(cycle); 1253 int save=debug_infolevel; 1254 if (debug_infolevel) 1255 --debug_infolevel; 1256 gen res(peval(pp,b,0)); 1257 debug_infolevel=save; 1258 return res; 1259 } 1260 peval_1(const polynome & p,const vecteur & v,const gen & mod)1261 polynome peval_1(const polynome & p,const vecteur &v,const gen & mod){ 1262 #if defined(NO_STDEXCEPT) && !defined(RTOS_THREADX) && !defined(VISUALC) && !defined(KHICAS) 1263 assert(p.dim==signed(v.size()+1)); 1264 #else 1265 if (p.dim!=signed(v.size()+1)) 1266 setsizeerr(gettext("peval_1")); 1267 #endif 1268 polynome res(1); 1269 index_t i(1); 1270 std::vector< monomial<gen> >::const_iterator it=p.coord.begin(); 1271 std::vector< monomial<gen> >::const_iterator itend=p.coord.end(); 1272 for (;it!=itend;){ 1273 i[0]=it->index.front(); 1274 polynome pactuel(Tnextcoeff<gen>(it,itend)); 1275 gen g(peval(pactuel,v,mod)); 1276 if ( (g.type==_POLY) && (g._POLYptr->dim==0) ) 1277 g=g._POLYptr->coord.empty()?0:g._POLYptr->coord.front().value; 1278 if (!is_zero(g)) 1279 res.coord.push_back(monomial<gen>(g,i)); 1280 } 1281 return res; 1282 } 1283 unmodularize(const vector<int> & a)1284 polynome unmodularize(const vector<int> & a){ 1285 if (a.empty()) 1286 return polynome(1); 1287 polynome res(1); 1288 vector< monomial<gen> > & v=res.coord; 1289 index_t i; 1290 int deg=int(a.size())-1; 1291 i.push_back(deg); 1292 vector<int>::const_iterator it=a.begin(); 1293 vector<int>::const_iterator itend=a.end(); 1294 for (;it!=itend;++it,--i[0]){ 1295 if (*it) 1296 v.push_back(monomial<gen>(*it,i)); 1297 } 1298 return res; 1299 } 1300 convert_from_truncate(const vector<T_unsigned<int,hashgcd_U>> & p,hashgcd_U var,polynome & P)1301 static bool convert_from_truncate(const vector< T_unsigned<int,hashgcd_U> > & p,hashgcd_U var,polynome & P){ 1302 P.dim=1; 1303 P.coord.clear(); 1304 vector< T_unsigned<int,hashgcd_U> >::const_iterator it=p.begin(),itend=p.end(); 1305 P.coord.reserve(itend-it); 1306 index_t i(1); 1307 for (;it!=itend;++it){ 1308 i.front()=it->u/var; 1309 P.coord.push_back(monomial<gen>(gen(it->g),i)); 1310 } 1311 return true; 1312 } 1313 1314 // return true if a good eval point has been found find_good_eval(const polynome & F,const polynome & G,polynome & Fb,polynome & Gb,vecteur & b,bool debuglog,const gen & mod)1315 bool find_good_eval(const polynome & F,const polynome & G,polynome & Fb,polynome & Gb,vecteur & b,bool debuglog,const gen & mod){ 1316 int Fdeg=int(F.lexsorted_degree()),Gdeg=int(G.lexsorted_degree()),nvars=int(b.size()); 1317 gen Fg,Gg; 1318 int essai=0; 1319 int dim=F.dim; 1320 if ( //false && 1321 mod.type==_INT_ && mod.val){ 1322 int modulo=mod.val; 1323 std::vector<hashgcd_U> vars(dim); 1324 vector< T_unsigned<int,hashgcd_U> > f,g,fb,gb; 1325 index_t d(dim); 1326 if (convert(F,G,d,vars,f,g,modulo)){ 1327 vector<int> bi(dim-1); 1328 vecteur2vector_int(b,modulo,bi); 1329 for (;;++essai){ 1330 if (modulo && essai>modulo) 1331 return false; 1332 peval_x2_xn<int,hashgcd_U>(f,bi,vars,fb,modulo); 1333 if (&F==&G) 1334 gb=fb; 1335 else 1336 peval_x2_xn(g,bi,vars,gb,modulo); 1337 if (!fb.empty() && !gb.empty() && int(fb.front().u/vars.front())==Fdeg && int(gb.front().u/vars.front())==Gdeg){ 1338 // convert back fb and gb and return true 1339 convert_from_truncate(fb,vars.front(),Fb); 1340 convert_from_truncate(gb,vars.front(),Gb); 1341 return true; 1342 } 1343 for (int i=0;i<dim-1;++i) 1344 bi[i]=std_rand() % modulo; 1345 } 1346 } 1347 } 1348 for (;;++essai){ 1349 if (!is_zero(mod) && essai>mod.val) 1350 return false; 1351 if (debuglog) 1352 CERR << "Find_good_eval " << CLOCK() << " " << b << '\n'; 1353 Fb=peval_1(F,b,mod); 1354 if (debuglog) 1355 CERR << "Fb= " << CLOCK() << " " << gen(Fb) << '\n'; 1356 if (&F==&G) 1357 Gb=Fb; 1358 else { 1359 Gb=peval_1(G,b,mod); 1360 } 1361 if (debuglog) 1362 CERR << "Gb= " << CLOCK() << " " << gen(Gb) << '\n'; 1363 if ( (Fb.lexsorted_degree()==Fdeg) && (Gb.lexsorted_degree()==Gdeg) ){ 1364 if (debuglog) 1365 CERR << "FOUND good eval" << CLOCK() << " " << b << '\n'; 1366 return true; 1367 } 1368 b=vranm(nvars,0,0); // find another random point 1369 } 1370 } 1371 1372 // It is probably required that 0 is a good evaluation point to 1373 // have an efficient algorithm 1374 // max_gcddeg is used when ezgcd was not successful to find 1375 // the gcd even with 2 evaluations leading to the same gcd degree 1376 // in this case ezgcd calls itself with a bound on the gcd degree 1377 // is_sqff is true if we know that F_orig or G_orig is squarefree 1378 // is_primitive is true if F_orig and G_orig is primitive ezgcd(const polynome & F_orig,const polynome & G_orig,polynome & GCD,bool is_sqff,bool is_primitive,int max_gcddeg,double maxop)1379 bool ezgcd(const polynome & F_orig,const polynome & G_orig,polynome & GCD,bool is_sqff,bool is_primitive,int max_gcddeg,double maxop){ 1380 if (debug_infolevel) 1381 CERR << "// Starting EZGCD dimension " << F_orig.dim << '\n'; 1382 if (F_orig.dim<2){ 1383 #ifdef NO_STDEXCEPT 1384 return false; 1385 #else 1386 setsizeerr(gettext("Args must be multivariate polynomials")); 1387 #endif 1388 } 1389 int Fdeg=F_orig.lexsorted_degree(),Gdeg=G_orig.lexsorted_degree(); 1390 polynome F(F_orig.dim),G(F_orig.dim),cF(F_orig.dim),cG(F_orig.dim),cFG(F_orig.dim); 1391 if (is_primitive){ 1392 cFG=polynome(monomial<gen>(plus_one,0,F_orig.dim)); 1393 cF=cFG; 1394 cG=cFG; 1395 F=F_orig; 1396 G=G_orig; 1397 } 1398 else { 1399 cF=Tlgcd(F_orig); 1400 cG=Tlgcd(G_orig); 1401 cFG=gcd(cF.trunc1(),cG.trunc1()).untrunc1(); 1402 F=F_orig/cF; 1403 G=G_orig/cG; 1404 } 1405 if (Tis_constant(F) || Tis_constant(G) ){ 1406 GCD=cFG; 1407 return true; 1408 } 1409 polynome lcF(Tfirstcoeff(F)),lcG(Tfirstcoeff(G)); 1410 double nop=double(lcF.coord.size())*double(F.coord.size())+double(lcG.coord.size())*double(G.coord.size()); 1411 if (maxop>0){ 1412 if (maxop<nop/10) 1413 return false; 1414 } 1415 vecteur b(F.dim-1); 1416 polynome Fb(1),Gb(1),Db(1); 1417 int old_gcddeg; 1418 for (;;){ 1419 if (debug_infolevel) 1420 CERR << "// Back to EZGCD dimension " << F_orig.dim << '\n'; 1421 find_good_eval(F,G,Fb,Gb,b); 1422 Db=gcd(Fb,Gb); 1423 old_gcddeg=Db.lexsorted_degree(); 1424 if (debug_infolevel) 1425 CERR << "// Eval at " << b << " gcd degree " << old_gcddeg << '\n'; 1426 if (!old_gcddeg){ 1427 GCD=cFG; 1428 return true; 1429 } 1430 if ( (!max_gcddeg) || (old_gcddeg<max_gcddeg) ) 1431 break; 1432 } 1433 polynome new_Fb(1),new_Gb(1),quo(F.dim),rem(F.dim); 1434 for (;;){ 1435 vecteur new_b(vranm(F.dim-1,0,0)); 1436 find_good_eval(F,G,new_Fb,new_Gb,new_b); 1437 if (b==new_b) 1438 continue; 1439 polynome new_Db(gcd(new_Fb,new_Gb)); 1440 int new_gcddeg=new_Db.lexsorted_degree(); 1441 if (debug_infolevel) 1442 CERR << "// Eval at " << new_b << " gcd degree " << new_gcddeg << '\n'; 1443 if (!new_gcddeg){ 1444 GCD=cFG; 1445 return true; 1446 } 1447 if (new_gcddeg>old_gcddeg) // bad evaluation point 1448 continue; 1449 if (new_gcddeg==old_gcddeg) // might be a good guess! 1450 break; 1451 old_gcddeg=new_gcddeg; 1452 Db=new_Db; 1453 Fb=new_Fb; 1454 Gb=new_Gb; 1455 b=new_b; 1456 } 1457 // Found two times the same degree, try to lift! 1458 if ( (Fdeg<=Gdeg) && (old_gcddeg==Fdeg) ){ 1459 if (G.TDivRem1(F,quo,rem) && rem.coord.empty()){ 1460 GCD= F*cFG; 1461 return true; 1462 } 1463 } 1464 if ( (Gdeg<Fdeg) && (old_gcddeg==Gdeg) ){ 1465 if (G.TDivRem1(F,quo,rem) && rem.coord.empty()){ 1466 GCD=G*cFG; 1467 return true; 1468 } 1469 } 1470 if (debug_infolevel) 1471 CERR << "// EZGCD degree " << old_gcddeg << '\n'; 1472 if ( (old_gcddeg==Fdeg) || (old_gcddeg==Gdeg) ) 1473 return false; 1474 // this algo is fast if 0 is a good eval & the degree of the gcd is small 1475 if (!is_zero(b)) 1476 return false; 1477 //if ( (old_gcddeg>4) && (old_gcddeg>Fdeg/4) && (old_gcddeg>Gdeg/4) ) 1478 // return false; 1479 polynome cofacteur(Fb/Db); 1480 if (Tis_constant(gcd(cofacteur,Db))){ 1481 // lift Fb/Db *Db, more precisely insure that lc of each factor 1482 // is lcF(b) 1483 gen lcFb(peval_back(lcF,b)); 1484 if (lcFb.type==_POLY) 1485 lcFb=lcFb._POLYptr->coord.front().value; 1486 Db=(lcFb*Db)/Db.coord.front().value; 1487 cofacteur=(lcFb*cofacteur)/cofacteur.coord.front().value; 1488 polynome liftF(F*lcF); 1489 polynome D(F_orig.dim),cofacteur_F(F_orig.dim),quo,rem; 1490 if (hensel_lift(liftF,lcF,cofacteur,Db,b,cofacteur_F,D,!Tis_constant(lcF),maxop) ){ 1491 D=D/Tlgcd(D); 1492 if (F.TDivRem1(D,quo,rem) && is_zero(rem) && G.TDivRem1(D,quo,rem) && is_zero(rem)){ 1493 GCD=D*cFG; 1494 return true; 1495 } 1496 } 1497 return false; 1498 } 1499 cofacteur=Gb/Db; 1500 if (Tis_constant(gcd(cofacteur,Db))){ 1501 // lift Gb/Db *Db, more precisely insure that lc of each factor 1502 // is lcG(b) 1503 gen lcGb(peval_back(lcG,b)); 1504 if (lcGb.type==_POLY) 1505 lcGb=lcGb._POLYptr->coord.front().value; 1506 Db=(lcGb*Db)/Db.coord.front().value; 1507 cofacteur=(lcGb*cofacteur)/cofacteur.coord.front().value; 1508 polynome liftG(G*lcG); 1509 polynome D(G_orig.dim),cofacteur_G(G_orig.dim),quo,rem; 1510 if (hensel_lift(liftG,lcG,cofacteur,Db,b,cofacteur_G,D,!Tis_constant(lcG),maxop) ){ 1511 D=D/Tlgcd(D); 1512 if (F.TDivRem1(D,quo,rem) && is_zero(rem) && G.TDivRem1(D,quo,rem) && is_zero(rem)){ 1513 GCD=D*cFG; 1514 return true; 1515 } 1516 } 1517 return false; 1518 } 1519 // FIXME find an integer j such that (F+jG)/D_b is coprime with D_b 1520 return false; 1521 } 1522 1523 // algorithm=0 for HEUGCD, 1 for PRS, 2 for EZGCD, 3 for MODGCD heugcd_psrgcd_ezgcd_modgcd(const gen & args,int algorithm,GIAC_CONTEXT)1524 static gen heugcd_psrgcd_ezgcd_modgcd(const gen & args,int algorithm,GIAC_CONTEXT){ 1525 vecteur & v=*args._VECTptr; 1526 gen p1(v[0]),p2(v[1]),n1,n2,d1,d2; 1527 vecteur lv; 1528 if ( (v.size()==3) && (v[2].type==_VECT) ) 1529 lv=*v[2]._VECTptr; 1530 lvar(p1,lv); 1531 lvar(p2,lv); 1532 p1=e2r(p1,lv,contextptr); 1533 fxnd(p1,n1,d1); 1534 p2=e2r(p2,lv,contextptr); 1535 fxnd(p2,n2,d2); 1536 gen res,np_simp,nq_simp,d_content; 1537 polynome p,q,p_gcd; 1538 if ( (n1.type!=_POLY) || (n2.type!=_POLY) ) 1539 res=gcd(n1,n2,contextptr); 1540 else { 1541 polynome pres; 1542 bool result=false; 1543 switch(algorithm){ 1544 case 0: 1545 p_gcd.dim=n1._POLYptr->dim; 1546 result=gcdheu(*n1._POLYptr,*n2._POLYptr,p,np_simp,q,nq_simp,p_gcd,d_content,true); 1547 pres=p_gcd*d_content; 1548 break; 1549 case 1: 1550 pres=gcdpsr(*n1._POLYptr,*n2._POLYptr); 1551 result=true; 1552 break; 1553 case 2: 1554 result=ezgcd(*n1._POLYptr,*n2._POLYptr,pres); 1555 break; 1556 case 3: 1557 result=gcd_modular_algo(*n1._POLYptr,*n2._POLYptr,pres,false); 1558 break; 1559 } 1560 if (result) 1561 res=pres; 1562 else 1563 return gensizeerr(gettext("GCD not successfull")); 1564 } 1565 return r2e(res,lv,contextptr); 1566 } 1567 _ezgcd(const gen & args,GIAC_CONTEXT)1568 gen _ezgcd(const gen & args,GIAC_CONTEXT){ 1569 if ( args.type==_STRNG && args.subtype==-1) return args; 1570 if ( (args.type!=_VECT) || (args._VECTptr->size()<2) ) 1571 return symbolic(at_ezgcd,args); 1572 return heugcd_psrgcd_ezgcd_modgcd(args,2,contextptr); 1573 } 1574 static const char _ezgcd_s []="ezgcd"; 1575 static define_unary_function_eval (__ezgcd,&_ezgcd,_ezgcd_s); 1576 define_unary_function_ptr5( at_ezgcd ,alias_at_ezgcd,&__ezgcd,0,true); 1577 _modgcd(const gen & args,GIAC_CONTEXT)1578 gen _modgcd(const gen & args,GIAC_CONTEXT){ 1579 if ( args.type==_STRNG && args.subtype==-1) return args; 1580 if ( (args.type!=_VECT) || (args._VECTptr->size()<2) ) 1581 return symbolic(at_modgcd,args); 1582 return heugcd_psrgcd_ezgcd_modgcd(args,3,contextptr); 1583 } 1584 static const char _modgcd_s []="modgcd"; 1585 static define_unary_function_eval (__modgcd,&_modgcd,_modgcd_s); 1586 define_unary_function_ptr5( at_modgcd ,alias_at_modgcd,&__modgcd,0,true); 1587 _heugcd(const gen & args,GIAC_CONTEXT)1588 gen _heugcd(const gen & args,GIAC_CONTEXT){ 1589 if ( args.type==_STRNG && args.subtype==-1) return args; 1590 if ( (args.type!=_VECT) || (args._VECTptr->size()<2) ) 1591 return symbolic(at_heugcd,args); 1592 return heugcd_psrgcd_ezgcd_modgcd(args,0,contextptr); 1593 } 1594 static const char _heugcd_s []="heugcd"; 1595 static define_unary_function_eval (__heugcd,&_heugcd,_heugcd_s); 1596 define_unary_function_ptr5( at_heugcd ,alias_at_heugcd,&__heugcd,0,true); 1597 _psrgcd(const gen & args,GIAC_CONTEXT)1598 gen _psrgcd(const gen & args,GIAC_CONTEXT){ 1599 if ( args.type==_STRNG && args.subtype==-1) return args; 1600 if ( (args.type!=_VECT) || (args._VECTptr->size()<2) ) 1601 return symbolic(at_psrgcd,args); 1602 return heugcd_psrgcd_ezgcd_modgcd(args,1,contextptr); 1603 } 1604 static const char _psrgcd_s []="psrgcd"; 1605 static define_unary_function_eval (__psrgcd,&_psrgcd,_psrgcd_s); 1606 define_unary_function_ptr5( at_psrgcd ,alias_at_psrgcd,&__psrgcd,0,true); 1607 1608 1609 #ifndef NO_NAMESPACE_GIAC 1610 } // namespace giac 1611 #endif // ndef NO_NAMESPACE_GIAC 1612