1 // -*- mode:C++ ; compile-command: "g++ -I.. -g -c lin.cc -DHAVE_CONFIG_H -DIN_GIAC" -*- 2 #include "giacPCH.h" 3 /* 4 * Copyright (C) 2000,2014 B. Parisse, Institut Fourier, 38402 St Martin d'Heres 5 * 6 * This program is free software; you can redistribute it and/or modify 7 * it under the terms of the GNU General Public License as published by 8 * the Free Software Foundation; either version 3 of the License, or 9 * (at your option) any later version. 10 * 11 * This program is distributed in the hope that it will be useful, 12 * but WITHOUT ANY WARRANTY; without even the implied warranty of 13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 14 * GNU General Public License for more details. 15 * 16 * You should have received a copy of the GNU General Public License 17 * along with this program. If not, see <http://www.gnu.org/licenses/>. 18 */ 19 using namespace std; 20 #include <stdexcept> 21 #include <cmath> 22 #include <cstdlib> 23 #include "sym2poly.h" 24 #include "usual.h" 25 #include "lin.h" 26 #include "subst.h" 27 #include "modpoly.h" 28 #include "prog.h" 29 #include "giacintl.h" 30 31 #ifndef NO_NAMESPACE_GIAC 32 namespace giac { 33 #endif // ndef NO_NAMESPACE_GIAC 34 35 // Should be rewritten with a map container for better efficiency! 36 contains(const gen & e,const unary_function_ptr & mys)37 bool contains(const gen & e,const unary_function_ptr & mys){ 38 if (e.type!=_SYMB) 39 return false; 40 if (e._SYMBptr->sommet==mys) 41 return true; 42 if (e._SYMBptr->feuille.type!=_VECT) 43 return contains(e._SYMBptr->feuille,mys); 44 vecteur::const_iterator it=e._SYMBptr->feuille._VECTptr->begin(),itend=e._SYMBptr->feuille._VECTptr->end(); 45 for (;it!=itend;++it) 46 if (contains(*it,mys)) 47 return true; 48 return false; 49 } 50 compress(vecteur & res,GIAC_CONTEXT)51 void compress(vecteur & res,GIAC_CONTEXT){ 52 if (res.size()==2) return; 53 vecteur v,w; 54 const_iterateur it=res.begin(),itend=res.end(); 55 v.reserve(itend-it); 56 w.reserve((itend-it)/2); 57 int pos; 58 for (;it!=itend;++it){ 59 pos=equalposcomp(w,*(it+1)); 60 if (pos){ 61 v[2*(pos-1)] = recursive_normal(v[2*(pos-1)] + *it,false,contextptr); 62 ++it; 63 } 64 else { 65 v.push_back(*it); 66 ++it; 67 w.push_back(*it); 68 v.push_back(*it); 69 } 70 } 71 swap(res,v); 72 } 73 74 // back conversion unlin(vecteur & v,GIAC_CONTEXT)75 gen unlin(vecteur & v,GIAC_CONTEXT){ 76 vecteur w; 77 gen coeff; 78 vecteur::const_iterator it=v.begin(),itend=v.end(); 79 for (;it!=itend;++it){ 80 coeff = *it; 81 ++it; 82 if (!is_zero(coeff)) 83 w.push_back(coeff*exp(*it,contextptr)); 84 } 85 if (w.empty()) 86 return 0; 87 return _plus(w,contextptr); 88 } 89 convolution(const gen & coeff,const gen & arg,const vecteur & w,vecteur & res,GIAC_CONTEXT)90 void convolution(const gen & coeff, const gen & arg,const vecteur & w,vecteur & res,GIAC_CONTEXT){ 91 vecteur::const_iterator it=w.begin(),itend=w.end(); 92 for (;it!=itend;++it){ 93 res.push_back(coeff*(*it)); 94 ++it; 95 res.push_back(recursive_normal(arg+(*it),false,contextptr)); 96 } 97 compress(res,contextptr); 98 } 99 convolution(const vecteur & v,const vecteur & w,vecteur & res,GIAC_CONTEXT)100 void convolution(const vecteur & v,const vecteur & w, vecteur & res,GIAC_CONTEXT){ 101 res.clear(); 102 res.reserve(res.size()+v.size()*w.size()/2); 103 vecteur::const_iterator it=v.begin(),itend=v.end(); 104 gen coeff; 105 for (;it!=itend;++it){ 106 coeff = *it; 107 ++it; 108 convolution(coeff,*it,w,res,contextptr); 109 } 110 } 111 convolutionpower(const vecteur & v,int k,vecteur & res,GIAC_CONTEXT)112 void convolutionpower(const vecteur & v,int k,vecteur & res,GIAC_CONTEXT){ 113 res.clear(); 114 // should be improved for efficiency! 115 if (k==0){ 116 res.push_back(1); 117 res.push_back(0); 118 return; 119 } 120 if (k==1){ 121 res=v; 122 return; 123 } 124 convolutionpower(v,k/2,res,contextptr); 125 vecteur tmp=res; 126 convolution(tmp,tmp,res,contextptr); 127 if (k%2){ 128 tmp=res; 129 convolution(tmp,v,res,contextptr); 130 } 131 } 132 133 // coeff & argument of exponential lin(const gen & e,vecteur & v,GIAC_CONTEXT)134 void lin(const gen & e,vecteur & v,GIAC_CONTEXT){ 135 if (e.type!=_SYMB){ 136 v.push_back(e); 137 v.push_back(0); 138 return ; // e*exp(0) 139 } 140 // e is symbolic, look for exp, cosh, sinh, +, *, neg and inv, ^ 141 unary_function_ptr s=e._SYMBptr->sommet; 142 if ((s==at_plus) && (e._SYMBptr->feuille.type==_VECT)){ 143 vecteur::const_iterator it=e._SYMBptr->feuille._VECTptr->begin(),itend=e._SYMBptr->feuille._VECTptr->end(); 144 for (;it!=itend;++it) 145 lin(*it,v,contextptr); 146 compress(v,contextptr); 147 return; 148 } 149 if (s==at_neg){ 150 vecteur tmp; 151 lin(e._SYMBptr->feuille,tmp,contextptr); 152 const_iterateur it=tmp.begin(),itend=tmp.end(); 153 for (;it!=itend;++it){ 154 v.push_back(-*it); 155 ++it; 156 v.push_back(*it); 157 } 158 return; 159 } 160 if (s==at_inv){ 161 vecteur w; 162 lin(e._SYMBptr->feuille,w,contextptr); 163 if (w.size()==2){ 164 v.push_back(inv(w[0],contextptr)); 165 v.push_back(-w[1]); 166 } 167 else { 168 gen coeff(unlin(w,contextptr)); 169 v.push_back(inv(coeff,contextptr)); 170 v.push_back(0); 171 } 172 return ; 173 } 174 if (s==at_prod){ 175 if (e._SYMBptr->feuille.type!=_VECT){ 176 lin(e._SYMBptr->feuille,v,contextptr); 177 return; 178 } 179 vecteur w; 180 vecteur::const_iterator it=e._SYMBptr->feuille._VECTptr->begin(),itend=e._SYMBptr->feuille._VECTptr->end(); 181 lin(*it,w,contextptr); 182 ++it; 183 for (;it!=itend;++it){ 184 vecteur v0; 185 lin(*it,v0,contextptr); 186 vecteur res; 187 convolution(w,v0,res,contextptr); 188 w=res; 189 } 190 v=mergevecteur(v,w); 191 return; 192 } 193 if (s==at_pow){ 194 vecteur::const_iterator it=e._SYMBptr->feuille._VECTptr->begin(); 195 vecteur w; 196 lin(*it,w,contextptr); 197 ++it; 198 if (w.size()==2){ 199 if ( is_zero(w[1]) && (w[0].type==_INT_) ){ 200 w[1]=ln(w[0],contextptr); 201 w[0]=plus_one; 202 } 203 v.push_back(pow(w[0],*it,contextptr)); 204 v.push_back(w[1]*(*it)); 205 return ; 206 } 207 if ((it->type==_INT_) && (it->val>=0)){ 208 vecteur z(w),tmp; 209 convolutionpower(z,it->val,tmp,contextptr); 210 v=mergevecteur(v,tmp); 211 compress(v,contextptr); 212 return ; 213 } 214 gen coeff=unlin(w,contextptr); 215 v.push_back(pow(coeff,*it,contextptr)); 216 v.push_back(0); 217 return ; 218 } 219 gen f=_lin(e._SYMBptr->feuille,contextptr); 220 if (s==at_exp){ 221 v.push_back(1); 222 v.push_back(f); 223 return ; // 1*exp(arg) 224 } 225 if (s==at_cosh || s==at_sinh){ 226 v.push_back(plus_one_half); 227 v.push_back(f); 228 v.push_back(s==at_cosh?plus_one_half:minus_one_half); 229 v.push_back(-f); 230 return ; // 1/2*exp(arg)+-1/2*exp(-arg) 231 } 232 v.push_back(symbolic(s,f)); 233 v.push_back(0); 234 } 235 symb_lin(const gen & a)236 symbolic symb_lin(const gen & a){ 237 return symbolic(at_lin,a); 238 } 239 240 // "unary" version _lin(const gen & args,GIAC_CONTEXT)241 gen _lin(const gen & args,GIAC_CONTEXT){ 242 if ( args.type==_STRNG && args.subtype==-1) return args; 243 gen var,res; 244 if (is_algebraic_program(args,var,res)) 245 return symbolic(at_program,makesequence(var,0,_lin(res,contextptr))); 246 if (is_equal(args)) 247 return apply_to_equal(args,_lin,contextptr); 248 vecteur v; 249 if (args.type!=_VECT){ 250 lin(args,v,contextptr); 251 return unlin(v,contextptr); 252 } 253 return apply(args,_lin,contextptr); 254 } 255 static const char _lin_s []="lin"; 256 static define_unary_function_eval (__lin,&_lin,_lin_s); 257 define_unary_function_ptr5( at_lin ,alias_at_lin,&__lin,0,true); 258 259 static const char _lineariser_s []="lineariser"; 260 static define_unary_function_eval (__lineariser,&_lin,_lineariser_s); 261 define_unary_function_ptr5( at_lineariser ,alias_at_lineariser,&__lineariser,0,true); 262 263 // back conversion tunlin(vecteur & v,GIAC_CONTEXT)264 gen tunlin(vecteur & v,GIAC_CONTEXT){ 265 vecteur w; 266 gen coeff; 267 vecteur::const_iterator it=v.begin(),itend=v.end(); 268 for (;it!=itend;++it){ 269 coeff = *it; 270 ++it; 271 coeff=coeff*(*it); 272 if (!is_zero(coeff)) 273 w.push_back(coeff); 274 } 275 if (w.empty()) 276 return 0; 277 if (w.size()==1) 278 return w.front(); 279 return symbolic(at_plus,gen(w,_SEQ__VECT)); 280 } 281 tadd(vecteur & res,const gen & coeff,const gen & angle,GIAC_CONTEXT)282 static void tadd(vecteur & res,const gen & coeff,const gen & angle,GIAC_CONTEXT){ 283 gen newangle=angle,newcoeff=coeff; 284 if ( (newangle.type==_SYMB) && (newangle._SYMBptr->sommet==at_neg)){ 285 newcoeff=-coeff; 286 newangle=-newangle; 287 } 288 if ( (newangle.type==_SYMB) && ( (newangle._SYMBptr->sommet==at_sin) || (newangle._SYMBptr->sommet==at_cos) ) ){ 289 res.push_back(newcoeff); 290 res.push_back(newangle); 291 } 292 else { 293 newcoeff=newcoeff*newangle; 294 if (!is_zero(newcoeff)){ 295 res.push_back(newcoeff); 296 res.push_back(plus_one); 297 } 298 } 299 } 300 tconvolution(const gen & coeff,const gen & arg,const vecteur & w,vecteur & res,GIAC_CONTEXT)301 void tconvolution(const gen & coeff, const gen & arg,const vecteur & w,vecteur & res,GIAC_CONTEXT){ 302 gen newcoeff,tmp; 303 if ((arg.type==_SYMB) && (arg._SYMBptr->sommet==at_cos)){ 304 vecteur::const_iterator it=w.begin(),itend=w.end(); 305 for (;it!=itend;++it){ 306 newcoeff=coeff*(*it); 307 ++it; 308 bool iscos=it->type==_SYMB && it->_SYMBptr->sommet==at_cos; 309 if ( (it->type==_SYMB) && (iscos || it->_SYMBptr->sommet==at_sin) ){ 310 newcoeff=recursive_normal(rdiv(newcoeff,plus_two,contextptr),false,contextptr); 311 gen tmp1(normal(it->_SYMBptr->feuille+arg._SYMBptr->feuille,false,contextptr)); 312 gen tmp2(normal(it->_SYMBptr->feuille-arg._SYMBptr->feuille,false,contextptr)); 313 tadd(res,newcoeff,iscos?cos(tmp2,contextptr):sin(tmp1,contextptr),contextptr); 314 tadd(res,newcoeff,iscos?cos(tmp1,contextptr):sin(tmp2,contextptr),contextptr); 315 continue; 316 } 317 res.push_back(recursive_normal(newcoeff*(*it),false,contextptr)); 318 res.push_back(arg); 319 } 320 compress(res,contextptr); 321 return; 322 } 323 if ((arg.type==_SYMB) && (arg._SYMBptr->sommet==at_sin)){ 324 vecteur::const_iterator it=w.begin(),itend=w.end(); 325 for (;it!=itend;++it){ 326 newcoeff=coeff*(*it); 327 ++it; 328 bool iscos=it->type==_SYMB && it->_SYMBptr->sommet==at_cos; 329 if ( (it->type==_SYMB) && (iscos || it->_SYMBptr->sommet==at_sin) ){ 330 newcoeff=recursive_normal(rdiv(newcoeff,plus_two,contextptr),false,contextptr); 331 gen tmp1(normal(arg._SYMBptr->feuille+it->_SYMBptr->feuille,false,contextptr)); 332 gen tmp2(normal(arg._SYMBptr->feuille-it->_SYMBptr->feuille,false,contextptr)); 333 if (iscos){ 334 tadd(res,newcoeff,sin(tmp1,contextptr),contextptr); 335 tadd(res,newcoeff,sin(tmp2,contextptr),contextptr); 336 } 337 else { 338 tadd(res,newcoeff,cos(tmp2,contextptr),contextptr); 339 tadd(res,-newcoeff,cos(tmp1,contextptr),contextptr); 340 } 341 continue; 342 } 343 res.push_back(recursive_normal(newcoeff*(*it),false,contextptr)); 344 res.push_back(arg); 345 } 346 compress(res,contextptr); 347 return; 348 } 349 const_iterateur it=w.begin(),itend=w.end(); 350 newcoeff=coeff*arg; 351 for (;it!=itend;++it){ 352 res.push_back(recursive_normal(*it*newcoeff,false,contextptr)); 353 ++it; 354 res.push_back(*it); 355 } 356 } 357 tconvolution(const vecteur & v,const vecteur & w,vecteur & res,GIAC_CONTEXT)358 void tconvolution(const vecteur & v,const vecteur & w, vecteur & res,GIAC_CONTEXT){ 359 res.clear(); 360 res.reserve(res.size()+v.size()*w.size()/2); 361 vecteur::const_iterator it=v.begin(),itend=v.end(); 362 gen coeff; 363 for (;it!=itend;++it){ 364 coeff = *it; 365 ++it; 366 tconvolution(coeff,*it,w,res,contextptr); 367 } 368 } 369 tconvolutionpower(const vecteur & v,int k,vecteur & res,GIAC_CONTEXT)370 void tconvolutionpower(const vecteur & v,int k,vecteur & res,GIAC_CONTEXT){ 371 res.clear(); 372 // should be improved for efficiency! 373 if (k==0){ 374 res.push_back(1); 375 res.push_back(1); 376 return; 377 } 378 if (k==1){ 379 res=v; 380 return; 381 } 382 tconvolutionpower(v,k/2,res,contextptr); 383 vecteur tmp=res; 384 tconvolution(tmp,tmp,res,contextptr); 385 if (k%2){ 386 tmp=res; 387 tconvolution(tmp,v,res,contextptr); 388 } 389 } 390 391 // coeff & argument of sin/cos tlin(const gen & e,vecteur & v,GIAC_CONTEXT)392 void tlin(const gen & e,vecteur & v,GIAC_CONTEXT){ 393 if (e.type!=_SYMB){ 394 v.push_back(e); 395 v.push_back(1); 396 return ; // e*1 397 } 398 // e is symbolic, look for cos, sin, +, *, neg and inv, ^ 399 unary_function_ptr s=e._SYMBptr->sommet; 400 if ( (s==at_cos) || (s==at_sin)){ 401 v.push_back(1); 402 v.push_back(e); 403 return ; 404 } 405 if ((s==at_plus) && (e._SYMBptr->feuille.type==_VECT)){ 406 vecteur::const_iterator it=e._SYMBptr->feuille._VECTptr->begin(),itend=e._SYMBptr->feuille._VECTptr->end(); 407 for (;it!=itend;++it) 408 tlin(*it,v,contextptr); 409 compress(v,contextptr); 410 return; 411 } 412 if (s==at_neg){ 413 vecteur w; 414 tlin(e._SYMBptr->feuille,w,contextptr); 415 const_iterateur it=w.begin(),itend=w.end(); 416 for (;it!=itend;++it){ 417 v.push_back(-*it); 418 ++it; 419 v.push_back(*it); 420 } 421 return; 422 } 423 if (s==at_prod){ 424 if (e._SYMBptr->feuille.type!=_VECT){ 425 tlin(e._SYMBptr->feuille,v,contextptr); 426 return; 427 } 428 vecteur::const_iterator it=e._SYMBptr->feuille._VECTptr->begin(),itend=e._SYMBptr->feuille._VECTptr->end(); 429 vecteur w; 430 tlin(*it,w,contextptr); 431 ++it; 432 for (;it!=itend;++it){ 433 vecteur v0; 434 tlin(*it,v0,contextptr); 435 vecteur res; 436 tconvolution(w,v0,res,contextptr); 437 w=res; 438 } 439 v=mergevecteur(v,w); 440 return; 441 } 442 if (s==at_pow){ 443 vecteur::const_iterator it=e._SYMBptr->feuille._VECTptr->begin(); 444 /* 445 if ( (v.size()==2) && ((it+1)->type==_INT_) && ((it+1)->val>=0) ){ 446 tlin(*it,v); 447 ++it; 448 return tpow(v,it->val); 449 } 450 */ 451 if (((it+1)->type==_INT_) && ((it+1)->val>=0)){ 452 vecteur w; 453 tlin(*it,w,contextptr); 454 vecteur z(w); 455 ++it; 456 tconvolutionpower(z,it->val,w,contextptr); 457 v=mergevecteur(v,w); 458 return ; 459 } 460 } 461 gen te=_tlin(e._SYMBptr->feuille,contextptr); 462 if (s==at_pow && te.type==_VECT && te._VECTptr->size()==2 && te._VECTptr->back()==plus_one_half) 463 v.push_back(sqrt(te._VECTptr->front(),contextptr)); 464 else 465 v.push_back(s(te,contextptr)); 466 v.push_back(1); 467 } 468 symb_tlin(const gen & a)469 symbolic symb_tlin(const gen & a){ 470 return symbolic(at_tlin,a); 471 } 472 473 // "unary" version _tlin(const gen & args,GIAC_CONTEXT)474 gen _tlin(const gen & args,GIAC_CONTEXT){ 475 if ( args.type==_STRNG && args.subtype==-1) return args; 476 gen var,res; 477 if (is_algebraic_program(args,var,res)) 478 return symbolic(at_program,makesequence(var,0,_tlin(res,contextptr))); 479 if (is_equal(args)) 480 return apply_to_equal(args,_tlin,contextptr); 481 vecteur v; 482 if (args.type!=_VECT){ 483 tlin(args,v,contextptr); 484 return tunlin(v,contextptr); 485 } 486 return apply(args,_tlin,contextptr); 487 } 488 static const char _tlin_s []="tlin"; 489 static define_unary_function_eval (__tlin,&_tlin,_tlin_s); 490 define_unary_function_ptr5( at_tlin ,alias_at_tlin,&__tlin,0,true); 491 492 static const char _lineariser_trigo_s []="lineariser_trigo"; 493 static define_unary_function_eval (__lineariser_trigo,&_tlin,_lineariser_trigo_s); 494 define_unary_function_ptr5( at_lineariser_trigo ,alias_at_lineariser_trigo,&__lineariser_trigo,0,true); 495 496 // Expand and texpand split(const gen & e,gen & coeff,gen & arg,GIAC_CONTEXT)497 static void split(const gen & e, gen & coeff, gen & arg,GIAC_CONTEXT){ 498 if (e.type==_INT_){ 499 coeff=e; 500 arg=plus_one; 501 } 502 if ( (e.type==_SYMB) && (e._SYMBptr->sommet==at_neg)){ 503 split(e._SYMBptr->feuille,coeff,arg,contextptr); 504 coeff=-coeff; 505 return; 506 } 507 if ( (e.type!=_SYMB) || (e._SYMBptr->sommet!=at_prod) ) { 508 coeff=plus_one; 509 arg=e; 510 return; 511 } 512 coeff = plus_one; 513 arg = plus_one; 514 const_iterateur it=e._SYMBptr->feuille._VECTptr->begin(),itend=e._SYMBptr->feuille._VECTptr->end(); 515 for (;it!=itend;++it){ 516 if (it->type==_INT_) 517 coeff = coeff * (*it); 518 else { 519 if ( (it->type==_SYMB) && (it->_SYMBptr->sommet==at_neg)){ 520 coeff = -coeff; 521 arg = arg * (it->_SYMBptr->feuille); 522 } 523 else 524 arg= arg * (*it); 525 } 526 } 527 if ( (coeff.type==_INT_) && (coeff.val<0) ){ 528 coeff=-coeff; 529 arg=-arg; 530 } 531 } 532 533 /* 534 gen sigma(const vecteur & v){ 535 // an "intelligent" version of return symbolic(at_plus,v); 536 // split each element of v as an integer coeff * an arg 537 vecteur vcoeff,varg,res; 538 int pos; 539 gen coeff,arg; 540 const_iterateur it=v.begin(),itend=v.end(); 541 vcoeff.reserve(itend-it); 542 varg.reserve(itend-it); 543 for (;it!=itend;++it){ 544 split(*it,coeff,arg); 545 pos=equalposcomp(varg,arg); 546 if (!pos){ 547 vcoeff.push_back(coeff); 548 varg.push_back(arg); 549 } 550 else 551 vcoeff[pos-1]=vcoeff[pos-1]+coeff; 552 } 553 it=vcoeff.begin(),itend=vcoeff.end(); 554 res.reserve(itend-it); 555 const_iterateur vargit=varg.begin(); 556 for (;it!=itend;++it,++vargit) 557 if (!is_zero(*it)) 558 res.push_back((*it)*(*vargit)); 559 if (res.empty()) 560 return zero; 561 if (res.size()>1) 562 return symbolic(at_plus,res); 563 return res.front(); 564 } 565 */ 566 prod_expand(const gen & a,const gen & b,GIAC_CONTEXT)567 gen prod_expand(const gen & a,const gen & b,GIAC_CONTEXT){ 568 bool a_is_plus= (a.type==_SYMB) && (a._SYMBptr->sommet==at_plus); 569 bool b_is_plus= (b.type==_SYMB) && (b._SYMBptr->sommet==at_plus); 570 if ( (!a_is_plus) && (!b_is_plus) ) 571 return a*b; 572 if (!a_is_plus) // distribute wrt b 573 return symbolic(at_plus,a*(*b._SYMBptr->feuille._VECTptr)); 574 if (!b_is_plus) 575 return symbolic(at_plus,b*(*a._SYMBptr->feuille._VECTptr)); 576 // distribute wrt a AND b 577 const_iterateur ita=a._SYMBptr->feuille._VECTptr->begin(),itaend=a._SYMBptr->feuille._VECTptr->end(); 578 const_iterateur itb=b._SYMBptr->feuille._VECTptr->begin(),itbend=b._SYMBptr->feuille._VECTptr->end(); 579 vecteur v; 580 v.reserve((itbend-itb)*(itaend-ita)); 581 for (;ita!=itaend;++ita){ 582 for (itb=b._SYMBptr->feuille._VECTptr->begin();itb!=itbend;++itb) 583 v.push_back( (*ita) * (*itb) ); 584 } 585 return symbolic(at_plus,gen(v,_SEQ__VECT)); 586 } 587 prod_expand(const const_iterateur it,const const_iterateur itend,GIAC_CONTEXT)588 static gen prod_expand(const const_iterateur it,const const_iterateur itend,GIAC_CONTEXT){ 589 int s=int(itend-it); 590 if (s==0) 591 return plus_one; 592 if (s==1) 593 return *it; 594 return _simplifier(prod_expand(prod_expand(it,it+s/2,contextptr),prod_expand(it+s/2,itend,contextptr),contextptr),contextptr); 595 } prod_expand(const gen & e_orig,GIAC_CONTEXT)596 static gen prod_expand(const gen & e_orig,GIAC_CONTEXT){ 597 gen e=aplatir_fois_plus(expand(e_orig,contextptr)); 598 if (e.type!=_VECT) 599 return e; 600 // look for sommet=at_plus inside e 601 return prod_expand(e._VECTptr->begin(),e._VECTptr->end(),contextptr); 602 } 603 prod_expand_nosimp(const const_iterateur it,const const_iterateur itend,GIAC_CONTEXT)604 static gen prod_expand_nosimp(const const_iterateur it,const const_iterateur itend,GIAC_CONTEXT){ 605 int s=int(itend-it); 606 if (s==0) 607 return plus_one; 608 if (s==1) 609 return *it; 610 return prod_expand(prod_expand_nosimp(it,it+s/2,contextptr),prod_expand_nosimp(it+s/2,itend,contextptr),contextptr); 611 } prod_expand_nosimp(const gen & e_orig,GIAC_CONTEXT)612 static gen prod_expand_nosimp(const gen & e_orig,GIAC_CONTEXT){ 613 gen e=aplatir_fois_plus(e_orig); 614 if (e.type!=_VECT) 615 return e; 616 // look for sommet=at_plus inside e 617 return prod_expand_nosimp(e._VECTptr->begin(),e._VECTptr->end(),contextptr); 618 } 619 pow_expand_add_res(vector<gen> & factn,int pos,int sumexpo,const gen & coeff,const vecteur & w,const gen & p,int k,int n,vecteur & res,GIAC_CONTEXT)620 static void pow_expand_add_res(vector<gen> & factn,int pos,int sumexpo,const gen & coeff,const vecteur & w,const gen & p, int k,int n,vecteur & res,GIAC_CONTEXT){ 621 if (sumexpo==k){ 622 // End recursion 623 res.push_back(coeff*p); 624 return; 625 } 626 if (pos==n-1){ 627 // End recursion 628 res.push_back(coeff/factn[k-sumexpo]*p*expand(pow(w[pos],k-sumexpo),contextptr)); 629 return; 630 } 631 for (int i=k-sumexpo;i>=0;--i){ 632 pow_expand_add_res(factn,pos+1,sumexpo+i,coeff/factn[i],w,expand(p*pow(w[pos],i),contextptr),k,n,res,contextptr); 633 } 634 } 635 expand_pow_expand(const gen & e,GIAC_CONTEXT)636 static gen expand_pow_expand(const gen & e,GIAC_CONTEXT){ 637 if (e.type!=_VECT || e._VECTptr->size()!=2) 638 return e; 639 vecteur & v=*e._VECTptr; 640 gen base=expand(v[0],contextptr); 641 gen exponent=expand(v[1],contextptr); 642 if (v[1].is_symb_of_sommet(at_plus) && v[1]._SYMBptr->feuille.type==_VECT){ 643 vecteur & w=*v[1]._SYMBptr->feuille._VECTptr; 644 const_iterateur it=w.begin(),itend=w.end(); 645 vecteur prodarg; 646 prodarg.reserve(itend-it); 647 for (;it!=itend;++it){ 648 prodarg.push_back(pow(base,*it,contextptr)); 649 } 650 return _prod(prodarg,contextptr); 651 } 652 if (v[1].type==_INT_ ){ 653 if (v[0].is_symb_of_sommet(at_prod)&& v[0]._SYMBptr->feuille.type==_VECT){ 654 vecteur w(*v[0]._SYMBptr->feuille._VECTptr); 655 iterateur it=w.begin(),itend=w.end(); 656 for (;it!=itend;++it){ 657 *it=pow(expand(*it,contextptr),v[1],contextptr); 658 } 659 return _prod(w,contextptr); 660 } 661 } 662 if (v[0].is_symb_of_sommet(at_plus) && v[0]._SYMBptr->feuille.type==_VECT && v[1].type==_INT_ && v[1].val>=0){ 663 int k=v[1].val; 664 if (!k) 665 return plus_one; 666 if (k==1) 667 return base; 668 vector<gen> factn(k+1); 669 factn[0]=1; 670 for (int i=1;i<=k;++i){ 671 factn[i]=i*factn[i-1]; 672 } 673 // (x1+...+xn)^k -> sum_{j1+...+jn=k} k!/(j1!j2!...jn!) x^j1 *... *x^jk 674 vecteur & w=*v[0]._SYMBptr->feuille._VECTptr; 675 int n=int(w.size()); 676 if (!n) 677 return gensizeerr(contextptr); 678 vecteur res; 679 gen p; 680 for (int i=k;i>=0;--i){ 681 p=expand(pow(w[0],i),contextptr); 682 pow_expand_add_res(factn,1,i,factn[k]/factn[i],w,p,k,n,res,contextptr); 683 } 684 return symbolic(at_plus,res); 685 } 686 if (v[0].is_symb_of_sommet(at_prod) && v[0]._SYMBptr->feuille.type==_VECT){ 687 const vecteur & vb=*v[0]._SYMBptr->feuille._VECTptr; 688 gen r1(1),r2(1); 689 for (int i=0;i<vb.size();++i){ 690 if (fastsign(vb[i],contextptr)==1) 691 r1=r1*pow(vb[i],exponent,contextptr); 692 else 693 r2=r2*vb[i]; 694 } 695 if (r1!=1) 696 return r1*pow(r2,exponent,contextptr); 697 } 698 return symb_pow(base,exponent); 699 } 700 expand_neg_expand(const gen & g_orig,GIAC_CONTEXT)701 static gen expand_neg_expand(const gen & g_orig,GIAC_CONTEXT){ 702 gen g=expand(g_orig,contextptr); 703 return -g; 704 } 705 tchebycheff(int n,bool first_kind)706 vecteur tchebycheff(int n,bool first_kind){ 707 vecteur v0,v1,vtmp; 708 if (first_kind) { 709 v0.push_back(1); 710 v1.push_back(1); 711 v1.push_back(0); 712 } 713 else 714 v1.push_back(1); 715 if (!n) 716 return v0; 717 if (n==1) 718 return v1; 719 for (--n;n;--n){ 720 multvecteur(2,v1,vtmp); 721 vtmp.push_back(0); 722 vtmp=vtmp-v0; 723 v0=v1; 724 v1=vtmp; 725 } 726 return v1; 727 } 728 exp_expand(const gen & e,GIAC_CONTEXT)729 gen exp_expand(const gen & e,GIAC_CONTEXT){ 730 if (e.type!=_SYMB) 731 return exp(e,contextptr); 732 if (e._SYMBptr->sommet==at_plus) 733 return symbolic(at_prod,apply(e._SYMBptr->feuille,exp_expand,contextptr)); 734 gen coeff,arg; 735 split(e,coeff,arg,contextptr); 736 return pow(exp(arg,contextptr),coeff,contextptr); 737 } 738 ln_expand0(const gen & e,GIAC_CONTEXT)739 static gen ln_expand0(const gen & e,GIAC_CONTEXT){ 740 if (e.type==_FRAC) 741 return ln(e._FRACptr->num,contextptr)-ln(e._FRACptr->den,contextptr); 742 if (e.type!=_SYMB) 743 return ln(e,contextptr); 744 if (e._SYMBptr->sommet==at_prod) 745 return _plus(apply(e._SYMBptr->feuille,ln_expand0,contextptr),contextptr);//symbolic(at_plus,apply(e._SYMBptr->feuille,ln_expand0,contextptr)); 746 if (e._SYMBptr->sommet==at_inv) 747 return -ln_expand0(e._SYMBptr->feuille,contextptr); 748 if (e._SYMBptr->sommet==at_pow){ 749 gen & tmp=e._SYMBptr->feuille; 750 if (tmp.type==_VECT && tmp._VECTptr->size()==2){ 751 gen base=tmp._VECTptr->front(),expo=tmp._VECTptr->back(); 752 if (!complex_mode(contextptr) && expo.type==_INT_ && expo.val%2==0) 753 base=abs(base,contextptr); 754 return expo*ln_expand0(base,contextptr); 755 } 756 } 757 return ln(e,contextptr); 758 } 759 ln_expand(const gen & e0,GIAC_CONTEXT)760 gen ln_expand(const gen & e0,GIAC_CONTEXT){ 761 gen e(factor(e0,false,contextptr)); 762 return ln_expand0(e,contextptr); 763 } 764 symhorner(const vecteur & v,const gen & e)765 gen symhorner(const vecteur & v, const gen & e){ 766 if (v.empty()) 767 return zero; 768 if (is_zero(e)) 769 return v.back(); 770 gen res=zero; 771 const_iterateur it=v.begin(),itend=v.end(); 772 int n=int(itend-it)-1; 773 for (;it!=itend;++it,--n) 774 res = res + (*it)*pow(e,n); 775 return res; 776 } 777 778 static gen cos_expand(const gen & e,GIAC_CONTEXT); sin_expand(const gen & e,GIAC_CONTEXT)779 static gen sin_expand(const gen & e,GIAC_CONTEXT){ 780 if (e.type!=_SYMB) 781 return sin(e,contextptr); 782 if (lidnt(e)==vecteur(1,cst_pi)){ 783 gen sine=sin(e,contextptr); 784 if (!contains(lidnt(sine),cst_pi)) 785 return sine; 786 } 787 if (e._SYMBptr->sommet==at_plus){ 788 vecteur v=*e._SYMBptr->feuille._VECTptr; 789 gen last=v.back(),first; 790 v.pop_back(); 791 if (v.size()==1) 792 first=v.front(); 793 else 794 first=symbolic(at_plus,v); 795 return sin_expand(first,contextptr)*cos_expand(last,contextptr)+cos_expand(first,contextptr)*sin_expand(last,contextptr); 796 } 797 if (e._SYMBptr->sommet==at_neg) 798 return -sin_expand(e._SYMBptr->feuille,contextptr); 799 gen coeff,arg; 800 split(e,coeff,arg,contextptr); 801 if (!is_one(coeff) && coeff.type==_INT_ && coeff.val<max_texpand_expansion_order) 802 return symhorner(tchebycheff(coeff.val,false),cos(arg,contextptr))*sin(arg,contextptr); 803 else 804 return sin(e,contextptr); 805 } 806 cos_expand(const gen & e,GIAC_CONTEXT)807 static gen cos_expand(const gen & e,GIAC_CONTEXT){ 808 if (e.type!=_SYMB) 809 return cos(e,contextptr); 810 if (lidnt(e)==vecteur(1,cst_pi)){ 811 gen cose=cos(e,contextptr); 812 if (!contains(lidnt(cose),cst_pi)) 813 return cose; 814 } 815 if (e._SYMBptr->sommet==at_plus){ 816 vecteur v=*e._SYMBptr->feuille._VECTptr; 817 gen last=v.back(),first; 818 v.pop_back(); 819 if (v.size()==1) 820 first=v.front(); 821 else 822 first=symbolic(at_plus,v); 823 return cos_expand(first,contextptr)*cos_expand(last,contextptr)-sin_expand(first,contextptr)*sin_expand(last,contextptr); 824 } 825 if (e._SYMBptr->sommet==at_neg) 826 return cos_expand(e._SYMBptr->feuille,contextptr); 827 gen coeff,arg; 828 split(e,coeff,arg,contextptr); 829 if (!is_one(coeff) && coeff.type==_INT_ && coeff.val<max_texpand_expansion_order) 830 return symhorner(tchebycheff(coeff.val,true),cos(arg,contextptr)); 831 else 832 return cos(e,contextptr); 833 } 834 sin2tancos(const gen & g,GIAC_CONTEXT)835 static gen sin2tancos(const gen & g,GIAC_CONTEXT){ 836 return symb_cos(g)*symb_tan(g); 837 } 838 even_pow_cos2tan(const gen & g,GIAC_CONTEXT)839 static gen even_pow_cos2tan(const gen & g,GIAC_CONTEXT){ 840 if ( (g.type!=_VECT) || (g._VECTptr->size()!=2) ) 841 return g; 842 gen a(g._VECTptr->front()),b(g._VECTptr->back()); 843 if ( (b.type!=_INT_) || (b.val%2) || (a.type!=_SYMB) || (a._SYMBptr->sommet!=at_cos) ) 844 return symbolic(at_pow,g); 845 int i=b.val/2; 846 return pow(plus_one+pow(symb_tan(a._SYMBptr->feuille),plus_two,contextptr),-i,contextptr); 847 } 848 tan_expand(const gen & e,GIAC_CONTEXT)849 static gen tan_expand(const gen & e,GIAC_CONTEXT){ 850 if (e.type!=_SYMB) 851 return tan(e,contextptr); 852 if (lidnt(e)==vecteur(1,cst_pi)){ 853 gen tane=tan(e,contextptr); 854 if (!contains(lidnt(tane),cst_pi)) 855 return tane; 856 } 857 if (e._SYMBptr->sommet==at_plus){ 858 vecteur v=*e._SYMBptr->feuille._VECTptr; 859 gen last=v.back(),first; 860 v.pop_back(); 861 if (v.size()==1) 862 first=v.front(); 863 else 864 first=symbolic(at_plus,v); 865 gen ta=tan_expand(first,contextptr); 866 gen tb=tan_expand(last,contextptr); 867 if (is_inf(ta)) 868 return -inv(tb,contextptr); 869 if (is_inf(tb)) 870 return -inv(ta,contextptr); 871 return rdiv(ta+tb,1-ta*tb,contextptr); 872 } 873 if (e._SYMBptr->sommet==at_neg) 874 return -tan_expand(e._SYMBptr->feuille,contextptr); 875 gen coeff,arg; 876 split(e,coeff,arg,contextptr); 877 if (!is_one(coeff) && coeff.type==_INT_ && coeff.val<max_texpand_expansion_order){ 878 gen g=rdiv(symhorner(tchebycheff(coeff.val,false),cos(arg,contextptr))*sin(arg,contextptr),symhorner(tchebycheff(coeff.val,true),cos(arg,contextptr)),contextptr); 879 vector<const unary_function_ptr *> v; 880 vector< gen_op_context > w; 881 v.push_back(at_sin); 882 w.push_back(&sin2tancos); 883 g=subst(g,v,w,false,contextptr); 884 v[0]=at_pow; 885 w[0]=(&even_pow_cos2tan); 886 g=subst(recursive_normal(g,false,contextptr),v,w,false,contextptr); 887 return recursive_normal(g,false,contextptr); 888 } 889 else 890 return tan(e,contextptr); 891 } 892 symb_texpand(const gen & a)893 symbolic symb_texpand(const gen & a){ 894 return symbolic(at_texpand,a); 895 } 896 897 // "unary" version _texpand(const gen & args,GIAC_CONTEXT)898 gen _texpand(const gen & args,GIAC_CONTEXT){ 899 if ( args.type==_STRNG && args.subtype==-1) return args; 900 gen var,res; 901 if (is_algebraic_program(args,var,res)) 902 return symbolic(at_program,makesequence(var,0,_texpand(res,contextptr))); 903 if (is_equal(args)) 904 return apply_to_equal(args,_texpand,contextptr); 905 vector<const unary_function_ptr *> v; 906 vector< gen_op_context > w; 907 v.push_back(at_exp); 908 w.push_back(&exp_expand); 909 v.push_back(at_ln); 910 w.push_back(&ln_expand); 911 v.push_back(at_prod); 912 w.push_back(&prod_expand); 913 v.push_back(at_sin); 914 w.push_back(&sin_expand); 915 v.push_back(at_cos); 916 w.push_back(&cos_expand); 917 v.push_back(at_tan); 918 w.push_back(&tan_expand); 919 return subst(args,v,w,false,contextptr); 920 } 921 static const char _texpand_s []="texpand"; 922 static define_unary_function_eval (__texpand,&_texpand,_texpand_s); 923 define_unary_function_ptr5( at_texpand ,alias_at_texpand,&__texpand,0,true); 924 925 static const char _developper_transcendant_s []="developper_transcendant"; 926 static define_unary_function_eval (__developper_transcendant,&_texpand,_developper_transcendant_s); 927 define_unary_function_ptr5( at_developper_transcendant ,alias_at_developper_transcendant,&__developper_transcendant,0,true); 928 expand(const gen & e,GIAC_CONTEXT)929 gen expand(const gen & e,GIAC_CONTEXT){ 930 if (is_equal(e)) 931 return apply_to_equal(e,expand,contextptr); 932 gen var,res; 933 if (e.type!=_VECT && is_algebraic_program(e,var,res)) 934 return symbolic(at_program,makesequence(var,0,expand(res,contextptr))); 935 if (e.type==_VECT && e.subtype==_SEQ__VECT && e._VECTptr->size()==2){ 936 gen last=e._VECTptr->back(); 937 if (last.type==_STRNG || last.type==_FUNC){ 938 vector<const unary_function_ptr *> v; 939 vector< gen_op_context > w; 940 if (contains(last,gen(at_prod)) || (last.type==_STRNG && !strcmp(last._STRNGptr->c_str(),"*"))){ // expand * with no further simplification 941 v.push_back(at_prod); 942 w.push_back(prod_expand_nosimp); 943 } 944 if (contains(last,gen(at_ln))){ 945 v.push_back(at_ln); 946 w.push_back(&ln_expand); 947 } 948 if (contains(last,gen(at_exp))){ 949 v.push_back(at_exp); 950 w.push_back(&exp_expand); 951 } 952 if (contains(last,gen(at_sin))){ 953 v.push_back(at_sin); 954 w.push_back(&sin_expand); 955 } 956 if (contains(last,gen(at_cos))){ 957 v.push_back(at_cos); 958 w.push_back(&cos_expand); 959 } 960 if (contains(last,gen(at_tan))){ 961 v.push_back(at_tan); 962 w.push_back(&tan_expand); 963 } 964 return subst(e._VECTptr->front(),v,w,false,contextptr); 965 } 966 } 967 vector<const unary_function_ptr *> v; 968 vector< gen_op_context > w; 969 v.push_back(at_prod); 970 v.push_back(at_pow); 971 v.push_back(at_neg); 972 w.push_back(&prod_expand); 973 w.push_back(&expand_pow_expand); 974 w.push_back(&expand_neg_expand); 975 return _simplifier(subst(e,v,w,false,contextptr),contextptr); 976 } 977 static const char _expand_s []="expand"; symb_expand(const gen & args)978 symbolic symb_expand(const gen & args){ 979 return symbolic(at_expand,args); 980 } 981 static define_unary_function_eval (__expand,&expand,_expand_s); 982 define_unary_function_ptr( at_expand ,alias_at_expand ,&__expand); 983 expexpand(const gen & e,GIAC_CONTEXT)984 gen expexpand(const gen & e,GIAC_CONTEXT){ 985 if (is_equal(e)) 986 return apply_to_equal(e,expexpand,contextptr); 987 gen var,res; 988 if (is_algebraic_program(e,var,res)) 989 return symbolic(at_program,makesequence(var,0,expexpand(res,contextptr))); 990 vector<const unary_function_ptr *> v(1,at_exp); 991 vector< gen_op_context > w(1,&exp_expand); 992 return subst(e,v,w,false,contextptr); 993 } 994 static const char _expexpand_s []="expexpand"; 995 static define_unary_function_eval (__expexpand,&expexpand,_expexpand_s); 996 define_unary_function_ptr5( at_expexpand ,alias_at_expexpand,&__expexpand,0,true); 997 lnexpand(const gen & e,GIAC_CONTEXT)998 gen lnexpand(const gen & e,GIAC_CONTEXT){ 999 if (is_equal(e)) 1000 return apply_to_equal(e,lnexpand,contextptr); 1001 gen var,res; 1002 if (is_algebraic_program(e,var,res)) 1003 return symbolic(at_program,makesequence(var,0,lnexpand(res,contextptr))); 1004 vector<const unary_function_ptr *> v(1,at_ln); 1005 vector< gen_op_context > w(1,&ln_expand); 1006 return subst(e,v,w,false,contextptr); 1007 } 1008 static const char _lnexpand_s []="lnexpand"; 1009 static define_unary_function_eval (__lnexpand,&lnexpand,_lnexpand_s); 1010 define_unary_function_ptr5( at_lnexpand ,alias_at_lnexpand,&__lnexpand,0,true); 1011 trigexpand(const gen & e,GIAC_CONTEXT)1012 gen trigexpand(const gen & e,GIAC_CONTEXT){ 1013 if (is_equal(e)) 1014 return apply_to_equal(e,trigexpand,contextptr); 1015 gen var,res; 1016 if (is_algebraic_program(e,var,res)) 1017 return symbolic(at_program,makesequence(var,0,trigexpand(res,contextptr))); 1018 vector<const unary_function_ptr *> v; 1019 vector< gen_op_context > w; 1020 v.push_back(at_sin); 1021 w.push_back(&sin_expand); 1022 v.push_back(at_cos); 1023 w.push_back(&cos_expand); 1024 v.push_back(at_tan); 1025 w.push_back(&tan_expand); 1026 v.push_back(at_prod); 1027 w.push_back(&prod_expand); 1028 return subst(e,v,w,false,contextptr); 1029 } 1030 static const char _trigexpand_s []="trigexpand"; 1031 static define_unary_function_eval (__trigexpand,&trigexpand,_trigexpand_s); 1032 define_unary_function_ptr5( at_trigexpand ,alias_at_trigexpand,&__trigexpand,0,true); 1033 1034 #ifndef NO_NAMESPACE_GIAC 1035 } // namespace giac 1036 #endif // ndef NO_NAMESPACE_GIAC 1037