1 // -*- mode:C++ ; compile-command: "g++-3.4 -DHAVE_CONFIG_H -I. -I.. -DIN_GIAC -g -c sym2poly.cc" -*- 2 #include "giacPCH.h" 3 /* 4 * This file implements several functions that work on univariate and 5 * multivariate polynomials and rational functions. 6 * These functions include polynomial quotient and remainder, GCD and LCM 7 * computation, factorization and rational function normalization. */ 8 9 /* 10 * Copyright (C) 2000,2014 B. Parisse, Institut Fourier, 38402 St Martin d'Heres 11 * 12 * This program is free software; you can redistribute it and/or modify 13 * it under the terms of the GNU General Public License as published by 14 * the Free Software Foundation; either version 3 of the License, or 15 * (at your option) any later version. 16 * 17 * This program is distributed in the hope that it will be useful, 18 * but WITHOUT ANY WARRANTY; without even the implied warranty of 19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 20 * GNU General Public License for more details. 21 * 22 * You should have received a copy of the GNU General Public License 23 * along with this program. If not, see <http://www.gnu.org/licenses/>. 24 */ 25 using namespace std; 26 #if !defined GIAC_HAS_STO_38 && !defined NSPIRE && !defined FXCG && !defined POCKETCAS 27 #include <fstream> 28 #endif 29 #include <string> 30 #include "gen.h" 31 #include "sym2poly.h" 32 #include "usual.h" 33 #include "unary.h" 34 #include "subst.h" 35 #include "modpoly.h" 36 #include "alg_ext.h" 37 #include "solve.h" 38 #include "input_parser.h" 39 #include "ezgcd.h" 40 #include "prog.h" 41 #include "ifactor.h" 42 #include "poly.h" 43 #include "plot.h" 44 #include "misc.h" 45 #include "giacintl.h" 46 #ifdef HAVE_UNISTD_H 47 #include <unistd.h> 48 #endif 49 50 #ifndef NO_NAMESPACE_GIAC 51 namespace giac { 52 #endif // ndef NO_NAMESPACE_GIAC 53 54 // "instantiate" debugging functions dbgprint(const polynome & p)55 void dbgprint(const polynome &p){ 56 p.dbgprint(); 57 } 58 dbgprint(const gen & e)59 void dbgprint(const gen & e){ 60 e.dbgprint(); 61 } 62 cdr_VECT(const vecteur & l)63 vecteur cdr_VECT(const vecteur & l){ 64 if (l.empty()) 65 return vecteur(l); 66 vecteur::const_iterator it=l.begin(),itend=l.end(); 67 vecteur res; 68 ++it; 69 for (;it!=itend;++it) 70 res.push_back(*it); 71 return vecteur(res); 72 } 73 74 //*************************** 75 // functions relative to lvar 76 //*************************** equalposcomp(const vecteur & l,const gen & e)77 int equalposcomp(const vecteur & l,const gen & e){ 78 int n=1; 79 for (vecteur::const_iterator it=l.begin();it!=l.end();++it){ 80 if ((*it)==e) 81 return(n); 82 else 83 n++; 84 } 85 return(0); 86 } 87 addtolvar(const gen & e,vecteur & l)88 void addtolvar(const gen & e, vecteur & l){ 89 if (equalposcomp(l,e)) 90 return; 91 l.push_back(e); 92 } 93 lvar_symbolic(const gen & g,vecteur & l)94 static void lvar_symbolic(const gen & g, vecteur &l){ 95 const symbolic & s = *g._SYMBptr; 96 if ( (s.sommet==at_plus) || (s.sommet==at_prod)){ 97 if (s.feuille.type!=_VECT){ 98 lvar(s.feuille,l); 99 return; 100 } 101 vecteur::iterator it=s.feuille._VECTptr->begin(), itend=s.feuille._VECTptr->end(); 102 for (;it!=itend;++it) 103 lvar(*it,l); 104 return; 105 } 106 if ( (s.sommet==at_neg) || (s.sommet==at_inv) ){ 107 lvar(s.feuille,l); 108 return; 109 } 110 if ( (s.sommet==at_pow) && ( (*s.feuille._VECTptr)[1].type==_INT_)) 111 lvar(s.feuille._VECTptr->front(),l); 112 else 113 addtolvar(g,l); 114 } 115 lvar(const vecteur & v,vecteur & l)116 void lvar(const vecteur & v,vecteur & l){ 117 vecteur::const_iterator it=v.begin(),itend=v.end(); 118 for (;it!=itend;++it) 119 lvar(*it,l); 120 } 121 lvar(const sparse_poly1 & p,vecteur & l)122 void lvar(const sparse_poly1 & p,vecteur & l){ 123 sparse_poly1::const_iterator it=p.begin(),itend=p.end(); 124 for (;it!=itend;++it){ 125 lvar(it->coeff,l); 126 // lvar(it->exponent,l); 127 } 128 } 129 lvar(const gen & e,vecteur & l)130 void lvar(const gen & e, vecteur & l) { 131 switch (e.type){ 132 case _INT_: case _DOUBLE_: case _ZINT: case _CPLX: case _POLY: case _EXT: case _USER: case _REAL: 133 return; 134 case _IDNT: 135 if (strcmp(e._IDNTptr->id_name,string_undef)) 136 addtolvar(e,l); 137 return ; 138 case _SYMB: 139 lvar_symbolic(e,l); 140 return ; 141 case _VECT: 142 lvar(*e._VECTptr,l); 143 return ; 144 case _FRAC: 145 lvar(e._FRACptr->num,l); 146 lvar(e._FRACptr->den,l); 147 return; 148 case _MOD: 149 lvar(*e._MODptr,l); 150 lvar(*(e._MODptr+1),l); 151 return; 152 case _SPOL1: 153 lvar(*e._SPOL1ptr,l); 154 return; 155 default: 156 return; 157 } 158 } 159 lvar(const gen & e)160 vecteur lvar(const gen & e){ 161 vecteur l; 162 lvar(e,l); 163 return l; 164 } 165 symb_lvar(const gen & e)166 gen symb_lvar(const gen & e){ 167 return symbolic(at_lvar,e); 168 } cklvar(const gen & e,GIAC_CONTEXT)169 gen cklvar(const gen & e,GIAC_CONTEXT){ 170 vecteur l; 171 lvar(e,l); 172 return l; 173 } 174 static const char _lvar_s []="lvar"; 175 static define_unary_function_eval (__lvar,&cklvar,_lvar_s); 176 define_unary_function_ptr5( at_lvar ,alias_at_lvar,&__lvar,0,true); 177 is_algebraic_EXTension(const gen & e)178 static bool is_algebraic_EXTension(const gen & e){ 179 if (e.type==_EXT) 180 return true; 181 if (e.type!=_SYMB) 182 return false; 183 if ( (e._SYMBptr->sommet==at_sqrt) || (e._SYMBptr->sommet==at_rootof) ) 184 return true; 185 if ( (e._SYMBptr->sommet==at_pow) && (e._SYMBptr->feuille._VECTptr->back().type==_FRAC) && (e._SYMBptr->feuille._VECTptr->back()._FRACptr->den.type==_INT_) &&(absint(e._SYMBptr->feuille._VECTptr->back()._FRACptr->den.val)<=MAX_ALG_EXT_ORDER_SIZE) ) 186 return true; 187 return false; 188 } 189 algebraic_argument(const gen & e)190 static gen algebraic_argument(const gen & e){ 191 if (e.type==_EXT) 192 return makevecteur(*e._EXTptr,*(e._EXTptr+1)); 193 if (e.type!=_SYMB) 194 return gensizeerr(gettext("sym2poly.cc/algebraic_argument")); 195 if ((e._SYMBptr->sommet==at_sqrt) || (e._SYMBptr->sommet==at_rootof) ) 196 return e._SYMBptr->feuille; 197 if ( (e._SYMBptr->sommet==at_pow) && (e._SYMBptr->feuille._VECTptr->back().type==_FRAC) && (e._SYMBptr->feuille._VECTptr->back()._FRACptr->den.type==_INT_) ) 198 return e._SYMBptr->feuille._VECTptr->front(); 199 return gentypeerr(gettext("algebraic_argument")); 200 } 201 equalposmat(const matrice & m,const gen & e,int & i,int & j)202 static bool equalposmat(const matrice & m,const gen & e,int &i,int & j){ 203 i=0; 204 const_iterateur it=m.begin(),itend=m.end(),jt,jtend; 205 for (;it!=itend;++it,++i){ 206 if (*it==e){ 207 j=-1; 208 return true; 209 } 210 else { 211 if (it->type!=_VECT){ 212 #ifdef NO_STDEXCEPT 213 return false; 214 #else 215 setsizeerr(gettext("sym2poly.cc/equalposmat")); 216 #endif 217 } 218 for (j=0,jt=it->_VECTptr->begin(),jtend=it->_VECTptr->end();jt!=jtend;++jt,++j) 219 if (*jt==e) 220 return true; 221 } 222 } 223 return false; 224 } 225 addfirstrow(const gen & e,matrice & m)226 static void addfirstrow(const gen & e,matrice & m){ 227 if (m.empty()){ 228 vecteur v(1,e); 229 m.push_back(v); 230 } 231 else { 232 if (m.front().type!=_VECT){ 233 #ifdef NO_STDEXCEPT 234 return; 235 #else 236 setsizeerr(gettext("sym2poly.cc/addfirstrow")); 237 #endif 238 } 239 vecteur v(*m.front()._VECTptr); 240 v.push_back(e); 241 m.front()=v; 242 } 243 } 244 ext_glue_matrices(const matrice & a,const matrice & b)245 static matrice ext_glue_matrices(const matrice & a,const matrice & b){ 246 if (a.size()>b.size()) 247 return ext_glue_matrices(b,a); 248 // Algorithm should be fixed in all generality 249 // the loop below fixes assume(a>0); normal(-sqrt(a)*sqrt(a*pi)*sqrt(pi)) 250 // a.size()<=b.size() 251 if (a.size()==b.size()){ 252 for (int i=0;i<a.size();++i){ 253 if (a[i].type==_VECT && b[i].type==_VECT){ 254 if (a[i]._VECTptr->size()<b[i]._VECTptr->size()) 255 break; 256 if (a[i]._VECTptr->size()>b[i]._VECTptr->size()) 257 return ext_glue_matrices(b,a); 258 } 259 } 260 } 261 if (b.empty() || a.empty() || (a==b)) 262 return b; 263 int i,j; 264 matrice res(b.begin()+1,b.end()); // all alg. extensions of b 265 // look in a for vars that are not inside res 266 matrice a_not_in_res; 267 const_iterateur it=a.begin(),itend=a.end(),jt,jtend; 268 for (;it!=itend;++it){ 269 vecteur temp; 270 for (jt=it->_VECTptr->begin(),jtend=it->_VECTptr->end();jt!=jtend;++jt){ 271 if (!equalposmat(res,*jt,i,j)) 272 temp.push_back(*jt); 273 } 274 if (it==a.begin() || (!temp.empty()) ) 275 a_not_in_res.push_back(temp); 276 } 277 // look in the first row of b for vars that are not inside a_not_in_res 278 jt=b.begin()->_VECTptr->begin(),jtend=b.begin()->_VECTptr->end(); 279 for (;jt!=jtend;++jt){ 280 if (!equalposmat(a_not_in_res,*jt,i,j)) 281 addfirstrow(*jt,a_not_in_res); 282 } 283 return mergevecteur(a_not_in_res,res); 284 } 285 alg_lvar(const sparse_poly1 & p,vecteur & l)286 void alg_lvar(const sparse_poly1 & p,vecteur & l){ 287 sparse_poly1::const_iterator it=p.begin(),itend=p.end(); 288 for (;it!=itend;++it){ 289 alg_lvar(it->coeff,l); 290 // lvar(it->exponent,l); 291 } 292 } 293 294 // return true if there is at least one algebraic extension alg_lvar(const gen & e,matrice & m)295 void alg_lvar(const gen & e,matrice & m){ 296 vecteur temp; 297 lvar(e,temp); 298 int i,j; 299 // For each variable of temp, 300 // if not alg var look if still inside m else add it to the first line 301 // else make a "merge" 302 const_iterateur it=temp.begin(),itend=temp.end(); 303 for (;it!=itend;++it){ 304 if ( !is_algebraic_EXTension(*it) ){ 305 if (!equalposmat(m,*it,i,j)){ 306 addfirstrow(*it,m); 307 } 308 } 309 else { // *it is an algebraic extension! 310 #if 1 311 matrice ext_mat; 312 vecteur v,vt; 313 ext_mat.push_back(v); 314 vt=alg_lvar(algebraic_argument(*it)); 315 int s=int(vt.size()); 316 if (s>1 || (s==1 && !vt.front()._VECTptr->empty()) ){ 317 ext_mat=mergevecteur(ext_mat,vt); 318 } 319 m=ext_glue_matrices(ext_mat,m); 320 #else 321 vecteur vt; 322 alg_lvar(algebraic_argument(*it),vt); 323 if (vt==m) return; 324 matrice ext_mat; 325 vecteur v; 326 ext_mat.push_back(v); 327 int s=int(vt.size()); 328 if (s>1 || (s==1 && !vt.front()._VECTptr->empty()) ) 329 ext_mat=mergevecteur(ext_mat,vt); 330 m=ext_glue_matrices(ext_mat,m); 331 #endif 332 } 333 } 334 } 335 alg_lvar_f(const gen & g,GIAC_CONTEXT)336 gen alg_lvar_f(const gen & g,GIAC_CONTEXT){ 337 vecteur w=lvar(g); 338 return symbolic(at_pow,makesequence(_product(w,contextptr),plus_one_half)); 339 } alg_lvar(const gen & e)340 vecteur alg_lvar(const gen & e){ 341 vecteur l; 342 l.push_back(l); // insure a null line inside the matrix of alg_lvar 343 if (has_op(e,*at_rootof)){ 344 vecteur w0=lvar(e),w; 345 for (int i=0;i<w0.size();++i){ 346 if (w0[i].is_symb_of_sommet(at_rootof)) 347 w.push_back(w0[i]); 348 } 349 vector<const unary_function_ptr *> vu; 350 vu.push_back(at_rootof); 351 vector <gen_op_context> vv; 352 vv.push_back(alg_lvar_f); // was _nop 353 // FIXME: arg of e is two vectors, if rootof is ^, this raises a warning 354 gen er=subst(w,vu,vv,false,context0); 355 if (er.type==_VECT) 356 er=subst(e,w,*er._VECTptr,false,context0); 357 else 358 er=e; 359 alg_lvar(er,l); 360 } 361 else 362 alg_lvar(e,l); 363 return l; 364 } 365 symb_algvar(const gen & e)366 gen symb_algvar(const gen & e){ 367 return symbolic(at_algvar,e); 368 } ckalgvar(const gen & e,GIAC_CONTEXT)369 gen ckalgvar(const gen & e,GIAC_CONTEXT){ 370 vecteur l; 371 alg_lvar(e,l); 372 return l; 373 } 374 static const char _algvar_s []="algvar"; 375 static define_unary_function_eval (__algvar,&ckalgvar,_algvar_s); 376 define_unary_function_ptr5( at_algvar ,alias_at_algvar,&__algvar,0,true); 377 378 379 //*********************************************** 380 // functions relative to fractions for efficiency 381 //*********************************************** 382 /* 383 static fraction fpow(const fraction & p,const gen & n){ 384 if (n.type!=_INT_) 385 setsizeerr(gettext("sym2poly.cc/fraction pow")); 386 return pow(p,n.val); 387 } 388 */ 389 ref_polynome2gen(ref_polynome * refptr)390 gen ref_polynome2gen(ref_polynome * refptr){ 391 if (refptr->t.coord.empty()){ 392 delete refptr; 393 return 0; 394 } 395 if (refptr->t.coord.front().index.is_zero() && is_atomic(refptr->t.coord.front().value) ){ 396 gen res=refptr->t.coord.front().value; 397 delete refptr; 398 return res; 399 } 400 return refptr; 401 } 402 monomial2gen(const monomial<gen> & m)403 gen monomial2gen(const monomial<gen> & m){ 404 if (m.index.is_zero() && is_atomic(m.value) ){ 405 return m.value; 406 } 407 return new ref_polynome(m); 408 } 409 simplify3(gen & n,gen & d)410 gen simplify3(gen & n,gen & d){ 411 if (is_one(n) || is_one(d)) 412 return plus_one; 413 if (n.type==_EXT){ 414 gen n_EXT=*n._EXTptr; 415 gen g=simplify(n_EXT,d); 416 if (!is_one(g)) n=algebraic_EXTension(n_EXT,*(n._EXTptr+1)); 417 return g; 418 } 419 if ((n.type==_POLY) && (d.type==_POLY)){ 420 ref_polynome * pptr = new ref_polynome(n._POLYptr->dim); 421 if ( 422 #ifndef SMARTPTR64 423 n.__POLYptr->ref_count==1 && d.__POLYptr->ref_count==1 424 #else 425 ((ref_polynome*)(* (ulonglong *) &n >> 16))->ref_count ==1 && 426 ((ref_polynome*)(* (ulonglong *) &d >> 16))->ref_count ==1 427 #endif 428 ){ 429 simplify(*n._POLYptr,*d._POLYptr,pptr->t); 430 return pptr; 431 } 432 polynome np(*n._POLYptr),dp(*d._POLYptr); 433 simplify(np,dp,pptr->t); 434 gen tmpmult(plus_one); 435 lcmdeno(pptr->t,tmpmult); 436 if (tmpmult.type==_POLY) 437 tmpmult=polynome(monomial<gen>(tmpmult,index_t(pptr->t.dim))); 438 gen g; 439 if (is_one(tmpmult)) 440 g=pptr; 441 else 442 g=fraction(tmpmult*pptr,tmpmult); // polynome(tmpmult,np.dim)); 443 // changed 4 nov. 2012 failure on f(x):=sqrt(x)/(x+1); F:=normal(f'(x)); 444 // Fix 9/2/2013: np and dp may have denominators if there are _EXT coeffs 445 tmpmult=1; 446 lcmdeno(np,tmpmult); 447 lcmdeno(dp,tmpmult); 448 n=np*tmpmult; 449 d=dp*tmpmult; 450 // moved 18/3/2014 for simplify(acos2atan(acos(sqrt(1-x^2))+acos(x))) 451 if (tmpmult.type==_POLY) tmpmult=polynome(monomial<gen>(tmpmult,index_t(pptr->t.dim))); 452 if (is_one(tmpmult)) 453 return g; 454 return fraction(g,tmpmult); // g*tmpmult; 455 } 456 if (n.type==_POLY) { 457 gen l; 458 vector< monomial<gen> > :: const_iterator it=n._POLYptr->coord.begin(),itend=n._POLYptr->coord.end(); 459 if (it<itend-1){ 460 l=gcd(it->value,(itend-1)->value,context0); 461 ++it; --itend; 462 } 463 for (;it!=itend;++it){ 464 l=gcd(l,it->value,context0); 465 if (is_one(l)) 466 return l; 467 } 468 gen g=simplify3(l,d); 469 if (is_one(g)) 470 return g; 471 polynome np(*n._POLYptr); 472 np=np/g; 473 n=np; 474 if (g.type>_DOUBLE_){ 475 // FIXME embedd g and d inside a polynomial like np was 476 polynome pg(np.dim); 477 pg.coord.push_back(monomial<gen>(g,np.dim)); 478 g=pg; 479 polynome pd(np.dim); 480 pd.coord.push_back(monomial<gen>(d,np.dim)); 481 d=pd; 482 } 483 return g; 484 } 485 if (d.type==_POLY){ 486 polynome np(n,d._POLYptr->dim),dp(*d._POLYptr); 487 polynome g(np.dim); 488 g=simplify(np,dp); 489 n=np; 490 d=dp; 491 return g; 492 } 493 gen g=gcd(n,d,context0); 494 n=n/g; // iquo(n,g); 495 d=d/g; // iquo(d,g); 496 return g; 497 } 498 has_EXT(const gen & g)499 static bool has_EXT(const gen & g){ 500 if (g.type==_EXT) 501 return true; 502 if (g.type!=_POLY) 503 return false; 504 polynome & p = *g._POLYptr; 505 vector< monomial<gen> >::const_iterator it=p.coord.begin(),itend=p.coord.end(); 506 for (;it!=itend;++it){ 507 if (has_EXT(it->value)) 508 return true; 509 } 510 return false; 511 } 512 _FRACadd(const gen & n1,const gen & d1,const gen & n2,const gen & d2,gen & num,gen & den)513 static void _FRACadd(const gen & n1, const gen & d1,const gen & n2, const gen & d2, gen & num, gen & den){ 514 // COUT << n1 << "/" << d1 << "+" << n2 << "/" << d2 << "="; 515 if (is_one(d1)){ 516 num=n1*d2+n2; 517 den=d2; 518 return; 519 } 520 if (is_one(d2)){ 521 num=n2*d1+n1; 522 den=d1; 523 // COUT << num << "/" << den << "\n"; 524 return; 525 } 526 // n1/d1+n2/d2 with g=gcd(d1,d2), d1=d1g*g, d2=d2g*g is 527 // (n1*d2g+n2*d1g)/g * 1/(d1g*d2g) 528 gen d1g(d1),d2g(d2),coeff(1); 529 den=simplify3(d1g,d2g); 530 if (den.type==_FRAC){ 531 coeff=den._FRACptr->den; 532 den=den._FRACptr->num; 533 } 534 #if 0 535 gen n1g(n1),n2g(n2); 536 num=simplify3(n1g,n2g); 537 if (num.type==_FRAC){ 538 den=den*num._FRACptr->den; 539 num=num._FRACptr->num; 540 } 541 n1g=(n1g*d2g+n2g*d1g)*coeff; 542 simplify3(n1g,den); 543 num=num*n1g; 544 den=den*d1g*d2g; 545 #else 546 num=(n1*d2g+n2*d1g)*coeff; 547 simplify3(num,den); 548 den=den*d1g*d2g; 549 #endif 550 } 551 _FRACmul(const gen & n1,const gen & d1,const gen & n2,const gen & d2,gen & num,gen & den)552 static void _FRACmul(const gen & n1, const gen & d1,const gen & n2, const gen & d2, gen & num, gen & den){ 553 // COUT << n1 << "/" << d1 << "*" << n2 << "/" << d2 << "="; 554 if (is_one(d1)){ 555 num=n1; 556 den=d2; 557 simplify3(num,den); 558 num=num*n2; 559 // COUT << num << "/" << den << "\n"; 560 return; 561 } 562 if (is_one(d2)){ 563 num=n2; 564 den=d1; 565 simplify3(num,den); 566 num=num*n1; 567 // COUT << num << "/" << den << "\n"; 568 return; 569 } 570 num=n1; 571 den=d2; 572 simplify3(num,den); 573 gen ntemp(n2),dtemp(d1); 574 simplify3(ntemp,dtemp); 575 num=num*ntemp; 576 den=den*dtemp; 577 // Further simplifications may occur with _EXT multiplications 578 if (has_EXT(ntemp)) 579 simplify3(num,den); 580 // COUT << num << "/" << den << "\n"; 581 } 582 583 //********************************** 584 // symbolic to tensor 585 //********************************** sym2radd(vecteur::const_iterator debut,vecteur::const_iterator fin,const gen & iext,const vecteur & l,const vecteur & lv,const vecteur & lvnum,const vecteur & lvden,int l_size,gen & num,gen & den,GIAC_CONTEXT,bool sequentially)586 static bool sym2radd (vecteur::const_iterator debut,vecteur::const_iterator fin,const gen & iext,const vecteur &l,const vecteur & lv, const vecteur & lvnum,const vecteur & lvden, int l_size, gen & num, gen & den,GIAC_CONTEXT,bool sequentially){ 587 bool totally_converted=true; 588 if (sequentially || fin-debut<4){ 589 gen n1,d1,n2,d2; 590 num=zero; 591 den=plus_one; 592 for (;debut!=fin;++debut){ 593 totally_converted=totally_converted && sym2r(*debut,iext,l,lv,lvnum,lvden,l_size,n1,d1,contextptr); 594 n2=num; 595 d2=den; 596 _FRACadd(n1,d1,n2,d2,num,den); 597 } 598 } 599 else { 600 vecteur::const_iterator milieu=debut+(fin-debut)/2; 601 gen n1,d1,n2,d2; 602 totally_converted=totally_converted && sym2radd(debut,milieu,iext,l,lv,lvnum,lvden,l_size,n1,d1,contextptr,sequentially); 603 totally_converted=totally_converted && sym2radd(milieu,fin,iext,l,lv,lvnum,lvden,l_size,n2,d2,contextptr,sequentially); 604 _FRACadd(n1,d1,n2,d2,num,den); 605 } 606 return totally_converted; 607 } 608 609 /* 610 static bool sym2rmulold (vecteur::const_iterator debut,vecteur::const_iterator fin,const vecteur &l, const vecteur & lv, const vecteur & lvnum,const vecteur & lvden,int l_size, gen & num, gen & den,GIAC_CONTEXT){ 611 bool totally_converted=true; 612 // First check for a "normal" monomial 613 gen coeff=plus_one; 614 if (!l.empty()){ 615 bool embedd = l.front().type==_VECT ; 616 vecteur l1; 617 if (embedd) 618 l1=*l.front()._VECTptr; 619 else 620 l1=l; 621 gen tmp; 622 int pui,pos; 623 index_t i(l_size); 624 for (;debut!=fin;++debut){ 625 if (debut->type<_IDNT) 626 coeff=coeff*(*debut); 627 else { 628 if (debut->is_symb_of_sommet(at_pow) && debut->_SYMBptr->feuille._VECTptr->back().type==_INT_){ 629 tmp=debut->_SYMBptr->feuille._VECTptr->front(); 630 pui=debut->_SYMBptr->feuille._VECTptr->back().val; 631 } 632 else { 633 tmp=*debut; 634 pui=1; 635 } 636 if ( tmp.type==_IDNT && (pos=equalposcomp(l1,tmp)) && pui>=0){ 637 i[pos-1] += pui; 638 } 639 else 640 break; 641 } 642 } 643 if (!is_zero(coeff)) 644 coeff=monomial2gen(monomial<gen>(coeff,i)); 645 } 646 if (fin-debut<4){ 647 gen n1,d1,n2,d2; 648 num=coeff; 649 den=plus_one; 650 for (;debut!=fin;++debut){ 651 totally_converted=totally_converted && sym2r(*debut,iext,l,lv,lvnum,lvden,l_size,n1,d1,contextptr); 652 n2=num; 653 d2=den; 654 _FRACmul(n1,d1,n2,d2,num,den); 655 } 656 } 657 else { 658 vecteur::const_iterator milieu=debut+(fin-debut)/2; 659 gen n1,d1,n2,d2; 660 totally_converted=totally_converted && sym2rmulold(debut,milieu,l,lv,lvnum,lvden,l_size,n1,d1,contextptr); 661 totally_converted=totally_converted && sym2rmulold(milieu,fin,l,lv,lvnum,lvden,l_size,n2,d2,contextptr); 662 _FRACmul(n1,d1,n2,d2,num,den); 663 simplify3(coeff,den); 664 num=num*coeff; 665 } 666 return totally_converted; 667 } 668 */ 669 670 sym2rmul(vecteur::const_iterator debut,vecteur::const_iterator fin,const gen & iext,const vecteur & l,const vecteur & lv,const vecteur & lvnum,const vecteur & lvden,int l_size,gen & num,gen & den,GIAC_CONTEXT)671 static bool sym2rmul(vecteur::const_iterator debut,vecteur::const_iterator fin,const gen & iext,const vecteur &l, const vecteur & lv, const vecteur & lvnum,const vecteur & lvden,int l_size, gen & num, gen & den,GIAC_CONTEXT){ 672 bool totally_converted=true; 673 // First check for a "normal" monomial 674 gen coeffnum=plus_one,coeffden=plus_one,coeffn=plus_one,coeffd=plus_one; 675 if (!l.empty()){ 676 bool embedd = l.front().type==_VECT ; 677 const vecteur & l1 = embedd?*l.front()._VECTptr:l; 678 gen tmp; 679 int pui,pos; 680 index_t inum(l_size),iden(l_size); 681 for (;debut!=fin;++debut){ 682 if (debut->type<_IDNT) 683 coeffn=coeffn*(*debut); 684 else { 685 if (debut->is_symb_of_sommet(at_inv)){ 686 tmp=debut->_SYMBptr->feuille; 687 if (tmp.type<_IDNT){ 688 coeffd=coeffd*tmp; 689 continue; 690 } 691 pui=-1; 692 } 693 else { 694 if (debut->is_symb_of_sommet(at_pow) && debut->_SYMBptr->feuille._VECTptr->back().type==_INT_){ 695 tmp=debut->_SYMBptr->feuille._VECTptr->front(); 696 pui=debut->_SYMBptr->feuille._VECTptr->back().val; 697 } 698 else { 699 tmp=*debut; 700 pui=1; 701 } 702 } 703 if (absint(pui)>=(1<<15)) 704 coeffn=gensizeerr(gettext("Polynomial exponent overflow.")); 705 if ( (tmp.type==_IDNT || tmp.type==_SYMB) && (pos=equalposcomp(l1,tmp)) ){ 706 if (pui>=0) 707 inum[pos-1] += pui; 708 else 709 iden[pos-1] -= pui; 710 } 711 else 712 break; 713 } 714 } 715 if (!is_exactly_zero(coeffn)){ 716 // simplify inum and iden 717 bool hasnum=false,hasden=false; 718 for (int i=0;i<l_size;i++){ 719 if (inum[i]<iden[i]){ 720 iden[i] -= inum[i]; 721 inum[i] = 0; 722 hasden=true; 723 } 724 else { 725 inum[i] -= iden[i]; 726 iden[i] = 0; 727 hasnum=true; 728 } 729 } 730 simplify3(coeffn,coeffd); 731 if (hasnum) 732 coeffnum=monomial2gen(monomial<gen>(coeffn,inum)); 733 else 734 coeffnum=coeffn; 735 if (hasden) 736 coeffden=monomial2gen(monomial<gen>(coeffd,iden)); 737 else 738 coeffden=coeffd; 739 } 740 else 741 coeffnum=0; 742 } 743 if (iext!=0){ 744 gen numr,numi; 745 reim(coeffden,numr,numi,contextptr); 746 if (!is_zero(numi)){ 747 coeffnum=coeffnum*(numr-cst_i*numi); 748 coeffden=numr*numr+numi*numi; 749 } 750 reim(coeffnum,numr,numi,contextptr); 751 if (!is_zero(numi)){ 752 coeffnum=numr+iext*numi; 753 fxnd(coeffnum,numr,numi); 754 coeffnum=numr; 755 coeffden=coeffden*numi; 756 } 757 } 758 if (fin-debut<4){ 759 num=coeffnum; 760 den=coeffden; 761 gen n1,d1,n2,d2; 762 for (;debut!=fin;++debut){ 763 totally_converted=totally_converted && sym2r(*debut,iext,l,lv,lvnum,lvden,l_size,n1,d1,contextptr); 764 n2=num; 765 d2=den; 766 _FRACmul(n1,d1,n2,d2,num,den); 767 } 768 } 769 else { 770 vecteur::const_iterator milieu=debut+(fin-debut)/2; 771 gen n1,d1,n2,d2,n3,d3; 772 totally_converted=totally_converted && sym2rmul(debut,milieu,iext,l,lv,lvnum,lvden,l_size,n1,d1,contextptr); 773 totally_converted=totally_converted && sym2rmul(milieu,fin,iext,l,lv,lvnum,lvden,l_size,n2,d2,contextptr); 774 _FRACmul(n1,d1,n2,d2,n3,d3); 775 _FRACmul(coeffnum,coeffden,n3,d3,num,den); 776 } 777 return totally_converted; 778 } 779 sym2rxrootnum(vecteur & v,const vecteur & lv,gen & num,gen & tmpden,GIAC_CONTEXT)780 static bool sym2rxrootnum(vecteur & v,const vecteur & lv,gen & num,gen & tmpden,GIAC_CONTEXT){ 781 // make the polynomial X^d - num 782 // check for irreducibility 783 polynome p(poly12polynome(v)); 784 polynome p_content(p.dim); 785 factorization f; 786 gen extra_div=1; 787 if (!factor(p,p_content,f,true,false,false,1,extra_div)){ 788 #ifndef NO_STDEXCEPT 789 setsizeerr(gettext("Can't check irreducibility extracting rootof")); 790 #endif 791 return false; 792 } 793 // now we choose the factor of lowest degree of the factorization 794 int lowest_degree=int(v.size()),deg; 795 factorization::const_iterator f_it,f_itend=f.end(); 796 // Rewrite num if it's an ext with i, because it is rewritten by factor 797 if (num.type==_EXT && has_i(num)){ 798 gen b=1; 799 for (f_it=f.begin();f_it!=f_itend;++f_it){ 800 b=b*f_it->fact.coord.back().value; 801 } 802 if (b.type==_EXT){ 803 gen q=evalf(b,1,contextptr)/evalf(num,1,contextptr); 804 gen qr=_round(q,contextptr); 805 if (is_zero(q-qr,contextptr)){ 806 num=b/qr; 807 if (num.type==_FRAC){ 808 // commented since x^d=num is solved using irr_p, translated to b below 809 // tmpden=tmpden*num._FRACptr->den; 810 num=num._FRACptr->num; 811 } 812 } 813 } 814 } 815 for (f_it=f.begin();f_it!=f_itend;++f_it){ 816 polynome irr_p(f_it->fact); 817 deg=irr_p.lexsorted_degree(); 818 if (!deg) 819 continue; 820 if (deg==1){ 821 v=polynome2poly1(irr_p); 822 lowest_degree=1; 823 num=rdiv(-v.back(),v.front(),contextptr); 824 // cerr << "xroot" << num << "\n"; 825 gen numlv=r2sym(num,lv,contextptr); 826 if (!lvar(evalf(numlv,1,contextptr)).empty()) 827 *logptr(contextptr) << gettext("Warning, checking for positivity of a root depending of parameters might return wrong sign: ")<< numlv << "\n"; 828 if (is_positive(numlv,contextptr)) 829 break; 830 } 831 if (deg==2 && deg==lowest_degree){ 832 vecteur tmp=polynome2poly1(irr_p); 833 gen delta=tmp[1]*tmp[1]-4*tmp[0]*tmp[2]; 834 delta=r2sym(delta,lv,contextptr); 835 if (is_positive(delta,contextptr)) 836 v=tmp; 837 } 838 if (deg>=lowest_degree) 839 continue; 840 v=polynome2poly1(irr_p); 841 lowest_degree=deg; 842 } 843 if (lowest_degree>1){ 844 // here we must check that num is not an extension!! 845 if (num.type==_EXT){ 846 gen a=*(num._EXTptr+1),b; 847 gen a__VECTg; 848 if (a.type==_VECT) 849 a__VECTg=a; 850 else { 851 if( a.type!=_EXT || (a._EXTptr+1)->type!=_VECT){ 852 #ifndef NO_STDEXCEPT 853 setsizeerr(gettext("sym2poly.cc/sym2rxroot")); 854 #endif 855 return false; 856 } 857 a__VECTg=*(a._EXTptr+1); 858 } 859 int k; 860 gen new_v=common_minimal_POLY(a__VECTg,v,a,b,k,contextptr); 861 if (is_undef(new_v)) 862 return false; 863 *(num._EXTptr+1)=a; 864 if (b.type==_FRAC){ 865 num=b._FRACptr->num; 866 tmpden=tmpden*b._FRACptr->den; 867 } 868 else 869 num=b; 870 } 871 else { 872 vecteur w(2); 873 w[0]=plus_one; 874 num=algebraic_EXTension(w,v); 875 } 876 } 877 return true; 878 } 879 fxnd(const gen & e,gen & num,gen & den)880 void fxnd(const gen & e,gen & num, gen & den){ 881 if (e.type==_FRAC){ 882 num=e._FRACptr->num; 883 den=e._FRACptr->den; 884 } 885 else { 886 num=e; 887 den=plus_one; 888 } 889 } 890 rsqff(const polynome & p)891 static factorization rsqff(const polynome & p){ 892 polynome s(lgcd(p)); 893 factorization f(sqff(p/s)); 894 if (p.dim==0){ 895 f.push_back(facteur<polynome>(p,1)); 896 return f; 897 } 898 if (p.dim==1){ 899 // adjust const coeff 900 gen p1=1; 901 factorization::iterator it=f.begin(),itend=f.end(); 902 for (;it!=itend;++it){ 903 p1 = p1*pow(it->fact.coord.front().value,it->mult); 904 } 905 p1=p.coord.front().value/p1; 906 if (is_positive(-p1,context0)){ // ok 907 for (it=f.begin();it!=itend;++it){ 908 if (it->mult%2){ 909 it->fact=-it->fact; 910 p1=-p1; 911 break; 912 } 913 } 914 } 915 if (!is_one(p1)) 916 f.push_back(facteur<polynome>(polynome(p1,1),1)); 917 return f; 918 } 919 factorization ff(rsqff(s.trunc1())); 920 factorization::const_iterator it=ff.begin(),itend=ff.end(); 921 for (;it!=itend;++it) 922 f.push_back(facteur<polynome>(it->fact.untrunc1(),it->mult)); 923 return f; 924 } 925 smaller_factor(vecteur & v)926 void smaller_factor(vecteur & v){ 927 // factor v, select smaller factor 928 polynome vp(1),vp_content; 929 vp=poly12polynome(v,1); 930 gen tmp(1); lcmdeno(vp,tmp); vp=tmp*vp; 931 factorization vf; 932 gen extra_div=1; 933 if (factor(vp,vp_content,vf,true,false,false,1,extra_div) && vf.size()>1){ 934 polynome2poly1(vf.front().fact,1,v); 935 for (unsigned i=1;i<vf.size();++i){ 936 vecteur vi; 937 polynome2poly1(vf[i].fact,1,vi); 938 if (vi.size()<v.size()) 939 v=vi; 940 } 941 } 942 } 943 sym2rxroot(gen & num,gen & den,int n,int d,const vecteur & l,GIAC_CONTEXT)944 static bool sym2rxroot(gen & num,gen & den,int n,int d,const vecteur & l,GIAC_CONTEXT){ 945 if (is_zero(num)) 946 return true; 947 if (den.type==_CPLX){ 948 gen denr,deni,numr,numi; 949 reim(den,denr,deni,contextptr); 950 den=denr*denr+deni*deni; 951 num=num*(denr-cst_i*deni); 952 reim(num,numr,numi,contextptr); 953 deni=gcd(gcd(numr,denr,contextptr),numi,contextptr); 954 num=num/deni; 955 den=den/deni; 956 } 957 if (is_positive(r2e(-den,l,contextptr),contextptr)){ 958 num=-num; 959 den=-den; 960 } 961 bool sign_changed=false; 962 if (d<0){ 963 n=-n; 964 d=-d; 965 } 966 if (n<0){ 967 if (num.type==_EXT){ 968 gen temp(inv_EXT(num)*den); 969 fxnd(temp,num,den); 970 } 971 else { 972 #if defined RTOS_THREADX || defined BESTA_OS || defined USTL 973 gen g=num; num=den; den=g; 974 #else 975 swap(num,den); 976 #endif 977 } 978 num=pow(num,-n)*pow(den,n*(1-d)); 979 den=pow(den,-n); 980 } 981 else { 982 num=pow(num,n)*pow(den,n*(d-1)); 983 den=pow(den,n); 984 } 985 /* d==2 test makes normal(sqrt(1-x)) -> i*sqrt(x-1) 986 commented!, should common_minimal_poly detect [1,0...0,+/-same_poly]? 987 if ( (d%2 || d=2) && ( 988 (num.type==_POLY && is_positive(evalf_double(-num._POLYptr->coord.front().value,2,0),0)) 989 ) ){ 990 num=-num; 991 sign_changed=true; 992 } 993 */ 994 // compute number of cst polynomial in num and keep track of dims 995 int embeddings=0; 996 vector<int> embeddings_s; 997 if (is_atomic(num)){ 998 const_iterateur it=l.begin(),itend=l.end(); 999 embeddings=int(itend-it); 1000 if ( 1001 0 && // disable for expressions with mixed rootof/fracpow? 1002 embeddings==1 && it->type==_VECT && it->_VECTptr->empty()) 1003 embeddings=0; 1004 else { 1005 for (int j=0;j<embeddings;++it,++j){ 1006 if (it->type!=_VECT){ 1007 string s("sym2rxroot error num="+num.print(contextptr)+" den="+den.print(contextptr)+" l="+gen(l).print(contextptr)); 1008 CERR << s << "\n"; 1009 #ifndef NO_STDEXCEPT 1010 setsizeerr(s); 1011 #endif 1012 return false; 1013 } 1014 embeddings_s.push_back(int(it->_VECTptr->size())); 1015 } 1016 } 1017 } 1018 else { 1019 for (;((num.type==_POLY) && (Tis_constant<gen>(*num._POLYptr)));++embeddings){ 1020 embeddings_s.push_back(num._POLYptr->dim); 1021 num=num._POLYptr->coord.front().value; 1022 } 1023 } 1024 if (num.type==_FRAC){ 1025 den=den*num._FRACptr->den; 1026 num=num._FRACptr->num*pow(num._FRACptr->den,d-1); 1027 } 1028 vecteur v(d+1,zero); 1029 gen tmpden=1; 1030 if (num.type==_POLY){ 1031 const factorization f=rsqff(*num._POLYptr); 1032 polynome S(plus_one,num._POLYptr->dim),D(plus_one,S.dim); 1033 factorization::const_iterator jt=f.begin(),jtend=f.end(); 1034 for (;jt!=jtend;++jt){ 1035 if (jt->mult%d) 1036 S=S*pow(jt->fact,jt->mult % d); 1037 D=D*pow(jt->fact,jt->mult/d); 1038 } 1039 // Check sign of D 1040 vecteur Dl(l); 1041 if (embeddings && Dl[embeddings].type==_VECT) 1042 Dl=*Dl[embeddings]._VECTptr; 1043 if (is_positive(r2e(-D,Dl,contextptr),contextptr)){ 1044 D=-D; 1045 if (d%2) 1046 S=-S; 1047 } 1048 gen cont=ppz(S); 1049 gen simpl,doubl; bool pos; 1050 zint2simpldoublpos(cont,simpl,doubl,pos,d,contextptr); 1051 if (!pos) simpl=-simpl; 1052 num=simpl*S; 1053 D=doubl*D; 1054 v.front()=plus_one; 1055 v.back()=-num; 1056 sym2rxrootnum(v,vecteur(l.begin()+embeddings,l.end()),num,tmpden,contextptr); 1057 if (num.type==_EXT && num._EXTptr->type==_VECT){ 1058 num=algebraic_EXTension(multvecteur(D,*(*num._EXTptr)._VECTptr),*(num._EXTptr+1)); 1059 } 1060 else { 1061 if (num.type==_POLY) 1062 num=D**num._POLYptr; 1063 else 1064 num=D*num; 1065 } 1066 } 1067 else { 1068 v.front()=plus_one; 1069 v.back()=-num; 1070 if (is_integer(num)){ 1071 gen simpl,doubl; 1072 bool pos; 1073 zint2simpldoublpos(num,simpl,doubl,pos,d,contextptr); 1074 v.back()=pos?-simpl:simpl; 1075 smaller_factor(v); 1076 if (is_minus_one(v.back())) 1077 num=doubl; 1078 else 1079 num=algebraic_EXTension(makevecteur(doubl,0),v); 1080 } 1081 else { 1082 bool neg=false; 1083 #if 1 1084 gen numd; 1085 if (has_evalf(num,numd,1,contextptr) && is_positive(-numd,contextptr)) 1086 neg=true; 1087 if (neg){ 1088 num=-num; 1089 v.back()=-num; 1090 } 1091 #endif 1092 sym2rxrootnum(v,vecteur(l.begin()+embeddings,l.end()),num,tmpden,contextptr); 1093 if (neg) 1094 num=cst_i*num; 1095 } 1096 } 1097 // end else num integer 1098 // and eventually we embedd num 1099 for (;embeddings;--embeddings){ 1100 num=polynome(num,embeddings_s.back()); 1101 tmpden=polynome(tmpden,embeddings_s.back()); 1102 embeddings_s.pop_back(); 1103 } 1104 if (sign_changed){ 1105 if (d==2) 1106 num=cst_i*num; 1107 else 1108 num=-num; 1109 } 1110 den=den*tmpden; 1111 return true; 1112 } 1113 sym2r(const symbolic & s,const gen & iext,const vecteur & l,const vecteur & lv,const vecteur & lvnum,const vecteur & lvden,int l_size,gen & num,gen & den,GIAC_CONTEXT)1114 static bool sym2r (const symbolic &s,const gen & iext,const vecteur &l, const vecteur & lv, const vecteur & lvnum,const vecteur & lvden, int l_size,gen & num,gen & den,GIAC_CONTEXT){ 1115 if (s.sommet==at_plus){ 1116 if (s.feuille.type!=_VECT){ 1117 return sym2r(s.feuille,iext,l,lv,lvnum,lvden,l_size,num,den,contextptr); 1118 } 1119 vecteur::iterator debut=s.feuille._VECTptr->begin(); 1120 vecteur::iterator fin=s.feuille._VECTptr->end(); 1121 bool sequentially=has_op(s.feuille,*at_inv) && (fin-debut<512); 1122 return sym2radd(debut,fin,iext,l,lv,lvnum,lvden,l_size,num,den,contextptr,sequentially); 1123 } 1124 if (s.sommet==at_prod){ 1125 vecteur::iterator debut=s.feuille._VECTptr->begin(); 1126 vecteur::iterator fin=s.feuille._VECTptr->end(); 1127 bool res=sym2rmul(debut,fin,iext,l,lv,lvnum,lvden,l_size,num,den,contextptr); 1128 /* 1129 gen numtmp,dentmp; 1130 sym2rmulold(debut,fin,l,lv,lvnum,lvden,l_size,numtmp,dentmp,contextptr); 1131 if (numtmp!=num || dentmp!=den) 1132 return false; 1133 */ 1134 return res; 1135 } 1136 if (s.sommet==at_neg){ 1137 bool totally_converted=sym2r(s.feuille,iext,l,lv,lvnum,lvden,l_size,num,den,contextptr); 1138 num=-num; 1139 return totally_converted; 1140 } 1141 if (s.sommet==at_inv){ 1142 bool totally_converted=sym2r(s.feuille,iext,l,lv,lvnum,lvden,l_size,num,den,contextptr); 1143 if (is_zero(num)){ 1144 num=unsigned_inf; 1145 den=plus_one; 1146 return true; 1147 } 1148 if (num.type==_EXT){ 1149 gen temp(inv_EXT(num)*den); 1150 fxnd(temp,num,den); 1151 } 1152 else { 1153 if ( (num.type==_POLY) && (num._POLYptr->dim==0) && (!num._POLYptr->coord.empty()) && (num._POLYptr->coord.front().value.type==_EXT) ){ 1154 gen temp(inv_EXT(num._POLYptr->coord.front().value)); 1155 if ( (den.type==_POLY) && (den._POLYptr->dim==0) && (!den._POLYptr->coord.empty()) ) 1156 temp=temp*den._POLYptr->coord.front().value; 1157 else 1158 temp=temp*den; 1159 gen tempnum,tempden; 1160 fxnd(temp,tempnum,tempden); 1161 polynome tmpnum(0),tmpden(0); 1162 tmpnum.coord.push_back(monomial<gen>(tempnum,index_t(0))); 1163 tmpden.coord.push_back(monomial<gen>(tempden,index_t(0))); 1164 num=tmpnum; 1165 if (tempden.type==_POLY) 1166 den=tmpden; 1167 else 1168 den=tempden; 1169 } 1170 else { 1171 #if defined RTOS_THREADX || defined BESTA_OS || defined USTL 1172 gen g=num; num=den; den=g; 1173 #else 1174 swap(num,den); 1175 #endif 1176 } 1177 } 1178 return totally_converted; 1179 } 1180 if (s.sommet==at_pow){ 1181 if ((*s.feuille._VECTptr)[1].type==_INT_) { 1182 bool totally_converted=sym2r(s.feuille._VECTptr->front(),iext,l,lv,lvnum,lvden,l_size,num,den,contextptr); 1183 int n=(*s.feuille._VECTptr)[1].val; 1184 if (absint(n)>=(1<<15)) 1185 num=gensizeerr(gettext("Polynomial exponent overflow.")); 1186 if (n<0){ 1187 if (num.type==_EXT){ 1188 gen temp(inv_EXT(num)*den); 1189 fxnd(temp,num,den); 1190 } 1191 else { 1192 #if defined RTOS_THREADX || defined BESTA_OS || defined USTL 1193 gen g=num; num=den; den=g; 1194 #else 1195 swap(num,den); 1196 #endif 1197 } 1198 num=pow(num,-n); 1199 den=pow(den,-n); 1200 } 1201 else { 1202 num=pow(num,n); 1203 den=pow(den,n); 1204 } 1205 if (num.type==_EXT) 1206 simplify3(num,den); 1207 // FIXME: embeddings like for rootof 1208 return totally_converted; 1209 } 1210 if ((*s.feuille._VECTptr)[1].type==_FRAC) { 1211 fraction f=*((*s.feuille._VECTptr)[1]._FRACptr); 1212 if ( (f.num.type==_INT_) && (f.den.type==_INT_)){ 1213 bool totally_converted=sym2r(s.feuille._VECTptr->front(),iext,l,lv,lvnum,lvden,l_size,num,den,contextptr); 1214 if (!sym2rxroot(num,den,f.num.val,f.den.val,l,contextptr)) 1215 return false; 1216 return totally_converted; 1217 } 1218 } 1219 } 1220 if (s.sommet==at_rootof){ 1221 bool totally_converted=sym2r(s.feuille._VECTptr->front(),iext,l,lv,lvnum,lvden,l_size,num,den,contextptr); 1222 gen pmin_num,pmin_den; 1223 totally_converted=totally_converted && sym2r(s.feuille._VECTptr->back(),iext,l,lv,lvnum,lvden,l_size,pmin_num,pmin_den,contextptr); 1224 if (!is_one(pmin_den)){ 1225 #ifndef NO_STDEXCEPT 1226 setsizeerr(gettext("Minimal poly. in rootof must be fraction free")); 1227 #endif 1228 return false; 1229 } 1230 int embeddings=0; 1231 vector<int> embeddings_s; 1232 // put out constant polynomials 1233 if (num.type!=_VECT || pmin_num.type!=_VECT) 1234 return totally_converted; 1235 vecteur vnum=*num._VECTptr,vpmin=*pmin_num._VECTptr; 1236 int s=int(l.size()); 1237 for (;embeddings<s;){ 1238 bool exitmainloop=false; 1239 iterateur it=vnum.begin(),itend=vnum.end(); 1240 int i=0; 1241 for (;;++it){ 1242 if (it==itend){ 1243 if (i==0){ 1244 ++i; 1245 it=vpmin.begin(); 1246 itend=vpmin.end(); 1247 } 1248 else 1249 break; 1250 } 1251 if (is_atomic(*it)) 1252 continue; 1253 if (it->type==_POLY && Tis_constant<gen>(*it->_POLYptr)) 1254 continue; 1255 exitmainloop=true; 1256 break; 1257 } 1258 if (exitmainloop) 1259 break; 1260 if (l[embeddings].type!=_VECT){ 1261 #ifndef NO_STDEXCEPT 1262 setsizeerr(gettext("sym2poly.cc/bool sym2r (const")); 1263 #endif 1264 return false; 1265 } 1266 embeddings_s.push_back(int(l[embeddings]._VECTptr->size())); 1267 ++embeddings; 1268 for (it=vnum.begin(),itend=vnum.end(),i=0;;++it){ 1269 if (it==itend){ 1270 if (i==0){ 1271 ++i; 1272 it=vpmin.begin(); 1273 itend=vpmin.end(); 1274 } 1275 else 1276 break; 1277 } 1278 if (it->type==_POLY) 1279 *it=it->_POLYptr->coord.front().value; 1280 } 1281 } 1282 num=algebraic_EXTension(vnum,vpmin); 1283 for (;embeddings;--embeddings){ 1284 num=polynome(num,embeddings_s.back()); 1285 embeddings_s.pop_back(); 1286 } 1287 return totally_converted; 1288 } 1289 num=s; 1290 den=plus_one; 1291 return false; 1292 } 1293 sym2r(const fraction & f,const gen & iext,const vecteur & l,const vecteur & lv,const vecteur & lvnum,const vecteur & lvden,int l_size,gen & num,gen & den,GIAC_CONTEXT)1294 static bool sym2r (const fraction &f,const gen &iext,const vecteur &l, const vecteur & lv, const vecteur & lvnum,const vecteur & lvden, int l_size,gen & num,gen & den,GIAC_CONTEXT){ 1295 gen dent,numt; 1296 bool totally_converted=sym2r(f.num,iext,l,lv,lvnum,lvden,l_size,num,dent,contextptr); 1297 totally_converted=totally_converted && sym2r(f.den,iext,l,lv,lvnum,lvden,l_size,den,numt,contextptr); 1298 num=num*numt; 1299 den=den*dent; 1300 return totally_converted; 1301 } 1302 sym2r(const vecteur & v,const gen & iext,const vecteur & l,const vecteur & lv,const vecteur & lvnum,const vecteur & lvden,int l_size,gen & num,gen & den,GIAC_CONTEXT)1303 bool sym2r(const vecteur &v,const gen & iext,const vecteur &l, const vecteur & lv, const vecteur & lvnum,const vecteur & lvden, int l_size,gen & num,gen & den,GIAC_CONTEXT){ 1304 den=plus_one; 1305 if (v.empty()){ 1306 num=zero; 1307 return true; 1308 } 1309 bool totally_converted=true; 1310 gen lcmdeno=plus_one; 1311 const_iterateur it=v.begin(),itend=v.end(); 1312 vecteur res,numv; 1313 res.reserve(2*(itend-it)); 1314 numv.reserve(itend-it); 1315 for (;it!=itend;++it){ 1316 totally_converted=totally_converted && sym2r(*it,iext,l,lv,lvnum,lvden,l_size,num,den,contextptr); 1317 lcmdeno = lcm(lcmdeno,den); 1318 res.push_back(num); 1319 res.push_back(den); 1320 } 1321 for (it=res.begin(),itend=res.end();it!=itend;){ 1322 num=*it; 1323 ++it; 1324 den=*it; 1325 ++it; 1326 numv.push_back(num*rdiv(lcmdeno,den,contextptr)); 1327 } 1328 den=lcmdeno; 1329 num=numv; 1330 return totally_converted; 1331 } 1332 sym2r(const vecteur & v,const vecteur & l,const vecteur & lv,const vecteur & lvnum,const vecteur & lvden,int l_size,gen & num,gen & den,GIAC_CONTEXT)1333 bool sym2r(const vecteur &v,const vecteur &l, const vecteur & lv, const vecteur & lvnum,const vecteur & lvden, int l_size,gen & num,gen & den,GIAC_CONTEXT){ 1334 return sym2r(v,1,l,lv,lvnum,lvden,l_size,num,den,contextptr); 1335 } 1336 sym2rmod(const gen * gptr,const gen & iext,const vecteur & l,const vecteur & lv,const vecteur & lvnum,const vecteur & lvden,int l_size,gen & num,gen & den,GIAC_CONTEXT)1337 static bool sym2rmod (const gen * gptr,const gen & iext,const vecteur &l, const vecteur & lv, const vecteur & lvnum,const vecteur & lvden, int l_size,gen & num,gen & den,GIAC_CONTEXT){ 1338 gen modnum,modden,modulo; 1339 bool totally_converted=sym2r(*(gptr+1),iext,l,lv,lvnum,lvden,l_size,modnum,modden,contextptr); 1340 modulo=rdiv(modnum,modden,contextptr); 1341 totally_converted=totally_converted && sym2r(*gptr,iext,l,lv,lvnum,lvden,l_size,num,den,contextptr); 1342 num=makemod(num,modulo); 1343 den=makemod(den,modulo); 1344 return totally_converted; 1345 } 1346 1347 // iext is an additional parameter (added 5 sept 2014) 1348 // value is 0 if i=sqrt(-1) is not in the algebraic number field from lvnum 1349 // otherwise it is an extension or fraction extension/integer. sym2r(const gen & e,const gen & iext,const vecteur & l,const vecteur & lv,const vecteur & lvnum,const vecteur & lvden,int l_size,gen & num,gen & den,GIAC_CONTEXT)1350 bool sym2r (const gen &e,const gen & iext,const vecteur &l, const vecteur & lv, const vecteur & lvnum,const vecteur & lvden, int l_size,gen & num,gen & den,GIAC_CONTEXT){ 1351 int n=0; 1352 switch (e.type ) { 1353 case _INT_: case _DOUBLE_: case _ZINT: case _REAL: case _FLOAT_: 1354 num=e; 1355 den=plus_one; 1356 return true; 1357 case _CPLX: 1358 if (iext==0){ 1359 num=e; 1360 den=plus_one; 1361 return true; 1362 } 1363 if (iext.type!=_FRAC){ 1364 num=*e._CPLXptr+iext*(*(e._CPLXptr+1)); 1365 den=plus_one; 1366 return true; 1367 } 1368 num=*(e._CPLXptr+1); 1369 den=iext._FRACptr->den; 1370 simplify3(num,den); 1371 num=*e._CPLXptr*den+iext._FRACptr->num*num; 1372 return true; 1373 case _IDNT: case _SYMB: 1374 if (e.is_symb_of_sommet(at_rootof) && e._SYMBptr->feuille.type==_VECT && e._SYMBptr->feuille._VECTptr->size()==2){ 1375 gen r=e._SYMBptr->feuille._VECTptr->back(); 1376 for (int i=0;i<lv.size() && i<lvnum.size();++i){ 1377 if (lv[i].is_symb_of_sommet(at_rootof) && lv[i]._SYMBptr->feuille.type==_VECT && lv[i]._SYMBptr->feuille._VECTptr->size()==2){ 1378 gen gl=lv[i]._SYMBptr->feuille._VECTptr->front(); 1379 if (gl.type==_VECT && gl._VECTptr->size()==2 && gl._VECTptr->front()==1 && gl._VECTptr->back()==0){ 1380 gl=lv[i]._SYMBptr->feuille._VECTptr->back(); 1381 if (r.type==_VECT && gl.type==_VECT && *gl._VECTptr==*r._VECTptr 1382 //&& lidnt(gl).empty() 1383 ){ 1384 n=i+1; // break; 1385 r=e._SYMBptr->feuille._VECTptr->front(); 1386 sym2r(r,iext,l,lv,lvnum,lvden,l_size,num,den,contextptr); 1387 r=num; 1388 if (is_one(den) && r.type==_VECT){ 1389 gen N=lvnum[n-1]; 1390 gen D=lvden[n-1]; 1391 hornerfrac(*r._VECTptr,N,D,num,den); 1392 return true; 1393 } 1394 } 1395 } 1396 } 1397 } 1398 if (n && n<=lvnum.size() && e._SYMBptr->feuille._VECTptr->front().type==_VECT){ 1399 gen N=lvnum[n-1]; 1400 gen D=lvden[n-1]; 1401 hornerfrac(*e._SYMBptr->feuille._VECTptr->front()._VECTptr,N,D,num,den); 1402 return true; 1403 } 1404 } 1405 n=equalposcomp(lv,e); 1406 if (n && (unsigned(n)<=lvnum.size())){ 1407 num=lvnum[n-1]; 1408 den=lvden[n-1]; 1409 return true; 1410 } 1411 if ((!l.empty()) && (l.front().type==_VECT) ){ 1412 int i,j; 1413 if (equalposmat(l,e,i,j)){ 1414 num=monomial2gen(monomial<gen>(gen(1),j+1,int(l[i]._VECTptr->size()))); 1415 for (int k=i-1;k>=0;--k) 1416 num=monomial2gen(monomial<gen>(num,int(l[k]._VECTptr->size()))); 1417 den=plus_one; 1418 return true; 1419 } 1420 } 1421 else { 1422 n=equalposcomp(l,e); 1423 if (n){ 1424 num=monomial2gen(monomial<gen>(gen(1),n,l_size)); 1425 den=plus_one; 1426 return true; 1427 } 1428 } 1429 if (e.type!=_SYMB){ 1430 num=e; 1431 den=plus_one; 1432 return true; 1433 } 1434 return sym2r(*e._SYMBptr,iext,l,lv,lvnum,lvden,l_size,num,den,contextptr); 1435 case _FRAC: 1436 return sym2r(*e._FRACptr,iext,l,lv,lvnum,lvden,l_size,num,den,contextptr); 1437 case _VECT: 1438 return sym2r(*e._VECTptr,iext,l,lv,lvnum,lvden,l_size,num,den,contextptr); 1439 case _POLY: case _EXT: 1440 if ((!l.empty()) && (l.front().type==_VECT) ){ 1441 num=e; 1442 for (int k=int(l.size())-1;k>=0;--k) // was l_size 1443 num=monomial2gen(monomial<gen>(num,int(l[k]._VECTptr->size()))); 1444 den=plus_one; 1445 } 1446 else { 1447 num=monomial2gen(monomial<gen>(e,l_size)); 1448 den=plus_one; 1449 } 1450 return true; 1451 case _MOD: 1452 return sym2rmod(e._MODptr,iext,l,lv,lvnum,lvden,l_size,num,den,contextptr); 1453 case _USER: 1454 num=e; 1455 den=plus_one; 1456 return true; 1457 default: 1458 #ifndef NO_STDEXCEPT 1459 settypeerr(gettext("Unable to handle type [sym2poly.cc]")+e.print(contextptr)); 1460 #endif 1461 num=gensizeerr(gettext("Unable to handle ")+e.print(contextptr)); 1462 den=plus_one; 1463 return 0; 1464 } 1465 return 0; 1466 } 1467 extract_ext(const gen & g,gen & extracted,vector<int> & embeddings)1468 static void extract_ext(const gen & g,gen & extracted,vector<int> & embeddings){ 1469 if (g.type==_EXT){ 1470 extracted= ext_reduce(g); 1471 return; 1472 } 1473 if (g.type==_POLY && !g._POLYptr->coord.empty() && Tis_constant<gen>(*g._POLYptr)){ 1474 embeddings.push_back(g._POLYptr->dim); 1475 extract_ext(g._POLYptr->coord.front().value,extracted,embeddings); 1476 } 1477 else 1478 extracted=0; 1479 } 1480 check_ext(const gen & oldext,const gen & curext)1481 static bool check_ext(const gen & oldext,const gen & curext){ 1482 if (oldext.type!=_VECT) 1483 return false; 1484 if (curext.type==_EXT) 1485 return oldext!=*(curext._EXTptr+1); 1486 if (curext.type==_FRAC) 1487 return check_ext(oldext,curext._FRACptr->num); 1488 return false; 1489 } 1490 1491 // rewrite extension inside a in terms of a common rootof 1492 // extensions in a must be at the top level 1493 // alg_extin is the current list of already translated algebraic extension 1494 // alg_extout is the current list of corresponding values 1495 // w.r.t the common minpoly 1496 // return a in terms of the common ext common_algext(const vecteur & a,const vecteur & l,GIAC_CONTEXT)1497 static vecteur common_algext(const vecteur & a,const vecteur & l,GIAC_CONTEXT){ 1498 const_iterateur it,itend=a.end(); 1499 vecteur alg_extin,alg_extoutnum,alg_extoutden; 1500 #if 0 // ndef GIAC_HAS_STO38 1501 // trying with rational univariate rep. 1502 for (it=a.begin();it!=itend;++it){ 1503 if (it->type==_EXT){ 1504 gen minpoly=*(it->_EXTptr+1); 1505 if (minpoly.type!=_VECT) 1506 break; 1507 unsigned i=0,s=minpoly._VECTptr->size(); 1508 for (;i<s;++i){ 1509 if (!is_integer((*minpoly._VECTptr)[i])) 1510 break; 1511 } 1512 if (i!=s) 1513 break; 1514 alg_extin.push_back(minpoly); 1515 } 1516 } 1517 if (it==itend){ 1518 // extract system of equations from alg_extin and call gbasis 1519 // problem: for e.g. sqrt(2),sqrt(3),sqrt(6) it will return a poly of deg. 8 1520 } 1521 alg_extin.clear(); alg_extoutnum.clear(); alg_extoutden.clear(); 1522 #endif 1523 for (it=a.begin();it!=itend;++it){ 1524 if (it->type==_EXT) 1525 break; 1526 } 1527 if (itend==it) 1528 return a; 1529 alg_extin.push_back(*(it->_EXTptr+1)); 1530 gen curext=*(it->_EXTptr+1); 1531 gen minpoly=curext; 1532 alg_extoutnum.push_back(makevecteur(1,0)); 1533 alg_extoutden.push_back(1); 1534 ++it; 1535 for (;it!=itend;++it){ 1536 if (it->type!=_EXT) 1537 continue; 1538 curext=*(it->_EXTptr+1); 1539 if (equalposcomp(alg_extin,curext)) 1540 continue; 1541 // make common extension of curext and minpoly 1542 gen oldminpoly=minpoly; 1543 gen oldcurext=curext; 1544 if (curext==oldminpoly){ 1545 alg_extin.push_back(oldcurext); 1546 alg_extoutnum.push_back(makevecteur(1,0)); 1547 alg_extoutden.push_back(1); 1548 continue; 1549 } 1550 else 1551 minpoly=common_EXT(curext,oldminpoly,&l,contextptr); 1552 if (minpoly.type!=_VECT) 1553 return vecteur(1,gensizeerr(contextptr)); 1554 if (minpoly._VECTptr->size()>unsigned(MAX_COMMON_ALG_EXT_ORDER_SIZE)) 1555 return vecteur(1,undef); 1556 // compute alg_extoutnum/den using newminpoly 1557 int s=int(alg_extin.size()); 1558 for (int i=0;i<s;++i){ 1559 if (alg_extoutnum[i].type!=_VECT) 1560 return vecteur(1,gensizeerr(contextptr)); 1561 if (oldminpoly.type==_FRAC){ 1562 gen resn,resd; 1563 hornerfrac(*alg_extoutnum[i]._VECTptr,oldminpoly._FRACptr->num,oldminpoly._FRACptr->den,resn,resd); 1564 alg_extoutden[i]=alg_extoutden[i]*resd; 1565 if (resn.type!=_EXT || *(resn._EXTptr+1)!=minpoly) 1566 return vecteur(1,gensizeerr(contextptr)); 1567 alg_extoutnum[i]=*resn._EXTptr; 1568 } 1569 else { 1570 gen tmp=horner(*alg_extoutnum[i]._VECTptr,oldminpoly); 1571 if (tmp.type!=_EXT || *(tmp._EXTptr+1)!=minpoly) 1572 return vecteur(1,gensizeerr(contextptr)); 1573 alg_extoutnum[i]=*tmp._EXTptr; 1574 } 1575 } 1576 // add oldcurext/curext to the list 1577 alg_extin.push_back(oldcurext); 1578 if (curext.type==_FRAC){ 1579 alg_extoutden.push_back(curext._FRACptr->den); 1580 curext=curext._FRACptr->num; 1581 } 1582 else 1583 alg_extoutden.push_back(1); 1584 if (curext.type!=_EXT || *(curext._EXTptr+1)!=minpoly) 1585 return vecteur(1,gensizeerr(contextptr)); 1586 alg_extoutnum.push_back(*curext._EXTptr); 1587 } 1588 vecteur res; 1589 for (it=a.begin();it!=itend;++it){ 1590 if (it->type!=_EXT) 1591 res.push_back(*it); 1592 int n=equalposcomp(alg_extin,*(it->_EXTptr+1)); 1593 if (!n) 1594 return vecteur(1,gensizeerr(contextptr)); 1595 --n; 1596 gen tmp=alg_extoutnum[n]; 1597 if (tmp.type==_VECT) 1598 tmp.subtype=_POLY1__VECT; 1599 if (!is_one(alg_extoutden[n])){ 1600 gen resn,resd; 1601 hornerfrac(*it->_EXTptr->_VECTptr,tmp,alg_extoutden[n],resn,resd); 1602 if (resn.type==_VECT && minpoly.type==_VECT && resn._VECTptr->size()>=minpoly._VECTptr->size()){ 1603 resn = *resn._VECTptr % *minpoly._VECTptr; 1604 gen tmp; 1605 lcmdeno_converted(*resn._VECTptr,tmp,contextptr); 1606 resd=resd*tmp; 1607 } 1608 res.push_back(algebraic_EXTension(resn,minpoly)/resd); 1609 } 1610 else { 1611 tmp=horner(*it->_EXTptr,tmp); 1612 if (tmp.type!=_VECT) 1613 res.push_back(tmp); 1614 else 1615 res.push_back(algebraic_EXTension(tmp,minpoly)); 1616 } 1617 } 1618 return res; 1619 } 1620 compute_lv_lvnum_lvden(const vecteur & l,vecteur & lv,vecteur & lvnum,vecteur & lvden,bool & totally_converted,int l_size,GIAC_CONTEXT)1621 static bool compute_lv_lvnum_lvden(const vecteur & l,vecteur & lv,vecteur & lvnum,vecteur & lvden,bool & totally_converted,int l_size,GIAC_CONTEXT){ 1622 gen num,den; 1623 // sort by ascending complexity 1624 iterateur it=lv.begin(),itend=lv.end(); 1625 lvnum.reserve(itend-it); 1626 lvden.reserve(itend-it); 1627 gen_sort_f(it,itend,symb_size_less); 1628 bool algext=false; 1629 // find num/den for each var of lv 1630 for (;it!=itend;++it){ 1631 algext = algext || (it->type==_SYMB && (it->_SYMBptr->sommet==at_pow || it->_SYMBptr->sommet==at_rootof)) ; 1632 totally_converted = totally_converted && sym2r(*it,0,l,lv,lvnum,lvden,l_size,num,den,contextptr); 1633 lvnum.push_back(num); 1634 lvden.push_back(den); 1635 } 1636 if (!algext || lv.size()==1 || l==lv) 1637 return true; 1638 it=lv.begin(); 1639 // create common extensions for rootofs 1640 // first find extensions and embeddings levels 1641 std::map< int,vecteur> extensions,newext; 1642 for (int i=0;it!=itend;++it,++i){ 1643 if (it->is_symb_of_sommet(at_pow) || it->is_symb_of_sommet(at_rootof)){ 1644 gen ext; 1645 vector<int> emb; 1646 extract_ext(lvnum[i],ext,emb); 1647 if (ext.type==_FRAC){ 1648 lvnum[i]=lvnum[i]*ext._FRACptr->den; 1649 lvden[i]=lvden[i]*ext._FRACptr->den; 1650 ext=ext._FRACptr->num; 1651 } 1652 if (!is_zero(ext)) 1653 extensions[int(emb.size())].push_back(ext); 1654 } 1655 } 1656 // now rewrite each element of the list at each embedding level 1657 std::map< int,vecteur >::iterator jt=extensions.begin(),jtend=extensions.end(); 1658 for (;jt!=jtend;++jt){ 1659 if (is_undef( (newext[jt->first]=common_algext(jt->second,l,contextptr)) )) 1660 return false; 1661 } 1662 //stores cur 1663 it=lv.begin(); 1664 for (int i=0;it!=itend;++it,++i){ 1665 if (it->is_symb_of_sommet(at_pow) || it->is_symb_of_sommet(at_rootof)){ 1666 gen ext; 1667 vector<int> emb; 1668 extract_ext(lvnum[i],ext,emb); 1669 if (is_zero(ext)) 1670 continue; 1671 int n=equalposcomp(extensions[int(emb.size())],ext); 1672 if (!n) 1673 return false; // setsizeerr(); 1674 gen tmp=newext[int(emb.size())],den=1; 1675 if (tmp.type!=_VECT || int(tmp._VECTptr->size())<n) 1676 return false; // setsizeerr(); 1677 tmp=(*tmp._VECTptr)[n-1]; 1678 if (tmp.type==_FRAC){ 1679 den=tmp._FRACptr->den; 1680 tmp=tmp._FRACptr->num; 1681 } 1682 for (;!emb.empty();){ 1683 tmp=polynome(tmp,emb.back()); 1684 den=polynome(den,emb.back()); 1685 emb.pop_back(); 1686 } 1687 lvden[i]=lvden[i]*den; 1688 lvnum[i]=tmp; 1689 } 1690 } 1691 return true; 1692 } 1693 in_find_extension(const gen & extension,vecteur & already,vector<int> & embed,gen & iext,GIAC_CONTEXT)1694 bool in_find_extension(const gen & extension,vecteur & already,vector<int> & embed,gen & iext,GIAC_CONTEXT){ 1695 iext=0; 1696 if (extension.type==_POLY && !extension._POLYptr->coord.empty()){ 1697 const polynome & p=*extension._POLYptr; 1698 embed.push_back(p.dim); 1699 for (unsigned i=0;i<p.coord.size();++i){ 1700 if (in_find_extension(p.coord[i].value,already,embed,iext,contextptr)) 1701 return true; 1702 } 1703 embed.pop_back(); 1704 return false; 1705 } 1706 if (extension.type!=_EXT) 1707 return false; 1708 gen Extension=*(extension._EXTptr+1); 1709 if (Extension.type!=_VECT || equalposcomp(already,Extension)) 1710 return false; 1711 bool done=false; 1712 for (unsigned i=0;i<Extension._VECTptr->size();++i){ 1713 gen c=(*Extension._VECTptr)[i]; 1714 if (!is_integer(c)){ 1715 done=true; 1716 // recurse 1717 if (in_find_extension(c,already,embed,iext,contextptr)) 1718 return true; 1719 } 1720 } 1721 if (done) 1722 return false; 1723 iext=makevecteur(1,0,1); 1724 gen currentext=Extension; 1725 common_EXT(iext,currentext,0,contextptr); 1726 if (currentext.type==_EXT) 1727 currentext=*(currentext._EXTptr+1); 1728 Extension=change_subtype(Extension,_POLY1__VECT); 1729 currentext=change_subtype(currentext,_POLY1__VECT); 1730 if (Extension==currentext) 1731 return true; 1732 already.push_back(Extension); 1733 return false; 1734 } 1735 1736 // remove embedded i in n clean_iext(gen & n,gen & d,const gen & iext,GIAC_CONTEXT)1737 bool clean_iext(gen & n,gen & d,const gen & iext,GIAC_CONTEXT){ 1738 if (iext==0) return true; 1739 if (n.type==_POLY){ 1740 polynome p=*n._POLYptr; 1741 gen iext_=iext; 1742 if (iext.type==_FRAC){ 1743 if (iext._FRACptr->num.type==_POLY){ 1744 if (iext._FRACptr->num._POLYptr->dim!=0 || iext._FRACptr->num._POLYptr->coord.empty()) 1745 return false; 1746 iext_=iext._FRACptr->num._POLYptr->coord.front().value/iext._FRACptr->den; 1747 } 1748 } 1749 else { 1750 if (iext.type==_POLY){ 1751 if (iext._POLYptr->dim!=0 || iext._POLYptr->coord.empty()) 1752 return false; 1753 iext_=iext._POLYptr->coord.front().value; 1754 } 1755 } 1756 unsigned i=0; 1757 for (i=0;i<p.coord.size();++i){ 1758 gen D=d; 1759 clean_iext(p.coord[i].value,D,iext_,contextptr); 1760 p.coord[i].value=p.coord[i].value/D; 1761 } 1762 d=1, 1763 lcmdeno(p,d); 1764 n=d*p; 1765 return true; 1766 } 1767 if (n.type==_EXT){ 1768 gen n1=*n._EXTptr,n2=*(n._EXTptr+1); 1769 if (n1.type==_VECT){ 1770 vecteur nv=*n1._VECTptr; 1771 if (has_i(nv)){ 1772 gen r=algebraic_EXTension(makevecteur(1,0),n2),R,I,res=0; 1773 for (unsigned i=0;i<nv.size();++i){ 1774 res = res*r; 1775 reim(nv[i],R,I,contextptr); 1776 res += R+I*iext; 1777 } 1778 if (res.type==_FRAC){ 1779 d=d*res._FRACptr->den; 1780 res=res._FRACptr->num; 1781 } 1782 n=res; 1783 } 1784 } 1785 return true; 1786 } 1787 if (n.type==_CPLX){ 1788 n=(*n._CPLXptr)+(*(n._CPLXptr+1))*iext; 1789 gen N,D; 1790 fxnd(n,N,D); 1791 n=N; 1792 d=d*D; 1793 return true; 1794 } 1795 return true; 1796 } 1797 clean_iext(vecteur & lvnum,vecteur & lvden,const gen & iext,GIAC_CONTEXT)1798 bool clean_iext(vecteur & lvnum,vecteur & lvden,const gen & iext,GIAC_CONTEXT){ 1799 if (iext==0) return true; 1800 for (unsigned i=0;i<lvnum.size();++i){ 1801 if (!clean_iext(lvnum[i],lvden[i],iext,contextptr)) 1802 return false; 1803 } 1804 return true; 1805 } 1806 find_iext(const gen & e,const vecteur & lvnum,GIAC_CONTEXT)1807 gen find_iext(const gen & e,const vecteur & lvnum,GIAC_CONTEXT){ 1808 gen iext(0); 1809 if (has_i(e)){ 1810 // check if i is inside algebraic extensions from lvnum 1811 vecteur already; unsigned k; 1812 for (k=0;k<lvnum.size();++k){ 1813 gen extension=lvnum[k]; 1814 vector<int> embed; 1815 if (in_find_extension(extension,already,embed,iext,contextptr)){ 1816 gen iextnum,iextden; 1817 fxnd(iext,iextnum,iextden); 1818 for (;!embed.empty();){ 1819 iextnum=polynome(iextnum,embed.back()); 1820 iextden=polynome(iextden,embed.back()); 1821 embed.pop_back(); 1822 } 1823 return iextnum/iextden; 1824 } 1825 } 1826 } 1827 return 0; 1828 } 1829 rootof_extract(const gen & e,GIAC_CONTEXT)1830 gen rootof_extract(const gen & e,GIAC_CONTEXT){ 1831 if (e.type==_VECT && e._VECTptr->size()==2) 1832 return e._VECTptr->front()*symb_rootof(gen(makevecteur(1,0),_POLY1__VECT),e._VECTptr->back(),contextptr); 1833 return symbolic(at_rootof,e); 1834 } 1835 lvar_rootof(const gen & e,const vecteur & l,vecteur & lv,GIAC_CONTEXT)1836 void lvar_rootof(const gen & e,const vecteur &l,vecteur & lv,GIAC_CONTEXT){ 1837 if (!l.empty() && l.front().type==_VECT && has_op(e,*at_rootof)){ 1838 vecteur lve(lvar(e)),lvr; 1839 for (int i=0;i<lve.size();++i){ 1840 if (lve[i].is_symb_of_sommet(at_rootof)){ 1841 lvr.push_back(lve[i]); 1842 continue; 1843 } 1844 } 1845 if (!lvr.empty()){ 1846 vector<const unary_function_ptr *> vu; 1847 vu.push_back(at_rootof); 1848 vector <gen_op_context> vv; 1849 vv.push_back(rootof_extract); 1850 gen er=subst(e,vu,vv,true,contextptr); 1851 lvar(er,lv); 1852 } 1853 else 1854 lvar(e,lv); 1855 } 1856 else 1857 lvar(e,lv); 1858 } 1859 sym2r(const gen & e,const vecteur & l,int l_size,gen & num,gen & den,GIAC_CONTEXT)1860 bool sym2r (const gen &e,const vecteur &l, int l_size,gen & num,gen & den,GIAC_CONTEXT){ 1861 if (e.type<_POLY){ 1862 num=e; 1863 den=plus_one; 1864 return true; 1865 } 1866 if ( (e.type==_POLY) || (e.type==_EXT)){ 1867 if ((!l.empty()) && (l.front().type==_VECT) ){ 1868 num=e; 1869 for (int k=int(l.size())-1;k>=0;--k) // was l.size() 1870 num=monomial2gen(monomial<gen>(num,int(l[k]._VECTptr->size()))); 1871 den=plus_one; 1872 } 1873 else { 1874 num=monomial2gen(monomial<gen>(e,l_size)); 1875 den=plus_one; 1876 } 1877 return true; 1878 } 1879 bool totally_converted=true; 1880 vecteur lv,lvnum,lvden; 1881 lvar_rootof(e,l,lv,context0); 1882 if (!compute_lv_lvnum_lvden(l,lv,lvnum,lvden,totally_converted,l_size,contextptr)){ 1883 num=undef; 1884 return false; 1885 } 1886 gen iext=find_iext(e,lvnum,contextptr); 1887 clean_iext(lvnum,lvden,iext,contextptr); 1888 totally_converted =totally_converted && sym2r(e,iext,l,lv,lvnum,lvden,l_size,num,den,contextptr); 1889 if (is_inf(num) && !is_inf(den)) den=1; 1890 if (is_inf(den) && !is_inf(num)) num=1; 1891 // If den is a _POLY, multiply den by the _EXT conjugate of it's lcoeff 1892 // FIXME this should be done recursively if the 1st coeff is a _POLY! 1893 if (den.type==_POLY && !den._POLYptr->coord.empty() && den._POLYptr->coord.front().value.type==_EXT){ 1894 gen a = ext_reduce(den._POLYptr->coord.front().value); 1895 if (is_undef(a)) return false; 1896 if (a.type == _EXT && a._EXTptr->type==_VECT){ 1897 vecteur u,v,d; 1898 egcd(*(a._EXTptr->_VECTptr),*((a._EXTptr+1)->_VECTptr),0,u,v,d); 1899 if (d.size()==1){ 1900 gen aconj=algebraic_EXTension(u,*(a._EXTptr+1)); 1901 aconj=polynome(aconj,int(den._POLYptr->coord.front().index.size())); 1902 num=aconj*num; 1903 den=aconj*den; 1904 simplify3(num,den); 1905 } 1906 } 1907 } 1908 return totally_converted; 1909 } 1910 sym2r(const gen & e,const vecteur & l,GIAC_CONTEXT)1911 fraction sym2r(const gen & e, const vecteur & l,GIAC_CONTEXT){ 1912 int l_size; 1913 if (!l.empty() && l.front().type==_VECT) 1914 l_size=int(l.front()._VECTptr->size()); 1915 else 1916 l_size=int(l.size()); 1917 gen num,den; 1918 bool ok=sym2r(e,l,l_size,num,den,contextptr); 1919 if (!ok){ 1920 num=string2gen("Error in normal",false); 1921 num.subtype=-1; //gensizeerr(contextptr); 1922 } 1923 if (is_positive(-den,contextptr)) 1924 return fraction(-num,-den); 1925 else 1926 return fraction(num,den); 1927 } 1928 1929 // fraction / x -> fraction of vecteur e2r(const gen & e,const gen & x,GIAC_CONTEXT)1930 gen e2r(const gen & e,const gen & x,GIAC_CONTEXT){ 1931 vecteur l(1,x); 1932 lvar(e,l); 1933 gen r=polynome2poly1(e2r(e,l,contextptr),1); 1934 return r2e(r,cdr_VECT(l),contextptr); 1935 } 1936 e2r(const gen & e,const vecteur & l,GIAC_CONTEXT)1937 gen e2r(const gen & e,const vecteur & l,GIAC_CONTEXT){ 1938 if (e.type!=_VECT) 1939 return sym2r(e,l,contextptr); 1940 bool totally_converted=true; 1941 int l_size; 1942 if (!l.empty() && l.front().type==_VECT) 1943 l_size=int(l.front()._VECTptr->size()); 1944 else 1945 l_size=int(l.size()); 1946 gen num,den; 1947 vecteur lv,lvnum,lvden; 1948 lvar_rootof(e,l,lv,contextptr); 1949 if (!compute_lv_lvnum_lvden(l,lv,lvnum,lvden,totally_converted,l_size,contextptr)) 1950 return gensizeerr(contextptr); 1951 gen iext=find_iext(e,lvnum,contextptr); 1952 clean_iext(lvnum,lvden,iext,contextptr); 1953 vecteur res; 1954 const_iterateur jt=e._VECTptr->begin(),jtend=e._VECTptr->end(); 1955 for (;jt!=jtend;++jt){ 1956 sym2r(*jt,iext,l,lv,lvnum,lvden,l_size,num,den,contextptr); 1957 res.push_back(num/den); 1958 } 1959 return gen(res,e.subtype); 1960 } 1961 1962 //********************************** 1963 // tensor to symbolic 1964 //********************************** 1965 niceprod(ref_vecteur * res,bool negatif)1966 static gen niceprod(ref_vecteur * res,bool negatif){ 1967 gen tmp; 1968 if (res->v.size()!=1) 1969 tmp=symbolic(at_prod,gen(res,_SEQ__VECT)); 1970 else { 1971 tmp=res->v.front(); 1972 delete_ref_vecteur(res); 1973 } 1974 return negatif?-tmp:tmp; 1975 } 1976 r2sym(const gen & e,const index_m & i,const vecteur & l,GIAC_CONTEXT)1977 gen r2sym(const gen & e,const index_m & i,const vecteur & l,GIAC_CONTEXT){ 1978 if (is_undef(e)) 1979 return e; 1980 if (i.size()!=l.size()) 1981 return gensizeerr(gettext("sym2poly/r2sym(const gen & e,const index_m & i,const vecteur & l)")); 1982 vecteur::const_iterator l_it=l.begin(); 1983 index_t::const_iterator it=i.begin(),itend=i.end(); 1984 ref_vecteur * res=new_ref_vecteur(0); 1985 res->v.reserve(itend-it+1); 1986 bool negatif=false; 1987 if (!is_one(e) || (e.type==_MOD) ){ 1988 if ( (e.type<=_REAL || e.type==_FRAC) && is_positive(-e,contextptr)){ 1989 negatif=true; 1990 if (!is_minus_one(e)) 1991 res->v.push_back(-e); 1992 } 1993 else 1994 res->v.push_back(e); 1995 } 1996 for (;it!=itend;++it,++l_it){ 1997 if ((*it)) 1998 res->v.push_back(pow(*l_it,*it,contextptr)); // change for normal(abs(z)^2), was pow(*l_it,*it) 1999 } 2000 if (res->v.empty()){ 2001 delete_ref_vecteur(res); 2002 return e; 2003 } 2004 return niceprod(res,negatif); 2005 } 2006 r2sym2(const polynome & p,const vecteur & l,GIAC_CONTEXT)2007 static gen r2sym2(const polynome & p, const vecteur & l,GIAC_CONTEXT){ 2008 monomial_v::const_iterator it=p.coord.begin(); 2009 monomial_v::const_iterator itend=p.coord.end(); 2010 if (itend-it==1) 2011 return r2sym(it->value,it->index,l,contextptr); 2012 ref_vecteur * res=new_ref_vecteur(0); 2013 res->v.reserve(itend-it); 2014 for (;it!=itend;++it){ 2015 res->v.push_back(r2sym(it->value,it->index,l,contextptr)) ; 2016 } 2017 if (res->v.size()==1){ 2018 gen tmp=res->v.front(); 2019 delete res; 2020 return tmp; 2021 } 2022 if (increasing_power(contextptr)) 2023 reverse(res->v.begin(),res->v.end()); 2024 return symbolic(at_plus,gen(res,_SEQ__VECT)); 2025 } 2026 r2sym(const polynome & p,const vecteur & l,GIAC_CONTEXT)2027 gen r2sym(const polynome & p, const vecteur & l,GIAC_CONTEXT){ 2028 if (p.coord.empty()) 2029 return 0; 2030 if (p.dim==0){ 2031 return p.constant_term(); 2032 } 2033 if (is_positive(-p.coord.front())) 2034 return -r2sym2(-p,l,contextptr); 2035 return r2sym2(p,l,contextptr); 2036 } 2037 r2sym(const fraction & f,const vecteur & l,GIAC_CONTEXT)2038 gen r2sym(const fraction & f, const vecteur & l,GIAC_CONTEXT){ 2039 if (f.den.type==_POLY && is_positive(-f.den._POLYptr->coord.front())){ 2040 return rdiv(r2sym(-f.num,l,contextptr),r2sym(-f.den,l,contextptr),contextptr); 2041 } 2042 // workaround for e.g. normal(sqrt(3)*sqrt(15)) 2043 gen fn=r2sym(f.num,l,contextptr); 2044 gen fd=r2sym(f.den,l,contextptr); 2045 if (fn.is_symb_of_sommet(at_prod) && fn._SYMBptr->feuille.type==_VECT && fd.type==_INT_ && absint(fd.val)>1){ 2046 gen c=1; 2047 vecteur v(*fn._SYMBptr->feuille._VECTptr); 2048 fn=1; 2049 for (unsigned i=0;i<v.size();++i){ 2050 if (v[i].type==_INT_){ 2051 c=c*v[i]; 2052 simplify(c,fd); 2053 } 2054 else 2055 fn=fn*v[i]; 2056 } 2057 fn=c*fn; 2058 } 2059 gen res=rdiv(fn,fd,contextptr); 2060 if (has_inf_or_undef(res)) 2061 return fraction(fn,fd); 2062 return res; 2063 } 2064 r2sym(const vecteur & v,const vecteur & l,GIAC_CONTEXT)2065 gen r2sym(const vecteur & v,const vecteur & l,GIAC_CONTEXT){ 2066 const_iterateur it=v.begin(),itend=v.end(); 2067 ref_vecteur * res=new_ref_vecteur(0); 2068 res->v.reserve(itend-it); 2069 for (;it!=itend;++it) 2070 res->v.push_back(r2sym(*it,l,contextptr)); 2071 return res; 2072 } 2073 2074 static gen r2sym(const gen & p, const const_iterateur & lt, const const_iterateur & ltend,GIAC_CONTEXT); 2075 solve_deg2(const gen & p,const gen & pmin,const const_iterateur & lt,const const_iterateur & ltend,GIAC_CONTEXT)2076 static gen solve_deg2(const gen & p,const gen & pmin,const const_iterateur & lt, const const_iterateur & ltend,GIAC_CONTEXT){ 2077 if ( p.type!=_VECT || pmin.type!=_VECT) 2078 return gensizeerr(gettext("sym2poly.cc/ckdeg2_rootof")); 2079 vecteur & w = *pmin._VECTptr; 2080 int s=int(w.size()); 2081 if (s!=3)// || is_greater(linfnorm(pmin,contextptr),(1<<30),contextptr)) 2082 return symb_rootof(p,r2sym(pmin,lt,ltend,contextptr),contextptr); 2083 if (p._VECTptr->size()!=2) 2084 return gensizeerr(gettext("sym2poly.cc/ckdeg2_rootof")); 2085 gen b_over_2=rdiv(w[1],plus_two,contextptr); 2086 if (is_zero(b_over_2)) 2087 return p._VECTptr->front()*sqrt(r2sym(-w.back()/w.front(),lt,ltend,contextptr),contextptr)+p._VECTptr->back(); 2088 gen x; 2089 if (b_over_2.type!=_FRAC){ 2090 gen a=r2sym(w.front(),lt,ltend,contextptr); 2091 gen minus_b_over_2=r2sym(-b_over_2,lt,ltend,contextptr); 2092 gen delta_prime=r2sym(pow(b_over_2,2,contextptr)-w.front()*w.back(),lt,ltend,contextptr); 2093 x=rdiv(minus_b_over_2+sqrt(delta_prime,contextptr),a,contextptr); 2094 } 2095 else { 2096 gen two_a=r2sym(plus_two*w.front(),lt,ltend,contextptr); 2097 gen minus_b=r2sym(-w[1],lt,ltend,contextptr); 2098 gen delta=r2sym(w[1]*w[1]-gen(4)*w.front()*w.back(),lt,ltend,contextptr); 2099 x=rdiv(minus_b+sqrt(delta,contextptr),two_a,contextptr); 2100 } 2101 return p._VECTptr->front()*x+p._VECTptr->back(); 2102 } 2103 2104 /* 2105 static gen ckdeg2_rootof(const gen & p,const gen & pmin,GIAC_CONTEXT){ 2106 if ( p.type!=_VECT || pmin.type!=_VECT) 2107 return gensizeerr(gettext("sym2poly.cc/ckdeg2_rootof")); 2108 if (pmin._VECTptr->size()!=3) 2109 return symb_rootof(p,pmin,contextptr); 2110 if (p._VECTptr->size()!=2) 2111 return gensizeerr(gettext("sym2poly.cc/ckdeg2_rootof")); 2112 vecteur v; 2113 identificateur x(" x"); 2114 in_solve(symb_horner(*pmin._VECTptr,x),x,v,1,contextptr); 2115 if (v.empty()) 2116 return gensizeerr(gettext("No root found for pmin")); 2117 return p._VECTptr->front()*v.front()+p._VECTptr->back(); 2118 } 2119 */ 2120 2121 // check that an _EXT e is root of a second order poly 2122 // return 0 if not, 2 if pmin(e) is 2nd order, 1 otherwise 2123 // FIXME: if pmin has integer coeffs as well as e 2124 // translate so that we just need to find the sign 2125 // then we are reduced to find the sign of P(alpha) where alpha is the 2126 // largest real root of a polynomial Q 2127 // Using VAS [csturm.cc: gen complexroot(const gen & g,bool complexe,GIAC_CONTEXT)] 2128 // find isolation interval for roots of Q, select the largest one, 2129 // find isolation interval for roots of P, if they intersect with previous one 2130 // improve precision by dichotomy until sign is resolved is_root_of_deg2(const gen & e,vecteur & v,GIAC_CONTEXT)2131 static int is_root_of_deg2(const gen & e,vecteur & v,GIAC_CONTEXT){ 2132 // find a,b,c such that a*e*e+b*e+c=0 2133 #ifdef DEBUG_SUPPORT 2134 if (e.type!=_EXT || e._EXTptr->type!=_VECT) 2135 setsizeerr(gettext("sym2poly.cc/is_root_of_deg2")); 2136 #endif // DEBUG_SUPPORT 2137 gen pmin=*(e._EXTptr+1),value; 2138 if (has_rootof_value(pmin,value,contextptr)) 2139 return 0; 2140 vecteur vpmin(min_pol(pmin)); 2141 if (!vpmin.empty() && is_undef(vpmin)) 2142 return 0; 2143 if (vpmin.size()==3) 2144 return 2; 2145 v.clear(); 2146 gen e_square(e*e); 2147 if (e_square.type!=_EXT){ // b=0, a=1, e*e=-c 2148 v.push_back(plus_one); 2149 v.push_back(zero); 2150 v.push_back(-e_square); 2151 return 1; 2152 } 2153 if (e_square._EXTptr->type!=_VECT) 2154 return 0; // setsizeerr(gettext("sym2poly.cc/is_root_of_deg2")); 2155 int s=int(e._EXTptr->_VECTptr->size()); 2156 int s2=int(e_square._EXTptr->_VECTptr->size()); 2157 if (s!=s2) 2158 return 0; 2159 gen b=-e_square._EXTptr->_VECTptr->front(); 2160 gen a=e._EXTptr->_VECTptr->front(); 2161 simplify(a,b); 2162 gen c=a*e_square+b*e; 2163 if (c.type==_EXT) 2164 return 0; 2165 v.push_back(a); 2166 v.push_back(b); 2167 v.push_back(-c); 2168 return 1; 2169 } 2170 randassume(const vecteur & boundaries,GIAC_CONTEXT)2171 static gen randassume(const vecteur & boundaries,GIAC_CONTEXT){ 2172 if (boundaries.empty()) 2173 return 0; 2174 if (boundaries[0].type==_VECT) 2175 return randassume(*boundaries[0]._VECTptr,contextptr); 2176 if (boundaries.size()<2) 2177 return 0; 2178 gen a=boundaries[0],b=boundaries[1]; 2179 if (a==minus_inf){ 2180 if (b==plus_inf) 2181 return _rand(makesequence(-100,100),contextptr); 2182 return _rand(makesequence(b-100,b-1),contextptr); 2183 } 2184 if (b==plus_inf) 2185 return _rand(makesequence(a+1,a+100),contextptr); 2186 return (b-a)*_rand(makesequence(0.0,1.0),contextptr)+a; 2187 } 2188 check_assume(vecteur & vzero,const vecteur & vassume,GIAC_CONTEXT)2189 static void check_assume(vecteur & vzero,const vecteur & vassume,GIAC_CONTEXT){ 2190 for (unsigned i=0;i<vzero.size();++i){ 2191 gen assi=vassume[i],tmp; 2192 if (assi.type==_IDNT && has_evalf(assi,tmp,1,contextptr)){ 2193 vzero[i]=tmp; 2194 continue; 2195 } 2196 if (assi.type==_VECT && assi.subtype==_ASSUME__VECT && assi._VECTptr->size()==3){ 2197 if((*assi._VECTptr)[1].type==_VECT && !(*assi._VECTptr)[1]._VECTptr->empty()){ 2198 vzero[i]=randassume(*(*assi._VECTptr)[1]._VECTptr,contextptr); 2199 } 2200 } 2201 } 2202 } 2203 minimal_polynomial(const gen & pp,bool minonly,GIAC_CONTEXT)2204 gen minimal_polynomial(const gen & pp,bool minonly,GIAC_CONTEXT){ 2205 gen f=*(pp._EXTptr+1); 2206 if (f.type!=_VECT) 2207 return undef; 2208 int d=int(f._VECTptr->size()); 2209 gen r=evalf_double(pp,1,contextptr); 2210 matrice m(d); 2211 m[0]=vecteur(d-1); 2212 m[0]._VECTptr->back()=1; 2213 gen ppi(1); 2214 for (int i=1;i<d;++i){ 2215 ppi=ppi*pp; 2216 if (ppi.type!=_EXT){ 2217 m[i]=vecteur(d-1); 2218 m[i]._VECTptr->back()=ppi; 2219 } 2220 else { 2221 m[i]=*ppi._EXTptr; 2222 if (m[i].type!=_VECT) 2223 return gensizeerr(contextptr); 2224 m[i]=*m[i]._VECTptr; 2225 if (int(m[i]._VECTptr->size())<d-1) 2226 *m[i]._VECTptr=mergevecteur(vecteur(d-1-m[i]._VECTptr->size()),*m[i]._VECTptr); 2227 } 2228 } // end for i 2229 if (!ckmatrix(m)) 2230 return gensizeerr(contextptr); 2231 m=mtran(m); // d-1 rows, d columns 2232 int st=step_infolevel(contextptr); 2233 step_infolevel(contextptr)=0; 2234 vecteur mk=mker(m,contextptr); 2235 step_infolevel(contextptr)=st; 2236 gen pm(mk.front()); // min poly= 1st kernel element 2237 if (pm.type==_VECT){ 2238 mk=*pm._VECTptr; 2239 reverse(mk.begin(),mk.end()); 2240 mk=trim(mk,0); 2241 if (minonly){ 2242 if (is_positive(-mk.front(),contextptr)) 2243 return gen(-mk,_POLY1__VECT); 2244 return gen(mk,_POLY1__VECT); 2245 } 2246 if ((d=int(mk.size()))<=5 && d!=4){ 2247 identificateur x(" x"); 2248 vecteur w; 2249 in_solve(symb_horner(mk,x),x,w,1,contextptr); 2250 if (r.type!=_DOUBLE_ && r.type!=_CPLX && !w.empty()) return w.front(); 2251 for (unsigned k=0;k<w.size();++k){ 2252 if (lop(w[k],at_rootof).empty()){ 2253 gen wkd=evalf_double(w[k],1,contextptr); 2254 if (wkd!=w[k] && is_greater(1e-6,abs(1-r/wkd,contextptr),contextptr)) 2255 return w[k]; 2256 } 2257 } 2258 } 2259 } 2260 return undef; 2261 } 2262 accurate_evalf_until(const gen & g_,GIAC_CONTEXT)2263 gen accurate_evalf_until(const gen & g_,GIAC_CONTEXT){ 2264 #ifdef HAVE_LIBMPFR 2265 gen g(g_); 2266 int n=32; 2267 gen cur=evalf_double(g,1,contextptr); 2268 if (g.type==_EXT){ 2269 gen p=*g._EXTptr; 2270 gen x=symb_rootof(makevecteur(1,0),*(g._EXTptr+1),contextptr); 2271 for (;n<=1024;n*=2){ 2272 gen tmp=_evalf(makesequence(x,n),contextptr); 2273 tmp=_horner(makesequence(p,tmp),contextptr); 2274 gen test=(1-cur/tmp); 2275 if (is_greater(1e-12,test,contextptr)) 2276 return cur; 2277 cur=tmp; 2278 } 2279 } 2280 for (;n<=1024;n*=2){ 2281 gen tmp=_evalf(makesequence(g,n),contextptr); 2282 gen test=(1-cur/tmp); 2283 if (is_greater(1e-12,test,contextptr)) 2284 return cur; 2285 cur=tmp; 2286 } 2287 return cur; 2288 #else 2289 return evalf_double(g_,1,contextptr); 2290 #endif 2291 } 2292 r2sym(const gen & p,const const_iterateur & lt,const const_iterateur & ltend,GIAC_CONTEXT)2293 static gen r2sym(const gen & p, const const_iterateur & lt, const const_iterateur & ltend,GIAC_CONTEXT){ 2294 // Note that if p.type==_FRAC && p._FRACptr->num.type==_EXT 2295 // it might require another simplification with the denom 2296 if (p.type==_FRAC){ 2297 gen res; 2298 if (p._FRACptr->den.type==_CPLX) 2299 res=rdiv(r2sym(p._FRACptr->num*conj(p._FRACptr->den,contextptr),lt,ltend,contextptr),r2sym(p._FRACptr->den*conj(p._FRACptr->den,contextptr),lt,ltend,contextptr),contextptr); 2300 else 2301 res=rdiv(r2sym(p._FRACptr->num,lt,ltend,contextptr),r2sym(p._FRACptr->den,lt,ltend,contextptr),contextptr); 2302 return (p._FRACptr->num.type==_EXT)?ratnormal(res,contextptr):res; 2303 } 2304 if (p.type==_MOD) 2305 return makemodquoted(r2sym(*p._MODptr,lt,ltend,contextptr),r2sym(*(p._MODptr+1),lt,ltend,contextptr)); 2306 if (p.type==_VECT){ 2307 const_iterateur it=p._VECTptr->begin(),itend=p._VECTptr->end(); 2308 vecteur res; 2309 res.reserve(itend-it); 2310 for (;it!=itend;++it) 2311 res.push_back(r2sym(*it,lt,ltend,contextptr)); 2312 return res; 2313 } 2314 if (p.type==_EXT){ 2315 gen pp=ext_reduce(p); 2316 if (pp.type!=_EXT){ 2317 if (is_undef(pp)) return pp; 2318 return r2sym(pp,lt,ltend,contextptr); 2319 } 2320 gen f=*(pp._EXTptr+1),fvalue; 2321 #if 1 2322 if ( ( ltend==lt || (ltend-lt==1 && lt->type==_VECT && lt->_VECTptr->empty()) ) 2323 && f.type==_VECT && !has_rootof_value(f,fvalue,contextptr) && !keep_algext(contextptr) ){ 2324 // univariate case 2325 // find minimal poly of the whole _EXT if extension degree is > 2 2326 int d=int(f._VECTptr->size()); 2327 // FIXME remove d<=10, requires better handling of rref with Gauss integers 2328 if (d>3 && d<=10){ 2329 gen res=minimal_polynomial(pp,false,contextptr); 2330 if (!is_undef(res)) 2331 return res; 2332 } // end if d>3 2333 } 2334 #endif 2335 if ((pp._EXTptr+1)->type!=_VECT) 2336 f=min_pol(*(pp._EXTptr+1)); // f is a _VECT 2337 if (is_undef(f)) 2338 return f; 2339 if (f._VECTptr->size()==3){ 2340 gen value; 2341 if (!has_rootof_value(f,value,contextptr)) 2342 return solve_deg2(r2sym(*(pp._EXTptr),lt,ltend,contextptr),f,lt,ltend,contextptr); 2343 } 2344 vecteur v; 2345 int t=0; 2346 fvalue=undef; 2347 if (!has_rootof_value(f,fvalue,contextptr)) 2348 t=is_root_of_deg2(pp,v,contextptr); 2349 // if (t==2) return solve_deg2(r2sym(*(pp._EXTptr),lt,ltend,contextptr),f,lt,ltend,contextptr); 2350 // if (t && f==v) t=0; 2351 if (t==1){ 2352 vecteur w; 2353 identificateur x(" x"); 2354 in_solve(symb_horner(*(r2sym(v,lt,ltend,contextptr)._VECTptr),x),x,w,1,contextptr); 2355 if (w.size()==1) 2356 return w.front(); 2357 #ifndef NO_STDEXCEPT 2358 try { 2359 #endif 2360 vecteur vinit(lt,ltend); 2361 if (lt!=ltend && vinit.front().type==_VECT) 2362 vinit=*vinit.front()._VECTptr; 2363 vecteur vassume(vinit); 2364 for (unsigned i=0;i<vassume.size();++i){ 2365 if (vinit[i].type==_IDNT) 2366 vinit[i]._IDNTptr->in_eval(0,vinit[i],vassume[i],contextptr); 2367 } 2368 gen vinitd; 2369 vecteur vzero=vecteur(vinit.size()); 2370 bool tst=has_evalf(vinit,vinitd,1,contextptr); 2371 if (tst) 2372 vzero=*vinitd._VECTptr; 2373 else 2374 check_assume(vzero,vassume,contextptr); 2375 vzero=subst(vzero,vinit,vzero,false,contextptr); 2376 gen tmp0=r2sym(*pp._EXTptr,lt,ltend,contextptr); 2377 gen tmp1=r2sym(f,lt,ltend,contextptr); 2378 for (int ntry=0;ntry<10;++ntry,tst=false){ 2379 gen tmp00=subst(tmp0,vinit,vzero,false,contextptr); 2380 gen tmp10=subst(tmp1,vinit,vzero,false,contextptr); 2381 if (!lvar(makevecteur(tmp00,tmp10)).empty()) 2382 break; 2383 if (tmp10.type==_VECT && _size(_gcd(makesequence(tmp10,derivative(*tmp10._VECTptr)),contextptr),contextptr)>1) 2384 tmp00=plus_one_half; // retry 2385 else { 2386 tmp00=algebraic_EXTension(tmp00,tmp10); 2387 tmp00=accurate_evalf_until(tmp00,contextptr); // tmp00.evalf_double(1,contextptr); 2388 } 2389 if (tmp00.type<=_CPLX){ 2390 gen tmp3=eval(subst(w.front(),vinit,vzero,false,contextptr),1,contextptr); 2391 tmp3=accurate_evalf_until(tmp3,contextptr); // tmp3.evalf_double(1,contextptr); 2392 gen tmp2=eval(subst(w.back(),vinit,vzero,false,contextptr),1,contextptr); 2393 tmp2=accurate_evalf_until(tmp2,contextptr);// tmp2.evalf_double(1,contextptr); 2394 if (tmp3.type<=_CPLX && tmp2.type<=_CPLX && tmp2!=tmp3){ 2395 if (!vzero.empty() && !tst) 2396 *logptr(contextptr) << gettext("Warning, choosing root of ") << f << gettext(" at parameters values ") << vzero << "\n"; 2397 if (is_greater(abs(tmp3-tmp00,contextptr),abs(tmp2-tmp00,contextptr),contextptr)) 2398 return w.back(); 2399 else 2400 return w.front(); 2401 } 2402 } 2403 // tmp2 and tmp3 are identical or tmp0 is not real, retry 2404 vzero=vranm(int(vzero.size()),0,contextptr); 2405 check_assume(vzero,vassume,contextptr); 2406 } 2407 #ifndef NO_STDEXCEPT 2408 } 2409 catch (std::runtime_error & e){ 2410 last_evaled_argptr(contextptr)=NULL; 2411 CERR << "sym2poly exception caught " << e.what() << "\n"; 2412 } 2413 #endif 2414 return w.front(); 2415 /* gen tmp0=algebraic_EXTension(r2sym(*pp._EXTptr,lt,ltend),r2sym(f,lt,ltend)).evalf_double(); 2416 if (tmp0.type==_DOUBLE_){ 2417 gen tmp1=evalf_double(w.front()),tmp2=evalf_double(w.back()); 2418 if (tmp1.type==_DOUBLE_ && tmp2.type==_DOUBLE_ && fabs((tmp2-tmp0)._DOUBLE_val)<fabs((tmp1-tmp0)._DOUBLE_val)) 2419 return w.back(); 2420 } 2421 return w.front(); */ 2422 } 2423 int s=int(f._VECTptr->size()); 2424 v=vecteur(s,zero); 2425 v.front()=plus_one; 2426 v.back()=f._VECTptr->back(); 2427 if (f==v){ 2428 if (is_undef(fvalue)) 2429 fvalue=pow(r2sym(-v.back(),lt,ltend,contextptr),fraction(1,s-1),contextptr); 2430 } 2431 else 2432 return symb_rootof(r2sym(*(pp._EXTptr),lt,ltend,contextptr),r2sym(f,lt,ltend,contextptr),contextptr); 2433 return symb_horner(*(r2sym(*(pp._EXTptr),lt,ltend,contextptr)._VECTptr),fvalue); 2434 } 2435 if (p.type==_SPOL1){ 2436 sparse_poly1 s=*p._SPOL1ptr; 2437 sparse_poly1::iterator it=s.begin(),itend=s.end(); 2438 vecteur l(lt,ltend); 2439 for (;it!=itend;++it){ 2440 it->coeff=r2sym(it->coeff,l,contextptr); 2441 } 2442 return s; 2443 } 2444 if ((p.type!=_POLY) || (lt==ltend)) 2445 return p; 2446 if (p._POLYptr->coord.empty()) 2447 return zero; 2448 if (p._POLYptr->dim==0){ 2449 return r2sym(p._POLYptr->coord.front().value,lt+1,ltend,contextptr); 2450 } 2451 if (is_positive(-p,contextptr) && !is_positive(p,contextptr)) 2452 return -r2sym(-p,lt,ltend,contextptr); 2453 monomial_v::const_iterator it=p._POLYptr->coord.begin(); 2454 monomial_v::const_iterator itend=p._POLYptr->coord.end(); 2455 if (itend==it) 2456 return zero; 2457 vecteur res; 2458 res.reserve(itend-it); 2459 for (;it!=itend;++it){ 2460 res.push_back(r2sym(r2sym(it->value,lt+1,ltend,contextptr),it->index,*(lt->_VECTptr),contextptr)) ; 2461 } 2462 if (res.size()==1) 2463 return res.front(); 2464 if (increasing_power(contextptr)) 2465 reverse(res.begin(),res.end()); 2466 return symbolic(at_plus,gen(res,_SEQ__VECT)); 2467 } 2468 r2sym(const gen & p,const vecteur & l,GIAC_CONTEXT)2469 gen r2sym(const gen & p, const vecteur & l,GIAC_CONTEXT){ 2470 if (p.type==_VECT){ 2471 gen res=r2sym(*p._VECTptr,l,contextptr); 2472 if (res.type==_VECT) 2473 res.subtype=p.subtype; 2474 return res; 2475 } 2476 if ( (!l.empty()) && (l.front().type==_VECT) ) 2477 return r2sym(p,l.begin(),l.end(),contextptr); 2478 if (p.type==_FRAC) 2479 return r2sym(*p._FRACptr,l,contextptr); 2480 if (p.type==_MOD) 2481 return makemodquoted(r2sym(*p._MODptr,l,contextptr),r2sym(*(p._MODptr+1),l,contextptr)); 2482 if (p.type==_EXT){ 2483 gen pp=ext_reduce(*(p._EXTptr),*(p._EXTptr+1)); 2484 if (pp.type!=_EXT){ 2485 if (is_undef(pp)) return pp; 2486 return r2sym(pp,l,contextptr); 2487 } 2488 gen f=*(pp._EXTptr+1); 2489 if ((pp._EXTptr+1)->type!=_VECT) 2490 f=min_pol(*(pp._EXTptr+1)); // f is a _VECT 2491 if (is_undef(f)) 2492 return f; 2493 vecteur v; 2494 int t=is_root_of_deg2(pp,v,contextptr); 2495 if (f._VECTptr->size()==3){ 2496 return solve_deg2(r2sym(*(pp._EXTptr),l,contextptr),f,l.begin(),l.end(),contextptr); 2497 } 2498 if (t==2){ 2499 // return ckdeg2_rootof(r2sym(*(pp._EXTptr),l,contextptr),r2sym(f,l,contextptr),contextptr); 2500 return solve_deg2(r2sym(*(pp._EXTptr),l,contextptr),f,l.begin(),l.end(),contextptr); 2501 } 2502 if (t==1){ 2503 vecteur w; 2504 identificateur x(" x"); 2505 in_solve(symb_horner(*(r2sym(v,l,contextptr)._VECTptr),x),x,w,1,contextptr); 2506 // find the nearest root from pp 2507 gen ww,ppp,mindist,curdist; 2508 if (has_evalf(w,ww,1,contextptr) && has_evalf(pp,ppp,1,contextptr)){ 2509 int pos=0; 2510 vecteur & www = *ww._VECTptr; 2511 mindist=abs(www[0]-ppp,contextptr); 2512 for (int i=1;i<(int)www.size();++i){ 2513 curdist=abs(www[i]-ppp,contextptr); 2514 if (is_strictly_greater(mindist,curdist,contextptr)){ 2515 pos=i; 2516 mindist=curdist; 2517 } 2518 } 2519 return w[pos]; 2520 } 2521 return w.front(); 2522 } 2523 int s=int(f._VECTptr->size()); 2524 v=vecteur(s,zero); 2525 v.front()=plus_one; 2526 v.back()=f._VECTptr->back(); 2527 gen theta; 2528 if (f==v) 2529 theta=pow(r2sym(-v.back(),l,contextptr),fraction(1,s-1),contextptr); 2530 else 2531 return symb_rootof(r2sym(*(pp._EXTptr),l,contextptr),r2sym(f,l,contextptr),contextptr); 2532 return symb_horner(*(r2sym(*(pp._EXTptr),l,contextptr)._VECTptr),theta); 2533 } 2534 if (p.type==_POLY) 2535 return r2sym(*p._POLYptr,l,contextptr); 2536 if (p.type==_SPOL1) 2537 return r2sym(p,l.begin(),l.end(),contextptr); 2538 return p; 2539 } 2540 r2e(const gen & p,const vecteur & l,GIAC_CONTEXT)2541 gen r2e(const gen & p, const vecteur & l,GIAC_CONTEXT){ 2542 return r2sym(p,l,contextptr); 2543 } 2544 2545 // alias vecteur -> polynome / x r2e(const gen & r,const gen & x,GIAC_CONTEXT)2546 gen r2e(const gen & r,const gen & x,GIAC_CONTEXT){ 2547 if (r.type==_FRAC) 2548 return fraction(r2e(r._FRACptr->num,x,contextptr),r2e(r._FRACptr->den,x,contextptr)); 2549 if (r.type==_VECT) 2550 return symb_horner(*r._VECTptr,x); 2551 else 2552 return r; 2553 } 2554 2555 // convert factorization to symbolic form r2sym(const factorization & vnum,const vecteur & l,GIAC_CONTEXT)2556 gen r2sym(const factorization & vnum,const vecteur & l,GIAC_CONTEXT){ 2557 gen resnum(1); 2558 factorization::const_iterator it=vnum.begin(),itend=vnum.end(); 2559 for (;it!=itend;++it){ 2560 // insure no embedded sqrt inside 2561 polynome p=it->fact; 2562 if (l.size()==1){ 2563 vector< monomial<gen> >::iterator pt=p.coord.begin(),ptend=p.coord.end(); 2564 vecteur vtmp(1,vecteur(0)); 2565 for (;pt!=ptend;++pt){ 2566 pt->value=r2sym(pt->value,vtmp,contextptr); 2567 } 2568 } 2569 gen tmp=r2sym(gen(p),l,contextptr); 2570 resnum=resnum*pow(tmp,it->mult); 2571 } 2572 return resnum; 2573 } 2574 2575 // convert pfde_VECT to symbolic form r2sym(const vector<pf<gen>> & pfde_VECT,const vecteur & l,GIAC_CONTEXT)2576 gen r2sym(const vector< pf<gen> > & pfde_VECT,const vecteur & l,GIAC_CONTEXT){ 2577 gen res(0); 2578 vector< pf<gen> >::const_iterator it=pfde_VECT.begin(),itend=pfde_VECT.end(); 2579 for (;it!=itend;++it){ 2580 res=res+rdiv(r2sym(gen(it->num),l,contextptr),(r2sym(gen(it->den/pow(it->fact,it->mult)),l,contextptr)*pow(r2sym(gen(it->fact),l,contextptr),it->mult)),contextptr); 2581 } 2582 return res; 2583 } 2584 2585 //***************************** 2586 /* Fonctions relatives to ex */ 2587 //***************************** 2588 2589 // special version of lop(.,at_pow) that removes integral powers lop_pow(const gen & g)2590 static vecteur lop_pow(const gen & g){ 2591 if (g.type==_SYMB) { 2592 if (g._SYMBptr->sommet==at_pow && g._SYMBptr->feuille._VECTptr->back().type!=_INT_) 2593 return vecteur(1,g); 2594 else { 2595 if (g._SYMBptr->sommet!=at_ln) 2596 return lop_pow(g._SYMBptr->feuille); 2597 } 2598 } 2599 if (g.type!=_VECT) 2600 return vecteur(0); 2601 vecteur res; 2602 const_iterateur it=g._VECTptr->begin(),itend=g._VECTptr->end(); 2603 for (;it!=itend;++it){ 2604 if (it->is_symb_of_sommet(at_pow) && it->_SYMBptr->feuille._VECTptr->back().type==_INT_) 2605 continue; 2606 res=mergeset(res,lop_pow(*it)); 2607 } 2608 return res; 2609 } 2610 2611 // sqrt(g), check whether g is a+b+/-2*sqrt(a*b) is_sqrtxy(const gen & g_,gen & res,GIAC_CONTEXT)2612 static bool is_sqrtxy(const gen & g_,gen & res,GIAC_CONTEXT){ 2613 gen g(g_); 2614 vecteur l=lop_pow(g); 2615 if (l.size()==2 && l[0]._SYMBptr->feuille[1]==plus_one_half && l[1]._SYMBptr->feuille[1]==plus_one_half){ 2616 gen c,d,e,f; 2617 if (!is_linear_wrt(g,l[0],c,d,contextptr)) 2618 return false; 2619 // c*l[0]+d 2620 if (!is_constant_wrt(d,l[1],contextptr)) 2621 return false; 2622 if (!is_linear_wrt(c,l[1],e,f,contextptr) || !is_zero(f,contextptr)) 2623 return false; 2624 // rewrite g 2625 l=vecteur(1,symb_pow(l[0]._SYMBptr->feuille[0]*l[1]._SYMBptr->feuille[0],plus_one_half)); 2626 g=d+e*l[0]; 2627 } 2628 if (l.size()!=1) // improve for sqrt(a)*sqrt(b) 2629 return false; 2630 gen e(l[0]._SYMBptr->feuille); 2631 if (e.type!=_VECT || e._VECTptr->size()!=2 || e._VECTptr->back()!=plus_one_half) 2632 return false; 2633 e=e._VECTptr->front(); 2634 gen c,d; 2635 if (!is_linear_wrt(g,l[0],c,d,contextptr)) 2636 return false; 2637 e=e*(c*c/gen(4)); 2638 gen cs=sign(c,contextptr); 2639 // d=a+b, e=a*b, factor x^2-d*x+e 2640 l=lidnt(makevecteur(d,e)); 2641 gen x; 2642 for (;;){ 2643 x=identificateur("x"+print_INT_(giac_rand(contextptr))); 2644 if (!equalposcomp(l,x)) 2645 break; 2646 } 2647 gen f=x*x-d*x+e; 2648 f=factor(f,x,false,contextptr); 2649 if (f.type!=_SYMB) 2650 return false; 2651 if (f._SYMBptr->sommet!=at_prod) 2652 return false; 2653 f=f._SYMBptr->feuille; 2654 if (f.type!=_VECT || f._VECTptr->size()!=2) 2655 return false; 2656 if (!is_linear_wrt(f._VECTptr->front(),x,c,e,contextptr) || is_zero(c,contextptr)) 2657 return false; 2658 gen a=-e/c; 2659 if (!is_linear_wrt(f._VECTptr->back(),x,c,e,contextptr) || is_zero(c,contextptr)) 2660 return false; 2661 gen b=-e/c; 2662 if (!is_greater(a,b,contextptr)) 2663 swapgen(a,b); 2664 res=sqrt(a,contextptr)+cs*sqrt(b,contextptr); 2665 if (!is_greater(a,b,contextptr)) 2666 *logptr(contextptr) << "Warning, assuming " << a-b << " is positive. Hint: run assume to make assumptions on a variable\n"; 2667 return true; 2668 } 2669 inner_sqrt(const gen & g)2670 static gen inner_sqrt(const gen & g){ 2671 //return g; // disabled 2672 if (g.type==_VECT){ 2673 vecteur v(*g._VECTptr); 2674 for (int i=0;i<v.size();++i) 2675 v[i]=inner_sqrt(v[i]); 2676 return gen(v,g.subtype); 2677 } 2678 if (g.type!=_SYMB) 2679 return g; 2680 gen f(inner_sqrt(g._SYMBptr->feuille)); 2681 if (g._SYMBptr->sommet!=at_prod || f.type!=_VECT) 2682 return symbolic(g._SYMBptr->sommet,f); 2683 const vecteur & v=*f._VECTptr; 2684 gen G=1; 2685 vecteur res; 2686 for (int i=0;i<v.size();++i){ 2687 if (v[i].is_symb_of_sommet(at_pow) && v[i]._SYMBptr->feuille[1]==plus_one_half){ 2688 gen vi=v[i]._SYMBptr->feuille[0]; 2689 if (vi.type<=_CPLX){ 2690 G=G*vi; 2691 continue; 2692 } 2693 vi=G*vi; 2694 G=1; 2695 res.push_back(symbolic(at_pow,makesequence(vi,plus_one_half))); 2696 } 2697 else 2698 res.push_back(v[i]); 2699 } 2700 if (!is_one(G)) 2701 res.push_back(symbolic(at_pow,makesequence(G,plus_one_half))); 2702 if (res.empty()) 2703 return 1; 2704 if (res.size()==1) 2705 return res.front(); 2706 return symbolic(at_prod,gen(res,_SEQ__VECT)); 2707 } 2708 in_normalize_sqrt(const gen & e,vecteur & L,bool keep_abs,GIAC_CONTEXT)2709 static gen in_normalize_sqrt(const gen & e,vecteur & L,bool keep_abs,GIAC_CONTEXT){ 2710 if (complex_mode(contextptr) || has_i(e)) 2711 return e; 2712 // remove multiple factors inside sqrt 2713 // vecteur l0=lop(e,at_pow),lin,lout; 2714 vecteur l0=lop_pow(e),lin,lout; 2715 const_iterateur it=l0.begin(),itend=l0.end(); 2716 for (;it!=itend;++it){ 2717 vecteur & arg=*it->_SYMBptr->feuille._VECTptr; 2718 gen g=arg[1],expnum,expden; 2719 if (g.type==_FRAC){ 2720 expnum=g._FRACptr->num; 2721 expden=g._FRACptr->den; 2722 } 2723 else { 2724 if ( (g.type!=_SYMB) || (g._SYMBptr->sommet!=at_prod) ) 2725 continue; 2726 gen & arg1=g._SYMBptr->feuille; 2727 if (arg1.type!=_VECT) 2728 continue; 2729 vecteur & v=*arg1._VECTptr; 2730 if ( (v.size()!=2) || (v[1].type!=_SYMB) || (v[1]._SYMBptr->sommet==at_inv) ) 2731 continue; 2732 expnum=v[0]; 2733 expden=v[1]._SYMBptr->feuille; 2734 } 2735 if ( (expden.type!=_INT_) || (expden.val!=2 ) ) 2736 continue; 2737 if (is_one(expnum) && is_integer(arg[0])) 2738 continue; 2739 // sqrt(arg[0]), we may check that arg[0] is a+b+/-2*sqrt(a*b) 2740 gen var,a,b,hyp; 2741 if (is_sqrtxy(arg[0],a,contextptr)){ 2742 lin.push_back(*it); 2743 lout.push_back(a); 2744 continue; 2745 } 2746 var=ggb_var(arg[0]); 2747 // if var is not assigned and arg[0] depends linearly on var, add an assumption 2748 if (complex_mode(contextptr)==false && var.type==_IDNT && var._IDNTptr->eval(1,var,contextptr)==var){ 2749 if (is_linear_wrt(arg[0],var,a,b,contextptr) && !is_zero(a)){ 2750 if (is_strictly_positive(a,contextptr)) 2751 hyp=symbolic(at_superieur_egal,makesequence(var,-b/a)); 2752 if (is_strictly_positive(-a,contextptr)) 2753 hyp=symbolic(at_inferieur_egal,makesequence(var,-b/a)); 2754 if (!is_zero(hyp)) 2755 L.push_back(hyp); 2756 } 2757 } 2758 vecteur lv(alg_lvar(arg[0])); 2759 gen num,den; 2760 a=e2r(arg[0],lv,contextptr); 2761 fxnd(a,num,den); 2762 gen nd=num*den,out(1); 2763 gen nover2=rdiv(expnum,plus_two,contextptr); 2764 if (nd.type==_EXT && nd._EXTptr->type==_VECT){ // extract content 2765 gen tmp=lgcd(*nd._EXTptr->_VECTptr); 2766 out=r2e(nd/tmp,lv,contextptr); 2767 // removed the ^ to avoid larger extensions e.g. for normal(sin(pi/5)) 2768 nd=tmp; 2769 } 2770 if (nd.type==_INT_ || nd.type==_ZINT){ 2771 gen simpl,doubl; 2772 bool pos; 2773 zint2simpldoublpos(nd,simpl,doubl,pos,2,contextptr); 2774 if (!pos) simpl=-simpl; 2775 lin.push_back(*it); 2776 gen tmparg=r2e(den,lv,contextptr); 2777 if (keep_abs) 2778 tmparg=abs(tmparg,contextptr); 2779 lout.push_back(pow(doubl/tmparg,expnum,contextptr)*pow(simpl*out,nover2,contextptr)); 2780 continue; 2781 } 2782 if (nd.type!=_POLY) 2783 continue; 2784 lin.push_back(*it); 2785 const factorization & f=rsqff(*nd._POLYptr); 2786 polynome s(plus_one,nd._POLYptr->dim),d(plus_one,s.dim); 2787 factorization::const_iterator jt=f.begin(),jtend=f.end(); 2788 for (;jt!=jtend;++jt){ 2789 if (jt->mult%2) 2790 s=s*jt->fact; 2791 d=d*pow(jt->fact,jt->mult/2); 2792 } 2793 // Extract integer content of s 2794 gen cont=Tppz<gen>(s); 2795 gen simpl,doubl; bool pos; 2796 zint2simpldoublpos(cont,simpl,doubl,pos,2,contextptr); 2797 if (!pos) simpl=-simpl; 2798 simpl=r2e(polynome(simpl,nd._POLYptr->dim),lv,contextptr); // if simpl is not an integer 2799 doubl=r2e(polynome(doubl,nd._POLYptr->dim),lv,contextptr); // if doubl is not an integer 2800 gen tmparg=r2e(den,lv,contextptr),tmpd=r2e(d,lv,contextptr); 2801 if (keep_abs){ 2802 tmparg=abs(tmparg,contextptr); 2803 tmpd=abs(tmpd,contextptr); 2804 } 2805 lout.push_back(pow(simpl*out,nover2,contextptr)*pow(doubl,expnum,contextptr)*pow(recursive_normal(r2e(s,lv,contextptr),contextptr),nover2,contextptr)*pow(tmpd,expnum,contextptr)*pow(tmparg,-expnum,contextptr)); 2806 } 2807 return subst(e,lin,lout,false,contextptr); 2808 } 2809 normalize_sqrt(const gen & e,GIAC_CONTEXT,bool keep_abs)2810 gen normalize_sqrt(const gen & e,GIAC_CONTEXT,bool keep_abs){ 2811 vecteur L; 2812 return in_normalize_sqrt(e,L,keep_abs,contextptr); 2813 } 2814 has_embedded_fractions(const gen & g)2815 static bool has_embedded_fractions(const gen & g){ 2816 if (g.type==_POLY){ 2817 vector< monomial<gen> >::const_iterator it=g._POLYptr->coord.begin(),itend=g._POLYptr->coord.end(); 2818 for (;it!=itend;++it){ 2819 if (has_embedded_fractions(it->value)) 2820 return true; 2821 } 2822 return false; 2823 } 2824 if (g.type==_VECT){ 2825 const_iterateur it=g._VECTptr->begin(),itend=g._VECTptr->end(); 2826 for (;it!=itend;++it){ 2827 if (has_embedded_fractions(*it)) 2828 return true; 2829 } 2830 return false; 2831 } 2832 if (g.type==_FRAC){ 2833 if (is_one(g._FRACptr->den)) 2834 return has_embedded_fractions(g._FRACptr->num); 2835 return true; 2836 } 2837 return false; 2838 } 2839 sort_func(const gen & a,const gen & b)2840 static bool sort_func(const gen & a,const gen & b){ 2841 if (a.type!=b.type) 2842 return a.type<b.type; 2843 if (a.type==_IDNT) 2844 return strcmp(a._IDNTptr->id_name,b._IDNTptr->id_name)<0; 2845 if (a.type==_SYMB){ 2846 int cmp=strcmp(a._SYMBptr->sommet.ptr()->s,b._SYMBptr->sommet.ptr()->s); 2847 if (cmp) return cmp<0; 2848 } 2849 return a.print(context0)<b.print(context0); 2850 } sort1(const vecteur & l)2851 static vecteur sort1(const vecteur & l){ 2852 vecteur res(l); 2853 gen_sort_f(res.begin(),res.end(),sort_func); 2854 return res; 2855 } sort0(vecteur & l)2856 static void sort0(vecteur & l){ 2857 // changed 9 dec 2019 for b:=rootof([[-1,-2,-4*a+3,16*a],[1,0,8*a-6,-32*a+8,16*a^2+24*a-3]])/(32*a^2-8*a)*atan(x/sqrt(-(-1+sqrt(1-4*a))/2))-rootof([[1,-2,4*a-3,16*a],[1,0,8*a-6,32*a-8,16*a^2+24*a-3]])/(32*a^2-8*a)*atan(x/sqrt(-(-1-sqrt(1-4*a))/2));normal(diff(b)); 2858 for (int i=1;i<l.size()-1;++i){ 2859 if (l[i].type==_VECT && l[i]._VECTptr->empty()){ 2860 l.erase(l.begin()+i); 2861 --i; 2862 } 2863 } 2864 iterateur it=l.begin(),itend=l.end(); 2865 for (;it!=itend;++it){ 2866 if (it->type==_VECT) 2867 *it=sort1(*it->_VECTptr); 2868 } 2869 } 2870 do_mult_i(const gen & g)2871 static bool do_mult_i(const gen & g){ 2872 if (g.type==_CPLX && is_zero(*g._CPLXptr) && is_positive(*(g._CPLXptr+1),context0)) 2873 return true; 2874 if (g.type==_POLY && do_mult_i(g._POLYptr->coord.front().value)) 2875 return true; 2876 if (g.type!=_EXT) 2877 return false; 2878 if (g._EXTptr->type!=_VECT || g._EXTptr->_VECTptr->empty()) 2879 return false; 2880 return do_mult_i(g._EXTptr->_VECTptr->front()); 2881 } 2882 purgeassumelist(const vecteur & L,GIAC_CONTEXT)2883 static void purgeassumelist(const vecteur & L,GIAC_CONTEXT){ 2884 for (unsigned k=0;k<L.size();++k){ 2885 gen Lk=ggb_var(L[k]); 2886 purgenoassume(Lk,contextptr); 2887 } 2888 } 2889 normal(const gen & e,bool distribute_div,GIAC_CONTEXT)2890 gen normal(const gen & e,bool distribute_div,GIAC_CONTEXT){ 2891 if (has_num_coeff(e)) 2892 return ratnormal(e,contextptr); 2893 // COUT << e << "\n"; 2894 if (e.type==_VECT){ 2895 vecteur res; 2896 for (const_iterateur it=e._VECTptr->begin();it!=e._VECTptr->end();++it) 2897 res.push_back(normal(*it,distribute_div,contextptr)); 2898 return gen(res,e.subtype); 2899 } 2900 if (e.type==_FRAC){ 2901 gen n=e._FRACptr->num; 2902 gen d=e._FRACptr->den; 2903 simplify(n,d); 2904 if (is_one(d)) 2905 return n; 2906 if (is_minus_one(d)) 2907 return -n; 2908 if (is_zero(d)){ 2909 if (is_zero(n)) 2910 return undef; 2911 else 2912 return unsigned_inf; 2913 } 2914 if (is_zero(n)) 2915 return zero; 2916 return fraction(n,d); 2917 } 2918 if (e.type==_MOD){ 2919 gen p=*e._MODptr; 2920 gen m=*(e._MODptr+1); 2921 vecteur l(lvar(m)); 2922 if (l.empty()){ 2923 l=lvar(p); 2924 p=e2r(p,l,contextptr); 2925 gen num,den; 2926 fxnd(p,num,den); 2927 num=makemod(num,m); 2928 if (!is_one(den)) 2929 den=makemod(den,m); 2930 p=rdiv(num,den,contextptr); 2931 return r2e(p,l,contextptr); 2932 } 2933 return _quorem(makesequence(p,m),contextptr)[1]; 2934 } 2935 if (e.type!=_SYMB) 2936 return e; 2937 if (is_equal(e)){ 2938 vecteur & v=*e._SYMBptr->feuille._VECTptr; 2939 return symbolic(at_equal,makesequence(normal(v.front(),distribute_div,contextptr),normal(v.back(),distribute_div,contextptr))); 2940 } 2941 if (is_inf(e) || is_undef(e) ) 2942 return e; 2943 gen ee,tmp; 2944 matrice l; 2945 fraction f(0); 2946 vecteur L; 2947 #ifndef NO_STDEXCEPT 2948 try { 2949 #endif 2950 ee=in_normalize_sqrt(e,L,true,contextptr); 2951 // put back integer content in sqrt, this may decrease the number 2952 // of extensions required, like for 2953 // assume(0<a<1/4);b:=regroup(int(1/(x^4+x^2+a)));normal(diff(b)); 2954 vecteur LL; gen EE=in_normalize_sqrt(inner_sqrt(e),LL,true,contextptr); 2955 if (lvar(EE).size()<lvar(ee).size()){ 2956 ee=EE; L=LL; 2957 } 2958 l=alg_lvar(ee); 2959 sort0(l); 2960 if (!L.empty() && debug_infolevel) 2961 *logptr(contextptr) << gettext("Making implicit assumption for sqrt argument ") << L << "\n"; 2962 if (contextptr && contextptr->previous) 2963 L.clear(); // otherwise ggbsort(x):=sort(x); r:=[sqrt(a)/a,1]; ggbsort(r); VARS(); keeps implicit assumption a>=0 globally 2964 for (unsigned k=0;k<L.size();++k) 2965 giac_assume(L[k],contextptr); 2966 tmp=e2r(ee,l,contextptr); 2967 if (is_undef(tmp)){ 2968 purgeassumelist(L,contextptr); 2969 #ifdef GIAC_GGB 2970 return undef; 2971 #endif 2972 if (calc_mode(contextptr)==1) 2973 return undef; 2974 *logptr(contextptr) << gettext("Unable to build a single algebraic extension for simplifying.\nTrying rational simplification only. This might return a wrong answer if simplifying 0/0!") << "\n"; 2975 l=lvar(ee); 2976 tmp=e2r(ee,l,contextptr); 2977 gen tmpf=evalf_double(ee-tmp,1,contextptr); 2978 if (tmpf.type==_DOUBLE_ && is_greater(tmpf,epsilon(contextptr),contextptr)) 2979 return undef; 2980 } 2981 #ifndef NO_STDEXCEPT 2982 } 2983 catch (std::runtime_error & err){ 2984 last_evaled_argptr(contextptr)=NULL; 2985 CERR << err.what() << "\n"; 2986 purgeassumelist(L,contextptr); 2987 return e; 2988 } 2989 #endif 2990 if (tmp.type==_FRAC){ 2991 f.num=tmp._FRACptr->num; 2992 f.den=tmp._FRACptr->den; 2993 } 2994 else 2995 f=tmp; 2996 if (f.den.type==_CPLX){ 2997 f.num=f.num*conj(f.den,contextptr); 2998 f.den=f.den*conj(f.den,contextptr); 2999 } 3000 if (do_mult_i(f.den)){ 3001 f.num = -cst_i*f.num; 3002 f.den = -cst_i*f.den; 3003 } 3004 if (do_mult_i(-f.den)){ 3005 f.num = cst_i*f.num; 3006 f.den = cst_i*f.den; 3007 } 3008 // search for embedded fractions 3009 if (has_embedded_fractions(f.num) || has_embedded_fractions(f.den)){ 3010 gen res=r2sym(f,l,contextptr); 3011 purgeassumelist(L,contextptr); 3012 return normal(res,distribute_div,contextptr); 3013 } 3014 if (distribute_div && f.num.type==_POLY && f.num._POLYptr->dim && f.den.type<_POLY){ 3015 gen res=r2sym(gen(*f.num._POLYptr/f.den),l,contextptr); 3016 purgeassumelist(L,contextptr); 3017 return res; 3018 } 3019 if (!distribute_div){ 3020 if (f.den.type==_POLY && is_positive(-f.den._POLYptr->coord.front())){ 3021 f.num=-f.num; 3022 f.den=-f.den; 3023 } 3024 if (f.num.type==_POLY && !f.num._POLYptr->coord.empty() && is_positive(-f.num._POLYptr->coord.front())){ 3025 f.num=-f.num; 3026 gen res=r2sym(f,l,contextptr); 3027 purgeassumelist(L,contextptr); 3028 return symbolic(at_neg,res); 3029 } 3030 } 3031 ee=r2sym(f,l,contextptr); 3032 purgeassumelist(L,contextptr); 3033 if (!is_one(f.den) && is_integer(f.den)){ 3034 ee=ratnormal(ratnormal(ee,contextptr),contextptr); // first ratnormal will expand sqrt()^ 3035 // second will remove them 3036 } 3037 return ee; // ratnormal(ee,contextptr); // for sqrt(1-a^2)/sqrt(1-a) 3038 } 3039 normal(const gen & e,GIAC_CONTEXT)3040 gen normal(const gen & e,GIAC_CONTEXT){ 3041 return normal(e,true,contextptr); 3042 } 3043 3044 // guess if g is a program/function: 1 func, 0 unknown, -1 not func is_program(const gen & g)3045 static int is_program(const gen & g){ 3046 #ifdef GIAC_HAS_STO_38 3047 if (g.type==_IDNT){ 3048 const char * idname=g._IDNTptr->id_name; 3049 if (strlen(idname)==2 && (idname[0]=='F' || idname[0]=='R' || idname[0]=='X' || idname[0]=='Y') && idname[1]>='0' && idname[1]<='9') 3050 return 1; 3051 else 3052 return -1; 3053 } 3054 #endif 3055 if (g.type==_FUNC) 3056 return 1; 3057 if (g.type==_VECT){ 3058 const_iterateur it=g._VECTptr->begin(),itend=g._VECTptr->end(); 3059 for (;it!=itend;++it){ 3060 if (it->type==_STRNG) 3061 return 0; 3062 int res=is_program(*it); 3063 if (res) 3064 return res; 3065 } 3066 return 0; 3067 } 3068 if (g.type!=_SYMB) 3069 return 0; 3070 if (g.is_symb_of_sommet(at_of) && g._SYMBptr->feuille.type==_VECT){ 3071 return is_program(g._SYMBptr->feuille._VECTptr->back()); 3072 } 3073 if (g.is_symb_of_sommet(at_program)) 3074 return 1; 3075 vecteur v(lvar(g)); 3076 if (v.size()==1 && v.front()==g){ 3077 if (g.type!=_SYMB) 3078 return 0; 3079 return is_program(g._SYMBptr->feuille); 3080 } 3081 for (unsigned i=0;i<v.size();++i){ 3082 int res=is_program(v[i]); 3083 if (res) 3084 return res; 3085 } 3086 return 0; 3087 } 3088 guess_program(gen & g,GIAC_CONTEXT)3089 bool guess_program(gen & g,GIAC_CONTEXT){ 3090 if (is_program(g)!=1) 3091 return false; 3092 g=eval(g,1,contextptr); 3093 return true; 3094 } 3095 is_algebraic_program(const gen & g,gen & f1,gen & f3)3096 bool is_algebraic_program(const gen & g,gen & f1,gen & f3){ 3097 if (g.type==_FUNC){ 3098 if (g==at_equal || g==at_superieur_strict || g==at_superieur_egal || g==at_inferieur_strict || g==at_inferieur_egal) 3099 return false; 3100 identificateur _tmpi(" x"); 3101 f1=_tmpi; 3102 f3=g(f1,context0); 3103 return true; 3104 } 3105 if (!g.is_symb_of_sommet(at_program) || g._SYMBptr->feuille.type!=_VECT || g._SYMBptr->feuille._VECTptr->size()!=3){ 3106 if (is_program(g)!=1) 3107 return false; 3108 identificateur _tmpi(" x"); 3109 f1=_tmpi; 3110 f3=g(f1,context0); 3111 return true; 3112 } 3113 f1=g._SYMBptr->feuille._VECTptr->front(); 3114 if (f1.type==_VECT && f1.subtype==_SEQ__VECT && f1._VECTptr->size()==1) 3115 f1=f1._VECTptr->front(); 3116 f3=g._SYMBptr->feuille._VECTptr->back(); 3117 bool res= ((f3.type==_SYMB || f3.type<=_IDNT || f3.type==_FRAC) && !f3.is_symb_of_sommet(at_bloc)); 3118 if (res) 3119 f3=eval(f3,1,context0); // eval operators like / 3120 return res; 3121 } has_algebraic_program(const gen & g)3122 bool has_algebraic_program(const gen & g){ 3123 if (g.type==_VECT){ 3124 vecteur & v=*g._VECTptr; 3125 for (unsigned i=0;i<v.size();++i){ 3126 if (has_algebraic_program(v[i])) 3127 return true; 3128 } 3129 return false; 3130 } 3131 if (g.type==_FUNC){ 3132 if (g==at_equal || g==at_superieur_strict || g==at_superieur_egal || g==at_inferieur_strict || g==at_inferieur_egal) 3133 return false; 3134 return true; 3135 } 3136 if (!g.is_symb_of_sommet(at_program) || g._SYMBptr->feuille.type!=_VECT || g._SYMBptr->feuille._VECTptr->size()!=3){ 3137 if (is_program(g)!=1) 3138 return false; 3139 return true; 3140 } 3141 gen f3=g._SYMBptr->feuille._VECTptr->back(); 3142 bool res= ((f3.type==_SYMB || f3.type<=_IDNT) && !f3.is_symb_of_sommet(at_bloc)); 3143 return res; 3144 } 3145 recursive_normal(const gen & e,bool distribute_div,GIAC_CONTEXT)3146 gen recursive_normal(const gen & e,bool distribute_div,GIAC_CONTEXT){ 3147 #ifdef TIMEOUT 3148 control_c(); 3149 #endif 3150 if (is_inequation(e) && e._SYMBptr->feuille.type==_VECT){ 3151 vecteur & v=*e._SYMBptr->feuille._VECTptr; 3152 unary_function_ptr u=e._SYMBptr->sommet; 3153 gen e2= v[0]-v[1]; 3154 e2=normal(e2,true,contextptr); 3155 gen c=_content(e2,contextptr); 3156 if (is_integer(c) || c.type==_FRAC){ 3157 e2=ratnormal(e2/c,contextptr); 3158 if (is_positive(-c,contextptr)) 3159 e2=-e2; 3160 } 3161 gen l=_lcoeff(e2,contextptr); 3162 if ( is_integer(l) || l.type==_FRAC){ 3163 e2=normal(e2/l,true,contextptr); 3164 if (e2.is_symb_of_sommet(at_plus) && e2._SYMBptr->feuille.type==_VECT && e2._SYMBptr->feuille._VECTptr->size()==2) 3165 e2=makesequence(e2._SYMBptr->feuille._VECTptr->front(),-e2._SYMBptr->feuille._VECTptr->back()); 3166 else 3167 e2=makesequence(e2,0); 3168 if (is_positive(-l,contextptr)){ 3169 if (u==at_superieur_strict) 3170 return symbolic(at_inferieur_strict,e2); 3171 if (u==at_superieur_egal) 3172 return symbolic(at_inferieur_egal,e2); 3173 if (u==at_inferieur_strict) 3174 return symbolic(at_superieur_strict,e2); 3175 if (u==at_inferieur_egal) 3176 return symbolic(at_superieur_egal,e2); 3177 return symbolic(u,e2); // should not happen 3178 } 3179 return symbolic(u,e2); 3180 } 3181 return symbolic(u,makesequence(e2,0)); 3182 } 3183 if (ctrl_c || interrupted) { 3184 gen res; 3185 interrupted = true; ctrl_c=false; 3186 gensizeerr(gettext("Stopped by user interruption."),res); 3187 return res; 3188 } 3189 if (e.type==_VECT) 3190 return apply(e,contextptr,(const gen_op_context) recursive_normal); 3191 gen e_copy(e); // was eval(e,1,contextptr)); 3192 //recursive BUG F(x):=int(1/sqrt(1+t^2),t,0,x);u(x):=exp(x); F(u(x)) 3193 if (e_copy.is_symb_of_sommet(at_pnt)) 3194 e_copy=e; 3195 if (e_copy.type==_FRAC) 3196 return e_copy._FRACptr->normal(); 3197 if (e_copy.type!=_SYMB && e_copy.type!=_MOD) 3198 return e_copy; 3199 if (is_inf(e_copy) || is_undef(e_copy) ) 3200 return e_copy; 3201 vecteur l=lvar(e_copy); 3202 vecteur l_subst(l); 3203 iterateur it=l_subst.begin(),itend=l_subst.end(); 3204 for (;it!=itend;++it){ 3205 if (it->type!=_SYMB) 3206 continue; 3207 if (it->_SYMBptr->sommet==at_when) 3208 continue; 3209 if (it->_SYMBptr->sommet!=at_pow){ 3210 gen tmp=it->_SYMBptr->feuille; 3211 if (!has_algebraic_program(tmp)) 3212 tmp=liste2symbolique(symbolique2liste(tmp,contextptr)); 3213 tmp=recursive_normal(tmp,false,contextptr); 3214 if (is_undef(tmp)) return e; 3215 *it=it->_SYMBptr->sommet(tmp,contextptr); 3216 continue; 3217 } 3218 if (it->_SYMBptr->feuille._VECTptr->front().type==_INT_ && it->_SYMBptr->feuille._VECTptr->back()==plus_one_half) 3219 continue; 3220 vecteur l=lvar(it->_SYMBptr->feuille._VECTptr->back()); 3221 int l_size; 3222 if (!l.empty() && l.front().type==_VECT) 3223 l_size=int(l.front()._VECTptr->size()); 3224 else 3225 l_size=int(l.size()); 3226 gen f,f_num,f_den; 3227 f=e2r(it->_SYMBptr->feuille._VECTptr->back(),l,contextptr); 3228 fxnd(f,f_num,f_den); 3229 gen num=r2sym(f_num,l,contextptr),den=r2sym(f_den,l,contextptr); 3230 gen base=recursive_normal(it->_SYMBptr->feuille._VECTptr->front(),false,contextptr); 3231 if ( is_zero(base) || num.type!=_INT_ || den.type!=_INT_ || den.val>MAX_COMMON_ALG_EXT_ORDER_SIZE) 3232 *it= pow(base,rdiv(num,den,contextptr),contextptr); 3233 else { 3234 if (den.val>1 && (is_integer(base) || base.type==_FRAC) && is_positive(base,contextptr)){ 3235 // should also handle more general base.type 3236 vecteur vtmp(den.val+1); 3237 vtmp[0]=1; 3238 vtmp.back()=-pow(base,num.val % den.val,contextptr); 3239 smaller_factor(vtmp); 3240 gen tmp; 3241 if (vtmp.size()>2) 3242 tmp=r2e(algebraic_EXTension(makevecteur(1,0),vtmp),vecteur(1,vecteur(0)),contextptr); 3243 else 3244 tmp=-vtmp.back()/vtmp.front(); 3245 *it=pow(base,num.val/den.val,contextptr)*tmp; 3246 } 3247 else 3248 *it= pow(base, num.val /den.val,contextptr) *pow(base,rdiv(num.val%den.val,den),contextptr); 3249 } 3250 } 3251 if (l!=l_subst) 3252 e_copy=subst(e_copy,l,l_subst,false,contextptr); 3253 // return global_eval(normal(e_copy),100); 3254 bool b=calc_mode(contextptr)==1 || abs_calc_mode(contextptr)==38; 3255 int ca=calc_mode(contextptr); 3256 calc_mode(0,contextptr); 3257 calc_mode(0,context0); 3258 int z=MPZ_MAXLOG2; 3259 MPZ_MAXLOG2=int(8e7); 3260 gen res=normal(e_copy,distribute_div,contextptr); 3261 MPZ_MAXLOG2=z; 3262 calc_mode(ca,context0); 3263 calc_mode(ca,contextptr); 3264 if ( b && !lop(res,at_rootof).empty()) 3265 res=simplifier(ratnormal(normalize_sqrt(e_copy,contextptr),contextptr),contextptr); 3266 return res; 3267 // removed eval since it eats neg(x-y) 3268 // eval(normal(e_copy,distribute_div),contextptr); 3269 } recursive_normal(const gen & e,GIAC_CONTEXT)3270 gen recursive_normal(const gen & e,GIAC_CONTEXT){ 3271 gen var,res; 3272 res=recursive_normal(e,true,contextptr); 3273 return res; 3274 } 3275 _recursive_normal(const gen & e,GIAC_CONTEXT)3276 gen _recursive_normal(const gen & e,GIAC_CONTEXT){ 3277 if (e.is_symb_of_sommet(at_unit)) 3278 return _usimplify(e,contextptr); 3279 gen var,res; 3280 if (is_equal(e)) 3281 return apply_to_equal(e,recursive_normal,contextptr); // symb_equal(_recursive_normal(equal2diff(e),contextptr),0); 3282 if (is_algebraic_program(e,var,res)) 3283 return symbolic(at_program,makesequence(var,0,recursive_normal(res,contextptr))); 3284 res=recursive_normal(e,true,contextptr); 3285 return res; 3286 } 3287 ratnormal(const gen & e,GIAC_CONTEXT)3288 gen ratnormal(const gen & e,GIAC_CONTEXT){ 3289 // COUT << e << "\n"; 3290 if (e.type==_VECT) 3291 return apply(e,ratnormal,contextptr); 3292 if (e.type==_FRAC){ 3293 gen n=e._FRACptr->num; 3294 gen d=e._FRACptr->den; 3295 simplify(n,d); 3296 if (is_one(d)) 3297 return n; 3298 if (is_minus_one(d)) 3299 return -n; 3300 if (is_zero(d)){ 3301 if (is_zero(n)) 3302 return undef; 3303 else 3304 return unsigned_inf; 3305 } 3306 if (is_zero(n)) 3307 return zero; 3308 return fraction(n,d); 3309 } 3310 if (e.type!=_SYMB && e.type!=_MOD) 3311 return e; 3312 if (is_inf(e) || is_undef(e) ) 3313 return e; 3314 matrice l; lvar(e,l); 3315 if (l.size()>1) l=sort1(l); 3316 gen fg=e2r(e,l,contextptr); // ok 3317 if (fg.type==_FRAC && fg._FRACptr->num.type==_FRAC){ 3318 fraction f(fg._FRACptr->num._FRACptr->num,fg._FRACptr->den*fg._FRACptr->num._FRACptr->den); 3319 f.normal(); 3320 return r2sym(f,l,contextptr); // ok 3321 } 3322 if (fg.type==_FRAC && fg._FRACptr->den.type==_CPLX){ 3323 gen tmp=conj(fg._FRACptr->den,contextptr); 3324 fg._FRACptr->num = fg._FRACptr->num * tmp; 3325 fg._FRACptr->den = fg._FRACptr->den * tmp; 3326 } 3327 return r2sym(fg,l,contextptr); 3328 } recursive_ratnormal(const gen & e,GIAC_CONTEXT)3329 gen recursive_ratnormal(const gen & e,GIAC_CONTEXT){ 3330 // COUT << e << "\n"; 3331 if (e.type==_VECT) 3332 return apply(e,recursive_ratnormal,contextptr); 3333 if (e.type==_FRAC){ 3334 gen n=e._FRACptr->num; 3335 gen d=e._FRACptr->den; 3336 simplify(n,d); 3337 if (is_one(d)) 3338 return n; 3339 if (is_minus_one(d)) 3340 return -n; 3341 if (is_zero(d)){ 3342 if (is_zero(n)) 3343 return undef; 3344 else 3345 return unsigned_inf; 3346 } 3347 if (is_zero(n)) 3348 return zero; 3349 return fraction(n,d); 3350 } 3351 if (e.type!=_SYMB && e.type!=_MOD) 3352 return e; 3353 if (is_inf(e) || is_undef(e) ) 3354 return e; 3355 matrice l; lvar(e,l); 3356 if (l.size()>1) l=sort1(l); 3357 gen fg=e2r(e,l,contextptr); 3358 for (unsigned i=0;i<l.size();++i){ 3359 if (l[i].type==_SYMB) 3360 l[i]=l[i]._SYMBptr->sommet(recursive_ratnormal(l[i]._SYMBptr->feuille,contextptr),contextptr); 3361 } 3362 if (fg.type==_FRAC && fg._FRACptr->num.type==_FRAC){ 3363 fraction f(fg._FRACptr->num._FRACptr->num,fg._FRACptr->den*fg._FRACptr->num._FRACptr->den); 3364 f.normal(); 3365 return ratnormal(r2sym(f,l,contextptr),contextptr); 3366 } 3367 if (fg.type==_FRAC && fg._FRACptr->den.type==_CPLX){ 3368 gen tmp=conj(fg._FRACptr->den,contextptr); 3369 fg._FRACptr->num = fg._FRACptr->num * tmp; 3370 fg._FRACptr->den = fg._FRACptr->den * tmp; 3371 } 3372 return ratnormal(r2sym(fg,l,contextptr),contextptr); 3373 } rationalgcd(const gen & a,const gen & b,GIAC_CONTEXT)3374 gen rationalgcd(const gen & a, const gen & b,GIAC_CONTEXT){ 3375 gen A,B,C,D; 3376 if (is_algebraic_program(a,A,B) && is_algebraic_program(b,C,D)){ 3377 if (A==C) 3378 return symbolic(at_program,makesequence(A,0,gcd(B,D,contextptr))); 3379 D=subst(D,C,A,false,contextptr); 3380 return symbolic(at_program,makesequence(A,0,gcd(B,D,contextptr))); 3381 } 3382 vecteur l(alg_lvar(a)); 3383 alg_lvar(b,l); 3384 fraction fa(e2r(a,l,contextptr)),fb(e2r(b,l,contextptr)); // ok 3385 if (debug_infolevel) 3386 CERR << "rational gcd begin " << CLOCK() << "\n"; 3387 if (!is_one(fa.den) || !is_one(fb.den)) 3388 CERR << "Warning gcd of fractions " << fa << " " << fb ; 3389 if (fa.num.type==_FRAC) 3390 fa.num=fa.num._FRACptr->num; 3391 if (fb.num.type==_FRAC) 3392 fb.num=fb.num._FRACptr->num; 3393 return r2sym(gcd(fa.num,fb.num,contextptr),l,contextptr); // ok 3394 } 3395 3396 static gen factor(const polynome & p,const vecteur &l,bool fixed_order,bool with_sqrt,gen divide_an_by,gen & extra_div,GIAC_CONTEXT); 3397 factor_multivar(const polynome & p,const vecteur & l,bool fixed_order,bool with_sqrt,gen divide_an_by,gen & extra_div,GIAC_CONTEXT)3398 static gen factor_multivar(const polynome & p,const vecteur &l,bool fixed_order,bool with_sqrt,gen divide_an_by,gen & extra_div,GIAC_CONTEXT){ 3399 polynome pp(p); 3400 int nvars=int(l.size()); 3401 if (l.front().type==_VECT) 3402 nvars=int(l.front()._VECTptr->size()); 3403 vector<int> deg(nvars); 3404 int mindeg=pp.lexsorted_degree(); 3405 int posmin=0; 3406 for (int i=1;i<nvars;i++){ 3407 int d=pp.degree(i); 3408 deg[i]=d; 3409 if (d<mindeg){ 3410 posmin=i; 3411 mindeg=d; 3412 } 3413 } 3414 // move posmin variable at the beginning for p and l 3415 vecteur lp; 3416 if (posmin && !fixed_order){ 3417 vecteur lptemp; 3418 pp.reorder(transposition(0,posmin,nvars)); 3419 vecteur::const_iterator it; 3420 if (l.front().type==_VECT){ 3421 it=l.front()._VECTptr->begin(); 3422 ++it; 3423 for (int i=1;i<nvars;i++,++it){ 3424 if (i==posmin){ 3425 if (lptemp.empty()) 3426 lptemp.push_back(*it); 3427 else 3428 lptemp.insert(lptemp.begin(),*it); 3429 lptemp.push_back(l.front()._VECTptr->front()); 3430 } 3431 else 3432 lptemp.push_back(*it); 3433 } 3434 lp=l; 3435 lp[0]=lptemp; 3436 } 3437 else { 3438 it=l.begin(); 3439 ++it; 3440 for (int i=1;i<nvars;i++,++it){ 3441 if (i==posmin){ 3442 if (lp.empty()) 3443 lp.push_back(*it); 3444 else 3445 lp.insert(lp.begin(),*it); 3446 lp.push_back(l.front()); 3447 } 3448 else 3449 lp.push_back(*it); 3450 } 3451 } 3452 } 3453 else 3454 lp=l; 3455 factorization v; 3456 polynome p_content(pp.dim); 3457 if (!factor(pp,p_content,v,false,with_sqrt,complex_mode(contextptr),divide_an_by,extra_div)) 3458 return gentypeerr(gettext("Not implemented, e.g. for multivariate mod/approx polynomials")); 3459 // factor p_content 3460 pp=p_content.trunc1(); 3461 vecteur ll; 3462 if (lp.front().type==_VECT){ 3463 ll=lp; 3464 ll[0]=cdr_VECT(*(ll[0]._VECTptr)); 3465 } 3466 else 3467 ll=cdr_VECT(lp); 3468 gen deno=1; lcmdeno(pp,deno); pp=deno*pp; deno=deno*extra_div; extra_div=1; 3469 gen tmp=factor(pp,ll,false,with_sqrt,1,extra_div,contextptr); 3470 tmp=tmp*r2sym(v,lp,contextptr)/r2sym(extra_div,ll,contextptr)/r2sym(deno,ll,contextptr); 3471 extra_div=1; 3472 return tmp; 3473 } 3474 factor(const polynome & p,const vecteur & l,bool fixed_order,bool with_sqrt,gen divide_an_by,gen & extra_div,GIAC_CONTEXT)3475 static gen factor(const polynome & p,const vecteur &l,bool fixed_order,bool with_sqrt,gen divide_an_by,gen & extra_div,GIAC_CONTEXT){ 3476 // find main var, make permutation on f.num, f.den and l 3477 // COUT << "Factor " << p << " " << l << "\n"; 3478 if (is_one(p)) 3479 return 1; 3480 if (l.empty() ) 3481 return r2sym(p,l,contextptr); 3482 if (!p.dim){ 3483 gen tmp; 3484 if (l.front().type==_VECT) 3485 tmp=r2sym(p,l.begin(),l.end(),contextptr); 3486 else 3487 tmp=r2sym(p,l,contextptr); 3488 return ratfactor(tmp,with_sqrt,contextptr); 3489 } 3490 if (p.dim>1) 3491 return factor_multivar(p,l,fixed_order,with_sqrt,divide_an_by,extra_div,contextptr); 3492 factorization v; 3493 polynome p_content(p.dim); 3494 if (!factor(p,p_content,v,false,with_sqrt,complex_mode(contextptr),divide_an_by,extra_div)) 3495 return gentypeerr(gettext("Not implemented, e.g. for multivariate mod/approx polynomials")); 3496 // factor p_content 3497 if (is_one(p_content)) 3498 return r2sym(v,l,contextptr)/r2sym(extra_div,l,contextptr); 3499 else { 3500 gen tmp(p_content); 3501 simplify(tmp,extra_div); 3502 if ( (tmp.type>_POLY || (tmp.type==_POLY && !tmp._POLYptr->coord.empty() && tmp._POLYptr->coord.front().value.type==_USER)) && 3503 !v.empty() && v.front().mult==1 && (v.size()==1 || v.front().fact.degree(0)==0) 3504 ){ // for GF factorization 3505 v.front().fact = tmp.type==_POLY?(*tmp._POLYptr*v.front().fact):(tmp*v.front().fact); 3506 return r2sym(v,l,contextptr)/r2sym(extra_div,l,contextptr); 3507 } 3508 polynome pp=p_content.trunc1(); 3509 vecteur ll; 3510 if (l.front().type==_VECT){ 3511 ll=l; 3512 ll[0]=cdr_VECT(*(ll[0]._VECTptr)); 3513 } 3514 else 3515 ll=cdr_VECT(l); 3516 tmp=factor(pp,ll,false,with_sqrt,1,extra_div,contextptr); //cerr << v << "\n"; 3517 return tmp*r2sym(v,l,contextptr)/r2sym(extra_div,l,contextptr); 3518 // return r2sym(tmp,l,contextptr)*r2sym(v,l,contextptr)/r2sym(extra_div,l,contextptr); 3519 } 3520 } 3521 var_factor(const gen & e,const vecteur & l,bool fixed_order,bool with_sqrt,const gen & divide_an_by,GIAC_CONTEXT)3522 static gen var_factor(const gen & e,const vecteur & l,bool fixed_order,bool with_sqrt,const gen & divide_an_by,GIAC_CONTEXT){ 3523 if (e.type!=_POLY) 3524 return r2sym(e/divide_an_by,l,contextptr); 3525 gen extra_div=1; 3526 return factor(*e._POLYptr,l,fixed_order,with_sqrt,divide_an_by,extra_div,contextptr); 3527 } 3528 ordered_factor(const gen & ee,vecteur & l,bool with_sqrt,GIAC_CONTEXT)3529 gen ordered_factor(const gen & ee,vecteur & l,bool with_sqrt,GIAC_CONTEXT){ 3530 gen e=normalize_sqrt(ee,contextptr); 3531 alg_lvar(e,l); 3532 gen f_num,f_den,f; 3533 f=e2r(e,l,contextptr); 3534 fxnd(f,f_num,f_den); 3535 return rdiv(var_factor(f_num,l,true,with_sqrt,1,contextptr),var_factor(f_den,l,true,with_sqrt,1,contextptr)); 3536 } 3537 factor(const gen & e,const identificateur & x,bool with_sqrt,GIAC_CONTEXT)3538 gen factor(const gen & e,const identificateur & x,bool with_sqrt,GIAC_CONTEXT){ 3539 if (e.type==_VECT){ 3540 vecteur w; 3541 vecteur::const_iterator it=e._VECTptr->begin(),itend=e._VECTptr->end(); 3542 for (;it!=itend;++it) 3543 w.push_back(factor(*it,x,with_sqrt,contextptr)); 3544 return w; 3545 } 3546 vecteur l(1,vecteur(1,x)); // insure x is the main var 3547 return ordered_factor(e,l,with_sqrt,contextptr); 3548 } 3549 in_factor(const gen & e_,bool with_sqrt,const gen & divide_an_by,GIAC_CONTEXT)3550 static gen in_factor(const gen & e_,bool with_sqrt,const gen & divide_an_by,GIAC_CONTEXT){ 3551 gen e(normalize_sqrt(e_,contextptr)); 3552 if (e.type==_VECT){ 3553 vecteur w; 3554 vecteur::const_iterator it=e._VECTptr->begin(),itend=e._VECTptr->end(); 3555 for (;it!=itend;++it) 3556 w.push_back(in_factor(*it,with_sqrt,divide_an_by,contextptr)); 3557 return w; 3558 } 3559 vecteur l; 3560 alg_lvar(e,l); 3561 if (!l.empty() && l.front().type==_VECT && l.front()._VECTptr->empty()) 3562 l=lvar(e); 3563 #ifdef POCKETCAS 3564 if (l.size()==1 && contains(l.front(),vx_var)){ // x first 3565 l=vecteur(1,vecteur(1,vx_var)); 3566 alg_lvar(e,l); 3567 } 3568 #endif 3569 /* if (calc_mode(contextptr)==1 && l.size()==1 && l.front().type==_VECT){ 3570 sort(l.front()._VECTptr->begin(),l.front()._VECTptr->end(),islesscomplexthanf); 3571 } */ 3572 gen f_num,f_den,f,dnum,dden(1); 3573 f=e2r(e,l,contextptr); 3574 fxnd(f,f_num,f_den); 3575 if (has_EXT(f_num)) 3576 simplify3(f_num,f_den); 3577 if (divide_an_by.type>=_IDNT){ 3578 gen divide=e2r(divide_an_by,l,contextptr); 3579 fxnd(divide,dnum,dden); 3580 if (dnum.type==_POLY){ 3581 if (!Tis_constant(*dnum._POLYptr)){ 3582 simplify3(f_num,dnum); 3583 } 3584 if (dnum.type==_POLY) 3585 dnum=dnum._POLYptr->coord.front().value; 3586 } 3587 if (dden.type==_POLY) 3588 dden=dden._POLYptr->coord.front().value; 3589 } 3590 else { 3591 dnum=divide_an_by; 3592 } 3593 if (dden.type==_CPLX){ 3594 gen tmp=conj(dden,contextptr); 3595 dnum=dnum*tmp; 3596 dden=dden*tmp; 3597 f_num=f_num*tmp; 3598 f_den=f_den*tmp; 3599 } 3600 if (dnum==-1){ 3601 f_num=-f_num; 3602 dnum=1; 3603 } 3604 gen N=var_factor(f_num,l,false,with_sqrt,dnum,contextptr); 3605 gen D=var_factor(f_den,l,false,with_sqrt,dden,contextptr); 3606 return rdiv(N,D); 3607 } 3608 factor(const gen & ee,bool with_sqrt,const gen & divide_an_by,GIAC_CONTEXT)3609 gen factor(const gen & ee,bool with_sqrt,const gen & divide_an_by,GIAC_CONTEXT){ 3610 if (xcas_mode(contextptr)==3 && is_integer(ee)) 3611 return _ifactor(ee,contextptr); 3612 gen e(ee); 3613 if (has_num_coeff(ee)) 3614 e=e.evalf(1,contextptr); 3615 else { 3616 vecteur L=lop_pow(e); 3617 L=lidnt(L); 3618 // make cfactor(-x/4*pi^2+i*sqrt(3+x)/2+x^2+3/4) work 3619 e=in_factor(e,with_sqrt && L.empty(),divide_an_by,contextptr); 3620 return e; 3621 } 3622 return in_factor(e,with_sqrt,divide_an_by,contextptr); 3623 } 3624 factor(const gen & ee,bool with_sqrt,GIAC_CONTEXT)3625 gen factor(const gen & ee,bool with_sqrt,GIAC_CONTEXT){ 3626 return factor(ee,with_sqrt,plus_one,contextptr); 3627 } 3628 ratfactor(const gen & ee,bool with_sqrt,GIAC_CONTEXT)3629 gen ratfactor(const gen & ee,bool with_sqrt,GIAC_CONTEXT){ 3630 gen e(normalize_sqrt(ee,contextptr)); 3631 if (has_num_coeff(ee)) 3632 e=e.evalf(1,contextptr); 3633 if (e.type==_VECT){ 3634 vecteur w; 3635 vecteur::const_iterator it=e._VECTptr->begin(),itend=e._VECTptr->end(); 3636 for (;it!=itend;++it) 3637 w.push_back(ratfactor(*it,with_sqrt,contextptr)); 3638 return w; 3639 } 3640 vecteur l; 3641 lvar(e,l); 3642 gen f_num,f_den,f; 3643 f=e2r(e,l,contextptr); 3644 fxnd(f,f_num,f_den); 3645 if (with_sqrt) 3646 l=vecteur(1,l); 3647 gen tmp=rdiv(var_factor(f_num,l,false,with_sqrt,1,contextptr),var_factor(f_den,l,false,with_sqrt,1,contextptr)); 3648 return tmp; 3649 } 3650 factor(const gen & ee,const gen & f,bool with_sqrt,GIAC_CONTEXT)3651 gen factor(const gen & ee,const gen & f,bool with_sqrt,GIAC_CONTEXT){ 3652 if (ee.type==_VECT){ 3653 vecteur & v=*ee._VECTptr; 3654 int s=int(v.size()); 3655 vecteur res(s); 3656 for (int i=0;i<s;++i) 3657 res[i]=factor(v[i],f,with_sqrt,contextptr); 3658 return res; 3659 } 3660 gen e(ee); 3661 if (has_num_coeff(ee)) 3662 e=e.evalf(1,contextptr); 3663 if (f.type==_IDNT) 3664 return factor(e,*f._IDNTptr,with_sqrt,contextptr); 3665 if (f.type==_VECT){ 3666 // should check that *f._VECTptr is made only of atomic _SYMB 3667 return ordered_factor(e,*f._VECTptr,with_sqrt,contextptr); 3668 } 3669 return gentypeerr(contextptr); 3670 } 3671 factorcollect(const gen & args,bool with_sqrt,GIAC_CONTEXT)3672 gen factorcollect(const gen & args,bool with_sqrt,GIAC_CONTEXT){ 3673 if (args.type!=_VECT) 3674 return factor(args,with_sqrt,contextptr); 3675 vecteur & v=*args._VECTptr; 3676 if (v.empty()) 3677 return gensizeerr(contextptr); 3678 if (v.size()==1) 3679 return vecteur(1,factor(v.front(),with_sqrt,contextptr)); 3680 if (args.subtype==_SEQ__VECT){ 3681 if (v.size()>2) 3682 toomanyargs("factor"); 3683 if (v.back()==at_sqrt) 3684 return factor(v.front(),true,contextptr); 3685 if (v.back().type!=_IDNT){ // FIXME could be improved! 3686 gen f=v.back(); 3687 if (v.back().type==_VECT) 3688 f=symbolic(at_prod,v.back()); 3689 gen res=factor(v.front()*f,with_sqrt,f,contextptr); 3690 return res; 3691 } 3692 return factor(v.front(),v.back(),with_sqrt,contextptr); 3693 } 3694 int s=int(v.size()); 3695 vecteur res(s); 3696 for (int i=0;i<s;++i) 3697 res[i]=factor(v[i],with_sqrt,contextptr); 3698 return res; 3699 } partfrac(const gen & e_,const vecteur & l,bool with_sqrt,GIAC_CONTEXT)3700 gen partfrac(const gen & e_,const vecteur & l,bool with_sqrt,GIAC_CONTEXT){ 3701 gen ext(1),e(e_); 3702 if (e.type==_VECT && e.subtype==_SEQ__VECT && e._VECTptr->size()==2){ 3703 ext=e._VECTptr->back(); 3704 e=e._VECTptr->front(); 3705 } 3706 if (e.type==_VECT){ 3707 vecteur w; 3708 vecteur::const_iterator it=e._VECTptr->begin(),itend=e._VECTptr->end(); 3709 for (;it!=itend;++it) 3710 w.push_back(partfrac(*it,l,with_sqrt,contextptr)); 3711 return w; 3712 } 3713 int l_size; 3714 gen xvar; 3715 if (!l.empty() && l.front().type==_VECT){ 3716 l_size=int(l.front()._VECTptr->size()); 3717 xvar=l.front(); 3718 } 3719 else { 3720 l_size=int(l.size()); 3721 xvar=l; 3722 } 3723 if (!l_size) 3724 return e; 3725 else 3726 xvar=xvar._VECTptr->front(); 3727 gen r=e2r(e,l,contextptr); 3728 gen r_num,r_den; 3729 fxnd(r,r_num,r_den); 3730 if (r_den.type!=_POLY){ 3731 if (r_num.type==_POLY) 3732 return rdiv(r2sym(r_num,l,contextptr),r2sym(r_den,l,contextptr)); 3733 else 3734 return e; 3735 } 3736 polynome f_den(*r_den._POLYptr),f_num(l_size); 3737 if (r_num.type==_POLY) 3738 f_num=*r_num._POLYptr; 3739 else 3740 f_num=polynome(r_num,l_size); 3741 factorization vden; 3742 polynome p_content(l_size); 3743 gen extra_div=1; 3744 if (ext!=1){ 3745 ext=e2r(ext,vecteur(1,vecteur(0)),contextptr); 3746 if (ext.type!=_EXT) 3747 ext=1; 3748 } 3749 factor(ext*f_den,p_content,vden,false,with_sqrt,complex_mode(contextptr),ext,extra_div); 3750 vector< pf<gen> > pfde_VECT; 3751 polynome ipnum(l_size),ipden(l_size); 3752 partfrac(f_num,f_den,vden,pfde_VECT,ipnum,ipden,!with_sqrt); 3753 gen res=rdiv(r2sym(gen(ipnum),l,contextptr),r2sym(gen(ipden),l,contextptr)); 3754 vector< pf<gen> > ::const_iterator it=pfde_VECT.begin(),itend=pfde_VECT.end(); 3755 for (;it!=itend;++it){ 3756 const pf<gen> & current=*it; 3757 gen reste(r2e(gen(current.num),l,contextptr)),deno(r2e(gen(current.fact),l,contextptr)); 3758 polynome p=it->mult==1?it->fact:pow(it->fact,it->mult),quo,rem; 3759 it->den.TDivRem(p,quo,rem,true); 3760 gen cur_deno(r2e(quo,l,contextptr)); 3761 if (current.mult==1){ 3762 // unitarize 3763 gen tmp(_lcoeff(makesequence(deno,xvar),contextptr)); 3764 reste=ratnormal(reste/tmp/cur_deno,contextptr); 3765 deno=recursive_normal(deno/tmp,contextptr); 3766 res += reste/deno; 3767 } 3768 else { 3769 for (int i=0;i<current.mult;++i){ 3770 gen tmp(_quorem(makesequence(reste,deno,xvar),contextptr)); 3771 if (tmp.type!=_VECT) 3772 return gensizeerr(contextptr); 3773 vecteur & vtmp=*tmp._VECTptr; 3774 reste=vtmp.front(); 3775 res=res+normal(vtmp.back()/cur_deno,contextptr)/pow(deno,current.mult-i); 3776 } 3777 } 3778 } 3779 return res; // +r2sym(pfde_VECT,l); 3780 } 3781 partfrac(const gen & e_,const identificateur & x,bool with_sqrt,GIAC_CONTEXT)3782 gen partfrac(const gen & e_,const identificateur & x,bool with_sqrt,GIAC_CONTEXT){ 3783 gen e=normalize_sqrt(e_,contextptr); 3784 vecteur l; 3785 l.push_back(x); // insure x is the main var 3786 l=vecteur(1,l); 3787 alg_lvar(e,l); 3788 return partfrac(e,l,with_sqrt,contextptr); 3789 } 3790 partfrac(const gen & e_,bool with_sqrt,GIAC_CONTEXT)3791 gen partfrac(const gen & e_,bool with_sqrt,GIAC_CONTEXT){ 3792 gen e=normalize_sqrt(e_,contextptr); 3793 vecteur l; 3794 alg_lvar(e,l); 3795 if (!l.empty() && l.front().type==_VECT && l.front()._VECTptr->empty()) 3796 return e_; 3797 if (l.size()==1 && contains(l.front(),vx_var)){ // x first 3798 l=vecteur(1,vecteur(1,vx_var)); 3799 alg_lvar(e,l); 3800 } 3801 return partfrac(e,l,with_sqrt,contextptr); 3802 } 3803 partfrac(const gen & e,const gen & f,bool with_sqrt,GIAC_CONTEXT)3804 gen partfrac(const gen & e,const gen & f,bool with_sqrt,GIAC_CONTEXT){ 3805 if (f.type==_IDNT) 3806 return partfrac(e,*f._IDNTptr,with_sqrt,contextptr); 3807 if (f.type==_VECT){ 3808 // should check that *f._VECTptr is made only of atomic _SYMB 3809 return partfrac(e,*f._VECTptr,with_sqrt,contextptr); 3810 } 3811 if (f.type==_SYMB) 3812 return partfrac(makesequence(e,f),with_sqrt,contextptr); 3813 return gentypeerr(contextptr); 3814 } 3815 3816 /* User functions */ symb2poly(const gen & fr,const gen & ba,GIAC_CONTEXT)3817 static gen symb2poly(const gen & fr,const gen & ba,GIAC_CONTEXT){ 3818 if (fr.type==_VECT) 3819 return apply1st(fr,ba,contextptr,symb2poly); 3820 if (ba.type==_VECT) 3821 return e2r(fr,*(ba._VECTptr),contextptr); 3822 if (ba.type!=_IDNT && ba.type!=_SYMB) 3823 return gensizeerr(contextptr); 3824 vecteur l(1,ba); 3825 lvar(fr,l); 3826 gen temp=e2r(fr,l,contextptr); 3827 l.erase(l.begin()); 3828 gen res; 3829 gen tmp2(polynome2poly1(temp,1)); 3830 //res=l.empty()?tmp2:r2e(tmp2,l,contextptr); 3831 res=l.empty()?tmp2:((tmp2.type==_FRAC && tmp2._FRACptr->den.type==_VECT && tmp2._FRACptr->den._VECTptr->size()>1)?gen(fraction(r2e(tmp2._FRACptr->num,l,contextptr),r2e(tmp2._FRACptr->den,l,contextptr))):r2e(tmp2,l,contextptr)); 3832 if (res.type==_FRAC && res._FRACptr->num.type==_VECT && res._FRACptr->den.type<_POLY){ 3833 res=inv(res._FRACptr->den,contextptr)*res._FRACptr->num; 3834 } 3835 if (res.type==_VECT && calc_mode(contextptr)!=1) 3836 res.subtype=_POLY1__VECT; 3837 return res; 3838 } _e2r(const gen & args,GIAC_CONTEXT)3839 gen _e2r(const gen & args,GIAC_CONTEXT){ 3840 if ( args.type==_STRNG && args.subtype==-1) return args; 3841 if (args.type!=_VECT) 3842 return _e2r(makesequence(args,vx_var),contextptr); 3843 vecteur & v=*args._VECTptr; 3844 int s=int(v.size()); 3845 if (s<2) 3846 return gendimerr(contextptr); 3847 gen res=v.front(); 3848 if (s==2 && v[1].type==_VECT) 3849 return symb2poly(res,v[1],contextptr); 3850 for (int i=1;i<s;++i){ 3851 res=symb2poly(res,v[i],contextptr); 3852 } 3853 if (res.type!=_VECT && res.type!=_FRAC && res.type!=_POLY) 3854 return gen(vecteur(1,res),calc_mode(contextptr)!=1?_POLY1__VECT:0); 3855 else 3856 return res; 3857 } 3858 3859 static const char _e2r_s []="e2r"; symb_e2r(const gen & arg1,const gen & arg2)3860 symbolic symb_e2r(const gen & arg1,const gen & arg2){ 3861 return symbolic(at_e2r,makesequence(arg1,arg2)); 3862 } 3863 static define_unary_function_eval (__e2r,&_e2r,_e2r_s); 3864 define_unary_function_ptr5( at_e2r ,alias_at_e2r,&__e2r,0,true); 3865 _r2e(const gen & args,GIAC_CONTEXT)3866 gen _r2e(const gen & args,GIAC_CONTEXT){ 3867 if ( args.type==_STRNG && args.subtype==-1) return args; 3868 if (args.type!=_VECT || args.subtype!=_SEQ__VECT) 3869 return _r2e(gen(makevecteur(args,vx_var),_SEQ__VECT),contextptr); 3870 vecteur & v=*args._VECTptr; 3871 int s=int(v.size()); 3872 if (s<2) 3873 return _r2e(gen(makevecteur(args,vx_var),_SEQ__VECT),contextptr); 3874 gen res=v[0]; 3875 for (int i=1;i<s;++i){ 3876 if (v[i].type==_VECT) 3877 res=r2e(res,*v[i]._VECTptr,contextptr); 3878 else { 3879 if (res.type==_VECT) 3880 res=horner(*res._VECTptr,v[i]); 3881 } 3882 } 3883 return res; 3884 } 3885 3886 static const char _r2e_s []="r2e"; symb_r2e(const gen & arg1,const gen & arg2)3887 symbolic symb_r2e(const gen & arg1,const gen & arg2){ 3888 return symbolic(at_r2e,makesequence(arg1,arg2)); 3889 } 3890 static define_unary_function_eval (__r2e,&_r2e,_r2e_s); 3891 define_unary_function_ptr5( at_r2e ,alias_at_r2e,&__r2e,0,true); 3892 3893 static const char _normal_s []="normal"; symb_normal(const gen & args)3894 symbolic symb_normal(const gen & args){ 3895 return symbolic(at_normal,args); 3896 } 3897 static define_unary_function_eval (__normal,(const gen_op_context)_recursive_normal,_normal_s); 3898 define_unary_function_ptr5( at_normal ,alias_at_normal,&__normal,0,true); 3899 3900 static const char _non_recursive_normal_s []="non_recursive_normal"; symb_non_recursive_normal(const gen & args)3901 symbolic symb_non_recursive_normal(const gen & args){ 3902 return symbolic(at_non_recursive_normal,args); 3903 } 3904 static define_unary_function_eval (__non_recursive_normal,(const gen_op_context)normal,_non_recursive_normal_s); 3905 define_unary_function_ptr5( at_non_recursive_normal ,alias_at_non_recursive_normal,&__non_recursive_normal,0,true); 3906 3907 3908 static const char _factor_s []="factor"; symb_factor(const gen & args)3909 symbolic symb_factor(const gen & args){ 3910 return symbolic(at_factor,args); 3911 } _factor(const gen & args,GIAC_CONTEXT)3912 gen _factor(const gen & args,GIAC_CONTEXT){ 3913 if ( args.type==_STRNG && args.subtype==-1) return args; 3914 if (is_integer(args)){ 3915 #ifdef KHICAS 3916 return _ifactor(args,contextptr); 3917 #else 3918 *logptr(contextptr) << "Run ifactor(" << args << ") for integer factorization." << "\n"; 3919 return args; 3920 #endif 3921 } 3922 if (args.is_symb_of_sommet(at_unit)) 3923 return mksa_reduce(args,contextptr); 3924 if (is_equal(args)) 3925 return apply_to_equal(args,_factor,contextptr); 3926 if (args.type==_VECT && args._VECTptr->size()==2 && is_equal(args._VECTptr->front())){ 3927 gen x=args._VECTptr->back(); 3928 gen a=_left(args._VECTptr->front(),contextptr); 3929 gen b=_right(args._VECTptr->front(),contextptr); 3930 return symb_equal(_factor(makesequence(a,x),contextptr),_factor(makesequence(b,x),contextptr)); 3931 } 3932 gen var,res; 3933 if (args.type!=_VECT && is_algebraic_program(args,var,res)) 3934 return symbolic(at_program,makesequence(var,0,_factor(res,contextptr))); 3935 if (xcas_mode(contextptr)==3) 3936 res=factorcollect(args,lvar(args).size()==1,contextptr); 3937 else 3938 res=factorcollect(args,withsqrt(contextptr),contextptr); 3939 return res; 3940 } 3941 static define_unary_function_eval (__factor,&_factor,_factor_s); 3942 define_unary_function_ptr5( at_factor ,alias_at_factor,&__factor,0,true); 3943 3944 static const char _collect_s []="collect"; _collect(const gen & args,GIAC_CONTEXT)3945 gen _collect(const gen & args,GIAC_CONTEXT){ 3946 if ( args.type==_STRNG && args.subtype==-1) return args; 3947 gen var,res; 3948 if (is_algebraic_program(args,var,res)) 3949 return symbolic(at_program,makesequence(var,0,_collect(res,contextptr))); 3950 if (is_equal(args)) 3951 return apply_to_equal(args,_collect,contextptr); 3952 if (args.type==_VECT && args.subtype==_SEQ__VECT && args._VECTptr->size()>=2&& args._VECTptr->front().type!=_VECT){ 3953 vecteur v(args._VECTptr->begin()+1,args._VECTptr->end()); 3954 res=_symb2poly(args,contextptr); 3955 if (res.type!=_FRAC){ 3956 res=_poly2symb(gen(mergevecteur(vecteur(1,res),v),_SEQ__VECT),contextptr); 3957 return res; 3958 } 3959 } 3960 res=factorcollect(args,false,contextptr); 3961 return res; 3962 } 3963 static define_unary_function_eval (__collect,&_collect,_collect_s); 3964 define_unary_function_ptr5( at_collect ,alias_at_collect,&__collect,0,true); 3965 3966 static const char _partfrac_s []="partfrac"; symb_partfrac(const gen & args)3967 symbolic symb_partfrac(const gen & args){ 3968 return symbolic(at_partfrac,args); 3969 } _partfrac(const gen & args_,GIAC_CONTEXT)3970 gen _partfrac(const gen & args_,GIAC_CONTEXT){ 3971 if (args_.type==_STRNG && args_.subtype==-1) return args_; 3972 gen args(args_),var,res; 3973 if (is_algebraic_program(args,var,res)) 3974 return symbolic(at_program,makesequence(var,0,_partfrac(res,contextptr))); 3975 if (is_equal(args)) 3976 return apply_to_equal(args,_partfrac,contextptr); 3977 args=exact(args,contextptr); 3978 if (args.type!=_VECT) 3979 return partfrac(args,withsqrt(contextptr),contextptr); 3980 if (args._VECTptr->size()>2) 3981 return gentoomanyargs(_partfrac_s); 3982 return partfrac(args._VECTptr->front(),args._VECTptr->back(),withsqrt(contextptr),contextptr); 3983 } 3984 static define_unary_function_eval (__partfrac,&_partfrac,_partfrac_s); 3985 define_unary_function_ptr5( at_partfrac ,alias_at_partfrac,&__partfrac,0,true); 3986 3987 static const char _resultant_s []="resultant"; symb_resultant(const gen & args)3988 symbolic symb_resultant(const gen & args){ 3989 return symbolic(at_resultant,args); 3990 } _resultant(const gen & args,GIAC_CONTEXT)3991 gen _resultant(const gen & args,GIAC_CONTEXT){ 3992 if ( args.type==_STRNG && args.subtype==-1) return args; 3993 if (args.type!=_VECT) 3994 return gentypeerr(contextptr); 3995 vecteur v =*args._VECTptr; 3996 int s=int(v.size()); 3997 if (s<2) 3998 toofewargs(_resultant_s); 3999 if (has_num_coeff(v)){ 4000 gen g=exact(args,contextptr); 4001 if (!has_num_coeff(g)){ 4002 g=_resultant(g,contextptr); 4003 return evalf(g,1,contextptr); 4004 } 4005 } 4006 if (v[0].type==_VECT && v[1].type==_VECT){ 4007 gen g1,g2; 4008 int t1=coefftype(*v[0]._VECTptr,g1),t2=coefftype(*v[1]._VECTptr,g2); 4009 double eps=epsilon(contextptr); 4010 if (t1==0 && t2==0){ 4011 if (eps>0) 4012 return mod_resultant(*v[0]._VECTptr,*v[1]._VECTptr,eps); 4013 gen res; 4014 subresultant(*v[0]._VECTptr,*v[1]._VECTptr,res); 4015 return res; 4016 } 4017 if (t1==_MOD || t2==_MOD){ 4018 gen m=t1==_MOD?*(g1._MODptr+1):*(g2._MODptr+1); 4019 modpoly A=unmod(*v[0]._VECTptr,m); 4020 modpoly B=unmod(*v[1]._VECTptr,m); 4021 gen res; 4022 if (ntlresultant(A,B,m,res)) 4023 return res; 4024 if (m.type==_INT_){ 4025 vector<int> a,b,tmp1,tmp2; 4026 vecteur2vector_int(A,m.val,a); 4027 vecteur2vector_int(B,m.val,b); 4028 return makemod(resultant_int(a,b,tmp1,tmp2,m.val),m); 4029 } 4030 #ifdef INT128 4031 if (m.type==_ZINT && sizeinbase2(m)<61){ 4032 longlong p=mpz_get_si(*m._ZINTptr); 4033 int n=giacmax(A.size(),B.size()); 4034 int l=sizeinbase2(n); 4035 if (((p-1)>>l)<<l==p-1){ 4036 vector<longlong> a,b,tmp1,tmp2; 4037 vecteur2vector_ll(A,p,a); 4038 vecteur2vector_ll(B,p,b); 4039 return makemod(resultantll(a,b,tmp1,tmp2,p),m); 4040 } 4041 } 4042 #endif 4043 return makemod(mod_resultant(A,B,eps),m); 4044 } 4045 // not very efficient... 4046 gen g(identificateur("tresultant")); 4047 v[0]=_poly2symb(makesequence(v[0],g),contextptr); 4048 v[1]=_poly2symb(makesequence(v[1],g),contextptr); 4049 v.insert(v.begin()+2,g); 4050 s++; 4051 } 4052 if (s==2){ 4053 if (v[0].type==_POLY && v[1].type==_POLY) 4054 return resultant(*v[0]._POLYptr,*v[1]._POLYptr); 4055 v.push_back(vx_var); 4056 } 4057 if (has_num_coeff(v)) 4058 return _det(_sylvester(args,contextptr),contextptr); 4059 if (v.back()==at_sylvester || v.back()==at_det) 4060 return _det(_sylvester(gen(vecteur(args._VECTptr->begin(),args._VECTptr->begin()+s-1),_SEQ__VECT),contextptr),contextptr); 4061 if (v.back()==at_lagrange && s>=3){ 4062 if (debug_infolevel) 4063 CERR << CLOCK()*1e-6 << " interp resultant begin " << "\n"; 4064 gen p=v[0]; 4065 gen q=v[1]; 4066 gen x=vx_var; 4067 if (s>3) 4068 x=v[2]; 4069 gen dbg; // dbg=_resultant(makesequence(p,q,x),contextptr); 4070 vecteur w(1,x); 4071 lvar(p,w); 4072 lvar(q,w); 4073 int m=_degree(makesequence(p,x),contextptr).val; 4074 int n=_degree(makesequence(q,x),contextptr).val; 4075 if (w.size()==1 && m+n>GIAC_PADIC){ 4076 gen S=_sylvester(makesequence(p,q,x),contextptr); 4077 return _det(S,contextptr); 4078 } 4079 if (w.size()>1){ 4080 gen y=w[1]; 4081 vecteur rargs(3,x); 4082 if (s>4 && equalposcomp(w,v[3])){ 4083 y=v[3]; 4084 if (s>5){ 4085 for (int i=4;i<s-1;++i) 4086 rargs.push_back(v[i]); 4087 rargs.push_back(at_lagrange); 4088 } 4089 else { 4090 if (m+n>=GIAC_PADIC) 4091 rargs.push_back(at_lagrange); 4092 } 4093 } 4094 else { 4095 if (m+n>=GIAC_PADIC) 4096 rargs.push_back(at_lagrange); 4097 } 4098 if (v[s-2]==at_det){ 4099 gen S=_sylvester(makesequence(p,q,x),contextptr); 4100 return _det(gen(makevecteur(S,at_lagrange),_SEQ__VECT),contextptr); 4101 } 4102 // bound on degree of resultant with respect to y 4103 int a=_total_degree(makesequence(p,makevecteur(x,y)),contextptr).val; 4104 int b=_total_degree(makesequence(q,makevecteur(x,y)),contextptr).val; 4105 // first estimate n*(a-m)+m*b 4106 int d1=n*(a-m)+m*b; 4107 // second estimate 4108 a=_degree(makesequence(p,y),contextptr).val; 4109 b=_degree(makesequence(q,y),contextptr).val; 4110 // a*n+b*m 4111 int d2=a*n+b*m; 4112 int d=giacmin(d1,d2); 4113 if (debug_infolevel) 4114 CERR << CLOCK()*1e-6 << " interp degree " << d << "\n"; 4115 vecteur X(d+1),Y(d+1); 4116 int j=-d/2; 4117 gen pl=_lcoeff(makesequence(p,x),contextptr); 4118 gen ql=_lcoeff(makesequence(q,x),contextptr); 4119 if (rargs.back()!=at_lagrange || w.size()==2){ 4120 // prepare p and q in polynomial format 4121 vecteur l; 4122 l.push_back(y); 4123 l.push_back(x); 4124 l=vecteur(1,l); 4125 alg_lvar(p,l); 4126 alg_lvar(q,l); 4127 gen f1,f1_num,f1_den,f2,f2_num,f2_den; 4128 f1=e2r(makevecteur(p,q),l,contextptr); 4129 f2=f1[1]; 4130 f1=f1[0]; 4131 fxnd(f1,f1_num,f1_den); 4132 fxnd(f2,f2_num,f2_den); 4133 if ( (f1_num.type==_POLY) && (f2_num.type==_POLY)){ 4134 const polynome & pp=*f1_num._POLYptr; 4135 const polynome & qp=*f2_num._POLYptr; 4136 gen coeff; 4137 if (!interpolable_resultant(pp,d,coeff,true,contextptr) || !interpolable_resultant(qp,d,coeff,true,contextptr)) 4138 return gensizeerr(gettext("Not enough elements in field to interpolate. Try in a field extension")); 4139 if (coeff.type==_USER) j=0; 4140 int dim=pp.dim; 4141 vecteur vp,vq,vp0,vq0; 4142 polynome2poly1(pp,1,vp); 4143 polynome2poly1(qp,1,vq); 4144 polynome pp0(pp); 4145 pp0.reorder(transposition(0,1,dim)); 4146 pp0=firstcoeff(pp0).trunc1(); 4147 polynome qp0(qp); 4148 qp0.reorder(transposition(0,1,dim)); 4149 qp0=firstcoeff(qp0).trunc1(); 4150 polynome2poly1(pp0,1,vp0); 4151 polynome2poly1(qp0,1,vq0); 4152 gen den=pow(r2sym(f1_den,l,contextptr),qp.lexsorted_degree(),contextptr)*pow(r2sym(f2_den,l,contextptr),pp.lexsorted_degree(),contextptr); 4153 vecteur tmpv1,tmpv2,S; 4154 for (int i=0;i<=d;++i,++j){ 4155 if (debug_infolevel && i%16==0) 4156 CERR << CLOCK()*1e-6 << " interp horner " << i << "\n"; 4157 gen xi; 4158 for (;;++j){ 4159 // find evaluation preserving degree in x 4160 if (0 && j==0) 4161 CERR << "j" << "\n"; 4162 xi=interpolate_xi(j,coeff); 4163 gen hp=horner(vp0,xi); 4164 gen hq=horner(vq0,xi); 4165 if (!is_zero(hp) && !is_zero(hq)){ 4166 break; 4167 } 4168 } 4169 X[i]=xi; 4170 gen gp=horner(vp,xi); 4171 gen gq=horner(vq,xi); 4172 if (debug_infolevel && i%16==0) 4173 CERR << CLOCK()*1e-6 << " interp resultant " << j << ", " << 100*double(i)/(d+1) << "% done" << "\n"; 4174 if (gp.type==_POLY && gq.type==_POLY){ 4175 if (0 && 4176 m>=GIAC_PADIC/2 && n>=GIAC_PADIC/2 4177 ) 4178 resultant_sylvester(*gp._POLYptr,*gq._POLYptr,tmpv1,tmpv2,S,Y[i]); 4179 else 4180 Y[i]=resultant(*gp._POLYptr,*gq._POLYptr); 4181 continue; 4182 } 4183 if (gp.type==_POLY){ 4184 Y[i]=pow(gq,gp._POLYptr->lexsorted_degree(),contextptr); 4185 continue; 4186 } 4187 if (gq.type==_POLY){ 4188 Y[i]=pow(gp,gq._POLYptr->lexsorted_degree(),contextptr); 4189 continue; 4190 } 4191 Y[i]=1; 4192 } 4193 if (debug_infolevel) 4194 CERR << CLOCK()*1e-6 << " interp dd " << "\n"; 4195 vecteur R=divided_differences(X,Y); 4196 if (debug_infolevel) 4197 CERR << CLOCK()*1e-6 << " interp build " << "\n"; 4198 modpoly resp(1,R[d]),tmp; // cst in y 4199 for (int i=d-1;i>=0;--i){ 4200 operator_times(resp,makevecteur(1,-X[i]),0,tmp); 4201 if (tmp.empty()) 4202 tmp=vecteur(1,R[i]); 4203 else 4204 tmp.back() += R[i]; 4205 tmp.swap(resp); 4206 } 4207 resp=trim(resp,0); 4208 vecteur & lf=*l.front()._VECTptr; 4209 lf.erase(lf.begin()); // remove y 4210 if (debug_infolevel) 4211 CERR << CLOCK()*1e-6 << " interp convert " << "\n"; 4212 gen resg=r2sym(resp,l,contextptr); 4213 if (resg.type==_VECT) 4214 resg=symb_horner(*resg._VECTptr,y); 4215 return resg/den; 4216 } 4217 } 4218 for (int i=0;i<=d;++i,++j){ 4219 if (debug_infolevel) 4220 CERR << CLOCK()*1e-6 << " interp resultant " << i << "\n"; 4221 for (;;++j){ 4222 // find evaluation preserving degree in x 4223 if (!is_zero(subst(pl,y,j,false,contextptr))&& !is_zero(subst(ql,y,j,false,contextptr))) 4224 break; 4225 } 4226 gen py=subst(p,y,j,false,contextptr); 4227 gen qy=subst(q,y,j,false,contextptr); 4228 X[i]=j; 4229 rargs[0]=py; 4230 rargs[1]=qy; 4231 Y[i]=_resultant(gen(rargs,_SEQ__VECT),contextptr); 4232 if (0){ 4233 gen dbgtst=subst(dbg,y,j,false,contextptr); 4234 if (!is_zero(ratnormal(dbgtst-Y[i],contextptr))) 4235 CERR << "err" << "\n"; 4236 } 4237 } 4238 if (debug_infolevel) 4239 CERR << CLOCK()*1e-6 << " interp interpolation begin " << "\n"; 4240 gen r=_lagrange(makesequence(X,Y,y),contextptr); 4241 if (debug_infolevel) 4242 CERR << CLOCK()*1e-6 << " interp interpolation end " << "\n"; 4243 return r; 4244 } 4245 } 4246 if (v.size()>3) 4247 return gentoomanyargs(_resultant_s); 4248 if (v.back().type==_MOD) 4249 v.back()=*v.back()._MODptr; 4250 if (v.back().is_symb_of_sommet(at_prod)){ 4251 const gen & f = v.back()._SYMBptr->feuille; 4252 if (f.type==_VECT && f._VECTptr->size()==2 && f._VECTptr->front().type==_MOD) 4253 v.back()=f._VECTptr->back(); 4254 } 4255 gen x=v.back(); 4256 vecteur l; 4257 l.push_back(x); 4258 l=vecteur(1,l); 4259 gen p1=v.front(),p2=v[1]; 4260 alg_lvar(p1,l); 4261 alg_lvar(p2,l); 4262 int l_size; 4263 if (!l.empty() && l.front().type==_VECT) 4264 l_size=int(l.front()._VECTptr->size()); 4265 else 4266 l_size=int(l.size()); 4267 gen f1,f1_num,f1_den,f2,f2_num,f2_den; 4268 f1=e2r(makevecteur(p1,p2),l,contextptr); 4269 f2=f1[1]; 4270 f1=f1[0]; 4271 fxnd(f1,f1_num,f1_den); 4272 fxnd(f2,f2_num,f2_den); 4273 if ( (f1_num.type==_POLY) && (f2_num.type==_POLY)) 4274 return r2sym(gen(resultant(*f1_num._POLYptr,*f2_num._POLYptr)),l,contextptr)/pow(r2sym(f1_den,l,contextptr),f2_num._POLYptr->lexsorted_degree(),contextptr)/pow(r2sym(f2_den,l,contextptr),f1_num._POLYptr->lexsorted_degree(),contextptr); 4275 if (is_zero(f1)) 4276 return f1; 4277 if (is_zero(f2)) 4278 return f2; 4279 return 1; 4280 } 4281 static define_unary_function_eval (__resultant,&_resultant,_resultant_s); 4282 define_unary_function_ptr5( at_resultant ,alias_at_resultant,&__resultant,0,true); 4283 4284 /* I/O on stream */ 4285 #ifdef NSPIRE 4286 template<class T> readargs_from_stream(nio::ios_base<T> & inf,vecteur & args,GIAC_CONTEXT)4287 void readargs_from_stream(nio::ios_base<T> & inf,vecteur & args,GIAC_CONTEXT) 4288 #else 4289 void readargs_from_stream(istream & inf,vecteur & args,GIAC_CONTEXT) 4290 #endif 4291 { 4292 string to_parse; 4293 char c; 4294 for (bool notbackslash=true;;){ 4295 inf.get(c); 4296 if (!inf 4297 #ifdef NSPIRE 4298 || c==char(EOF) 4299 #endif 4300 ) 4301 break; 4302 if ( notbackslash || (c!='\n') ) 4303 to_parse +=c; 4304 else 4305 to_parse = to_parse.substr(0,to_parse.size()-1); 4306 notbackslash= (c!='\\'); 4307 } 4308 gen e(to_parse,contextptr); 4309 if (e.type==_VECT) 4310 args=*e._VECTptr; 4311 else 4312 args=makevecteur(e); 4313 } 4314 4315 #ifdef NSPIRE 4316 template<class T> read1arg_from_stream(nio::ios_base<T> & inf,GIAC_CONTEXT)4317 gen read1arg_from_stream(nio::ios_base<T> & inf,GIAC_CONTEXT) 4318 #else 4319 gen read1arg_from_stream(istream & inf,GIAC_CONTEXT) 4320 #endif 4321 { 4322 string to_parse; 4323 char c; 4324 for (bool notbackslash=true;;){ 4325 inf.get(c); 4326 if (!inf) 4327 break; 4328 if ( notbackslash || (c!='\n') ) 4329 to_parse +=c; 4330 else 4331 to_parse = to_parse.substr(0,to_parse.size()-1); 4332 notbackslash= (c!='\\'); 4333 } 4334 return gen(to_parse,contextptr); 4335 } 4336 readargs(int ARGC,char * ARGV[],vecteur & args,GIAC_CONTEXT)4337 void readargs(int ARGC, char *ARGV[],vecteur & args,GIAC_CONTEXT){ 4338 #ifdef FXCG 4339 return; 4340 #else 4341 // first initialize random generator for factorizations 4342 srand(0); 4343 //srand(time(NULL)); 4344 string s; 4345 #ifdef NSPIRE 4346 if (ARGC==0){ // fake, just to instantiate for file 4347 file bidon("bidon","r"); 4348 readargs_from_stream(bidon,args,contextptr); 4349 } 4350 #endif 4351 if (ARGC==1) 4352 readargs_from_stream(CIN,args,contextptr); 4353 else { 4354 if (ARGC==2) { 4355 if (!secure_run 4356 #if !defined(__MINGW_H) && !defined(NSPIRE) 4357 && access(ARGV[1],R_OK)==0 4358 #endif 4359 ){ 4360 #if !defined GIAC_HAS_STO_38 && !defined NSPIRE && !defined FXCG && !defined POCKETCAS 4361 ifstream inf(ARGV[1]); 4362 readargs_from_stream(inf,args,contextptr); 4363 #endif 4364 } 4365 else { 4366 s=string(ARGV[1]); 4367 gen e(s,contextptr); 4368 args.push_back(e); 4369 } 4370 } 4371 else { 4372 vecteur v; 4373 for (int i=1;i<ARGC;i++){ 4374 s=string(ARGV[i]); 4375 gen e(s,contextptr); 4376 v.push_back(e); 4377 } 4378 args.push_back(v); 4379 } 4380 } 4381 if (args.empty()) 4382 args.push_back(gentypeerr(contextptr)); 4383 #endif 4384 } 4385 4386 #if 0 4387 /* Functions below are not used anymore */ 4388 // rewrite embedded fraction inside g as num/den 4389 static void reduce_alg_ext(const gen & g,gen & num,gen & den){ 4390 num=g; 4391 den=plus_one; 4392 int embeddings=0; 4393 vector<int> embeddings_s; 4394 for (;((num.type==_POLY) && (Tis_constant<gen>(*num._POLYptr)));++embeddings){ 4395 embeddings_s.push_back(num._POLYptr->dim); 4396 num=num._POLYptr->coord.front().value; 4397 } 4398 if (num.type==_EXT) 4399 num=ext_reduce(num); 4400 if (is_undef(num)) return; 4401 if (num.type==_FRAC){ 4402 den=num._FRACptr->den; 4403 num=num._FRACptr->num; 4404 } 4405 /* else commented since ext_reduce above might reduce g 4406 else { 4407 num=g; 4408 return; 4409 } */ 4410 for (;embeddings;--embeddings){ 4411 num=polynome(num,embeddings_s.back()); 4412 den=polynome(den,embeddings_s.back()); 4413 embeddings_s.pop_back(); 4414 } 4415 } 4416 #endif 4417 4418 // check in g all ext1 ext and rewrite them with ext2 remove_ext_copies(gen & g,const gen & ext1,const gen & ext2)4419 static void remove_ext_copies(gen & g,const gen & ext1,const gen & ext2){ 4420 if (g.type==_POLY){ 4421 vector<monomial<gen> >::iterator it=g._POLYptr->coord.begin(),itend=g._POLYptr->coord.end(); 4422 for (;it!=itend;++it) 4423 remove_ext_copies(it->value,ext1,ext2); 4424 return; 4425 } 4426 if (g.type==_VECT){ 4427 iterateur it=g._VECTptr->begin(),itend=g._VECTptr->end(); 4428 for (;it!=itend;++it) 4429 remove_ext_copies(*it,ext1,ext2); 4430 } 4431 if (g.type==_FRAC) 4432 remove_ext_copies(g._FRACptr->num,ext1,ext2); 4433 if (g.type==_EXT){ 4434 gen & ext=*(g._EXTptr+1); 4435 if (ext==ext1) 4436 ext=ext2; 4437 else 4438 remove_ext_copies(ext,ext1,ext2); 4439 } 4440 } 4441 4442 #ifndef NO_NAMESPACE_GIAC 4443 } // namespace giac 4444 #endif // ndef NO_NAMESPACE_GIAC 4445 4446