1 // -*- mode:C++ ; compile-command: "g++-3.4 -I.. -DHAVE_CONFIG_H -DIN_GIAC -g -c subst.cc" -*- 2 #include "giacPCH.h" 3 4 /* 5 * Copyright (C) 2000,14 B. Parisse, Institut Fourier, 38402 St Martin d'Heres 6 * 7 * This program is free software; you can redistribute it and/or modify 8 * it under the terms of the GNU General Public License as published by 9 * the Free Software Foundation; either version 3 of the License, or 10 * (at your option) any later version. 11 * 12 * This program is distributed in the hope that it will be useful, 13 * but WITHOUT ANY WARRANTY; without even the implied warranty of 14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 15 * GNU General Public License for more details. 16 * 17 * You should have received a copy of the GNU General Public License 18 * along with this program. If not, see <http://www.gnu.org/licenses/>. 19 */ 20 using namespace std; 21 #include <stdexcept> 22 #include <cmath> 23 #include "subst.h" 24 #include "symbolic.h" 25 #include "sym2poly.h" 26 #include "unary.h" 27 #include "usual.h" 28 #include "derive.h" 29 #include "intg.h" 30 #include "prog.h" 31 #include "lin.h" 32 #include "solve.h" 33 #include "plot.h" 34 #include "modpoly.h" 35 #include "maple.h" 36 #include "ti89.h" 37 #include "alg_ext.h" 38 #include "giacintl.h" 39 40 #ifndef NO_NAMESPACE_GIAC 41 namespace giac { 42 #endif // ndef NO_NAMESPACE_GIAC 43 checkanglemode(GIAC_CONTEXT)44 gen checkanglemode(GIAC_CONTEXT){ 45 if (!angle_radian(contextptr)) 46 //grad 47 return gensizeerr(gettext("This function works only in radian mode")); 48 return 0; 49 } 50 51 // bool quote_subst=false; 52 /* 53 static vector <const unary_function_ptr *> merge(const vector <const unary_function_ptr *>& v,const vector <const unary_function_ptr *> & w){ 54 vector <const unary_function_ptr *> res(v); 55 vector <const unary_function_ptr *>::const_iterator it=w.begin(),itend=w.end(); 56 for (;it!=itend;++it){ 57 res.push_back(*it); 58 } 59 return res; 60 } 61 62 static vector <gen_op> merge(const vector <gen_op>& v,const vector <gen_op> & w){ 63 vector <gen_op> res(v); 64 vector <gen_op>::const_iterator it=w.begin(),itend=w.end(); 65 for (;it!=itend;++it){ 66 res.push_back(*it); 67 } 68 return res; 69 } 70 */ 71 72 // sin, cos, tan in terms of tan(x/2) sin2tan2(const gen & e,GIAC_CONTEXT)73 gen sin2tan2(const gen & e,GIAC_CONTEXT){ 74 gen a=symb_tan(rdiv(e,plus_two,contextptr)); 75 return rdiv(plus_two*a,pow(a,2)+1,contextptr); 76 } 77 cos2tan2(const gen & e,GIAC_CONTEXT)78 gen cos2tan2(const gen & e,GIAC_CONTEXT){ 79 gen a=symb_tan(rdiv(e,plus_two,contextptr)); 80 return rdiv(1-pow(a,2),pow(a,2)+1,contextptr); 81 } 82 tan2tan2(const gen & e,GIAC_CONTEXT)83 gen tan2tan2(const gen & e,GIAC_CONTEXT){ 84 gen a=symb_tan(rdiv(e,plus_two,contextptr)); 85 return rdiv(plus_two*a,1-pow(a,2),contextptr); 86 } 87 88 // hyperbolic trig to exp sinh2exp(const gen & e,GIAC_CONTEXT)89 gen sinh2exp(const gen & e,GIAC_CONTEXT){ 90 gen a=exp(e,contextptr); 91 return rdiv(a-inv(a,contextptr),plus_two,contextptr); 92 } 93 cosh2exp(const gen & e,GIAC_CONTEXT)94 gen cosh2exp(const gen & e,GIAC_CONTEXT){ 95 gen a=exp(e,contextptr); 96 return rdiv(a+inv(a,contextptr),plus_two,contextptr); 97 } 98 tanh2exp(const gen & e,GIAC_CONTEXT)99 gen tanh2exp(const gen & e,GIAC_CONTEXT){ 100 gen a=exp(2*e,contextptr);//pow(exp(e,contextptr),2); 101 return rdiv(a-plus_one,a+plus_one,contextptr); 102 } 103 104 // trig to exp degtorad(const gen & g,GIAC_CONTEXT)105 gen degtorad(const gen & g,GIAC_CONTEXT){ 106 if (angle_radian(contextptr)) 107 return g; 108 return g*deg2rad_e; 109 } radtodeg(const gen & g,GIAC_CONTEXT)110 gen radtodeg(const gen & g,GIAC_CONTEXT){ 111 if (angle_radian(contextptr)) 112 return g; 113 return g*gen(180)/cst_pi; 114 } 115 //grad (next few commands angletorad(const gen & g,GIAC_CONTEXT)116 gen angletorad(const gen & g,GIAC_CONTEXT){ 117 if (angle_radian(contextptr)) 118 return g; 119 else if(angle_degree(contextptr)) 120 return g*deg2rad_e; 121 else 122 return g*grad2rad_e; 123 } radtoangle(const gen & g,GIAC_CONTEXT)124 gen radtoangle(const gen & g,GIAC_CONTEXT){ 125 if(angle_radian(contextptr)) 126 return g; 127 else if(angle_degree(contextptr)) 128 return g*gen(180)/cst_pi; 129 else 130 return g*gen(200)/cst_pi; 131 } sin2exp(const gen & e,GIAC_CONTEXT)132 gen sin2exp(const gen & e,GIAC_CONTEXT){ 133 gen a=exp(cst_i*angletorad(e,contextptr),contextptr); 134 return rdiv(a-inv(a,contextptr),plus_two*cst_i,contextptr); 135 } cos2exp(const gen & e,GIAC_CONTEXT)136 gen cos2exp(const gen & e,GIAC_CONTEXT){ 137 gen a=exp(cst_i*angletorad(e,contextptr),contextptr); 138 return rdiv(a+inv(a,contextptr),plus_two,contextptr); 139 } tan2exp(const gen & e,GIAC_CONTEXT)140 gen tan2exp(const gen & e,GIAC_CONTEXT){ 141 gen a=pow(exp(cst_i*angletorad(e,contextptr),contextptr),2); 142 return rdiv(a-plus_one,cst_i*(a+plus_one),contextptr); 143 } 144 exp2sincos(const gen & e,GIAC_CONTEXT)145 gen exp2sincos(const gen & e,GIAC_CONTEXT){ 146 gen a=re(e,contextptr),b=im(e,contextptr); 147 return exp(a,contextptr)*(cos(radtoangle(b,contextptr),contextptr)+cst_i*sin(radtoangle(b,contextptr),contextptr)); 148 } 149 tantosincos(const gen & e,GIAC_CONTEXT)150 gen tantosincos(const gen & e,GIAC_CONTEXT){ 151 return rdiv(symb_sin(e),symb_cos(e),contextptr); 152 } 153 tantosincos2(const gen & e,GIAC_CONTEXT)154 gen tantosincos2(const gen & e,GIAC_CONTEXT){ 155 gen e2=ratnormal(2*e,contextptr); 156 return rdiv(symb_sin(e2),(1+symb_cos(e2)),contextptr); 157 } 158 tantocossin2(const gen & e,GIAC_CONTEXT)159 gen tantocossin2(const gen & e,GIAC_CONTEXT){ 160 gen e2=ratnormal(2*e,contextptr); 161 return rdiv((1-symb_cos(e2)),symb_sin(e2),contextptr); 162 } 163 asintoacos(const gen & e,GIAC_CONTEXT)164 gen asintoacos(const gen & e,GIAC_CONTEXT){ 165 if (angle_radian(contextptr)) 166 return cst_pi_over_2-acos(e,contextptr); 167 else if(angle_degree(contextptr)) 168 return 90-acos(e,contextptr); 169 //grad 170 else 171 return 100 - acos(e, contextptr); 172 } 173 acostoasin(const gen & e,GIAC_CONTEXT)174 gen acostoasin(const gen & e,GIAC_CONTEXT){ 175 if (angle_radian(contextptr)) 176 return cst_pi_over_2-asin(e,contextptr); 177 else if(angle_degree(contextptr)) 178 return 90-asin(e,0); 179 //grad 180 else 181 return 90-asin(e,0); 182 } 183 asintoatan(const gen & e,GIAC_CONTEXT)184 gen asintoatan(const gen & e,GIAC_CONTEXT){ 185 return symb_atan(rdiv(e,sqrt(1-pow(e,plus_two,contextptr),contextptr),contextptr)); 186 } 187 atantoasin(const gen & e,GIAC_CONTEXT)188 gen atantoasin(const gen & e,GIAC_CONTEXT){ 189 return symb_asin(rdiv(e,sqrt(1+pow(e,plus_two,contextptr),contextptr),contextptr)); 190 } 191 acostoatan(const gen & e,GIAC_CONTEXT)192 gen acostoatan(const gen & e,GIAC_CONTEXT){ 193 return cst_pi_over_2-asintoatan(e,contextptr); 194 } 195 atantoacos(const gen & e,GIAC_CONTEXT)196 gen atantoacos(const gen & e,GIAC_CONTEXT){ 197 return _asin2acos(atantoasin(e,contextptr),contextptr); 198 } 199 asin2ln(const gen & g_orig,GIAC_CONTEXT)200 gen asin2ln(const gen & g_orig,GIAC_CONTEXT){ 201 //grad 202 gen g=angletorad(g_orig,contextptr); 203 return cst_i*ln(g+sqrt(pow(g,2)-1,contextptr),contextptr)+cst_pi_over_2; 204 } 205 acos2ln(const gen & g_orig,GIAC_CONTEXT)206 gen acos2ln(const gen & g_orig,GIAC_CONTEXT){ 207 //grad 208 gen g=angletorad(g_orig,contextptr); 209 return -cst_i*ln(g+sqrt(pow(g,2)-1,contextptr),contextptr); 210 } 211 atan2ln(const gen & g_orig,GIAC_CONTEXT)212 gen atan2ln(const gen & g_orig,GIAC_CONTEXT){ 213 //grad 214 gen g=angletorad(g_orig,contextptr); 215 return rdiv(cst_i*ln(rdiv(cst_i+g,cst_i-g),contextptr),plus_two,contextptr); 216 } 217 218 // g=[base, exponant] trigcospow(const gen & g0,GIAC_CONTEXT)219 gen trigcospow(const gen & g0,GIAC_CONTEXT){ 220 gen g(g0); 221 if (g.type!=_VECT) 222 return gensizeerr(contextptr); 223 g.subtype=_SEQ__VECT; 224 vecteur & v=*g._VECTptr; 225 gen & base=v.front(); 226 gen & expo=v.back(); 227 if ( (base.type!=_SYMB) || (expo.type!=_INT_) ) 228 return symbolic(at_pow,g); 229 gen tmpcos=symb_cos(base._SYMBptr->feuille); 230 int ediv=expo.val/2,emod=expo.val%2; 231 if (base._SYMBptr->sommet==at_sin) 232 return pow(1-pow(tmpcos,2),ediv)*pow(base,emod); 233 if (base._SYMBptr->sommet==at_tan){ 234 gen tmp=pow(tmpcos,2); 235 tmp=rdiv(plus_one,tmp,contextptr)-plus_one; // 1/ cos^2 -1 236 return pow(tmp,ediv)*pow(rdiv(symb_sin(base._SYMBptr->feuille),tmpcos,contextptr),emod); 237 } 238 return symbolic(at_pow,g); 239 } 240 trigsinpow(const gen & g0,GIAC_CONTEXT)241 gen trigsinpow(const gen & g0,GIAC_CONTEXT){ 242 gen g(g0); 243 if (g.type!=_VECT) 244 return gensizeerr(contextptr); 245 g.subtype=_SEQ__VECT; 246 vecteur & v=*g._VECTptr; 247 gen & base=v.front(); 248 gen & expo=v.back(); 249 if ( (base.type!=_SYMB) || (expo.type!=_INT_) ) 250 return symbolic(at_pow,g); 251 gen tmpsin=symb_sin(base._SYMBptr->feuille); 252 int ediv=expo.val/2,emod=expo.val%2; 253 if (base._SYMBptr->sommet==at_cos) 254 return pow(1-pow(tmpsin,2),ediv)*pow(base,emod); 255 if (base._SYMBptr->sommet==at_tan){ 256 gen tmp=pow(tmpsin,2); 257 tmp=rdiv(tmp,plus_one-tmp,contextptr); // sin^2/ (1-sin^2) 258 return pow(tmp,ediv)*pow(rdiv(tmpsin,symb_cos(base._SYMBptr->feuille),contextptr),emod); 259 } 260 return symbolic(at_pow,g); 261 } 262 trigtanpow(const gen & g0,GIAC_CONTEXT)263 gen trigtanpow(const gen & g0,GIAC_CONTEXT){ 264 gen g(g0); 265 if (g.type!=_VECT) 266 return gensizeerr(contextptr); 267 g.subtype=_SEQ__VECT; 268 vecteur & v=*g._VECTptr; 269 gen & base=v.front(); 270 gen & expo=v.back(); 271 if ( (base.type!=_SYMB) || (expo.type!=_INT_) ) 272 return symbolic(at_pow,g); 273 gen tmptan=symb_tan(base._SYMBptr->feuille); 274 int ediv=expo.val/2,emod=expo.val%2; 275 if (base._SYMBptr->sommet==at_cos) 276 return pow(1+pow(tmptan,2),-ediv)*pow(base,emod); 277 if (base._SYMBptr->sommet==at_sin){ 278 gen tmp=pow(tmptan,2); 279 tmp=rdiv(tmp,plus_one+tmp,contextptr); // tan^2/ (1+tan^2) 280 return pow(tmp,ediv)*pow(tmptan*symb_cos(base._SYMBptr->feuille),emod); 281 } 282 return symbolic(at_pow,g); 283 } 284 pownegtoinvpow(const gen & g0,GIAC_CONTEXT)285 gen pownegtoinvpow(const gen & g0,GIAC_CONTEXT){ 286 gen g(g0); 287 if (g.type!=_VECT) 288 return gensizeerr(contextptr); 289 g.subtype=_SEQ__VECT; 290 vecteur & v(*g._VECTptr); 291 if (v.size()!=2) 292 return gensizeerr(contextptr); 293 if (v[1].type!=_SYMB) 294 return symbolic(at_pow,g); 295 unary_function_ptr &u=v[1]._SYMBptr->sommet; 296 if (u==at_neg) 297 return inv(powtopowexpand(makevecteur(v[0],v[1]._SYMBptr->feuille),contextptr),contextptr); 298 return symbolic(at_pow,g); 299 } 300 powtopowexpand(const gen & g0,GIAC_CONTEXT)301 gen powtopowexpand(const gen & g0,GIAC_CONTEXT){ 302 gen g(g0); 303 if (g.type!=_VECT) 304 return gensizeerr(contextptr); 305 g.subtype=_SEQ__VECT; 306 vecteur & v(*g._VECTptr); 307 if (v.size()!=2) 308 return gensizeerr(contextptr); 309 if (v[1].type!=_SYMB) 310 return symbolic(at_pow,g); 311 unary_function_ptr &u=v[1]._SYMBptr->sommet; 312 if (u==at_neg) 313 return inv(powtopowexpand(makevecteur(v[0],v[1]._SYMBptr->feuille),contextptr),contextptr); 314 if ( (v[1]._SYMBptr->feuille.type!=_VECT) || ((u!=at_plus) && (u!=at_prod)) ) 315 return symbolic(at_pow,g); 316 vecteur & w=*v[1]._SYMBptr->feuille._VECTptr; 317 const_iterateur it=w.begin(),itend=w.end(); 318 if (u==at_plus){ 319 gen res(plus_one); 320 for (;it!=itend;++it) 321 res=res*powtopowexpand(makevecteur(v[0],*it),contextptr); 322 return res; 323 } 324 if (u==at_prod){ 325 if (w.size()!=2) 326 return symbolic(at_pow,g); 327 if (w[0].type==_INT_) 328 return pow(powtopowexpand(makevecteur(v[0],w[1]),contextptr),w[0],contextptr); 329 if (w[1].type==_INT_) 330 return pow(powtopowexpand(makevecteur(v[0],w[0]),contextptr),w[1],contextptr); 331 } 332 return symbolic(at_pow,g); 333 } 334 exptopower(const gen & g,GIAC_CONTEXT)335 gen exptopower(const gen & g,GIAC_CONTEXT){ 336 if (is_zero(g)) return 1; 337 gen a,ar,ai,b; 338 if (has_i(g) && !complex_mode(contextptr) && contains(g,cst_pi) && is_linear_wrt(g,cst_pi,a,b,contextptr) && !is_zero(a)){ 339 // exp(pi*a+b) 340 reim(a,ar,ai,contextptr); 341 // checked added for ai otherwise expre:=2*pi*sin(a*pi)/(1-cos(2*a*pi));simplify(expre); is wrong 342 if (is_zero(ar) && is_assumed_integer(ai,contextptr)) 343 return exptopower(b,contextptr)*pow(-1,ai,contextptr); 344 } 345 vecteur l(lop(g,at_ln)); 346 if (l.size()!=1) 347 return symbolic(at_exp,g); 348 identificateur tmp(" x"); 349 gen gg=subst(g,l,vecteur(1,tmp),false,contextptr); 350 if (!is_linear_wrt(gg,tmp,a,b,contextptr) || has_i(a)) 351 return symbolic(at_exp,g); 352 return exp(b,contextptr)*pow(l[0]._SYMBptr->feuille,a,contextptr); 353 } 354 355 // One substitution of an identifier subst_vecteur(const vecteur & v,const gen & i,const gen & newi,vecteur & w,bool quotesubst,GIAC_CONTEXT)356 static void subst_vecteur(const vecteur & v,const gen & i,const gen & newi,vecteur & w,bool quotesubst,GIAC_CONTEXT){ 357 if (&v==&w){ 358 vecteur::iterator it=w.begin(),itend=w.end(); 359 for (;it!=itend;++it) 360 *it=subst(*it,i,newi,quotesubst,contextptr); 361 } 362 else { 363 w.reserve(v.size()); 364 vecteur::const_iterator it=v.begin(),itend=v.end(); 365 for (;it!=itend;++it) 366 w.push_back(subst(*it,i,newi,quotesubst,contextptr)); 367 } 368 } 369 subst(const vecteur & v,const gen & i,const gen & newi,bool quotesubst,GIAC_CONTEXT)370 vecteur subst(const vecteur & v,const gen & i,const gen & newi,bool quotesubst,GIAC_CONTEXT){ 371 vecteur w; 372 subst_vecteur(v,i,newi,w,quotesubst,contextptr); 373 return w; 374 } 375 376 static bool has_subst(const gen & e,const gen & i,const gen & newi,gen & newe,bool quotesubst,GIAC_CONTEXT); 377 has_subst_vect(const gen & e,const gen & i,const gen & newi,gen & newe,bool quotesubst,GIAC_CONTEXT)378 static bool has_subst_vect(const gen & e,const gen & i,const gen & newi,gen & newe,bool quotesubst,GIAC_CONTEXT){ 379 const vecteur & v =*e._VECTptr; 380 const_iterateur it=v.begin(),itend=v.end(); 381 for (;it!=itend;++it){ 382 if (has_subst(*it,i,newi,newe,quotesubst,contextptr)) 383 break; 384 } 385 if (it==itend) 386 return false; 387 gen res(new_ref_vecteur(*e._VECTptr),e.subtype); 388 vecteur & w =*res._VECTptr; 389 iterateur jt=w.begin()+(it-v.begin()); 390 *jt=newe; 391 newe=res; 392 for (++jt,++it;it!=itend;++it,++jt){ 393 if (!has_subst(*it,i,newi,*jt,quotesubst,contextptr)) 394 *jt=*it; 395 } 396 return true; 397 } 398 subst_integrate(const gen & e,const gen & i,const gen & newi,bool quotesubst,int intg,GIAC_CONTEXT)399 static gen subst_integrate(const gen & e,const gen & i,const gen & newi,bool quotesubst,int intg,GIAC_CONTEXT){ 400 vecteur v=*e._SYMBptr->feuille._VECTptr; 401 int s=int(v.size()); 402 if ( (s>1) && (v[1]==i) ){ 403 vecteur l(*_lname(newi,contextptr)._VECTptr); 404 if (!l.empty() && l.front()==cst_pi) 405 l.erase(l.begin()); 406 if (l.empty()) 407 return gensizeerr(contextptr); 408 v[1]=l.front(); 409 v=subst(v,i,newi,quotesubst,contextptr); 410 gen der=derive(newi,l.front(),contextptr); 411 if (intg==0 && der!=1) 412 return gensizeerr("Unable to handle sum change of variables"); 413 else 414 v[0]=v[0]*der; 415 if (is_undef(v[0])) 416 return v[0]; 417 if (s>2){ 418 identificateur t(" tsubst"); 419 vecteur w=solve(newi-t,*l.front()._IDNTptr,1,contextptr); 420 if (w.empty()) 421 return gensizeerr(gettext("Unable to solve")); 422 if (s>3){ 423 bool ordonne=is_greater(v[3],v[2],contextptr); 424 v[2]=limit(w.front(),t,v[2],ordonne?1:-1,contextptr); 425 v[3]=limit(w.front(),t,v[3],ordonne?-1:1,contextptr); 426 } 427 else 428 v[2]=limit(w.front(),t,v[2],0,contextptr); 429 } 430 return symbolic((intg==0?at_sum:at_integrate),gen(v,_SEQ__VECT)); 431 } 432 v=subst(v,i,newi,quotesubst,contextptr); 433 if (intg && s>=3 && v[2]==v[3]){ 434 // *logptr(contextptr) << "Warning, assuming that " << v[0] << " is regular at " << v[2] << '\n'; 435 return 0; 436 } 437 return symbolic((intg==0?at_sum:at_integrate),gen(v,_SEQ__VECT)); 438 } 439 subst_derive(const gen & e,const gen & i,const gen & newi,bool quotesubst,GIAC_CONTEXT)440 static gen subst_derive(const gen & e,const gen & i,const gen & newi,bool quotesubst,GIAC_CONTEXT){ 441 if (series_flags(contextptr) & (1<<7)) 442 return symbolic(at_of,makesequence(e,symb_equal(i,newi))); 443 vecteur v=*e._SYMBptr->feuille._VECTptr; 444 int s=int(v.size()); 445 if ( (s==2) && (v[1]==i) ){ 446 vecteur l(*_lname(newi,contextptr)._VECTptr); 447 if (l.empty()) 448 return gensizeerr(contextptr); 449 v[1]=l.front(); 450 v=subst(v,i,newi,quotesubst,contextptr); 451 return rdiv(symbolic(at_derive,gen(v,_SEQ__VECT)),derive(newi,l.front(),contextptr),contextptr); 452 } 453 // Warning: ? return e for desolve (is_linear_diffeq) 454 return symbolic(at_derive,subst(gen(v,_SEQ__VECT),i,newi,quotesubst,contextptr)); 455 // return e; 456 } 457 has_subst(const gen & e,const gen & i,const gen & newi,gen & newe,bool quotesubst,GIAC_CONTEXT)458 static bool has_subst(const gen & e,const gen & i,const gen & newi,gen & newe,bool quotesubst,GIAC_CONTEXT){ 459 switch (e.type){ 460 case _INT_: case _ZINT: case _CPLX: case _DOUBLE_: case _REAL: case _STRNG: case _MOD: case _SPOL1: case _USER: 461 return false; 462 case _IDNT: case _FUNC: 463 if (e==i){ 464 newe=newi; 465 return true; 466 } 467 else 468 return false; 469 case _SYMB: 470 if (e==i){ 471 newe=newi; 472 return true; 473 } 474 if (i.type==_FUNC && e._SYMBptr->sommet==i){ 475 if (!has_subst(e._SYMBptr->feuille,i,newi,newe,quotesubst,contextptr)) 476 newe=e._SYMBptr->feuille; 477 newe=newi(newe,contextptr); 478 return true; 479 } 480 if ( e._SYMBptr->sommet==at_pow && i.type==_SYMB && i._SYMBptr->sommet==at_exp && (*(e._SYMBptr->feuille._VECTptr))[1]*ln((*(e._SYMBptr->feuille._VECTptr))[0],contextptr) == i._SYMBptr->feuille ) { 481 CERR << e << "=" << i << '\n'; 482 newe=newi; 483 return true; 484 } 485 if ( (e._SYMBptr->sommet==at_integrate || !strcmp(e._SYMBptr->sommet.ptr()->s,"integration") || !strcmp(e._SYMBptr->sommet.ptr()->s,"int")) && (e._SYMBptr->feuille.type==_VECT) && (i.type==_IDNT) ){ 486 newe=subst_integrate(e,i,newi,quotesubst,1,contextptr); 487 return true; 488 } 489 if ( e._SYMBptr->sommet==at_sum && e._SYMBptr->feuille.type==_VECT && (i.type==_IDNT) ){ 490 newe=subst_integrate(e,i,newi,quotesubst,0,contextptr); 491 return true; 492 } 493 if ( (e._SYMBptr->sommet==at_derive) && (e._SYMBptr->feuille.type==_VECT) &&(i.type==_IDNT) ){ 494 newe=subst_derive(e,i,newi,quotesubst,contextptr); 495 return true; 496 } 497 if (has_subst(e._SYMBptr->feuille,i,newi,newe,quotesubst,contextptr)){ 498 if (quotesubst || e._SYMBptr->sommet.quoted()) 499 newe=symbolic(e._SYMBptr->sommet,newe); 500 else 501 newe=e._SYMBptr->sommet(newe,contextptr); 502 return true; 503 } 504 else 505 return false; 506 case _VECT: 507 return has_subst_vect(e,i,newi,newe,quotesubst,contextptr); 508 case _FRAC: 509 if (e==i){ 510 newe=newi; 511 return true; 512 } 513 newe=fraction(subst(e._FRACptr->num,i,newi,quotesubst,contextptr),subst(e._FRACptr->den,i,newi,quotesubst,contextptr)); 514 return true; 515 default: 516 #ifndef NO_STDEXCEPT 517 settypeerr(contextptr); 518 #endif 519 return false; 520 } 521 } 522 subst(const gen & e,const gen & i,const gen & newi,bool quotesubst,GIAC_CONTEXT)523 gen subst(const gen & e,const gen & i,const gen & newi,bool quotesubst,GIAC_CONTEXT){ 524 if (is_inequation(newi) || newi.is_symb_of_sommet(at_and) || newi.is_symb_of_sommet(at_ou)) 525 return gensizeerr(contextptr); 526 if (i.type==_VECT){ 527 if (newi.type!=_VECT || i._VECTptr->size()!=newi._VECTptr->size()){ 528 #ifndef NO_STDEXCEPT 529 setdimerr(contextptr); 530 #endif 531 return e; 532 } 533 return subst(e,*i._VECTptr,*newi._VECTptr,quotesubst,contextptr); 534 } 535 if (i.type!=_IDNT && i.type!=_SYMB && i.type!=_FUNC) 536 *logptr(contextptr) << gettext("Warning, replacing ") << i << gettext(" by ") << newi << gettext(", a substitution variable should perhaps be purged.") << '\n'; 537 gen res; 538 if (has_subst(e,i,newi,res,quotesubst,contextptr)) 539 return res; 540 else 541 return e; 542 } 543 quotesubst(const gen & e,const gen & i,const gen & newi,GIAC_CONTEXT)544 gen quotesubst(const gen & e,const gen & i,const gen & newi,GIAC_CONTEXT){ 545 return subst(e,i,newi,true,contextptr); 546 } 547 548 #ifdef GIAC_MULTISUBST 549 static int multisubst(const gen & e,vecteur & res,const gen & x,const vecteur & xval,GIAC_CONTEXT); 550 multisubst_frac(const gen & e,vecteur & res,const gen & x,const vecteur & xval,GIAC_CONTEXT)551 static int multisubst_frac(const gen & e,vecteur & res,const gen & x,const vecteur & xval,GIAC_CONTEXT){ 552 vecteur resn,resd; 553 const gen & N=e._FRACptr->num; 554 const gen & D=e._FRACptr->den; 555 int mn=multisubst_frac(N,resn,x,xval,contextptr); 556 int md=multisubst_frac(D,resd,x,xval,contextptr); 557 if (!mn && !md) 558 return 0; 559 int n=xval.size(); 560 res=xval; 561 if (mn){ 562 if (md){ 563 for (int i=0;i<n;++i) 564 res[i]=resn[i]/resd[i]; 565 } 566 else { 567 for (int i=0;i<n;++i) 568 res[i]=resn[i]/D; 569 } 570 } 571 else { 572 for (int i=0;i<n;++i) 573 res[i]=N/resd[i]; 574 } 575 return 1; 576 } 577 multisubst_vect(const gen & e,vecteur & res,const gen & x,const vecteur & xval,GIAC_CONTEXT)578 static int multisubst_vect(const gen & e,vecteur & res,const gen & x,const vecteur & xval,GIAC_CONTEXT){ 579 const vecteur & v = *e._VECTptr; 580 int s=int(v.size()); 581 std::vector<vecteur> vi(s) ; 582 std::vector<int> mi(s); 583 int m=0; 584 for (int i=0;i<s;++i){ 585 if (mi[i]=multisubst(v[i],vi[i],x,xval,contextptr)) 586 m=1; 587 } 588 if (!m) 589 return 0; 590 // build res 591 int n=xval.size(); 592 res=vecteur(n); 593 iterateur it=res.begin(),itend=res.end(); 594 for (;it!=itend;++it){ 595 #ifdef SMARTPTR64 596 *it=vecteur(s); 597 #else 598 it->type=_VECT; 599 it->__VECTptr=new_ref_vecteur(s); 600 #endif 601 } 602 for (int i=0;i<s;++i){ 603 // compute the i-th col of res 604 if (mi[i]){ 605 vecteur & vii=vi[i]; 606 for (int j=0;j<n;j++){ 607 vecteur & resv=*res[j]._VECTptr; 608 resv[i]=vii[j]; 609 } 610 } 611 else { 612 for (int j=0;j<n;j++){ 613 vecteur & resv=*res[j]._VECTptr; 614 resv[i]=v[i]; 615 } 616 } 617 } 618 return 1; 619 } 620 multisubst_symb(const gen & e,vecteur & res,const gen & x,const vecteur & xval,GIAC_CONTEXT)621 static int multisubst_symb(const gen & e,vecteur & res,const gen & x,const vecteur & xval,GIAC_CONTEXT){ 622 int m=multisubst(e._SYMBptr->feuille,res,x,xval,contextptr); 623 // FIXME integrate and derive 624 if (!m) 625 return 0; 626 iterateur it=res.begin(),itend=res.end(); 627 for (;it!=itend;++it){ 628 *it=e._SYMBptr->sommet.quoted()?symbolic(e._SYMBptr->sommet,*it):e._SYMBptr->sommet(*it,contextptr); 629 } 630 return 1; 631 } 632 633 // returns 1 if e evals to several values 634 // or 0 if they are all same = x multisubst(const gen & e,vecteur & res,const gen & x,const vecteur & xval,GIAC_CONTEXT)635 static int multisubst(const gen & e,vecteur & res,const gen & x,const vecteur & xval,GIAC_CONTEXT){ 636 switch (e.type){ 637 case _INT_: case _ZINT: case _CPLX: case _DOUBLE_: case _REAL: case _STRNG: case _MOD: case _SPOL1: case _USER: 638 return 0; 639 case _IDNT: case _FUNC: 640 if (e==x){ 641 res=xval; 642 return 1; 643 } 644 else 645 return 0; 646 case _SYMB: 647 if (e==x){ 648 res=xval; 649 return 1; 650 } 651 return multisubst_symb(e,res,x,xval,contextptr); 652 case _VECT: 653 return multisubst_vect(e,res,x,xval,contextptr); 654 case _FRAC: 655 if (e==x){ 656 res=xval; 657 return 1; 658 } 659 return multisubst_frac(e,res,x,xval,contextptr); 660 default: 661 return gentypeerr(contextptr); 662 } 663 return 0; 664 } 665 multisubst(const gen & g,const gen & x,const vecteur & xval,GIAC_CONTEXT)666 static vecteur multisubst(const gen & g,const gen & x,const vecteur & xval,GIAC_CONTEXT){ 667 int s=xval.size(); 668 vecteur res; 669 if (multisubst(g,res,x,xval,contextptr)) 670 return res; 671 else 672 return vecteur(s,g); 673 } 674 _multisubst(const gen & args,GIAC_CONTEXT)675 static gen _multisubst(const gen & args,GIAC_CONTEXT){ 676 if ( args.type==_STRNG && args.subtype==-1) return args; 677 if (args.type!=_VECT || args._VECTptr->size()!=3) 678 return gendimerr(contextptr); 679 const vecteur & v=*args._VECTptr; 680 if (v[2].type!=_VECT) 681 return gensizeerr(contextptr); 682 if (v[2]._VECTptr->empty()) 683 return v[2]; 684 return multisubst(v[0],v[1],*v[2]._VECTptr,contextptr); 685 } 686 static const char _multisubst_s []="multisubst"; 687 static define_unary_function_eval (__multisubst,&_multisubst,_multisubst_s); 688 define_unary_function_ptr5( at_multisubst ,alias_at_multisubst,&__multisubst,0,true); 689 #endif // GIAC_MULTISUBST 690 subst(const sparse_poly1 & v,const gen & i,const gen & newi,bool quotesubst,GIAC_CONTEXT)691 sparse_poly1 subst(const sparse_poly1 & v,const gen & i,const gen & newi,bool quotesubst,GIAC_CONTEXT){ 692 sparse_poly1 p; 693 sparse_poly1::const_iterator it=v.begin(),itend=v.end(); 694 p.reserve(itend-it); 695 gen e; 696 for (;it!=itend;++it){ 697 e=recursive_normal(subst(it->coeff,i,newi,quotesubst,contextptr),contextptr); 698 if (!is_zero(e)) 699 p.push_back(monome(e,it->exponent)); 700 } 701 return p; 702 } 703 sortsubst(const vecteur & v,const vecteur & i,const vecteur & newi,bool quotesubst,GIAC_CONTEXT)704 vecteur sortsubst(const vecteur & v,const vecteur & i,const vecteur & newi,bool quotesubst,GIAC_CONTEXT){ 705 if (i.empty()) 706 return v; 707 if (v==i) 708 return newi; 709 vecteur w; 710 w.reserve(v.size()); 711 vecteur::const_iterator it=v.begin(),itend=v.end(); 712 for (;it!=itend;++it) 713 w.push_back(sortsubst(*it,i,newi,quotesubst,contextptr)); 714 return w; 715 } 716 717 // used for sorting first_sort(const gen & a,const gen & b)718 static bool first_sort(const gen & a,const gen & b){ 719 return islesscomplexthanf(a[0],b[0]); // FIXME GIAC_CONTEXT 720 } 721 722 // g known to be in interval findpos(const_iterateur it,const_iterateur itend,const gen & g)723 static int findpos(const_iterateur it,const_iterateur itend,const gen & g){ 724 const gen & itg=*it; 725 if (g==itg) 726 return 1; 727 int n=int(itend-it); 728 if (n<2) 729 return 0; 730 n /= 2; 731 const_iterateur itmid=it+n; 732 const gen & itgmid = *itmid; 733 if (islesscomplexthanf(g,itgmid)) 734 return findpos(it,itmid,g); 735 else { 736 int f=findpos(itmid,itend,g); 737 return f?n+f:0; 738 } 739 } 740 findpos(const vecteur & v,const gen & g)741 int findpos(const vecteur & v,const gen & g){ 742 const_iterateur it=v.begin(),itend=v.end(); 743 if (it==itend) 744 return 0; 745 if (g==*it) 746 return 1; 747 if (g==*(itend-1)) 748 return int(itend-it); 749 if (itend-it<=2) 750 return 0; 751 if (islesscomplexthanf(g,*it) || islesscomplexthanf(*(itend-1),g)) 752 return 0; 753 return findpos(it,itend,g); 754 } 755 756 // Multiple substitutions sort2(vecteur & i,vecteur & newi,GIAC_CONTEXT)757 static void sort2(vecteur & i,vecteur & newi,GIAC_CONTEXT){ 758 for (unsigned k=0;k<i.size();++k){ 759 if (i[k].type!=_IDNT && i[k].type!=_SYMB && i[k].type!=_FUNC && !is_zero(i[k]-newi[k])) 760 *logptr(contextptr) << gettext("Warning, replacing ") << i[k] << gettext(" by ") << newi[k] << gettext(", a substitution variable should perhaps be purged.") << '\n'; 761 } 762 int is=int(i.size()); 763 if (is<2) 764 return; 765 if (is==2 && is==int(newi.size()) ){ 766 if (islesscomplexthanf(i[0],i[1])) 767 return; 768 swapgen(i[0],i[1]); 769 swapgen(newi[0],newi[1]); 770 return; 771 } 772 // set same size, required for mrv substition in series.cc 773 for (int j=int(newi.size());j<is;++j) 774 newi.push_back(i[j]); 775 matrice atrier=mtran(makevecteur(i,newi)); 776 gen_sort_f(atrier.begin(),atrier.end(),first_sort); 777 atrier=mtran(atrier); 778 i=*atrier[0]._VECTptr; 779 newi=*atrier[1]._VECTptr; 780 } 781 subst(const vecteur & v,const vecteur & i,const vecteur & newi,bool quotesubst,GIAC_CONTEXT)782 vecteur subst(const vecteur & v,const vecteur & i,const vecteur & newi,bool quotesubst,GIAC_CONTEXT){ 783 vecteur sorti(i),sortnewi(newi); 784 sort2(sorti,sortnewi,contextptr); 785 return sortsubst(v,sorti,sortnewi,quotesubst,contextptr); 786 } 787 sortsubst(const gen & e,const vecteur & i,const vecteur & newi,bool quotesubst,GIAC_CONTEXT)788 gen sortsubst(const gen & e,const vecteur & i,const vecteur & newi,bool quotesubst,GIAC_CONTEXT){ 789 if (i.empty()) 790 return e; 791 int pos; 792 switch (e.type){ 793 case _INT_: case _ZINT: case _DOUBLE_: case _CPLX: case _REAL: 794 return e; 795 case _IDNT: 796 pos=findpos(i,e); 797 if (pos) 798 return newi[pos-1]; 799 else 800 return e; 801 case _SYMB: 802 pos=findpos(i,e); 803 if (pos) 804 return newi[pos-1]; 805 if (!quotesubst && abs_calc_mode(contextptr)!=38 && e._SYMBptr->sommet==at_pow){ 806 pos=findpos(i,exp((*(e._SYMBptr->feuille._VECTptr))[1]*ln((*(e._SYMBptr->feuille._VECTptr))[0],contextptr),contextptr) ); 807 if (pos) 808 return newi[pos-1]; 809 } 810 if (e._SYMBptr->feuille.type==_VECT){ 811 gen ef(sortsubst(*e._SYMBptr->feuille._VECTptr,i,newi,quotesubst,contextptr)); 812 ef.subtype=e._SYMBptr->feuille.subtype; 813 if (quotesubst || e._SYMBptr->sommet.quoted() 814 // || e._SYMBptr->sommet==at_pow 815 || (e._SYMBptr->sommet==at_pow && ef.type==_VECT && ef._VECTptr->size()==2 && ef._VECTptr->front().type>_POLY && !is_zero(ef._VECTptr->front())) 816 ) 817 return symbolic(e._SYMBptr->sommet,ef); 818 else 819 return e._SYMBptr->sommet(ef,contextptr); 820 } 821 else { 822 gen ef=sortsubst(e._SYMBptr->feuille,i,newi,quotesubst,contextptr); 823 if (quotesubst || e._SYMBptr->sommet.quoted()) 824 return symbolic(e._SYMBptr->sommet,ef); 825 else { 826 // code below would allow assume(z,complex);subst(re(z),[re,im],[x->(x+conj(x))/2,x->(x-conj(x))/2/i]); 827 // if (int pos=equalposcomp(i,e._SYMBptr->sommet)) return newi[pos-1](ef,contextptr); 828 return e._SYMBptr->sommet(ef,contextptr); 829 } 830 } 831 case _VECT: 832 return gen(sortsubst(*e._VECTptr,i,newi,quotesubst,contextptr),e.subtype); 833 case _FRAC: 834 pos=findpos(i,e); 835 if (pos) 836 return newi[pos-1]; 837 return fraction(sortsubst(e._FRACptr->num,i,newi,quotesubst,contextptr),sortsubst(e._FRACptr->den,i,newi,quotesubst,contextptr)); 838 default: 839 pos=findpos(i,e); 840 if (pos) 841 return newi[pos-1]; 842 return e; 843 } 844 } 845 subst(const gen & e,const vecteur & i,const vecteur & newi,bool quotesubst,GIAC_CONTEXT)846 gen subst(const gen & e,const vecteur & i,const vecteur & newi,bool quotesubst,GIAC_CONTEXT){ 847 if (i.empty()) return e; 848 vecteur sorti(i),sortnewi(newi); 849 sort2(sorti,sortnewi,contextptr); 850 return sortsubst(e,sorti,sortnewi,quotesubst,contextptr); 851 } 852 subst(const gen & e,const vector<const unary_function_ptr * > & v,const vector<gen_op_context> & w,bool quotesubst,GIAC_CONTEXT)853 gen subst(const gen & e,const vector<const unary_function_ptr *> & v,const vector< gen_op_context > & w,bool quotesubst,GIAC_CONTEXT){ 854 if (v.empty()) 855 return e; 856 if (e.type==_VECT){ 857 const_iterateur it=e._VECTptr->begin(),itend=e._VECTptr->end(); 858 vecteur res; 859 res.reserve(itend-it); 860 for (;it!=itend;++it) 861 res.push_back(subst(*it,v,w,quotesubst,contextptr)); 862 return gen(res,e.subtype); 863 } 864 if (e.type!=_SYMB) 865 return e; 866 if (e._SYMBptr->sommet==at_entry || e._SYMBptr->sommet==at_ans) 867 return gensizeerr(contextptr); 868 gen arg=subst(e._SYMBptr->feuille,v,w,quotesubst,contextptr); 869 int n=equalposcomp(v,&e._SYMBptr->sommet); 870 if (!n){ 871 if (quotesubst){ 872 gen res=symbolic(e._SYMBptr->sommet,arg); 873 res.subtype=e.subtype; 874 return res; 875 } 876 return e._SYMBptr->sommet(arg,contextptr); 877 } 878 gen tmp=w[n-1](arg,contextptr); 879 return tmp; 880 } 881 882 // Quick check if e contains some ptr of v has_op_list(const gen & e,const unary_function_ptr * v)883 bool has_op_list(const gen & e,const unary_function_ptr * v){ 884 if (e.type==_SYMB){ 885 if (equalposcomp(v,&e._SYMBptr->sommet)) 886 return true; 887 return has_op_list(e._SYMBptr->feuille,v); 888 } 889 if (e.type==_VECT){ 890 const_iterateur it=e._VECTptr->begin(),itend=e._VECTptr->end(); 891 for (;it!=itend;++it){ 892 if (has_op_list(*it,v)) 893 return true; 894 } 895 } 896 return false; 897 } 898 899 // Quick check if e contains v has_op(const gen & e,const unary_function_ptr & u)900 bool has_op(const gen & e,const unary_function_ptr & u){ 901 if (e.type==_SYMB){ 902 if (u==e._SYMBptr->sommet) 903 return true; 904 return has_op(e._SYMBptr->feuille,u); 905 } 906 if (e.type==_VECT){ 907 const_iterateur it=e._VECTptr->begin(),itend=e._VECTptr->end(); 908 for (;it!=itend;++it){ 909 if (has_op(*it,u)) 910 return true; 911 } 912 } 913 return false; 914 } 915 subst(const gen & e,const unary_function_ptr * v,const gen_op_context * w,bool quotesubst,GIAC_CONTEXT,bool recursive_nonrat)916 gen subst(const gen & e,const unary_function_ptr * v,const gen_op_context * w,bool quotesubst,GIAC_CONTEXT,bool recursive_nonrat){ 917 if (*v==0 || !has_op_list(e,v)) 918 return e; 919 if (e.type==_VECT){ 920 const_iterateur it=e._VECTptr->begin(),itend=e._VECTptr->end(); 921 vecteur res; 922 res.reserve(itend-it); 923 for (;it!=itend;++it) 924 res.push_back(subst(*it,v,w,quotesubst,contextptr,recursive_nonrat)); 925 return gen(res,e.subtype); 926 } 927 if (e.type!=_SYMB) 928 return e; 929 // recursive call only for rational operators 930 gen arg=e._SYMBptr->feuille; int pos=-1; 931 if (recursive_nonrat || ((pos= equalposcomp(analytic_sommets,e._SYMBptr->sommet)) && (pos<=4 || (pos==5 && is_integer(e._SYMBptr->feuille[1]))))) 932 arg=subst(arg,v,w,quotesubst,contextptr,recursive_nonrat); 933 int n=equalposcomp(v,e._SYMBptr->sommet); 934 if (!n){ 935 if (quotesubst){ 936 gen res=symbolic(e._SYMBptr->sommet,arg); 937 res.subtype=e.subtype; 938 return res; 939 } 940 return e._SYMBptr->sommet(arg,contextptr); 941 } 942 gen tmp=w[n-1](arg,contextptr); 943 return tmp; 944 } 945 subst(const gen & e,const vector<const unary_function_ptr * > & v,const vector<gen_op> & w,bool quotesubst,GIAC_CONTEXT)946 gen subst(const gen & e,const vector<const unary_function_ptr *> & v,const vector< gen_op > & w,bool quotesubst,GIAC_CONTEXT){ 947 if (v.empty()) 948 return e; 949 if (e.type==_VECT){ 950 const_iterateur it=e._VECTptr->begin(),itend=e._VECTptr->end(); 951 vecteur res; 952 res.reserve(itend-it); 953 for (;it!=itend;++it) 954 res.push_back(subst(*it,v,w,quotesubst,contextptr)); 955 return gen(res,e.subtype); 956 } 957 if (e.type!=_SYMB) 958 return e; 959 gen arg=subst(e._SYMBptr->feuille,v,w,quotesubst,contextptr); 960 int n=equalposcomp(v,&e._SYMBptr->sommet); 961 if (!n){ 962 if (quotesubst){ 963 gen res=symbolic(e._SYMBptr->sommet,arg); 964 res.subtype=e.subtype; 965 return res; 966 } 967 return e._SYMBptr->sommet(arg,contextptr); 968 } 969 gen tmp=w[n-1](arg); 970 return tmp; 971 } 972 halftan(const gen & e,GIAC_CONTEXT)973 gen halftan(const gen & e,GIAC_CONTEXT){ 974 return subst(e,sincostan_tab,halftan_tab,false,contextptr); 975 } 976 _halftan(const gen & args,GIAC_CONTEXT)977 gen _halftan(const gen & args,GIAC_CONTEXT){ 978 if ( args.type==_STRNG && args.subtype==-1) return args; 979 gen var,res; 980 if (is_algebraic_program(args,var,res)) 981 return symbolic(at_program,makesequence(var,0,_halftan(res,contextptr))); 982 if (is_equal(args)) 983 return apply_to_equal(args,_halftan,contextptr); 984 return halftan(args,contextptr); 985 } 986 static const char _halftan_s []="halftan"; 987 static define_unary_function_eval (__halftan,&_halftan,_halftan_s); 988 define_unary_function_ptr5( at_halftan ,alias_at_halftan,&__halftan,0,true); 989 990 // sin(x)=-cos(x+pi/2) shift_sin(const gen & x,GIAC_CONTEXT)991 static gen shift_sin(const gen & x,GIAC_CONTEXT){ 992 return -symb_cos(ratnormal(x+cst_pi_over_2,contextptr)); 993 } 994 // cos(x)=sin(x+pi/2) shift_cos(const gen & x,GIAC_CONTEXT)995 static gen shift_cos(const gen & x,GIAC_CONTEXT){ 996 return symb_sin(ratnormal(x+cst_pi_over_2,contextptr)); 997 } 998 // tan(x+pi/2)=-1/tan(x) shift_tan(const gen & x,GIAC_CONTEXT)999 static gen shift_tan(const gen & x,GIAC_CONTEXT){ 1000 return minus_one/symb_tan(ratnormal(x+cst_pi_over_2,contextptr)); 1001 } 1002 const gen_op_context shift_phase_v_tab[]={shift_sin,shift_cos,shift_tan}; shift_phase(const gen & e,GIAC_CONTEXT)1003 gen shift_phase(const gen & e,GIAC_CONTEXT){ 1004 return subst(e,sincostan_tab,shift_phase_v_tab,false,contextptr); 1005 } 1006 _shift_phase(const gen & args,GIAC_CONTEXT)1007 gen _shift_phase(const gen & args,GIAC_CONTEXT){ 1008 if ( args.type==_STRNG && args.subtype==-1) return args; 1009 if (is_equal(args)) 1010 return apply_to_equal(args,_shift_phase,contextptr); 1011 return shift_phase(args,contextptr); 1012 } 1013 static const char _shift_phase_s []="shift_phase"; 1014 static define_unary_function_eval (__shift_phase,&_shift_phase,_shift_phase_s); 1015 define_unary_function_ptr5( at_shift_phase ,alias_at_shift_phase,&__shift_phase,0,true); 1016 inv_test_exp(const gen & e,GIAC_CONTEXT)1017 gen inv_test_exp(const gen & e,GIAC_CONTEXT){ 1018 if ( (e.type==_SYMB) && (e._SYMBptr->sommet==at_exp)) 1019 return symbolic(at_exp,-e._SYMBptr->feuille); 1020 return inv(e,contextptr); 1021 } 1022 rewrite_hyper(const gen & e,GIAC_CONTEXT)1023 gen rewrite_hyper(const gen & e,GIAC_CONTEXT){ 1024 /* 1025 vector<const unary_function_ptr *> vu; 1026 vu.push_back(at_sinh); 1027 vu.push_back(at_cosh); 1028 vu.push_back(at_tanh); 1029 vu.push_back(at_inv); 1030 vector <gen_op_context> vv(hyp2exp_tab,hyp2exp_tab+3); 1031 vv.push_back(inv_test_exp); 1032 return subst(e,vu,vv,false,contextptr); 1033 */ 1034 return subst(e,sinhcoshtanhinv_tab,hypinv2exp_tab,false,contextptr); 1035 } 1036 hyp2exp(const gen & e,GIAC_CONTEXT)1037 gen hyp2exp(const gen & e,GIAC_CONTEXT){ 1038 return subst(e,sinhcoshtanh_tab,hyp2exp_tab,false,contextptr); 1039 } _hyp2exp(const gen & args,GIAC_CONTEXT)1040 gen _hyp2exp(const gen & args,GIAC_CONTEXT){ 1041 if ( args.type==_STRNG && args.subtype==-1) return args; 1042 gen var,res; 1043 if (is_algebraic_program(args,var,res)) 1044 return symbolic(at_program,makesequence(var,0,_hyp2exp(res,contextptr))); 1045 if (is_equal(args)) 1046 return apply_to_equal(args,_hyp2exp,contextptr); 1047 return hyp2exp(args,contextptr); 1048 } 1049 static const char _hyp2exp_s []="hyp2exp"; 1050 static define_unary_function_eval (__hyp2exp,&_hyp2exp,_hyp2exp_s); 1051 define_unary_function_ptr5( at_hyp2exp ,alias_at_hyp2exp,&__hyp2exp,0,true); 1052 sincos(const gen & e,GIAC_CONTEXT)1053 gen sincos(const gen & e,GIAC_CONTEXT){ 1054 //grad 1055 if (angle_radian(contextptr)){ 1056 gen tmp=subst(e,tan_tab,tan2sincos_tab,true,contextptr); 1057 tmp=_pow2exp(tmp,contextptr); 1058 tmp=subst(tmp,exp_tab,exp2sincos_tab,false,contextptr); 1059 tmp=_exp2pow(tmp,contextptr); 1060 return tmp; 1061 } 1062 else 1063 return e; 1064 } _sincos(const gen & args,GIAC_CONTEXT)1065 gen _sincos(const gen & args,GIAC_CONTEXT){ 1066 if ( args.type==_STRNG && args.subtype==-1) return args; 1067 gen var,res; 1068 if (is_algebraic_program(args,var,res)) 1069 return symbolic(at_program,makesequence(var,0,_sincos(res,contextptr))); 1070 if (is_equal(args)) 1071 return apply_to_equal(args,_sincos,contextptr); 1072 return sincos(args,contextptr); 1073 } 1074 static const char _sincos_s []="sincos"; 1075 static define_unary_function_eval (__sincos,&_sincos,_sincos_s); 1076 define_unary_function_ptr5( at_sincos ,alias_at_sincos,&__sincos,0,true); 1077 trig2exp(const gen & e,GIAC_CONTEXT)1078 gen trig2exp(const gen & e,GIAC_CONTEXT){ 1079 //grad 1080 if (angle_radian(contextptr)) 1081 return subst(e,sincostan_tab,trig2exp_tab,false,contextptr); 1082 else 1083 return e; 1084 } _trig2exp(const gen & args,GIAC_CONTEXT)1085 gen _trig2exp(const gen & args,GIAC_CONTEXT){ 1086 if ( args.type==_STRNG && args.subtype==-1) return args; 1087 gen var,res; 1088 if (is_algebraic_program(args,var,res)) 1089 return symbolic(at_program,makesequence(var,0,_trig2exp(res,contextptr))); 1090 if (is_equal(args)) 1091 return apply_to_equal(args,_trig2exp,contextptr); 1092 return trig2exp(args,contextptr); 1093 } 1094 static const char _trig2exp_s []="trig2exp"; 1095 static define_unary_function_eval (__trig2exp,&_trig2exp,_trig2exp_s); 1096 define_unary_function_ptr5( at_trig2exp ,alias_at_trig2exp,&__trig2exp,0,true); 1097 halftan_hyp2exp(const gen & e,GIAC_CONTEXT)1098 gen halftan_hyp2exp(const gen & e,GIAC_CONTEXT){ 1099 return subst(e,sincostansinhcoshtanh_tab,halftan_hyp2exp_tab,false,contextptr); 1100 } 1101 _halftan_hyp2exp(const gen & args,GIAC_CONTEXT)1102 gen _halftan_hyp2exp(const gen & args,GIAC_CONTEXT){ 1103 if ( args.type==_STRNG && args.subtype==-1) return args; 1104 gen var,res; 1105 if (is_algebraic_program(args,var,res)) 1106 return symbolic(at_program,makesequence(var,0,_halftan_hyp2exp(res,contextptr))); 1107 if (is_equal(args)) 1108 return apply_to_equal(args,_halftan_hyp2exp,contextptr); 1109 return halftan_hyp2exp(args,contextptr); 1110 } 1111 static const char _halftan_hyp2exp_s []="halftan_hyp2exp"; 1112 static define_unary_function_eval (__halftan_hyp2exp,&_halftan_hyp2exp,_halftan_hyp2exp_s); 1113 define_unary_function_ptr5( at_halftan_hyp2exp ,alias_at_halftan_hyp2exp,&__halftan_hyp2exp,0,true); 1114 asin2acos(const gen & e,GIAC_CONTEXT)1115 gen asin2acos(const gen & e,GIAC_CONTEXT){ 1116 return subst(e,asin_tab,asin2acos_tab,false,contextptr); 1117 } _asin2acos(const gen & args,GIAC_CONTEXT)1118 gen _asin2acos(const gen & args,GIAC_CONTEXT){ 1119 if ( args.type==_STRNG && args.subtype==-1) return args; 1120 gen var,res; 1121 if (is_algebraic_program(args,var,res)) 1122 return symbolic(at_program,makesequence(var,0,_asin2acos(res,contextptr))); 1123 if (is_equal(args)) 1124 return apply_to_equal(args,_asin2acos,contextptr); 1125 return asin2acos(args,contextptr); 1126 } 1127 static const char _asin2acos_s []="asin2acos"; 1128 static define_unary_function_eval (__asin2acos,&_asin2acos,_asin2acos_s); 1129 define_unary_function_ptr5( at_asin2acos ,alias_at_asin2acos,&__asin2acos,0,true); 1130 asin2atan(const gen & e,GIAC_CONTEXT)1131 gen asin2atan(const gen & e,GIAC_CONTEXT){ 1132 return subst(e,asin_tab,asin2atan_tab,false,contextptr); 1133 } _asin2atan(const gen & args,GIAC_CONTEXT)1134 gen _asin2atan(const gen & args,GIAC_CONTEXT){ 1135 if ( args.type==_STRNG && args.subtype==-1) return args; 1136 gen var,res; 1137 if (is_algebraic_program(args,var,res)) 1138 return symbolic(at_program,makesequence(var,0,_asin2atan(res,contextptr))); 1139 if (is_equal(args)) 1140 return apply_to_equal(args,_asin2atan,contextptr); 1141 return asin2atan(args,contextptr); 1142 } 1143 static const char _asin2atan_s []="asin2atan"; 1144 static define_unary_function_eval (__asin2atan,&_asin2atan,_asin2atan_s); 1145 define_unary_function_ptr5( at_asin2atan ,alias_at_asin2atan,&__asin2atan,0,true); 1146 acos2asin(const gen & e,GIAC_CONTEXT)1147 gen acos2asin(const gen & e,GIAC_CONTEXT){ 1148 return subst(e,acos_tab,acos2asin_tab,false,contextptr); 1149 } _acos2asin(const gen & args,GIAC_CONTEXT)1150 gen _acos2asin(const gen & args,GIAC_CONTEXT){ 1151 if ( args.type==_STRNG && args.subtype==-1) return args; 1152 gen var,res; 1153 if (is_algebraic_program(args,var,res)) 1154 return symbolic(at_program,makesequence(var,0,_acos2asin(res,contextptr))); 1155 if (is_equal(args)) 1156 return apply_to_equal(args,acos2asin,contextptr); 1157 return acos2asin(args,contextptr); 1158 } 1159 static const char _acos2asin_s []="acos2asin"; 1160 static define_unary_function_eval (__acos2asin,&_acos2asin,_acos2asin_s); 1161 define_unary_function_ptr5( at_acos2asin ,alias_at_acos2asin,&__acos2asin,0,true); 1162 acos2atan(const gen & e,GIAC_CONTEXT)1163 gen acos2atan(const gen & e,GIAC_CONTEXT){ 1164 return subst(e,acos_tab,acos2atan_tab,false,contextptr); 1165 } _acos2atan(const gen & args,GIAC_CONTEXT)1166 gen _acos2atan(const gen & args,GIAC_CONTEXT){ 1167 if ( args.type==_STRNG && args.subtype==-1) return args; 1168 gen var,res; 1169 if (is_algebraic_program(args,var,res)) 1170 return symbolic(at_program,makesequence(var,0,_acos2atan(res,contextptr))); 1171 if (is_equal(args)) 1172 return apply_to_equal(args,_acos2atan,contextptr); 1173 return acos2atan(args,contextptr); 1174 } 1175 static const char _acos2atan_s []="acos2atan"; 1176 static define_unary_function_eval (__acos2atan,&_acos2atan,_acos2atan_s); 1177 define_unary_function_ptr5( at_acos2atan ,alias_at_acos2atan,&__acos2atan,0,true); 1178 atan2asin(const gen & e,GIAC_CONTEXT)1179 gen atan2asin(const gen & e,GIAC_CONTEXT){ 1180 return subst(e,atan_tab,atan2asin_tab,false,contextptr); 1181 } _atan2asin(const gen & args,GIAC_CONTEXT)1182 gen _atan2asin(const gen & args,GIAC_CONTEXT){ 1183 if ( args.type==_STRNG && args.subtype==-1) return args; 1184 gen var,res; 1185 if (is_algebraic_program(args,var,res)) 1186 return symbolic(at_program,makesequence(var,0,_atan2asin(res,contextptr))); 1187 if (is_equal(args)) 1188 return apply_to_equal(args,_atan2asin,contextptr); 1189 return atan2asin(args,contextptr); 1190 } 1191 static const char _atan2asin_s []="atan2asin"; 1192 static define_unary_function_eval (__atan2asin,&_atan2asin,_atan2asin_s); 1193 define_unary_function_ptr5( at_atan2asin ,alias_at_atan2asin,&__atan2asin,0,true); 1194 atan2acos(const gen & e,GIAC_CONTEXT)1195 gen atan2acos(const gen & e,GIAC_CONTEXT){ 1196 return subst(e,atan_tab,atan2acos_tab,false,contextptr); 1197 } _atan2acos(const gen & args,GIAC_CONTEXT)1198 gen _atan2acos(const gen & args,GIAC_CONTEXT){ 1199 if ( args.type==_STRNG && args.subtype==-1) return args; 1200 gen var,res; 1201 if (is_algebraic_program(args,var,res)) 1202 return symbolic(at_program,makesequence(var,0,_atan2acos(res,contextptr))); 1203 if (is_equal(args)) 1204 return apply_to_equal(args,_atan2acos,contextptr); 1205 return atan2acos(args,contextptr); 1206 } 1207 static const char _atan2acos_s []="atan2acos"; 1208 static define_unary_function_eval (__atan2acos,&_atan2acos,_atan2acos_s); 1209 define_unary_function_ptr5( at_atan2acos ,alias_at_atan2acos,&__atan2acos,0,true); 1210 atrig2ln(const gen & e,GIAC_CONTEXT)1211 gen atrig2ln(const gen & e,GIAC_CONTEXT){ 1212 if (angle_radian(contextptr)) 1213 return subst(e,asinacosatan_tab,atrig2ln_tab,false,contextptr); 1214 else 1215 return e; 1216 } _atrig2ln(const gen & args,GIAC_CONTEXT)1217 gen _atrig2ln(const gen & args,GIAC_CONTEXT){ 1218 if ( args.type==_STRNG && args.subtype==-1) return args; 1219 gen var,res; 1220 if (is_algebraic_program(args,var,res)) 1221 return symbolic(at_program,makesequence(var,0,_atrig2ln(res,contextptr))); 1222 if (is_equal(args)) 1223 return apply_to_equal(args,_atrig2ln,contextptr); 1224 return atrig2ln(args,contextptr); 1225 } 1226 static const char _atrig2ln_s []="atrig2ln"; 1227 static define_unary_function_eval (__atrig2ln,&_atrig2ln,_atrig2ln_s); 1228 define_unary_function_ptr5( at_atrig2ln ,alias_at_atrig2ln,&__atrig2ln,0,true); 1229 is_rational(const gen & g)1230 bool is_rational(const gen & g){ 1231 if (is_integer(g)) 1232 return true; 1233 if (g.type!=_FRAC) 1234 return false; 1235 return is_integer(g._FRACptr->num) && is_integer(g._FRACptr->den); 1236 } 1237 // if g is a symbolic depending linearly and rationnaly on a ln, 1238 // factors out this term before taking the exp 1239 // N.B. this should probably also extract constants otherwise tsimplify(exp(x+1)+exp(x-1)) is left as is rewrite_strong_exp(const gen & g_orig,GIAC_CONTEXT)1240 static gen rewrite_strong_exp(const gen & g_orig,GIAC_CONTEXT){ 1241 if (g_orig.type!=_SYMB) 1242 return exp(g_orig,contextptr); 1243 gen g(g_orig),res(plus_one); 1244 vecteur v(lop(lvar(g),at_ln)); 1245 v.insert(v.begin(),cst_pi); 1246 int s=int(v.size()); 1247 identificateur t(" t"); 1248 for (int i=0;i<s;++i){ 1249 gen gt=quotesubst(g,v[i],t,contextptr); 1250 gen dg=normal(subst(derive(gt,t,contextptr),t,zero,false,contextptr),contextptr); 1251 if (is_undef(dg)) 1252 return dg; 1253 gen gdg=g-dg*v[i]; 1254 if (!i) 1255 dg=dg/cst_i; 1256 if (is_zero(dg)) 1257 continue; 1258 if (!is_rational(dg) && !is_assumed_integer(dg,contextptr)) 1259 continue; 1260 g=gdg; 1261 if (!i) 1262 res=res*exp(cst_i*cst_pi*dg,contextptr); 1263 else { 1264 if (dg.is_symb_of_sommet(at_neg)) 1265 res=res/pow(v[i]._SYMBptr->feuille,dg._SYMBptr->feuille,contextptr); 1266 else 1267 res=res*pow(v[i]._SYMBptr->feuille,dg,contextptr); 1268 } 1269 } 1270 return res*exp(normal(g,contextptr),contextptr); 1271 } 1272 1273 // After extracting the cst coeff of g 1274 // if g is a linear combination of the components of wrt, 1275 // this will return the coeffs of the linear comb 1276 // otherwise it will add g at the end of wrt and return [0...0 1] 1277 // It assumes that g and the coeff of wrt are multivariate rat. fractions as_linear_combination(const gen & g,vecteur & wrt,GIAC_CONTEXT)1278 vecteur as_linear_combination(const gen & g,vecteur & wrt,GIAC_CONTEXT){ 1279 vecteur v(wrt); 1280 v.push_back(g); 1281 gen d; 1282 lcmdeno_converted(v,d,contextptr); 1283 gen cl=v.back(); 1284 int n=int(wrt.size()); 1285 vecteur res(n),mat; 1286 polynome p; 1287 if (cl.type<_POLY){ 1288 // code added 26 may 2015 for simplify(exp(1/2+t/2)-exp(t/2)*exp(1/2)); 1289 for (unsigned i=0;i<v.size();++i){ 1290 if (v[i].type==_POLY){ 1291 p=polynome(monomial<gen>(cl,v[i]._POLYptr->dim)); 1292 break; 1293 } 1294 } 1295 } 1296 if (p.coord.empty() && cl.type!=_POLY){ // search for a non poly in wrt 1297 gen gg; 1298 int i; 1299 for (i=0;i<n;++i){ 1300 if (v[i].type!=_POLY){ 1301 gg=rdiv(g,wrt[i],contextptr); 1302 if (gg.type==_INT_) 1303 break; 1304 if ( (gg.type==_FRAC) && (gg._FRACptr->num.type==_INT_) && (gg._FRACptr->den.type==_INT_) ) 1305 break; 1306 } 1307 } 1308 if (i==n){ 1309 wrt.push_back(g); 1310 res.push_back(plus_one); 1311 } 1312 else 1313 res[i]=gg; 1314 return res; 1315 } 1316 // we are now back to express cl as linear comb with integer coeff 1317 // but here in multivariate polynomials 1318 // we build now a linear system and solve it 1319 if (p.coord.empty()) 1320 p=*cl._POLYptr; 1321 vector< monomial<gen> >::const_iterator it=p.coord.begin(),itend=p.coord.end(); 1322 for (;it!=itend;++it){ 1323 const index_m & i=it->index; 1324 vecteur tmp(n+1); 1325 tmp.back()=it->value; 1326 for (int j=0,k;j<n;++j){ 1327 gen & gg=v[j]; 1328 if (gg.type==_POLY){ 1329 if ( (k=gg._POLYptr->position(i)) >=0 ) 1330 tmp[j]=gg._POLYptr->coord[k].value; 1331 } 1332 else { 1333 if (i.is_zero()) 1334 tmp[j]=gg; 1335 } 1336 } 1337 mat.push_back(tmp); 1338 } 1339 // add monomials that are not inside the last polynomial 1340 polynome fait(p); 1341 for (int j=0,k;j<n;++j){ 1342 gen & gg=v[j]; 1343 if (gg.type!=_POLY) 1344 continue; 1345 polynome & pj =*gg._POLYptr; 1346 vector< monomial<gen> >::const_iterator it=pj.coord.begin(),itend=pj.coord.end(); 1347 for (;it!=itend;++it){ 1348 const index_m & i=it->index; 1349 k=fait.position(i); 1350 if (k>=0 ) // this monomial is already in the system 1351 continue; 1352 fait=fait+monomial<gen>(plus_one,i); 1353 // make a new line in the system for this new monomial 1354 vecteur tmp(n+1); 1355 tmp.back()=zero; 1356 for (int J=0,K;J<n;++J){ 1357 gen & GG=v[J]; 1358 if (GG.type==_POLY){ 1359 if ( (K=GG._POLYptr->position(i)) >=0 ) 1360 tmp[J]=GG._POLYptr->coord[K].value; 1361 } 1362 else { 1363 if (i.is_zero()) 1364 tmp[J]=GG; 1365 } 1366 } 1367 mat.push_back(tmp); 1368 } 1369 } 1370 // take re and im since coeffs are reals 1371 // FIXME: should have a better algo for rational coeffs 1372 mat=mergevecteur(*re(mat,contextptr)._VECTptr,*im(mat,contextptr)._VECTptr); 1373 // find solutions in kernel of mat 1374 matrice noyau=mker(mat,contextptr); 1375 // try each row of noyau that has a non-zero last coeff 1376 // since the first row coeff is -1 all coeff must be rationals 1377 int ns=int(noyau.size()); 1378 for (int i=0;i<ns;++i){ 1379 vecteur & w=*noyau[i]._VECTptr; 1380 if (is_zero(w.back())) 1381 continue; 1382 gen gg; 1383 lcmdeno(w,gg,contextptr); // lcmdeno_converted? 1384 // check that all coeff are integers 1385 const_iterateur jt=w.begin(),jtend=w.end(); 1386 for (;jt!=jtend;++jt){ 1387 if (jt->type!=_INT_) 1388 break; 1389 } 1390 if (jt!=jtend) 1391 continue; 1392 // additional check: w.v=0 1393 gen check=dotvecteur(w,v); 1394 if (!is_zero(check)) 1395 continue; 1396 gg=w.back(); 1397 w.pop_back(); 1398 res=divvecteur(w,-gg); 1399 return res; 1400 } 1401 // no solution found 1402 wrt.push_back(g); 1403 res.push_back(plus_one); 1404 return res; 1405 } 1406 is_unit(const gen & g)1407 bool is_unit(const gen & g){ 1408 if (g.type==_INT_) 1409 return (g.val==1) || (g.val==-1); 1410 if (g.type==_ZINT) 1411 return (g==1) || (g==-1); 1412 if (g.type==_CPLX) 1413 return (g._CPLXptr->type==_INT_) && (g._CPLXptr->val==0) && ( (g._CPLXptr+1)->type==_INT_) && ( ( (g._CPLXptr+1)->val==1) || ( (g._CPLXptr+1)->val==-1) ); 1414 if (g.type==_POLY) 1415 return (g._POLYptr->coord.size()==1) && (g._POLYptr->coord.front().index.is_zero()) && is_unit(g._POLYptr->coord.front().value); 1416 return false; 1417 } 1418 is_algebraic_extension(const gen & g)1419 bool is_algebraic_extension(const gen & g){ 1420 if (g.type==_EXT) 1421 return true; 1422 if (g.type==_POLY && g._POLYptr->dim==0 && !g._POLYptr->coord.empty()) 1423 return is_algebraic_extension(g._POLYptr->coord.front().value); 1424 return false; 1425 } 1426 1427 // a is an extension, search if can be expressed as a product of powers of ext ext_relation(const gen & a0,const vecteur & ext,const vecteur & vars,vecteur & coeffs,GIAC_CONTEXT)1428 static bool ext_relation(const gen & a0,const vecteur & ext,const vecteur & vars,vecteur & coeffs,GIAC_CONTEXT){ 1429 if (ext.empty()) 1430 return false; 1431 gen a=r2e(a0,vars,contextptr); 1432 vecteur w=*r2e(ext,vars,contextptr)._VECTptr; 1433 w.push_back(a); 1434 vecteur l(lidnt(w)); 1435 if (!l.empty()){ // check for random values of the variables 1436 vecteur lval=vranm(int(l.size()),0,0); 1437 // should be improved to take care of assumptions and avoid bad eval points 1438 w=subst(w,l,lval,false,contextptr); 1439 } 1440 int s=int(w.size()); 1441 matrice test(s); 1442 gen precision=1<<30; // 2^30 1443 for (int i=0;i<s;i++){ 1444 vecteur tmp(s+2); 1445 tmp[i]=1; 1446 tmp[s]=_floor(precision*evalf(re(ln(w[i],contextptr),contextptr),1,contextptr),contextptr); 1447 tmp[s+1]=_floor(precision*evalf(im(ln(w[i],contextptr),contextptr),1,contextptr),contextptr); 1448 test[i]=tmp; 1449 } 1450 matrice S,A,L,O; 1451 S=lll(test,L,O,A,contextptr); 1452 if (is_undef(S)) 1453 return false; 1454 // if the first line of A has a small norm then it is the line of coeff 1455 coeffs=*A[0]._VECTptr; 1456 if (is_greater(linfnorm(coeffs,contextptr),20,contextptr)) 1457 return false; 1458 // found a probable relation, check it 1459 if (int(coeffs.size())!=s+2) 1460 return false; 1461 coeffs.pop_back(); 1462 coeffs.pop_back(); 1463 gen num(1),den(1); 1464 for (int i=0;i<s;++i){ 1465 if (is_positive(coeffs[i],contextptr)) 1466 num=num*pow((i==s-1?a0:ext[i]),coeffs[i],contextptr); 1467 else 1468 den=den*pow((i==s-1?a0:ext[i]),-coeffs[i],contextptr); 1469 } 1470 return is_zero(num-den); 1471 } 1472 1473 // Complete a list of prime together factors primeargs such that 1474 // a is a product of these factors decompose(const gen & a_orig,vecteur & primeargs,vecteur & extargs,int begin,const vecteur & vars,GIAC_CONTEXT)1475 static void decompose(const gen &a_orig,vecteur & primeargs,vecteur & extargs,int begin,const vecteur & vars,GIAC_CONTEXT){ 1476 if (is_unit(a_orig)) 1477 return; 1478 if (is_algebraic_extension(a_orig)){ 1479 vecteur coeffs; 1480 if (!ext_relation(a_orig,extargs,vars,coeffs,contextptr)) 1481 extargs.push_back(a_orig); 1482 return; 1483 } 1484 if (primeargs.empty()){ 1485 primeargs.push_back(a_orig); 1486 return; 1487 } 1488 gen a(a_orig); 1489 for (int i=begin;i<signed(primeargs.size());){ 1490 gen g=simplify3(a,primeargs[i]); 1491 if (g.type==_FRAC){ 1492 a=a*g; 1493 primeargs[i]=primeargs[i]*g; 1494 ++i; continue; 1495 } 1496 if (is_strictly_positive(r2e(-g,vars,contextptr),contextptr)){ 1497 g=-g; a=-a; primeargs[i]=-primeargs[i]; 1498 } 1499 if (is_unit(g)){ // prime together -> go to the next element of primeargs 1500 ++i; 1501 continue; 1502 } 1503 if (is_unit(primeargs[i])){ 1504 primeargs[i]=g; 1505 continue; 1506 } 1507 decompose(g,primeargs,extargs,i,vars,contextptr); 1508 } 1509 if (!is_unit(a)) 1510 primeargs.push_back(a); 1511 } 1512 branch_evalf(const gen & g,GIAC_CONTEXT)1513 static gen branch_evalf(const gen & g,GIAC_CONTEXT){ 1514 if (is_undef(g)) 1515 return g; 1516 vecteur v(*_lname(evalf(g,1,contextptr),contextptr)._VECTptr); 1517 gen gg(g); 1518 int s=int(v.size()); 1519 for (int i=0;i<s;++i){ 1520 vecteur w; 1521 gen point=0; 1522 int direction=1; 1523 find_range(v[i],w,contextptr); 1524 if (w.size()==1 && w.front().type==_VECT){ 1525 w=*w.front()._VECTptr; 1526 if (w.size()==2){ 1527 gen l=w.front(),m=w.back(); 1528 #if 1 1529 if (l==minus_inf && m==plus_inf) 1530 point=0; 1531 else 1532 point=ratnormal((l+m)/2,contextptr); 1533 if (!is_inf(point) && !is_undef(point)){ 1534 *logptr(contextptr) << gettext("Simplification assuming ") << v[i] << " near " << point << '\n'; 1535 point=subst(gg,*v[i]._IDNTptr,point,false,contextptr); 1536 if (!is_inf(point) && !is_undef(point)){ 1537 return evalf(point,1,contextptr); 1538 } 1539 } 1540 #endif 1541 if (!is_inf(m)){ 1542 point=m; 1543 direction=-1; 1544 } 1545 else { 1546 if (!is_inf(l)) 1547 point=l; 1548 } 1549 } 1550 } 1551 if (!is_inf(point) && !is_undef(point)) 1552 *logptr(contextptr) << gettext("Simplification assuming ") << v[i] << " near " << point << (direction==1?"+":"-") << '\n'; 1553 #ifdef NO_STDEXCEPT 1554 gg=limit(gg,*v[i]._IDNTptr,point,direction,contextptr); 1555 #ifdef TIMEOUT 1556 control_c(); 1557 #endif 1558 if (ctrl_c || interrupted) 1559 return gensizeerr(contextptr); 1560 if (is_undef(gg)) 1561 gg=0; 1562 #else 1563 try { 1564 gg=limit(gg,*v[i]._IDNTptr,point,direction,contextptr); 1565 } catch (std::runtime_error &){ 1566 last_evaled_argptr(contextptr)=NULL; 1567 gg=0; 1568 } 1569 #endif 1570 } 1571 return evalf(gg,1,contextptr); 1572 } expanded_ln(const gen & a_orig,const vecteur & primeargs,const vecteur & extargs,const vecteur & lnprimeargs,const vecteur & lnextargs,const vecteur & vars,GIAC_CONTEXT)1573 static gen expanded_ln(const gen & a_orig,const vecteur & primeargs,const vecteur & extargs,const vecteur & lnprimeargs,const vecteur & lnextargs,const vecteur & vars,GIAC_CONTEXT){ 1574 #ifdef TIMEOUT 1575 control_c(); 1576 #endif 1577 if (ctrl_c || interrupted) 1578 return gensizeerr(contextptr); 1579 gen res,gg,a=a_orig; 1580 int k; 1581 if (is_algebraic_extension(a)){ 1582 vecteur coeffs; 1583 if (ext_relation(a,extargs,vars,coeffs,contextptr)){ 1584 int s=int(coeffs.size())-1; 1585 for (int i=0;i<s;++i){ 1586 res = res + coeffs[i]*lnextargs[i]; 1587 } 1588 res=-res/coeffs[s]; 1589 } 1590 else 1591 res=ln(r2e(a,vars,contextptr),contextptr); 1592 return res; 1593 } 1594 int p=int(primeargs.size()); 1595 for (int j=0;j<p;++j){ 1596 for (k=0;;++k){ 1597 gen gp=primeargs[j]; 1598 gg=simplify3(a,gp); 1599 if (is_unit(gg)) 1600 break; 1601 } 1602 if (k) 1603 res=res+gen(k)*lnprimeargs[j]; 1604 res=res+ln(r2e(gg,vars,contextptr),contextptr); 1605 } 1606 res=res+ln(r2e(a,vars,contextptr),contextptr); 1607 gg=re(branch_evalf(rdiv(ln(r2e(a_orig,vars,contextptr),contextptr)-res,cst_pi_over_2*cst_i,contextptr),contextptr),contextptr); 1608 if (gg.type==_DOUBLE_) 1609 res=res+cst_pi_over_2*cst_i*gen(int(std::floor(gg._DOUBLE_val+0.5))); 1610 return res; 1611 } 1612 gen2poly(const gen & g,int s)1613 polynome gen2poly(const gen & g,int s){ 1614 if (g.type==_POLY) 1615 return *g._POLYptr; 1616 else 1617 return polynome(g,s); 1618 } 1619 // If g has exponential variables, factors out powers of 1620 // these exponential variables and return the ln of g simplified 1621 // plus the arg of the exponential vars simplifylnarg(const gen & g,GIAC_CONTEXT)1622 static gen simplifylnarg(const gen & g,GIAC_CONTEXT){ 1623 gen res; 1624 vecteur l(lvar(g)); 1625 int s=int(l.size()); 1626 fraction f(sym2r(g,l,contextptr)); 1627 gen nf,df; 1628 fxnd(f,nf,df); 1629 polynome n(gen2poly(nf,s)),d(gen2poly(df,s)); 1630 // for every exp variable inside l, look if n or d has a valuation 1631 // with respect to this degree 1632 const_iterateur it=l.begin(),itend=l.end(); 1633 for (int i=0;it!=itend;++it,++i){ 1634 if (!it->is_symb_of_sommet(at_exp)) 1635 continue; 1636 index_t shift_index(s); 1637 int nv=n.valuation(i); 1638 shift_index[i]=-nv; 1639 n=n.shift(shift_index); 1640 int dv=d.valuation(i); 1641 shift_index[i]=-dv; 1642 d=d.shift(shift_index); 1643 int v=nv-dv; 1644 res=res+v*it->_SYMBptr->feuille; 1645 } 1646 res=res+ln(r2sym(n,l,contextptr)/r2sym(d,l,contextptr),contextptr); 1647 return res; 1648 } simplifylnexp(const gen & g,GIAC_CONTEXT)1649 static gen simplifylnexp(const gen & g,GIAC_CONTEXT){ 1650 vecteur l(lop(g,at_ln)); 1651 int s=int(l.size()); 1652 if (!s) 1653 return g; 1654 vecteur lsub; 1655 for (int i=0;i<s;++i) 1656 lsub.push_back(simplifylnarg(l[i]._SYMBptr->feuille,contextptr)); 1657 return quotesubst(g,l,lsub,contextptr); 1658 } 1659 rlvar(const gen & e,vecteur & res,bool alg)1660 static void rlvar(const gen &e,vecteur & res,bool alg){ 1661 vecteur v; 1662 if (alg){ 1663 vecteur tmp=alg_lvar(e); 1664 // make one list from the matrix 1665 const_iterateur it=tmp.begin(),itend=tmp.end(); 1666 for (;it!=itend;++it){ 1667 if (it->type==_VECT){ 1668 const_iterateur jt=it->_VECTptr->begin(),jtend=it->_VECTptr->end(); 1669 for (;jt!=jtend;++jt){ 1670 if (!equalposcomp(v,*jt)) 1671 v.push_back(*jt); 1672 } 1673 } 1674 } 1675 } 1676 else 1677 v=lvar(e); 1678 vecteur::const_iterator it=v.begin(),itend=v.end(); 1679 for (;it!=itend;++it){ 1680 if (!equalposcomp(res,*it)){ 1681 // recursive call 1682 res.push_back(*it); 1683 if (it->type==_SYMB) { 1684 rlvar(it->_SYMBptr->feuille,res,alg); 1685 if (it->_SYMBptr->sommet==at_pow) 1686 rlvar(symbolic(at_ln,(*(it->_SYMBptr->feuille._VECTptr))[0]),res,alg); 1687 } 1688 } 1689 } 1690 } 1691 1692 // recursive (alg) list of var rlvar(const gen & e,bool alg)1693 vecteur rlvar(const gen &e,bool alg){ 1694 vecteur res; 1695 rlvar(e,res,alg); 1696 gen_sort_f(res.begin(),res.end(),symb_size_less); 1697 return res; 1698 } 1699 1700 _pow2exp(const gen & e,GIAC_CONTEXT)1701 gen _pow2exp(const gen & e,GIAC_CONTEXT){ 1702 if ( e.type==_STRNG && e.subtype==-1) return e; 1703 if (e.type==_VECT){ 1704 const_iterateur it=e._VECTptr->begin(),itend=e._VECTptr->end(); 1705 vecteur v; 1706 v.reserve(itend-it); 1707 for (;it!=itend;++it) 1708 v.push_back(_pow2exp(*it,contextptr)); 1709 return gen(v,e.subtype); 1710 } 1711 if (e.type!=_SYMB) 1712 return e; 1713 gen var,res; 1714 if (is_algebraic_program(e,var,res)) 1715 return symbolic(at_program,makesequence(var,0,_pow2exp(res,contextptr))); 1716 if ( (e._SYMBptr->sommet==at_pow || e._SYMBptr->sommet==at_surd || e._SYMBptr->sommet==at_NTHROOT) && e._SYMBptr->feuille.type==_VECT && e._SYMBptr->feuille._VECTptr->size()==2){ 1717 vecteur v=*e._SYMBptr->feuille._VECTptr; 1718 if (e._SYMBptr->sommet==at_NTHROOT) 1719 swapgen(v[0],v[1]); 1720 /* 1721 if (v[0].type==_VECT) 1722 return gensizeerr(gettext("Conversion of ^ of list/matrices to exp/ln not allowed. For symbolic power of square matrices, try matpow instead of ^")); 1723 */ 1724 if (v[1].type!=_INT_ && v[1].type!=_FRAC){ 1725 gen tmp=-v[0]; 1726 gen tmp1=(e._SYMBptr->sommet==at_surd || e._SYMBptr->sommet==at_NTHROOT)?inv(v[1],contextptr):v[1]; 1727 tmp1=_pow2exp(tmp1,contextptr); 1728 if (is_strictly_positive(tmp,contextptr)) 1729 return exp(tmp1*_pow2exp(ln(tmp,contextptr),contextptr),contextptr)*symb_exp(v[1]*cst_i*cst_pi); 1730 else 1731 return exp(tmp1*_pow2exp(ln(v[0],contextptr),contextptr),contextptr); 1732 } 1733 } 1734 return e._SYMBptr->sommet(_pow2exp(e._SYMBptr->feuille,contextptr),contextptr); 1735 } 1736 static const char _pow2exp_s []="pow2exp"; 1737 static define_unary_function_eval (__pow2exp,&_pow2exp,_pow2exp_s); 1738 define_unary_function_ptr5( at_pow2exp ,alias_at_pow2exp,&__pow2exp,0,true); 1739 pow2expln(const gen & e,const identificateur & x,GIAC_CONTEXT)1740 gen pow2expln(const gen & e,const identificateur & x,GIAC_CONTEXT){ 1741 if (e.type==_VECT){ 1742 const_iterateur it=e._VECTptr->begin(),itend=e._VECTptr->end(); 1743 vecteur v; 1744 v.reserve(itend-it); 1745 for (;it!=itend;++it) 1746 v.push_back(pow2expln(*it,x,contextptr)); 1747 return v; 1748 } 1749 if (e.type!=_SYMB) 1750 return e; 1751 if (e._SYMBptr->feuille.type==_VECT){ 1752 vecteur & v=*e._SYMBptr->feuille._VECTptr; 1753 if ((e._SYMBptr->sommet==at_pow) && ( contains(v[1],x) ||(v[1].type!=_INT_ && contains(v[0],x) ) ) ) 1754 return symb_exp(pow2expln(v[1],x,contextptr)*symb_ln(pow2expln(v[0],x,contextptr))); 1755 } 1756 return e._SYMBptr->sommet(pow2expln(e._SYMBptr->feuille,x,contextptr),contextptr); 1757 } 1758 pow2expln(const gen & e,GIAC_CONTEXT)1759 gen pow2expln(const gen & e,GIAC_CONTEXT){ 1760 if (e.type==_VECT) 1761 return apply(e,pow2expln,contextptr); 1762 if (e.type!=_SYMB) 1763 return e; 1764 if (e._SYMBptr->feuille.type==_VECT){ 1765 vecteur & v=*e._SYMBptr->feuille._VECTptr; 1766 if (e._SYMBptr->sommet==at_pow && v[1].type!=_INT_ && !(v[1].type==_FRAC && is_integer(v[0]))) 1767 return symb_exp(pow2expln(v[1],contextptr)*symb_ln(pow2expln(v[0],contextptr))); 1768 } 1769 return e._SYMBptr->sommet(pow2expln(e._SYMBptr->feuille,contextptr),contextptr); 1770 } 1771 simplifyfactorial(const gen & g,GIAC_CONTEXT)1772 gen simplifyfactorial(const gen & g,GIAC_CONTEXT){ 1773 vecteur l(lop(g,at_factorial)); 1774 int s=int(l.size()); 1775 if (s<2) 1776 return g; 1777 // look if difference of arguments in l are integers 1778 vecteur newl(s); 1779 // newl will contain couples of arg of factorials and integers shift 1780 newl[0]=makevecteur(l[0]._SYMBptr->feuille,zero); 1781 for (int i=1;i<s;++i){ 1782 // look in newl for an arg of factorial with integer difference 1783 gen current=l[i]._SYMBptr->feuille; 1784 newl[i]=makevecteur(current,zero); 1785 for (int j=0;j<i;++j){ 1786 if (!is_zero(newl[j][1])) 1787 continue; 1788 gen base=newl[j][0]; 1789 gen difference=current-base; 1790 difference=simplify(difference,contextptr); // recursive call 1791 if (difference.type!=_INT_) 1792 continue; 1793 int k=difference.val; 1794 if (k>=0){ 1795 newl[i]=makevecteur(base,difference); 1796 break; 1797 } 1798 // current is the new factorial basis, change base in newl 1799 newl[i]=makevecteur(current,zero); 1800 for (k=0;k<i;++k){ 1801 vecteur & n=*newl[k]._VECTptr; 1802 if (n[0]==base){ 1803 n[0]=current; 1804 n[1]=n[1]-difference; 1805 } 1806 } 1807 break; 1808 } 1809 } 1810 // now replace in newl the couple by a factorial*a product 1811 for (int i=0;i<s;++i){ 1812 gen base=newl[i][0]; 1813 int shift=newl[i][1].val; 1814 gen product=plus_one; 1815 for (int j=shift;j>0;--j){ 1816 product=product*(base+j); 1817 } 1818 newl[i]=product*symbolic(at_factorial,base); 1819 } 1820 return subst(g,l,newl,false,contextptr); 1821 } 1822 simplifypsi(const gen & g,GIAC_CONTEXT)1823 gen simplifypsi(const gen & g,GIAC_CONTEXT){ 1824 vecteur l(lop(g,at_Psi)); 1825 int s=int(l.size()); 1826 if (s<2) 1827 return g; 1828 // look if difference of arguments in l are integers 1829 vecteur newl(s); 1830 // newl will contain couples of arg of Psi and integers shift 1831 newl[0]=makevecteur(l[0]._SYMBptr->feuille,zero); 1832 for (int i=1;i<s;++i){ 1833 // look in newl for an arg of Psi with integer difference 1834 gen current=l[i]._SYMBptr->feuille; 1835 newl[i]=makevecteur(current,zero); 1836 for (int j=0;j<i;++j){ 1837 if (!is_zero(newl[j][1])) 1838 continue; 1839 gen base=newl[j][0]; 1840 gen difference=current-base; 1841 difference=simplify(difference,contextptr); // recursive call 1842 if (difference.type==_VECT && difference._VECTptr->size()==2 && difference._VECTptr->back()==0) 1843 difference=difference._VECTptr->front(); 1844 if (difference.type!=_INT_) 1845 continue; 1846 int k=difference.val; 1847 if (k>=0){ 1848 newl[i]=makevecteur(base,difference); 1849 break; 1850 } 1851 // current is the new Psi basis, change base in newl 1852 newl[i]=makevecteur(current,zero); 1853 for (k=0;k<i;++k){ 1854 vecteur & n=*newl[k]._VECTptr; 1855 if (n[0]==base){ 1856 n[0]=current; 1857 n[1]=n[1]-difference; 1858 } 1859 } 1860 break; 1861 } 1862 } 1863 // now replace in newl the couple by a factorial*a product 1864 for (int i=0;i<s;++i){ 1865 gen base=newl[i][0]; 1866 gen expo=-1; gen base0=base; 1867 if (base.type==_VECT && base._VECTptr->size()==2){ 1868 expo=-1-base._VECTptr->back(); 1869 base0=base._VECTptr->front(); 1870 } 1871 int shift=newl[i][1].val; 1872 gen somme=0; 1873 for (int j=shift-1;j>=0;--j){ 1874 somme += pow(base0+j,expo,contextptr); 1875 } 1876 newl[i]=somme+symbolic(at_Psi,base); 1877 } 1878 return subst(g,l,newl,false,contextptr); 1879 } 1880 gammatofactorial(const gen & x,GIAC_CONTEXT)1881 gen gammatofactorial(const gen & x,GIAC_CONTEXT){ 1882 if (x.is_symb_of_sommet(at_plus) && x._SYMBptr->feuille.type==_VECT){ 1883 vecteur v = *x._SYMBptr->feuille._VECTptr; 1884 if (v.size()==2 && v.back()==plus_one) 1885 return symbolic(at_factorial,v.front()); 1886 if (v.size()>2 && v.back()==plus_one){ 1887 v.pop_back(); 1888 return symbolic(at_factorial,symbolic(at_plus,gen(v,_SEQ__VECT))); 1889 } 1890 } 1891 if (x.type==_VECT) 1892 return symbolic(at_Gamma,x); 1893 return symbolic(at_factorial,x-1); 1894 } 1895 gamma2factorial(const gen & g,GIAC_CONTEXT)1896 gen gamma2factorial(const gen & g,GIAC_CONTEXT){ 1897 return subst(g,gamma_tab,gamma2factorial_tab,false,contextptr); 1898 } 1899 factorialtogamma(const gen & g,GIAC_CONTEXT)1900 gen factorialtogamma(const gen & g,GIAC_CONTEXT){ 1901 return symbolic(at_Gamma,g+1); 1902 } 1903 factorial2gamma(const gen & g,GIAC_CONTEXT)1904 gen factorial2gamma(const gen & g,GIAC_CONTEXT){ 1905 return subst(g,factorial_tab,factorial2gamma_tab,false,contextptr); 1906 } 1907 tsimplify_common(const gen & e,GIAC_CONTEXT)1908 gen tsimplify_common(const gen & e,GIAC_CONTEXT){ 1909 gen g=pow2expln(e,contextptr); 1910 g=gamma2factorial(g,contextptr); 1911 g=simplifyfactorial(g,contextptr); 1912 g=simplifypsi(g,contextptr); 1913 // analyse of args of ln 1914 g=simplifylnexp(g,contextptr); 1915 vecteur l(lop(g,at_ln)); 1916 int s=int(l.size()); 1917 if (s>1){ 1918 vecteur argln(s); 1919 for (int i=0;i<s;++i) 1920 argln[i]=tsimplify_common(l[i]._SYMBptr->feuille,contextptr); // was tsimplify, but replacing sin/cos in arg of ln is bad 1921 // check for common factors between args 1922 // now we make a vector of prime together values and rewrite 1923 // each arg as a product of these values + a multiple of 2*pi*i 1924 // arg <--> [ powers multiple_2*pi ] 1925 // Initialization 1926 vecteur vars(alg_lvar(argln)); 1927 // FIXME alg_lvar for alg ext but requires decompose to work! 1928 for (int i=0;i<s;++i) 1929 argln[i]=e2r(argln[i],vars,contextptr); 1930 // extensions must be treated separately 1931 gen num,den; 1932 vecteur primeargs,extargs; 1933 for (int i=0;i<s;++i){ 1934 // check for common factors between argln[i] and primeargs 1935 fxnd(argln[i],num,den); 1936 decompose(num,primeargs,extargs,0,vars,contextptr); 1937 decompose(den,primeargs,extargs,0,vars,contextptr); 1938 } 1939 // now rewrite ln[argln[i]] as a sum of ln[primeargs[]] 1940 // int p=primeargs.size(); 1941 vecteur lnprimeargs(*apply(r2e(primeargs,vars,contextptr),at_ln,contextptr)._VECTptr); 1942 vecteur lnextargs(*apply(r2e(extargs,vars,contextptr),at_ln,contextptr)._VECTptr); 1943 vecteur chk(lidnt(lop(lnextargs,at_rootof))); 1944 if (!chk.empty()) 1945 return e; 1946 vecteur newln(s); 1947 for (int i=0;i<s;++i){ 1948 #ifdef TIMEOUT 1949 control_c(); 1950 #endif 1951 if (ctrl_c || interrupted) 1952 return gensizeerr(contextptr); 1953 fxnd(argln[i],num,den); 1954 newln[i]=expanded_ln(num,primeargs,extargs,lnprimeargs,lnextargs,vars,contextptr)-expanded_ln(den,primeargs,extargs,lnprimeargs,lnextargs,vars,contextptr); 1955 #ifdef TIMEOUT 1956 control_c(); 1957 #endif 1958 if (ctrl_c || interrupted) 1959 return gensizeerr(contextptr); 1960 gen gg=evalf_double(re(branch_evalf(rdiv(l[i]-newln[i],cst_two_pi*cst_i,contextptr),contextptr),contextptr),0,contextptr); 1961 #ifdef TIMEOUT 1962 control_c(); 1963 #endif 1964 if (ctrl_c || interrupted) 1965 return gensizeerr(contextptr); 1966 if (gg.type==_DOUBLE_) 1967 newln[i]=newln[i]+cst_two_pi*cst_i*gen(int(std::floor(gg._DOUBLE_val+0.5))); 1968 } 1969 g=subst(g,l,newln,false,contextptr); 1970 } 1971 // analyse of arguments of exp: 1972 l=lop(g,at_exp); 1973 s=int(l.size()); 1974 if (!s) 1975 return recursive_normal(g,contextptr); 1976 // recursively simplify inside exp 1977 vecteur newl(s); // vector of args of the exponential 1978 for (int i=0;i<s;++i) 1979 newl[i]=tsimplify_common(l[i]._SYMBptr->feuille,contextptr); // was tsimplify, but replacing sin/cos in terms of complex exp is bad 1980 // check for linear relations with rational coefficients between args 1981 // add i*pi and ln to the linear relations checking 1982 // First convert everything to multivariate fractions 1983 vecteur vars(1,cst_pi); 1984 lvar(newl,vars); 1985 vecteur ln_vars(lop(newl,at_ln)); 1986 vecteur independant(*e2r(ln_vars,vars,contextptr)._VECTptr); 1987 int n_ln=int(independant.size()); 1988 independant.push_back(e2r(cst_i*cst_pi,vars,contextptr)); 1989 matrice m; 1990 int st=step_infolevel(contextptr); 1991 step_infolevel(0,contextptr); 1992 for (int i=0;i<s;++i){ 1993 m.push_back(as_linear_combination(e2r(newl[i],vars,contextptr),independant,contextptr)); 1994 } 1995 step_infolevel(st,contextptr); 1996 // we have a matrix of rational numbers and a vector "independant" 1997 // so that newl[i]=i-th line of the matrix dot independant 1998 // For each column of the matrix we take the lcm of the denominators 1999 // and multiply the column by the lcm, dividing the corresponding 2000 // coeff of independant 2001 // then we use exp[r2e(independant,vars)]^[i-th line] for exp[newl[i]] 2002 // we do the substitution l by exp[newl] in g 2003 // and we return normal(g) 2004 // First make m a rectangular array 2005 int c=int(m.back()._VECTptr->size()),r=int(m.size()); 2006 for (int i=0;i<r;++i){ 2007 int ms=int(m[i]._VECTptr->size()); 2008 if (ms<c) 2009 m[i]=mergevecteur(*m[i]._VECTptr,vecteur(c-ms,zero)); 2010 } 2011 // lcm 2012 matrice mt=mtran(m); 2013 gen lcmg; 2014 for (int i=0;i<c;++i){ 2015 lcmdeno_converted(*mt[i]._VECTptr,lcmg,contextptr); 2016 // check for arg of the form ln()/deno 2017 if (i<n_ln) 2018 independant[i]=pow(ln_vars[i]._SYMBptr->feuille,inv(lcmg,contextptr),contextptr); 2019 else 2020 independant[i]=rewrite_strong_exp(r2e(rdiv(independant[i],lcmg,contextptr),vars,contextptr),contextptr); 2021 } 2022 m=mtran(mt); 2023 // compute exp[newl] 2024 for (int i=0;i<s;++i){ 2025 vecteur & ligne=*m[i]._VECTptr; 2026 gen res(plus_one); 2027 for (int j=0;j<c;++j){ 2028 res=res*pow(independant[j],ligne[j],contextptr); 2029 } 2030 newl[i]=res; 2031 } 2032 g=subst(g,l,newl,false,contextptr); 2033 g=normal(g,contextptr); 2034 return g;// ratnormal(g,contextptr); 2035 } tsimplify(const gen & e,GIAC_CONTEXT)2036 gen tsimplify(const gen & e,GIAC_CONTEXT){ 2037 // replace trig/inv trig expressions with exp/ln 2038 gen g=trig2exp(e,contextptr); 2039 g=atrig2ln(g,contextptr); 2040 return tsimplify_common(g,contextptr); 2041 } tsimplify_noexpln(const gen & e,int s1,int s2,GIAC_CONTEXT)2042 gen tsimplify_noexpln(const gen & e,int s1,int s2,GIAC_CONTEXT){ 2043 int te=taille(e,65536); 2044 gen g=e; 2045 if (s1>1 && angle_radian(contextptr)){ 2046 #ifdef FXCG 2047 vecteur v1(loptab(e,sincostan_tab)); 2048 if (!v1[0].is_symb_of_sommet(at_tan) || !v1[1].is_symb_of_sommet(at_tan)) 2049 #endif 2050 g=subst(g,sincostan_tab,trig2exp_tab,false,contextptr,false); // g=trig2exp(e,contextptr); 2051 } 2052 if (s2>1 && angle_radian(contextptr)) 2053 g=subst(g,asinacosatan_tab,atrig2ln_tab,false,contextptr,false);//g=atrig2ln(g,contextptr); 2054 bool b=complex_mode(contextptr); 2055 complex_mode(true,contextptr); 2056 g=tsimplify_common(g,contextptr); 2057 complex_mode(b,contextptr); 2058 int tg=taille(g,8*te); // since sin/cos are coded as exp(i*...) 2059 if (tg>=8*te){ 2060 g=gamma2factorial(e,contextptr); 2061 g=simplifyfactorial(g,contextptr); 2062 return g; 2063 } 2064 return g; 2065 } _tsimplify(const gen & args,GIAC_CONTEXT)2066 gen _tsimplify(const gen & args,GIAC_CONTEXT){ 2067 if ( args.type==_STRNG && args.subtype==-1) return args; 2068 gen var,res; 2069 if (is_algebraic_program(args,var,res)) 2070 return symbolic(at_program,makesequence(var,0,_tsimplify(res,contextptr))); 2071 if (is_equal(args)) 2072 return apply_to_equal(args,_tsimplify,contextptr); 2073 return tsimplify(args,contextptr); 2074 } 2075 static const char _tsimplify_s []="tsimplify"; 2076 static define_unary_function_eval (__tsimplify,&_tsimplify,_tsimplify_s); 2077 define_unary_function_ptr5( at_tsimplify ,alias_at_tsimplify,&__tsimplify,0,true); 2078 2079 // find symbolic vars in g that have u has sommet lop(const gen & g,const unary_function_ptr & u)2080 vecteur lop(const gen & g,const unary_function_ptr & u){ 2081 if (!has_op(g,u)) 2082 return vecteur(0); 2083 if (g.type==_SYMB) { 2084 vecteur vrec=lop(g._SYMBptr->feuille,u); 2085 if (g._SYMBptr->sommet==u) 2086 vrec.push_back(g); 2087 return vrec; 2088 } 2089 if (g.type!=_VECT) 2090 return vecteur(0); 2091 vecteur res; 2092 const_iterateur it=g._VECTptr->begin(),itend=g._VECTptr->end(); 2093 for (;it!=itend;++it) 2094 res=mergeset(res,lop(*it,u)); 2095 return res; 2096 } 2097 2098 // find symbolic vars in g that have u has sommet lop(const gen & g,const unary_function_ptr * u)2099 vecteur lop(const gen & g,const unary_function_ptr * u){ 2100 if (!has_op(g,*u)) 2101 return vecteur(0); 2102 if (g.type==_SYMB) { 2103 vecteur vrec=lop(g._SYMBptr->feuille,u); 2104 if (g._SYMBptr->sommet==u) 2105 vrec.push_back(g); 2106 return vrec; 2107 } 2108 if (g.type!=_VECT) 2109 return vecteur(0); 2110 vecteur res; 2111 const_iterateur it=g._VECTptr->begin(),itend=g._VECTptr->end(); 2112 for (;it!=itend;++it) 2113 res=mergeset(res,lop(*it,u)); 2114 return res; 2115 } 2116 2117 // find symbolic vars in g that have u has sommet loptab(const gen & g,const unary_function_ptr * v)2118 vecteur loptab(const gen & g,const unary_function_ptr * v){ 2119 if (g.type==_SYMB) { 2120 if (equalposcomp(v,g._SYMBptr->sommet)) 2121 return vecteur(1,g); 2122 else 2123 return loptab(g._SYMBptr->feuille,v); 2124 } 2125 if (g.type!=_VECT) 2126 return vecteur(0); 2127 vecteur res; 2128 const_iterateur it=g._VECTptr->begin(),itend=g._VECTptr->end(); 2129 for (;it!=itend;++it) 2130 res=mergeset(res,loptab(*it,v)); 2131 return res; 2132 } 2133 2134 // find symbolic vars in g that have u has sommet lop(const gen & g,const vector<const unary_function_ptr * > & v)2135 vecteur lop(const gen & g,const vector<const unary_function_ptr *> & v){ 2136 if (g.type==_SYMB) { 2137 if (equalposcomp(v,&g._SYMBptr->sommet)) 2138 return vecteur(1,g); 2139 else 2140 return lop(g._SYMBptr->feuille,v); 2141 } 2142 if (g.type!=_VECT) 2143 return vecteur(0); 2144 vecteur res; 2145 const_iterateur it=g._VECTptr->begin(),itend=g._VECTptr->end(); 2146 for (;it!=itend;++it) 2147 res=mergeset(res,lop(*it,v)); 2148 return res; 2149 } 2150 expln2trig(const gen & g,GIAC_CONTEXT)2151 gen expln2trig(const gen & g,GIAC_CONTEXT){ 2152 if (g.type==_VECT) 2153 return apply(g,contextptr,expln2trig); 2154 if (g.type!=_SYMB) 2155 return g; 2156 if (g._SYMBptr->sommet==at_inv){ 2157 // inv[*]=*[inv] 2158 // inv[exp]=exp[neg] 2159 const gen & f=g._SYMBptr->feuille; 2160 if (f.type==_SYMB){ 2161 const unary_function_ptr & s=f._SYMBptr->sommet; 2162 const gen & ff=f._SYMBptr->feuille; 2163 if (s==at_exp) 2164 return expln2trig(symbolic(at_exp,-ff),contextptr); 2165 if (s==at_prod) 2166 return _prod(expln2trig(inv(ff,contextptr),contextptr),contextptr); 2167 if (s==at_pow) 2168 return pow(expln2trig(inv(ff._VECTptr->front(),contextptr),contextptr),ff._VECTptr->back(),contextptr); 2169 } 2170 // otherwise multiply by the conjugate 2171 gen tmp=expln2trig(g._SYMBptr->feuille,contextptr); 2172 gen retmp=re(tmp,contextptr); 2173 gen imtmp=im(tmp,contextptr); 2174 return (retmp-cst_i*imtmp)*inv(pow(retmp,2)+pow(imtmp,2),contextptr); 2175 } 2176 if (g._SYMBptr->sommet==at_exp) 2177 return sincos(g,contextptr); 2178 gen f=expln2trig(g._SYMBptr->feuille,contextptr); 2179 if (g._SYMBptr->sommet!=at_plus && g._SYMBptr->sommet!=at_prod && g._SYMBptr->sommet!=at_inv && g._SYMBptr->sommet!=at_pow && g._SYMBptr->sommet!=at_neg) 2180 f=normal(f,contextptr); 2181 if (g._SYMBptr->sommet==at_ln){ 2182 gen gre=simplify(re(f,contextptr),contextptr); 2183 gen gim=simplify(im(f,contextptr),contextptr); 2184 if (is_zero(gre)) 2185 return ln(pow(gim,2),contextptr)/2+sign(gim,contextptr)*cst_i*cst_pi_over_2; 2186 if (is_zero(gim)){ 2187 if (complex_mode(contextptr)) 2188 return rdiv(ln(pow(gre,2),contextptr),plus_two,contextptr)+cst_i*(plus_one-sign(gre,contextptr))*cst_pi_over_2; 2189 else 2190 return ln(gre,contextptr); 2191 } 2192 return rdiv(ln(pow(gre,2)+pow(gim,2),contextptr),plus_two,contextptr)+cst_i*(atan(gim/gre,contextptr)+sign(gim,contextptr)*(plus_one-sign(gre,contextptr))*cst_pi_over_2); 2193 } 2194 return symbolic(g._SYMBptr->sommet,f); 2195 } 2196 gen_feuille(const gen & g)2197 gen gen_feuille(const gen & g){ 2198 if (g.type==_SYMB) 2199 return g._SYMBptr->feuille; 2200 return g; 2201 } 2202 cossinexp2rootof(const gen & e,GIAC_CONTEXT,int maxalg)2203 gen cossinexp2rootof(const gen & e,GIAC_CONTEXT,int maxalg){ 2204 gen f=trig2exp(e,contextptr); 2205 vecteur l(lvar(f)),l1,l2; 2206 int n,d,q,r; 2207 for (unsigned i=0;i<l.size();++i){ 2208 if (l[i].is_symb_of_sommet(at_exp)){ 2209 gen li=l[i]._SYMBptr->feuille; 2210 if (contains(li,cst_pi)){ 2211 li=ratnormal(li/cst_pi/cst_i,contextptr); 2212 if (li.type==_FRAC && li._FRACptr->num.type==_INT_ && li._FRACptr->den.type==_INT_){ 2213 n=li._FRACptr->num.val; 2214 d=li._FRACptr->den.val; 2215 if (d<0){ n=-n; d=-d; } 2216 q=n/d; 2217 r=n%d; 2218 if (q%2) 2219 q=-1; 2220 else 2221 q=1; 2222 if (r<0) r += 2*d; 2223 // exp(r*i*pi/d) -> use rootof([1,..,0],cyclotomic(2*d)) 2224 vecteur vr(r+1); 2225 vr[0]=1; 2226 vecteur vc(cyclotomic(2*d)); 2227 if (vc.size()>(maxalg?maxalg:MAX_ALG_EXT_ORDER_SIZE)) 2228 return e; 2229 vr = vr % vc; 2230 if (!is_undef(vc)){ 2231 l1.push_back(l[i]); 2232 l2.push_back(q*symb_rootof(vr,vc,contextptr)); 2233 } 2234 } 2235 } 2236 } 2237 } 2238 if (l1.empty()) 2239 return e; 2240 f=subst(f,l1,l2,false,contextptr); 2241 return f; 2242 } 2243 powneg2invpow(const gen & e,GIAC_CONTEXT)2244 gen powneg2invpow(const gen & e,GIAC_CONTEXT){ 2245 return subst(e,pow_tab,powneg2invpow_tab,false,contextptr); 2246 } 2247 ataninvtoatan(const gen & g,GIAC_CONTEXT)2248 gen ataninvtoatan(const gen &g,GIAC_CONTEXT){ 2249 if (g.is_symb_of_sommet(at_inv)) 2250 return cst_pi_over_2*sign(g._SYMBptr->feuille,contextptr)-atan(g._SYMBptr->feuille,contextptr); 2251 return symb_atan(g); 2252 } 2253 ataninv2atan(const gen & g,GIAC_CONTEXT)2254 gen ataninv2atan(const gen & g,GIAC_CONTEXT){ 2255 const vector< const unary_function_ptr *> atan_v(1,at_atan); 2256 const vector< gen_op_context > ataninv2atan_v(1,ataninvtoatan); 2257 return subst(g,atan_v,ataninv2atan_v,false,contextptr); 2258 } in_cklin(const gen & tmp)2259 bool in_cklin(const gen & tmp){ 2260 if (tmp.is_symb_of_sommet(at_neg)) 2261 return in_cklin(tmp._SYMBptr->feuille); 2262 if (tmp.is_symb_of_sommet(at_exp) && has_i(tmp)) 2263 return true; 2264 if (tmp.is_symb_of_sommet(at_pow)) 2265 return in_cklin(tmp._SYMBptr->feuille[0]); 2266 if (tmp.is_symb_of_sommet(at_prod)){ 2267 gen t=tmp._SYMBptr->feuille; 2268 if (t.type==_VECT){ 2269 const_iterateur it=t._VECTptr->begin(),itend=t._VECTptr->end(); 2270 for (;it!=itend;++it) 2271 if (in_cklin(*it)) 2272 return true; 2273 } 2274 } 2275 return false; 2276 } 2277 2278 // ck if g has a denominator with exponentials, if so linearize cklin(const gen & g,GIAC_CONTEXT)2279 gen cklin(const gen & g,GIAC_CONTEXT){ 2280 vecteur gn,gd; 2281 prod2frac(g,gn,gd); 2282 if (gd.empty()) 2283 return g; 2284 for (unsigned i=0;i<gd.size();++i){ 2285 gen tmp=simplifier(gd[i],contextptr); 2286 if (in_cklin(tmp)) 2287 return _lin(g,contextptr); 2288 } 2289 return g; 2290 } 2291 simplify(const gen & e_orig,GIAC_CONTEXT)2292 gen simplify(const gen & e_orig,GIAC_CONTEXT){ 2293 if (e_orig.type<=_POLY || is_inf(e_orig) || has_num_coeff(e_orig)) 2294 return e_orig; 2295 gen e=simplifier(e_orig,contextptr); 2296 if (e.type==_FRAC) 2297 return _evalc(e_orig,contextptr); 2298 vecteur vsign=lop(e,at_sign); 2299 vecteur vabs=lop(e,at_abs),vs1,vs2; 2300 for (int i=0;i<int(vabs.size());++i){ 2301 vabs[i]=vabs[i]._SYMBptr->feuille; 2302 } 2303 for (int i=0;i<int(vsign.size());++i){ 2304 gen arg=vsign[i]._SYMBptr->feuille; 2305 if (equalposcomp(vabs,arg)){ 2306 vs1.push_back(symbolic(at_sign,arg)); 2307 vs2.push_back(symbolic(at_abs,arg)/arg); 2308 } 2309 } 2310 if (!vs1.empty()) 2311 e=subst(e,vs1,vs2,false,contextptr); 2312 // ratnormal added for E:=2*exp(t/25)/(19+exp(t/25)); F:=simplifier(int(E,t)); 2313 // M:=(1/50)*int(E,t,50,100); simplify(M) 2314 vecteur lnv=lop(e,at_ln); 2315 if (!lnv.empty()){ 2316 gen trye=lncollect(ratnormal(e_orig,contextptr),contextptr); 2317 if (lop(trye,at_ln).size()<lnv.size()) 2318 e=trye; 2319 } 2320 if (!lop(e,at_exp).empty()){ 2321 #ifdef NO_STDEXCEPT 2322 gen et=ratnormal(_texpand(e,contextptr),contextptr); 2323 if (!is_undef(et) && lvar(et).size()<lvar(e).size()) 2324 e=et; 2325 #else 2326 try { 2327 gen et=ratnormal(_texpand(e,contextptr),contextptr); 2328 if (lvar(et).size()<lvar(e).size()) 2329 e=et; 2330 } catch (std::runtime_error & err){ 2331 last_evaled_argptr(contextptr)=NULL; 2332 } 2333 #endif 2334 e=_exp2pow(e,contextptr); 2335 } 2336 if (!lop(e,at_pow).empty()) 2337 e=powneg2invpow(e,contextptr); 2338 if (!lop(e,at_atan).empty()) 2339 e=ataninv2atan(e,contextptr); 2340 if (contains(e,cst_pi)){ 2341 e=cossinexp2rootof(e,contextptr); 2342 e=normal(e,contextptr); 2343 } 2344 if (e.type==_SYMB && e._SYMBptr->feuille.type!=_VECT){ 2345 if (e._SYMBptr->sommet!=at_inv && e._SYMBptr->sommet!=at_neg) 2346 return e._SYMBptr->sommet(simplify(e._SYMBptr->feuille,contextptr),contextptr); 2347 } 2348 vabs=lop(e,at_abs); 2349 vecteur vabs2(vabs); 2350 iterateur it=vabs2.begin(),itend=vabs2.end(); 2351 for (;it!=itend;++it){ 2352 *it=symbolic(at_abs,factor(it->_SYMBptr->feuille,false,contextptr)); 2353 } 2354 e=quotesubst(e,vabs,vabs2,contextptr); 2355 vabs=lop(lvar(e),at_pow); 2356 vabs2=vabs; 2357 it=vabs2.begin(); itend=vabs2.end(); 2358 for (;it!=itend;++it){ 2359 gen tmp=it->_SYMBptr->feuille; 2360 if (tmp.type==_VECT && tmp._VECTptr->size()==2){ 2361 gen tmp1=simplify(tmp._VECTptr->front(),contextptr); 2362 tmp1=pow(tmp1,tmp._VECTptr->back(),contextptr); 2363 gen tmp2=normal(tmp1,contextptr); 2364 *it=tmp2; 2365 } 2366 } 2367 // try to rewrite powers with less indep. vars 2368 vecteur bases,bases2; 2369 if (vabs2.size()>1){ 2370 #ifdef NO_STDEXCEPT 2371 vecteur vabs2tmp=*tsimplify_common(vabs2,contextptr)._VECTptr; 2372 if (is_undef(vabs2tmp)){ 2373 *logptr(contextptr) << vabs2tmp << '\n'; 2374 return e_orig; 2375 } 2376 // check for rootof? 2377 vabs2=vabs2tmp; 2378 #else 2379 vabs2=*tsimplify_common(vabs2,contextptr)._VECTptr; 2380 #endif 2381 if (1){ 2382 int S=int(vabs2.size()); 2383 vector<int> base(S),expo(S); 2384 for (int i=0;i<int(S);++i){ 2385 gen g=vabs2[i]; 2386 if (!g.is_symb_of_sommet(at_pow)){ 2387 bases.push_back(g); 2388 base[i]=int(bases.size()-1); 2389 expo[i]=1; 2390 continue; 2391 } 2392 gen b=g[1],e=g[2]; 2393 if (!lop(b,at_pow).empty()){ 2394 bases.push_back(g); 2395 base[i]=int(bases.size()-1); 2396 expo[i]=1; 2397 continue; 2398 } 2399 if (e.type!=_FRAC || e._FRACptr->num!=1 || e._FRACptr->den.type!=_INT_){ 2400 bases.push_back(g); 2401 base[i]=int(bases.size()-1); 2402 expo[i]=1; 2403 continue; 2404 } 2405 int p=equalposcomp(bases,b); 2406 if (p==0){ 2407 bases.push_back(b); 2408 base[i]=int(bases.size()-1); 2409 expo[i]=e._FRACptr->den.val; 2410 continue; 2411 } 2412 base[i]=p-1; 2413 expo[i]=e._FRACptr->den.val; 2414 } 2415 // find lcm of exponent of indices i having the same base 2416 vector<int> lcms(bases.size(),1); 2417 for (int i=0;i<S;++i){ 2418 int p=base[i]; 2419 lcms[p]=(lcms[p]*long(expo[i]))/gcd(lcms[p],expo[i]); 2420 } 2421 for (int p=0;p<int(bases.size());++p){ 2422 bases[p]=symb_pow(bases[p],fraction(1,lcms[p])); 2423 bases2.push_back(identificateur(" simplify_pow"+print_INT_(p))); 2424 } 2425 // rewrite vabs2 2426 for (int i=0;i<S;++i){ 2427 int p=base[i]; 2428 if (lcms[p]==1) 2429 continue; 2430 gen basesp1=bases[p][1]; 2431 if (basesp1.type<_IDNT) 2432 continue; 2433 // bases[p]^expo[i] is rewritten as (bases[p]^1/lcms[p])^(1/(lcms[p]/expo[i])) 2434 vabs2[i]=symb_pow(bases2[p],lcms[p]/expo[i]); 2435 vabs.push_back(basesp1); 2436 vabs2.push_back(symb_pow(bases2[p],lcms[p])); 2437 } 2438 } 2439 } 2440 e=quotesubst(e,vabs,vabs2,contextptr); 2441 e=quotesubst(e,vabs,vabs2,contextptr); // second replacement because vabs2 might contain expression in vabs 2442 vecteur ve(lop(e,at_exp)); 2443 ve=lop(ve,at_atan); 2444 if (!ve.empty()) 2445 e=_exp2trig(e,contextptr); // exp(i*atan()) 2446 e=recursive_normal(e,contextptr); 2447 if (is_undef(e)) return e; 2448 if (!bases.empty()) 2449 e=quotesubst(e,bases2,bases,contextptr); 2450 // Don't touch fractional powers and absolute value anymore 2451 vabs2=lvar(e); 2452 vabs.clear(); 2453 for (unsigned int i=0;i<vabs2.size();++i){ 2454 if (vabs2[i].is_symb_of_sommet(at_abs) || vabs2[i].is_symb_of_sommet(at_pow) || vabs2[i].is_symb_of_sommet(at_integrate)){ 2455 vabs.push_back(vabs2[i]); 2456 continue; 2457 } 2458 if (vabs2[i].is_symb_of_sommet(at_exp)){ 2459 // detect sin/cos/tan/exp/ inside an exp 2460 gen tmp=_lin(_trig2exp(lvar(vabs2[i]._SYMBptr->feuille),contextptr),contextptr); 2461 if (!lop(tmp,at_exp).empty()) 2462 vabs.push_back(vabs2[i]); 2463 } 2464 } 2465 if (!vabs.empty() && debug_infolevel) 2466 *logptr(contextptr) << gettext("simplify preserving ") << vabs << '\n'; 2467 int s=int(vabs.size()); 2468 vabs2=vecteur(s); 2469 for (int i=0;i<s;++i){ 2470 // vabs2[i]=identificateur(" "+vabs[i].print(contextptr)); 2471 vabs2[i]=identificateur(" simplify_"+print_INT_(i)); 2472 } 2473 gen esave=e; 2474 // check for the presence of trig/atrig functions 2475 vecteur v1(loptab(e,sincostan_tab)); 2476 int s1=int(v1.size()),s2=int(loptab(e,asinacosatan_tab).size()); 2477 if ( s1 &&s2 ){ 2478 // retry with e texpanded 2479 gen e2=_texpand(v1,contextptr); 2480 if (e2.type==_VECT){ 2481 vecteur v2(loptab(e2,asinacosatan_tab)); 2482 if (int(v2.size())<s2){ 2483 e2=simplify(e2,contextptr); 2484 if (e2.type==_VECT && e2!=v1){ 2485 e=subst(e,v1,e2,false,contextptr); 2486 return simplify(e,contextptr); 2487 } 2488 } 2489 } 2490 } 2491 if (s1>1){ 2492 // retry with trigtan/trigcos/trigsin/halftan 2493 gen e2=recursive_normal(_trigtan(e,contextptr),contextptr); 2494 if (int(loptab(e2,sincostan_tab).size())<s1) 2495 return simplify(e2,contextptr); 2496 e2=recursive_normal(_trigcos(e,contextptr),contextptr); 2497 if (int(loptab(e2,sincostan_tab).size())<s1) 2498 return simplify(e2,contextptr); 2499 e2=recursive_normal(_trigsin(e,contextptr),contextptr); 2500 if (int(loptab(e2,sincostan_tab).size())<s1) 2501 return simplify(e2,contextptr); 2502 e2=_halftan(v1,contextptr); 2503 if (e2.type==_VECT){ 2504 vecteur w1(loptab(e2,sincostan_tab)); 2505 if (w1.size()<v1.size()){ 2506 e=subst(e,v1,e2,false,contextptr); 2507 return simplify(e,contextptr); 2508 } 2509 } 2510 } 2511 e=quotesubst(e,vabs,vabs2,contextptr); 2512 #if defined NO_STDEXCEPT // GIAC_HAS_STO_38 2513 { 2514 vecteur v2=loptab(e,asinacosatan_tab); 2515 if (!v2.empty()){ 2516 unsigned count=0; 2517 for (unsigned i=0;i<v2.size();++i){ 2518 if (v2[i].is_symb_of_sommet(at_asin)||v2[i].is_symb_of_sommet(at_acos)) 2519 count+=int(lidnt(v2[i]).size()); 2520 } 2521 if (count>1) 2522 return quotesubst(e,vabs2,vabs,contextptr); 2523 } 2524 } 2525 #endif 2526 gen g=tsimplify_noexpln(e,s1,s2,contextptr); 2527 gen glin=cklin(g,contextptr); 2528 bool glinb=glin!=g; 2529 g=glin; 2530 g=_exp2pow(g,contextptr); 2531 // if s2==0 and s1>1 and only tan, should return trigtan(exp2trig(g))? 2532 if (s1<=1 && s2<= 1){ 2533 g=quotesubst(g,vabs2,vabs,contextptr); 2534 return ratnormal(g,contextptr);//glinb?g:ratnormal(g,contextptr); 2535 } 2536 int te=taille(e,RAND_MAX); 2537 int tg=taille(g,10*te); 2538 if (tg>=10*te) 2539 return esave; 2540 // convert back to trig and atrig functions 2541 g=expln2trig(g,contextptr); 2542 if (!complex_mode(contextptr) && !has_i(g)){ 2543 if (s1){ 2544 if (v1.front().is_symb_of_sommet(at_sin)){ 2545 g=trigsin(g,contextptr); 2546 g=recursive_normal(g,contextptr); 2547 return quotesubst(g,vabs2,vabs,contextptr); 2548 } 2549 } 2550 #ifdef FXCG 2551 if (s1!=2 || !v1[0].is_symb_of_sommet(at_tan) || !v1[1].is_symb_of_sommet(at_tan)) 2552 #endif 2553 g=recursive_normal(trigcos(g,contextptr),contextptr); 2554 return quotesubst(g,vabs2,vabs,contextptr); 2555 } 2556 gen reg,img; 2557 reim(g,reg,img,contextptr); 2558 reg=recursive_normal(re(g,contextptr),contextptr); 2559 img=recursive_normal(im(g,contextptr),contextptr); 2560 if (s1){ 2561 gen g1=normal(trigcos(reg,contextptr),contextptr)+cst_i*normal(trigcos(img,contextptr),contextptr); 2562 gen g2=normal(trigsin(reg,contextptr),contextptr)+cst_i*normal(trigsin(img,contextptr),contextptr); 2563 g1=quotesubst(g1,vabs2,vabs,contextptr); 2564 g2=quotesubst(g2,vabs2,vabs,contextptr); 2565 int g1s=int(lvar(g1).size()), g2s=int(lvar(g2).size()); 2566 if (g1s!=g2s) 2567 return g1s<g2s?g1:g2; 2568 g1s=taille(g1,RAND_MAX),g2s=taille(g2,RAND_MAX); 2569 if (g1s!=g2s) 2570 return g1s<g2s?g1:g2; 2571 if (v1.front().is_symb_of_sommet(at_sin)) 2572 return g2; 2573 return g1; 2574 } 2575 g=normal(trigcos(reg,contextptr),contextptr)+cst_i*normal(trigcos(img,contextptr),contextptr); 2576 g=quotesubst(g,vabs2,vabs,contextptr); 2577 return g; 2578 } 2579 static const char _expln2trig_s []="expln2trig"; 2580 static define_unary_function_eval (__expln2trig,&expln2trig,_expln2trig_s); 2581 define_unary_function_ptr5( at_expln2trig ,alias_at_expln2trig,&__expln2trig,0,true); 2582 _simplify(const gen & args,GIAC_CONTEXT)2583 gen _simplify(const gen & args,GIAC_CONTEXT){ 2584 if ( args.type==_STRNG && args.subtype==-1) return args; 2585 gen var,res; 2586 if (is_algebraic_program(args,var,res)) 2587 return symbolic(at_program,makesequence(var,0,_simplify(res,contextptr))); 2588 if (args.type==_VECT){ 2589 vecteur & v =*args._VECTptr; 2590 int vs=int(v.size()); 2591 if ( (vs==2 || vs==3) && args.subtype==_SEQ__VECT && args[1].type==_VECT && !ckmatrix(args) && !ckmatrix(args._VECTptr->back())){ 2592 // simplify with side relations 2593 #ifdef NO_STDEXCEPT 2594 return _greduce(args,contextptr); 2595 #else 2596 try { 2597 return _greduce(args,contextptr); 2598 } 2599 catch(std::runtime_error & e){ 2600 last_evaled_argptr(contextptr)=NULL; 2601 *logptr(contextptr) << e.what() << '\n'; 2602 } 2603 #endif 2604 } 2605 return apply(args,_simplify,contextptr); 2606 } 2607 if (is_equal(args)) 2608 return apply_to_equal(args,_simplify,contextptr); 2609 int st=step_infolevel(contextptr); 2610 step_infolevel(0,contextptr); 2611 int c=calc_mode(contextptr); 2612 calc_mode(0,contextptr); 2613 vecteur sub1,sub2; 2614 surd2pow(args,sub1,sub2,contextptr); 2615 if (!sub2.empty()) 2616 sub2=subst(sub2,at_abs,at_nop,false,contextptr); 2617 res=args; 2618 if (!sub1.empty()) 2619 res=subst(res,sub1,sub2,false,contextptr); 2620 #ifdef NO_STDEXCEPT 2621 res=simplify(res,contextptr); 2622 #else 2623 try { 2624 res=simplify(res,contextptr); 2625 } catch (...){ 2626 last_evaled_argptr(contextptr)=NULL; 2627 } 2628 #endif 2629 if (!sub1.empty()) 2630 res=subst(res,sub2,sub1,false,contextptr); 2631 step_infolevel(st,contextptr); 2632 calc_mode(c,contextptr); 2633 if ( (c==1 || c==-38 || c==38) && !lop(res,at_rootof).empty()) 2634 res=ratnormal(args,contextptr); 2635 return res; 2636 } 2637 static const char _simplify_s []="simplify"; 2638 static define_unary_function_eval (__simplify,&_simplify,_simplify_s); 2639 define_unary_function_ptr5( at_simplify ,alias_at_simplify,&__simplify,0,true); 2640 trigcos(const gen & e,GIAC_CONTEXT)2641 gen trigcos(const gen & e,GIAC_CONTEXT){ 2642 return subst(simplifier(e,contextptr),pow_tab,trigcos_tab,false,contextptr); 2643 } _trigcos(const gen & args,GIAC_CONTEXT)2644 gen _trigcos(const gen & args,GIAC_CONTEXT){ 2645 if ( args.type==_STRNG && args.subtype==-1) return args; 2646 gen var,res; 2647 if (is_algebraic_program(args,var,res)) 2648 return symbolic(at_program,makesequence(var,0,_trigcos(res,contextptr))); 2649 if (is_equal(args)) 2650 return apply_to_equal(args,_trigcos,contextptr); 2651 gen g=ratnormal(_tan2sincos(args,contextptr),contextptr); 2652 return normal(trigcos(g,contextptr),contextptr); 2653 } 2654 static const char _trigcos_s []="trigcos"; 2655 static define_unary_function_eval (__trigcos,&_trigcos,_trigcos_s); 2656 define_unary_function_ptr5( at_trigcos ,alias_at_trigcos,&__trigcos,0,true); 2657 trigsin(const gen & e,GIAC_CONTEXT)2658 gen trigsin(const gen & e,GIAC_CONTEXT){ 2659 return subst(simplifier(e,contextptr),pow_tab,trigsin_tab,false,contextptr); 2660 } _trigsin(const gen & args,GIAC_CONTEXT)2661 gen _trigsin(const gen & args,GIAC_CONTEXT){ 2662 if ( args.type==_STRNG && args.subtype==-1) return args; 2663 gen var,res; 2664 if (is_algebraic_program(args,var,res)) 2665 return symbolic(at_program,makesequence(var,0,_trigsin(res,contextptr))); 2666 if (is_equal(args)) 2667 return apply_to_equal(args,_trigsin,contextptr); 2668 gen g=ratnormal(_tan2sincos(args,contextptr),contextptr); 2669 return normal(trigsin(g,contextptr),contextptr); 2670 } 2671 static const char _trigsin_s []="trigsin"; 2672 static define_unary_function_eval (__trigsin,&_trigsin,_trigsin_s); 2673 define_unary_function_ptr5( at_trigsin ,alias_at_trigsin,&__trigsin,0,true); 2674 trigtan(const gen & e,GIAC_CONTEXT)2675 gen trigtan(const gen & e,GIAC_CONTEXT){ 2676 return subst(simplifier(e,contextptr),pow_tab,trigtan_tab,false,contextptr); 2677 } 2678 gen _sin2costan(const gen & args,GIAC_CONTEXT); _trigtan(const gen & args,GIAC_CONTEXT)2679 gen _trigtan(const gen & args,GIAC_CONTEXT){ 2680 if ( args.type==_STRNG && args.subtype==-1) return args; 2681 gen var,res; 2682 if (is_algebraic_program(args,var,res)) 2683 return symbolic(at_program,makesequence(var,0,_trigtan(res,contextptr))); 2684 if (is_equal(args)) 2685 return apply_to_equal(args,_trigtan,contextptr); 2686 gen g=ratnormal(_sin2costan(args,contextptr),contextptr); 2687 return normal(trigtan(g,contextptr),contextptr); 2688 } 2689 static const char _trigtan_s []="trigtan"; 2690 static define_unary_function_eval (__trigtan,&_trigtan,_trigtan_s); 2691 define_unary_function_ptr5( at_trigtan ,alias_at_trigtan,&__trigtan,0,true); 2692 tan2sincos(const gen & e,GIAC_CONTEXT)2693 gen tan2sincos(const gen & e,GIAC_CONTEXT){ 2694 return subst(e,tan_tab,tan2sincos_tab,false,contextptr); 2695 } _tan2sincos(const gen & args,GIAC_CONTEXT)2696 gen _tan2sincos(const gen & args,GIAC_CONTEXT){ 2697 if ( args.type==_STRNG && args.subtype==-1) return args; 2698 gen var,res; 2699 if (is_algebraic_program(args,var,res)) 2700 return symbolic(at_program,makesequence(var,0,_tan2sincos(res,contextptr))); 2701 if (is_equal(args)) 2702 return apply_to_equal(args,_tan2sincos,contextptr); 2703 return tan2sincos(args,contextptr); 2704 } 2705 static const char _tan2sincos_s []="tan2sincos"; 2706 static define_unary_function_eval (__tan2sincos,&_tan2sincos,_tan2sincos_s); 2707 define_unary_function_ptr5( at_tan2sincos ,alias_at_tan2sincos,&__tan2sincos,0,true); 2708 sintocostan(const gen & e,GIAC_CONTEXT)2709 static gen sintocostan(const gen & e,GIAC_CONTEXT){ 2710 return tan(e,contextptr)*cos(e,contextptr); 2711 } sin2_costan(const gen & e,GIAC_CONTEXT)2712 gen sin2_costan(const gen & e,GIAC_CONTEXT){ 2713 vector< gen_op_context > sin2costan_v(1,sintocostan); 2714 vector<const unary_function_ptr *> sin_v(1,at_sin); 2715 return subst(e,sin_v,sin2costan_v,false,contextptr); 2716 } _sin2costan(const gen & args,GIAC_CONTEXT)2717 gen _sin2costan(const gen & args,GIAC_CONTEXT){ 2718 if ( args.type==_STRNG && args.subtype==-1) return args; 2719 gen var,res; 2720 if (is_algebraic_program(args,var,res)) 2721 return symbolic(at_program,makesequence(var,0,_sin2costan(res,contextptr))); 2722 if (is_equal(args)) 2723 return apply_to_equal(args,_sin2costan,contextptr); 2724 return sin2_costan(args,contextptr); 2725 } 2726 static const char _sin2costan_s []="sin2costan"; 2727 static define_unary_function_eval (__sin2costan,&_sin2costan,_sin2costan_s); 2728 define_unary_function_ptr5( at_sin2costan ,alias_at_sin2costan,&__sin2costan,0,true); 2729 costosintan(const gen & e,GIAC_CONTEXT)2730 static gen costosintan(const gen & e,GIAC_CONTEXT){ 2731 return sin(e,contextptr)/tan(e,contextptr); 2732 } cos2sintan(const gen & e,GIAC_CONTEXT)2733 gen cos2sintan(const gen & e,GIAC_CONTEXT){ 2734 vector< gen_op_context > cos2sintan_v(1,costosintan); 2735 vector<const unary_function_ptr *> cos_v(1,at_cos); 2736 return subst(e,cos_v,cos2sintan_v,false,contextptr); 2737 } _cos2sintan(const gen & args,GIAC_CONTEXT)2738 gen _cos2sintan(const gen & args,GIAC_CONTEXT){ 2739 if ( args.type==_STRNG && args.subtype==-1) return args; 2740 gen var,res; 2741 if (is_algebraic_program(args,var,res)) 2742 return symbolic(at_program,makesequence(var,0,_cos2sintan(res,contextptr))); 2743 if (is_equal(args)) 2744 return apply_to_equal(args,_cos2sintan,contextptr); 2745 return cos2sintan(args,contextptr); 2746 } 2747 static const char _cos2sintan_s []="cos2sintan"; 2748 static define_unary_function_eval (__cos2sintan,&_cos2sintan,_cos2sintan_s); 2749 define_unary_function_ptr5( at_cos2sintan ,alias_at_cos2sintan,&__cos2sintan,0,true); 2750 tan2sincos2(const gen & e,GIAC_CONTEXT)2751 gen tan2sincos2(const gen & e,GIAC_CONTEXT){ 2752 return subst(e,tan_tab,tan2sincos2_tab,false,contextptr); 2753 } _tan2sincos2(const gen & args,GIAC_CONTEXT)2754 gen _tan2sincos2(const gen & args,GIAC_CONTEXT){ 2755 if ( args.type==_STRNG && args.subtype==-1) return args; 2756 gen var,res; 2757 if (is_algebraic_program(args,var,res)) 2758 return symbolic(at_program,makesequence(var,0,_tan2sincos2(res,contextptr))); 2759 if (is_equal(args)) 2760 return apply_to_equal(args,_tan2sincos2,contextptr); 2761 return tan2sincos2(args,contextptr); 2762 } 2763 static const char _tan2sincos2_s []="tan2sincos2"; 2764 static define_unary_function_eval (__tan2sincos2,&_tan2sincos2,_tan2sincos2_s); 2765 define_unary_function_ptr5( at_tan2sincos2 ,alias_at_tan2sincos2,&__tan2sincos2,0,true); 2766 tan2cossin2(const gen & e,GIAC_CONTEXT)2767 gen tan2cossin2(const gen & e,GIAC_CONTEXT){ 2768 return subst(e,tan_tab,tan2cossin2_tab,false,contextptr); 2769 } _tan2cossin2(const gen & args,GIAC_CONTEXT)2770 gen _tan2cossin2(const gen & args,GIAC_CONTEXT){ 2771 if ( args.type==_STRNG && args.subtype==-1) return args; 2772 gen var,res; 2773 if (is_algebraic_program(args,var,res)) 2774 return symbolic(at_program,makesequence(var,0,_tan2cossin2(res,contextptr))); 2775 if (is_equal(args)) 2776 return apply_to_equal(args,_tan2cossin2,contextptr); 2777 return tan2cossin2(args,contextptr); 2778 } 2779 static const char _tan2cossin2_s []="tan2cossin2"; 2780 static define_unary_function_eval (__tan2cossin2,&_tan2cossin2,_tan2cossin2_s); 2781 define_unary_function_ptr5( at_tan2cossin2 ,alias_at_tan2cossin2,&__tan2cossin2,0,true); 2782 tcollect(const gen & args,bool output_cos,GIAC_CONTEXT)2783 gen tcollect(const gen & args,bool output_cos,GIAC_CONTEXT){ 2784 gen errcode=checkanglemode(contextptr); 2785 if (is_undef(errcode)) return errcode; 2786 vecteur v; 2787 tlin(args,v,contextptr); 2788 // v= [coeff, sin/cos/1] 2789 int s=int(v.size()); 2790 vector<int> deja; 2791 gen res,argu; 2792 for (int i=1;i<s;i+=2){ 2793 if (equalposcomp(deja,i)) 2794 continue; 2795 if (v[i].type!=_SYMB){ 2796 res = res + v[i-1]*v[i]; 2797 continue; 2798 } 2799 argu=v[i]._SYMBptr->feuille; 2800 int j=i+2; 2801 for (;j<s;j+=2){ 2802 if ( (v[j].type==_SYMB) && (v[j]._SYMBptr->feuille==argu) ) 2803 break; 2804 } 2805 if (j>=s) 2806 res = res + v[i-1]*v[i]; 2807 else { // found 2 trig terms with the same arg, combine 2808 deja.push_back(j); 2809 gen coscoeff,sincoeff; 2810 if (v[i]._SYMBptr->sommet==at_cos){ 2811 coscoeff=v[i-1]; 2812 sincoeff=v[j-1]; 2813 } 2814 else { 2815 coscoeff=v[j-1]; 2816 sincoeff=v[i-1]; 2817 } 2818 gen newcoeff=normal(sqrt(pow(coscoeff,2)+pow(sincoeff,2),contextptr),contextptr); 2819 if (output_cos){ 2820 gen angle=normal(arg(coscoeff+cst_i*sincoeff,contextptr),contextptr); 2821 res=res+newcoeff*symbolic(at_cos,argu-angle); 2822 } 2823 else { 2824 gen angle=normal(arg(sincoeff+cst_i*coscoeff,contextptr),contextptr); 2825 res=res+newcoeff*symbolic(at_sin,argu+angle); 2826 } 2827 } 2828 } 2829 return res; 2830 } tcollect(const gen & args,GIAC_CONTEXT)2831 gen tcollect(const gen & args,GIAC_CONTEXT){ 2832 return tcollect(args,true,contextptr); 2833 } _tcollect(const gen & args,GIAC_CONTEXT)2834 gen _tcollect(const gen & args,GIAC_CONTEXT){ 2835 if ( args.type==_STRNG && args.subtype==-1) return args; 2836 gen var,res; 2837 if (is_algebraic_program(args,var,res)) 2838 return symbolic(at_program,makesequence(var,0,_tcollect(res,contextptr))); 2839 if (is_equal(args)) 2840 return apply_to_equal(args,_tcollect,contextptr); 2841 return apply(args,contextptr,tcollect); 2842 } 2843 static const char _tcollect_s []="tcollect"; 2844 static define_unary_function_eval (__tcollect,&_tcollect,_tcollect_s); 2845 define_unary_function_ptr5( at_tcollect ,alias_at_tcollect,&__tcollect,0,true); 2846 tcollectsin(const gen & args,GIAC_CONTEXT)2847 gen tcollectsin(const gen & args,GIAC_CONTEXT){ 2848 return tcollect(args,false,contextptr); 2849 } _tcollectsin(const gen & args,GIAC_CONTEXT)2850 gen _tcollectsin(const gen & args,GIAC_CONTEXT){ 2851 if ( args.type==_STRNG && args.subtype==-1) return args; 2852 gen var,res; 2853 if (is_algebraic_program(args,var,res)) 2854 return symbolic(at_program,makesequence(var,0,_tcollectsin(res,contextptr))); 2855 if (is_equal(args)) 2856 return apply_to_equal(args,_tcollectsin,contextptr); 2857 return apply(args,contextptr,tcollectsin); 2858 } 2859 static const char _tcollectsin_s []="tcollectsin"; 2860 static define_unary_function_eval (__tcollectsin,&_tcollectsin,_tcollectsin_s); 2861 define_unary_function_ptr5( at_tcollectsin ,alias_at_tcollectsin,&__tcollectsin,0,true); 2862 2863 static const char _rassembler_trigo_s []="rassembler_trigo"; 2864 static define_unary_function_eval (__rassembler_trigo,&_tcollect,_rassembler_trigo_s); 2865 define_unary_function_ptr5( at_rassembler_trigo ,alias_at_rassembler_trigo,&__rassembler_trigo,0,true); 2866 postlncollect(vecteur & rescoeff,vecteur & resargln,vecteur & res,GIAC_CONTEXT)2867 static void postlncollect(vecteur & rescoeff,vecteur & resargln,vecteur & res,GIAC_CONTEXT){ 2868 gen d; 2869 lcmdeno(rescoeff,d,contextptr); 2870 const_iterateur it=rescoeff.begin(),itend=rescoeff.end(); 2871 const_iterateur jt=resargln.begin(); 2872 gen lnint(plus_one); 2873 res.reserve(2*(itend-it)); 2874 for (;it!=itend;++it,++jt){ 2875 if (is_zero(*it)) 2876 continue; 2877 if (it->type==_INT_ && absint(it->val)<MAX_ALG_EXT_ORDER_SIZE) 2878 lnint=lnint*pow(*jt,it->val); 2879 else { 2880 res.push_back(*it/d); 2881 res.push_back(*jt); 2882 } 2883 } 2884 if (!is_one(lnint)){ 2885 res.push_back(inv(d,contextptr)); 2886 res.push_back(lnint); 2887 } 2888 } 2889 2890 // -> [ coeffln argln ...] inlncollect(const gen & args,GIAC_CONTEXT)2891 static vecteur inlncollect(const gen & args,GIAC_CONTEXT){ 2892 if (args.type!=_SYMB) 2893 return makevecteur(args,symbolic(at_exp,1)); 2894 const unary_function_ptr & u=args._SYMBptr->sommet; 2895 gen g=args._SYMBptr->feuille; 2896 if (u==at_ln) 2897 return makevecteur(1,g); 2898 if (u==at_plus){ 2899 if (g.type!=_VECT) 2900 return inlncollect(g,contextptr); 2901 vecteur rescoeff,resargln; 2902 const_iterateur it=g._VECTptr->begin(),itend=g._VECTptr->end(),jt,jtend; 2903 for (;it!=itend;++it){ 2904 vecteur tmp(inlncollect(*it,contextptr)); 2905 jt=tmp.begin(); 2906 jtend=tmp.end(); 2907 for (;jt!=jtend;jt+=2){ 2908 int i; 2909 if ( (i=equalposcomp(resargln,*(jt+1))) ){ 2910 rescoeff[i-1]=normal(rescoeff[i-1]+*jt,contextptr); 2911 } 2912 else { 2913 rescoeff.push_back(*jt); 2914 resargln.push_back(*(jt+1)); 2915 } 2916 } 2917 } 2918 vecteur res; 2919 postlncollect(rescoeff,resargln,res,contextptr); 2920 return res; 2921 } 2922 if (u==at_neg){ 2923 vecteur tmp(inlncollect(g,contextptr)); 2924 iterateur it=tmp.begin(),itend=tmp.end(); 2925 for (;it!=itend;it+=2){ 2926 *it=-*it; 2927 } 2928 return tmp; 2929 } 2930 if (u==at_prod){ 2931 if (g.type!=_VECT) 2932 return inlncollect(g,contextptr); 2933 // is there a ln in g? 2934 const_iterateur it=g._VECTptr->begin(),it0=g._VECTptr->begin(),itend=g._VECTptr->end(); 2935 for (;it!=itend;++it){ 2936 if ( (it->type==_SYMB) && (it->_SYMBptr->sommet==at_ln) ) 2937 break; 2938 } 2939 if (it!=itend){ 2940 vecteur v(mergevecteur(vecteur(it0,it),vecteur(it+1,itend))); 2941 gen tmp=_prod(v,contextptr); 2942 vecteur vcoeff(1,tmp),varg(1,it->_SYMBptr->feuille),res; 2943 postlncollect(vcoeff,varg,res,contextptr); 2944 return res; 2945 } 2946 } 2947 return makevecteur((u)(_lncollect(g,contextptr),contextptr), 2948 symbolic(at_exp,1)); 2949 } lncollect(const gen & args,GIAC_CONTEXT)2950 gen lncollect(const gen & args,GIAC_CONTEXT){ 2951 vecteur v(inlncollect(args,contextptr)); 2952 gen res; 2953 const_iterateur it=v.begin(),itend=v.end(); 2954 for (;it!=itend;it+=2){ 2955 res = res + (*it) * ln (*(it+1),contextptr); 2956 } 2957 return res; 2958 } _lncollect(const gen & args,GIAC_CONTEXT)2959 gen _lncollect(const gen & args,GIAC_CONTEXT){ 2960 if ( args.type==_STRNG && args.subtype==-1) return args; 2961 gen var,res; 2962 if (is_algebraic_program(args,var,res)) 2963 return symbolic(at_program,makesequence(var,0,_lncollect(res,contextptr))); 2964 if (is_equal(args)) 2965 return apply_to_equal(args,_lncollect,contextptr); 2966 return apply(args,lncollect,contextptr); 2967 } 2968 static const char _lncollect_s []="lncollect"; 2969 static define_unary_function_eval (__lncollect,&_lncollect,_lncollect_s); 2970 define_unary_function_ptr5( at_lncollect ,alias_at_lncollect,&__lncollect,0,true); 2971 powexpand(const gen & e,GIAC_CONTEXT)2972 gen powexpand(const gen & e,GIAC_CONTEXT){ 2973 return subst(e,pow_tab,powexpand_tab,false,contextptr); 2974 } _powexpand(const gen & args,GIAC_CONTEXT)2975 gen _powexpand(const gen & args,GIAC_CONTEXT){ 2976 if ( args.type==_STRNG && args.subtype==-1) return args; 2977 gen var,res; 2978 if (is_algebraic_program(args,var,res)) 2979 return symbolic(at_program,makesequence(var,0,_powexpand(res,contextptr))); 2980 if (is_equal(args)) 2981 return apply_to_equal(args,_powexpand,contextptr); 2982 return apply(args,powexpand,contextptr); 2983 } 2984 static const char _powexpand_s []="powexpand"; 2985 static define_unary_function_eval (__powexpand,&_powexpand,_powexpand_s); 2986 define_unary_function_ptr5( at_powexpand ,alias_at_powexpand,&__powexpand,0,true); 2987 exp2pow(const gen & e,GIAC_CONTEXT)2988 gen exp2pow(const gen & e,GIAC_CONTEXT){ 2989 return subst(e,exp_tab,exp2power_tab,false,contextptr); 2990 } _exp2pow(const gen & args,GIAC_CONTEXT)2991 gen _exp2pow(const gen & args,GIAC_CONTEXT){ 2992 if ( args.type==_STRNG && args.subtype==-1) return args; 2993 gen var,res; 2994 if (is_algebraic_program(args,var,res)) 2995 return symbolic(at_program,makesequence(var,0,_exp2pow(res,contextptr))); 2996 if (is_equal(args)) 2997 return apply_to_equal(args,_exp2pow,contextptr); 2998 return apply(args,exp2pow,contextptr); 2999 } 3000 static const char _exp2pow_s []="exp2pow"; 3001 static define_unary_function_eval (__exp2pow,&_exp2pow,_exp2pow_s); 3002 define_unary_function_ptr5( at_exp2pow ,alias_at_exp2pow,&__exp2pow,0,true); 3003 factor_xn(const gen & args,const gen & x,GIAC_CONTEXT)3004 gen factor_xn(const gen & args,const gen & x,GIAC_CONTEXT){ 3005 vecteur l(1,x); 3006 lvar(args,l); 3007 gen temp=e2r(args,l,contextptr); 3008 gen n,d; 3009 fxnd(temp,n,d); 3010 l.erase(l.begin()); 3011 vecteur nv(gen2vecteur(r2e(polynome2poly1(n,1),l,contextptr))); 3012 vecteur dv(gen2vecteur(r2e(polynome2poly1(d,1),l,contextptr))); 3013 int ns=int(nv.size()),ds=int(dv.size()); 3014 int n1=ns-ds; 3015 return pow(x,n1)*symb_horner(nv,x,ns-1)/symb_horner(dv,x,ds-1); 3016 } factor_xn(const gen & args,GIAC_CONTEXT)3017 gen factor_xn(const gen & args,GIAC_CONTEXT){ 3018 return factor_xn(args,vx_var,contextptr); 3019 } _factor_xn(const gen & args,GIAC_CONTEXT)3020 gen _factor_xn(const gen & args,GIAC_CONTEXT){ 3021 if ( args.type==_STRNG && args.subtype==-1) return args; 3022 gen var,res; 3023 if (is_algebraic_program(args,var,res)) 3024 return symbolic(at_program,makesequence(var,0,_factor_xn(res,contextptr))); 3025 if (is_equal(args)) 3026 return apply_to_equal(args,_factor_xn,contextptr); 3027 return apply(args,factor_xn,contextptr); 3028 } 3029 static const char _factor_xn_s []="factor_xn"; 3030 static define_unary_function_eval (__factor_xn,&_factor_xn,_factor_xn_s); 3031 define_unary_function_ptr5( at_factor_xn ,alias_at_factor_xn,&__factor_xn,0,true); 3032 3033 static const char _developper_s []="developper"; 3034 static define_unary_function_eval (__developper,&expand,_developper_s); 3035 define_unary_function_ptr5( at_developper ,alias_at_developper,&__developper,0,true); 3036 3037 static const char _factoriser_s []="factoriser"; 3038 static define_unary_function_eval (__factoriser,&_factor,_factoriser_s); 3039 define_unary_function_ptr5( at_factoriser ,alias_at_factoriser,&__factoriser,0,true); 3040 3041 static const char _factoriser_xn_s []="factoriser_xn"; 3042 static define_unary_function_eval (__factoriser_xn,&_factor_xn,_factoriser_xn_s); 3043 define_unary_function_ptr5( at_factoriser_xn ,alias_at_factoriser_xn,&__factoriser_xn,0,true); 3044 3045 static const char _resoudre_s []="resoudre"; 3046 static define_unary_function_eval_quoted (__resoudre,&_solve,_resoudre_s); 3047 define_unary_function_ptr5( at_resoudre ,alias_at_resoudre,&__resoudre,_QUOTE_ARGUMENTS,true); 3048 3049 static const char _substituer_s []="substituer"; 3050 static define_unary_function_eval (__substituer,&_subst,_substituer_s); 3051 define_unary_function_ptr5( at_substituer ,alias_at_substituer,&__substituer,0,true); 3052 3053 static const char _deriver_s []="deriver"; 3054 static define_unary_function_eval (__deriver,&_derive,_deriver_s); 3055 define_unary_function_ptr5( at_deriver ,alias_at_deriver,&__deriver,0,true); 3056 3057 static const char _integrer_s []="integrer"; 3058 static define_unary_function_eval (__integrer,&_integrate,_integrer_s); 3059 define_unary_function_ptr5( at_integrer ,alias_at_integrer,&__integrer,0,true); 3060 3061 static const char _limite_s []="limite"; 3062 static define_unary_function_eval (__limite,&_limit,_limite_s); 3063 define_unary_function_ptr5( at_limite ,alias_at_limite,&__limite,0,true); 3064 find_conjugates(const gen & g,vecteur & v_in,vecteur & v_out)3065 static void find_conjugates(const gen & g,vecteur & v_in,vecteur & v_out){ 3066 v_in.clear(); 3067 v_out.clear(); 3068 vecteur v1(lop(g,at_pow)); 3069 int n=int(v1.size()); 3070 for (int i=0;i<n;++i){ 3071 gen & tmp =v1[i]._SYMBptr->feuille; 3072 if (tmp.type!=_VECT) 3073 continue; 3074 vecteur & v=*tmp._VECTptr; 3075 if (v.size()==2 && v.back().type==_FRAC && v.back()._FRACptr->den==plus_two){ 3076 v_in.push_back(v1[i]); 3077 v_out.push_back(-v1[i]); 3078 break ; // change only 1 sqrt 3079 } 3080 } 3081 } 3082 _mult_conjugate(const gen & g0,GIAC_CONTEXT)3083 gen _mult_conjugate(const gen & g0,GIAC_CONTEXT){ 3084 if ( g0.type==_STRNG && g0.subtype==-1) return g0; 3085 gen g(g0); 3086 bool deno=true; 3087 if (g.type==_VECT && g._VECTptr->size()==2){ 3088 if (g._VECTptr->back()==at_numer) 3089 deno=false; 3090 g=g._VECTptr->front(); 3091 } 3092 gen n,d; 3093 vecteur v_in,v_out; 3094 prod2frac(g,v_in,v_out); 3095 n=_prod(v_in,contextptr); 3096 d=_prod(v_out,contextptr); 3097 find_conjugates(d,v_in,v_out); 3098 // search sqrt in d 3099 if (!deno || v_in.empty()){ 3100 find_conjugates(n,v_in,v_out); 3101 gen mult=subst(n,v_in,v_out,false,contextptr); 3102 n=n*mult; // n=simplify(n*mult); 3103 d=d*mult; 3104 } 3105 else { 3106 gen mult=subst(d,v_in,v_out,false,contextptr); 3107 d=d*mult; // d=simplify(d*mult); 3108 n=n*mult; 3109 } 3110 return n/d; 3111 } 3112 static const char _mult_conjugate_s []="mult_conjugate"; 3113 static define_unary_function_eval (__mult_conjugate,&_mult_conjugate,_mult_conjugate_s); 3114 define_unary_function_ptr5( at_mult_conjugate ,alias_at_mult_conjugate,&__mult_conjugate,0,true); 3115 _mult_c_conjugate(const gen & g0,GIAC_CONTEXT)3116 gen _mult_c_conjugate(const gen & g0,GIAC_CONTEXT){ 3117 if ( g0.type==_STRNG && g0.subtype==-1) return g0; 3118 gen g(g0); 3119 bool deno=true; 3120 if (g.type==_VECT && g._VECTptr->size()==2){ 3121 if (g._VECTptr->back()==at_numer) 3122 deno=false; 3123 g=g._VECTptr->front(); 3124 } 3125 gen n,d; 3126 vecteur v_in,v_out; 3127 prod2frac(g,v_in,v_out); 3128 n=_prod(v_in,contextptr); 3129 d=_prod(v_out,contextptr); 3130 if (deno && !is_zero(im(d,contextptr))){ 3131 gen mult=conj(d,contextptr); 3132 if (n==1) 3133 n=mult; 3134 else 3135 n=symbolic(at_prod,gen(makevecteur(n,mult),_SEQ__VECT)); 3136 if (d==1) 3137 d=mult; 3138 else 3139 d=symbolic(at_prod,gen(makevecteur(d,mult),_SEQ__VECT)); 3140 } 3141 else { 3142 if (!is_zero(im(n,contextptr))){ 3143 gen mult=conj(n,contextptr); 3144 if (n==1) 3145 n=mult; 3146 else 3147 n=symbolic(at_prod,gen(makevecteur(n,mult),_SEQ__VECT)); 3148 if (d==1) 3149 d=mult; 3150 else 3151 d=symbolic(at_prod,gen(makevecteur(d,mult),_SEQ__VECT)); 3152 } 3153 } 3154 return n/d; 3155 } 3156 static const char _mult_c_conjugate_s []="mult_c_conjugate"; 3157 static define_unary_function_eval_quoted (__mult_c_conjugate,&_mult_c_conjugate,_mult_c_conjugate_s); 3158 define_unary_function_ptr5( at_mult_c_conjugate ,alias_at_mult_c_conjugate,&__mult_c_conjugate,_QUOTE_ARGUMENTS,true); 3159 3160 static const char _multiplier_conjugue_s []="multiplier_conjugue"; 3161 static define_unary_function_eval (__multiplier_conjugue,&_mult_conjugate,_multiplier_conjugue_s); 3162 define_unary_function_ptr5( at_multiplier_conjugue ,alias_at_multiplier_conjugue,&__multiplier_conjugue,0,true); 3163 3164 static const char _multiplier_conjugue_complexe_s []="multiplier_conjugue_complexe"; 3165 static define_unary_function_eval (__multiplier_conjugue_complexe,&_mult_c_conjugate,_multiplier_conjugue_complexe_s); 3166 define_unary_function_ptr5( at_multiplier_conjugue_complexe ,alias_at_multiplier_conjugue_complexe,&__multiplier_conjugue_complexe,0,true); 3167 _combine(const gen & args,const context * contextptr)3168 gen _combine(const gen & args,const context * contextptr){ 3169 if ( args.type==_STRNG && args.subtype==-1) return args; 3170 if (args.type!=_VECT) 3171 return gensizeerr(contextptr); 3172 vecteur & v=*args._VECTptr; 3173 int s=int(v.size()); 3174 if (s<2) 3175 return gensizeerr(contextptr); 3176 gen & f=v[s-1]; 3177 gen g=v.front(); 3178 if (f.type==_FUNC){ 3179 if (f==at_sincos || f==at_sin || f==at_cos) 3180 return tcollect(g,contextptr); 3181 if (f==at_exp) 3182 return _lin(g,contextptr); 3183 if (f==at_ln) 3184 return lncollect(g,contextptr); 3185 return gensizeerr(contextptr); 3186 } 3187 if (f.type==_INT_ && f.val>=0) { 3188 int i=f.val; 3189 if (f.subtype==_INT_MAPLECONVERSION){ 3190 switch (i){ 3191 case _TRIG: 3192 return tcollect(g,contextptr); 3193 case _EXPLN: 3194 return _lin(lncollect(g,contextptr),contextptr); 3195 default: 3196 return gensizeerr(contextptr); 3197 } 3198 } 3199 } 3200 return gensizeerr(contextptr); 3201 } 3202 static const char _combine_s []="combine"; 3203 static define_unary_function_eval (__combine,&_combine,_combine_s); 3204 define_unary_function_ptr5( at_combine ,alias_at_combine,&__combine,0,true); 3205 rectangular2polar(const gen & g,const context * contextptr)3206 static gen rectangular2polar(const gen & g,const context * contextptr){ 3207 gen args=remove_at_pnt(g); 3208 gen module=abs(args,contextptr),argument=arg(args,contextptr); 3209 module=normal(module,contextptr); 3210 if (is_zero(argument)) 3211 return module; 3212 gen res=module*symbolic(at_exp,cst_i*argument); 3213 if (calc_mode(contextptr)==1) 3214 return symb_quote(res); 3215 return res; 3216 } 3217 _rectangular2polar(const gen & args,const context * contextptr)3218 gen _rectangular2polar(const gen & args,const context * contextptr){ 3219 if ( args.type==_STRNG && args.subtype==-1) return args; 3220 return apply(args,rectangular2polar,contextptr); 3221 } 3222 static const char _rectangular2polar_s []="rectangular2polar"; 3223 static define_unary_function_eval (__rectangular2polar,&_rectangular2polar,_rectangular2polar_s); 3224 define_unary_function_ptr5( at_rectangular2polar ,alias_at_rectangular2polar,&__rectangular2polar,0,true); 3225 polar2rectangular(const gen & g,const context * contextptr)3226 static gen polar2rectangular(const gen & g,const context * contextptr){ 3227 gen args=remove_at_pnt(g); 3228 gen reel=re(args,contextptr), imag=im(args,contextptr); 3229 return reel+cst_i*imag; 3230 } 3231 _polar2rectangular(const gen & args,const context * contextptr)3232 gen _polar2rectangular(const gen & args,const context * contextptr){ 3233 if ( args.type==_STRNG && args.subtype==-1) return args; 3234 return apply(args,polar2rectangular,contextptr); 3235 } 3236 static const char _polar2rectangular_s []="polar2rectangular"; 3237 static define_unary_function_eval (__polar2rectangular,&_polar2rectangular,_polar2rectangular_s); 3238 define_unary_function_ptr5( at_polar2rectangular ,alias_at_polar2rectangular,&__polar2rectangular,0,true); 3239 3240 // x,y,z->rho,theta in [-pi/2..pi/2],psi 3241 // x=rho*cos(theta)*cos(phi), y=rho*cos(theta)*sin(phi), z=rho*sin(theta) 3242 // rho=sqrt(x^2+y^2+z^2) rectangular2spherical(const gen & g,const context * contextptr)3243 static gen rectangular2spherical(const gen & g,const context * contextptr){ 3244 if (g.type!=_VECT || g._VECTptr->size()!=3) 3245 return gensizeerr(contextptr); 3246 vecteur & v =*g._VECTptr; 3247 gen x =v[0],y=v[1],z=v[2]; 3248 gen rho2=x*x+y*y+z*z; 3249 gen rho=sqrt(rho2,contextptr); 3250 gen theta=asin(z/rho,contextptr); 3251 gen phi=arg(x+cst_i*y,contextptr); 3252 return makevecteur(rho,theta,phi); 3253 } 3254 _rectangular2spherical(const gen & args,const context * contextptr)3255 gen _rectangular2spherical(const gen & args,const context * contextptr){ 3256 if ( args.type==_STRNG && args.subtype==-1) return args; 3257 if (args.type==_VECT && args._VECTptr->size()==3 && args._VECTptr->front().type!=_VECT) 3258 return rectangular2spherical(args,contextptr); 3259 return apply(args,rectangular2spherical,contextptr); 3260 } 3261 static const char _rectangular2spherical_s []="rectangular2spherical"; 3262 static define_unary_function_eval (__rectangular2spherical,&_rectangular2spherical,_rectangular2spherical_s); 3263 define_unary_function_ptr5( at_rectangular2spherical ,alias_at_rectangular2spherical,&__rectangular2spherical,0,true); 3264 spherical2rectangular(const gen & g,const context * contextptr)3265 static gen spherical2rectangular(const gen & g,const context * contextptr){ 3266 if (g.type!=_VECT || g._VECTptr->size()!=3) 3267 return gensizeerr(contextptr); 3268 vecteur & v =*g._VECTptr; 3269 gen rho=v[0],theta=v[1],phi=v[2]; 3270 gen rhocos=rho*cos(theta,contextptr); 3271 return makevecteur(rhocos*cos(phi,contextptr),rhocos*sin(phi,contextptr),rho*sin(theta,contextptr)); 3272 } 3273 _spherical2rectangular(const gen & args,const context * contextptr)3274 gen _spherical2rectangular(const gen & args,const context * contextptr){ 3275 if ( args.type==_STRNG && args.subtype==-1) return args; 3276 if (args.type==_VECT && args._VECTptr->size()==3 && args._VECTptr->front().type!=_VECT) 3277 return spherical2rectangular(args,contextptr); 3278 return apply(args,spherical2rectangular,contextptr); 3279 } 3280 static const char _spherical2rectangular_s []="spherical2rectangular"; 3281 static define_unary_function_eval (__spherical2rectangular,&_spherical2rectangular,_spherical2rectangular_s); 3282 define_unary_function_ptr5( at_spherical2rectangular ,alias_at_spherical2rectangular,&__spherical2rectangular,0,true); 3283 heavisidetosign(const gen & args,GIAC_CONTEXT)3284 static gen heavisidetosign(const gen & args,GIAC_CONTEXT){ 3285 return (1+sign(args,contextptr))/2; 3286 } 3287 const gen_op_context heaviside2sign_tab[]={heavisidetosign,0}; 3288 Heavisidetosign(const gen & args,GIAC_CONTEXT)3289 gen Heavisidetosign(const gen & args,GIAC_CONTEXT){ 3290 return subst(args,Heaviside_tab,heaviside2sign_tab,false,contextptr); 3291 // return subst(args,vector<const unary_function_ptr *>(1,at_Heaviside), vector< gen_op_context >(1,heavisidetosign),false,contextptr); 3292 } _Heavisidetosign(const gen & args,GIAC_CONTEXT)3293 gen _Heavisidetosign(const gen & args,GIAC_CONTEXT){ 3294 if ( args.type==_STRNG && args.subtype==-1) return args; 3295 if (is_equal(args)) 3296 return apply_to_equal(args,Heavisidetosign,contextptr); 3297 return apply(args,Heavisidetosign,contextptr); 3298 } 3299 heavisidetopiecewise(const gen & args,GIAC_CONTEXT)3300 static gen heavisidetopiecewise(const gen & args,GIAC_CONTEXT){ 3301 return symbolic(at_piecewise,makesequence(symb_superieur_strict(args,0),1,0)); 3302 } 3303 const gen_op_context heaviside2piecewise_tab[]={heavisidetopiecewise,0}; 3304 Heavisidetopiecewise(const gen & args,GIAC_CONTEXT)3305 gen Heavisidetopiecewise(const gen & args,GIAC_CONTEXT){ 3306 return subst(args,Heaviside_tab,heaviside2piecewise_tab,false,contextptr); 3307 // return subst(args,vector<const unary_function_ptr *>(1,at_Heaviside), vector< gen_op_context >(1,heavisidetopiecewise),false,contextptr); 3308 } _Heavisidetopiecewise(const gen & args,GIAC_CONTEXT)3309 gen _Heavisidetopiecewise(const gen & args,GIAC_CONTEXT){ 3310 if ( args.type==_STRNG && args.subtype==-1) return args; 3311 if (is_equal(args)) 3312 return apply_to_equal(args,Heavisidetopiecewise,contextptr); 3313 return apply(args,Heavisidetopiecewise,contextptr); 3314 } 3315 3316 #if defined FXCG || !defined USE_GMP_REPLACEMENTS 3317 // find simplest between some trig simplifications, by Luka Marohnić _trigsimplify(const gen & g,GIAC_CONTEXT)3318 gen _trigsimplify(const gen & g,GIAC_CONTEXT) { 3319 if (g.type==_STRNG && g.subtype==-1) return g; 3320 vecteur can(1,_simplify(g,contextptr)); 3321 can.push_back(_texpand(can.back(),contextptr)); 3322 can.push_back(_tcollect(can.back(),contextptr)); 3323 for (int i=1;i<3;++i) { 3324 can.push_back(_trigtan(can[i],contextptr)); 3325 can.push_back(_trigsin(can[i],contextptr)); 3326 can.push_back(_trigcos(can[i],contextptr)); 3327 can.push_back(_tlin(can[i],contextptr)); 3328 } 3329 int n=can.size(); 3330 for (int i=0;i<n;++i) { 3331 can.push_back(_tcollect(can[i],contextptr)); 3332 } 3333 n=can.size(); 3334 for (int i=0;i<n;++i) { 3335 can.push_back(_trigtan(can[i],contextptr)); 3336 } 3337 gen simplest=g; 3338 int len=taille(g,0); 3339 for (const_iterateur it=can.begin();it!=can.end();++it) { 3340 int c=taille(*it,len); 3341 if (c<len) { 3342 simplest=*it; 3343 len=c; 3344 } 3345 } 3346 return simplest; 3347 } 3348 static const char _trigsimplify_s []="trigsimplify"; 3349 static define_unary_function_eval (__trigsimplify,&_trigsimplify,_trigsimplify_s); 3350 define_unary_function_ptr5(at_trigsimplify,alias_at_trigsimplify,&__trigsimplify,0,true) 3351 #endif 3352 3353 #ifndef NO_NAMESPACE_GIAC 3354 } // namespace giac 3355 #endif // ndef NO_NAMESPACE_GIAC 3356