1 // -*- mode:C++ ; compile-command: "g++-3.4 -I.. -g -c series.cc -DIN_GIAC -DHAVE_CONFIG_H " -*- 2 #include "giacPCH.h" 3 /* 4 * Copyright (C) 2000,14 B. Parisse, Institut Fourier, 38402 St Martin d'Heres 5 * 6 * This program is free software; you can redistribute it and/or modify 7 * it under the terms of the GNU General Public License as published by 8 * the Free Software Foundation; either version 3 of the License, or 9 * (at your option) any later version. 10 * 11 * This program is distributed in the hope that it will be useful, 12 * but WITHOUT ANY WARRANTY; without even the implied warranty of 13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 14 * GNU General Public License for more details. 15 * 16 * You should have received a copy of the GNU General Public License 17 * along with this program. If not, see <http://www.gnu.org/licenses/>. 18 */ 19 using namespace std; 20 #include <stdexcept> 21 #include <cmath> 22 #include "derive.h" 23 #include "subst.h" 24 #include "series.h" 25 #include "symbolic.h" 26 #include "unary.h" 27 #include "usual.h" 28 #include "poly.h" 29 #include "sym2poly.h" 30 #include "tex.h" 31 #include "prog.h" 32 #include "misc.h" 33 #include "intg.h" 34 #include "maple.h" 35 #include "lin.h" 36 #include "plot.h" 37 #include "giacintl.h" 38 39 #ifndef NO_NAMESPACE_GIAC 40 namespace giac { 41 #endif // ndef NO_NAMESPACE_GIAC 42 43 static int mrv_begin_order=2; 44 taylor_(const gen & f_x,const gen & x,const gen & lim_point,int ordre,vecteur & v,GIAC_CONTEXT)45 static bool taylor_(const gen & f_x,const gen & x,const gen & lim_point,int ordre,vecteur & v,GIAC_CONTEXT){ 46 gen current_derf(f_x),value,factorielle(1); 47 for (int i=0;;++i){ 48 value=subst(current_derf,x,lim_point,false,contextptr); 49 if (is_undef(value)){ 50 // if (x.type==_IDNT) value=limit(current_derf,*x._IDNTptr,lim_point,0,contextptr); 51 if (is_undef(value)) 52 return false; 53 } 54 v.push_back(ratnormal(rdiv(value,factorielle,contextptr),contextptr)); 55 if (i==ordre){ 56 v.push_back(undef); 57 return true; 58 } 59 factorielle = factorielle * gen(i+1); 60 current_derf=ratnormal(derive(current_derf,x,contextptr),contextptr); 61 if (is_undef(current_derf)) 62 return false; 63 } 64 v.dbgprint(); 65 return false; 66 } 67 taylor(const gen & f_x,const gen & x,const gen & lim_point,int ordre,vecteur & v,GIAC_CONTEXT)68 bool taylor(const gen & f_x,const gen & x,const gen & lim_point,int ordre,vecteur & v,GIAC_CONTEXT){ 69 int i=series_flags(contextptr); 70 series_flags(contextptr)=series_flags(contextptr) | (1<<7) ; 71 bool b=taylor_(f_x,x,lim_point,ordre,v,contextptr); 72 series_flags(i,contextptr); 73 return b; 74 } 75 76 // direction is always ignored for taylor, but might not 77 // for generic series_expansion 78 // shift coeff =0 for taylor taylor(const gen & lim_point,int ordre,const unary_function_ptr & f,int direction,gen & shift_coeff,GIAC_CONTEXT)79 gen taylor(const gen & lim_point,int ordre,const unary_function_ptr & f,int direction,gen & shift_coeff,GIAC_CONTEXT){ 80 // Special handling for sin/cos expansion inside limit 81 if ( is_inf(lim_point) && ( (f==at_cos) || (f==at_sin) ) ){ 82 gen g=bounded_function(contextptr); 83 /* 84 int i=sincosinf.size(); 85 sincosinf.push_back(gen(" sincosinf"+print_INT_(i))); 86 gen g=sincosinf.back(); 87 if (!g._IDNTptr->value){ 88 vecteur minusone_one(2); 89 minusone_one[0]=minus_one; 90 minusone_one[1]=plus_one; 91 gen v(vecteur(1,gen(minusone_one,_LINE__VECT))); 92 gen d(_DOUBLE_); 93 d.subtype=_INT_TYPE; 94 gen aa(makevecteur(d,v,vecteur(0)),_ASSUME__VECT); 95 g._IDNTptr->value=new gen(aa); 96 } 97 */ 98 return vecteur(1,g); 99 } 100 // if preprocessing is needed for f, series_expansion for ordre==-1 should 101 // push back in a global vector f and it's substitution 102 if (ordre<0) 103 return 0; 104 shift_coeff=0; 105 if (is_undef(lim_point) || is_inf(lim_point)){ 106 invalidserieserr(gettext("non tractable function ")+(f.ptr()->print(contextptr)+(" at "+lim_point.print(contextptr)))); 107 return undef; 108 } 109 identificateur x(" "); 110 vecteur v; 111 gen fx=f(x,contextptr); 112 if (taylor(fx,x,lim_point,ordre,v,contextptr)) 113 return v; 114 else 115 return undef; 116 } 117 porder(const sparse_poly1 & a)118 gen porder(const sparse_poly1 & a){ 119 if (a.empty()) 120 return plus_inf; 121 sparse_poly1::const_iterator a_end=a.end()-1; 122 if (is_undef(a_end->coeff)) 123 return a_end->exponent; 124 else 125 return plus_inf; 126 } 127 sparse_poly12vecteur(const sparse_poly1 & p,vecteur & v,int & shift)128 bool sparse_poly12vecteur(const sparse_poly1 & p,vecteur & v,int & shift){ 129 sparse_poly1::const_iterator it=p.begin(),itend=p.end(); 130 v.clear(); 131 if (p.empty()) 132 return true; 133 if (p.back().exponent.type!=_INT_) 134 return false; 135 int n1=p.front().exponent.val,n2=p.back().exponent.val; 136 if (n1>n2 || (n2-n1)+1<0) // if n==RAND_MAX, n+1<0 137 return false; 138 if (n1<0) 139 shift=n1; 140 else 141 shift=n1=0; 142 v.resize(n2-n1+1); 143 for (;it!=itend;++it){ 144 if (it->exponent.type!=_INT_) 145 return false; 146 int m=it->exponent.val; 147 if (m<n1 || m>n2) 148 return false; 149 v[m-n1]=it->coeff; 150 } 151 reverse(v.begin(),v.end()); 152 return true; 153 } 154 vecteur2sparse_poly1(const vecteur & v,sparse_poly1 & p)155 void vecteur2sparse_poly1(const vecteur & v,sparse_poly1 & p){ 156 p.clear(); 157 vecteur::const_iterator it=v.begin(),itend=v.end(); 158 p.reserve(itend-it); 159 for (int i=0;it!=itend;++i,++it){ 160 if (!is_zero(*it)) 161 p.push_back(monome(*it,i)); 162 } 163 } 164 gen2spol1(const gen & g)165 sparse_poly1 gen2spol1(const gen &g){ 166 if (g.type!=_VECT) 167 return sparse_poly1(1,monome(g,0)); 168 sparse_poly1 p; 169 vecteur2sparse_poly1(*g._VECTptr,p); 170 return p; 171 } 172 vecteur2sparse_poly1(const vecteur & v)173 sparse_poly1 vecteur2sparse_poly1(const vecteur & v){ 174 sparse_poly1 p; 175 vecteur2sparse_poly1(v,p); 176 return p; 177 } 178 spol12gen(const sparse_poly1 & p,GIAC_CONTEXT)179 gen spol12gen(const sparse_poly1 & p,GIAC_CONTEXT){ 180 string t; 181 t = t+series_variable_name(contextptr); 182 identificateur tt(t); 183 gen T(tt),remains; 184 gen g=sparse_poly12gen(p,T,remains,false); 185 if (!is_zero(remains)) 186 g += remains*order_size(T,contextptr); 187 return g; 188 } 189 spol12gen(const gen & coeff,const gen & x)190 static gen spol12gen(const gen & coeff,const gen & x){ 191 if (coeff.type==_VECT){ 192 vecteur v=*coeff._VECTptr; 193 int s=int(v.size()); 194 for (int i=0;i<s;++i){ 195 v[i]=spol12gen(v[i],x); 196 } 197 return gen(v,coeff.subtype); 198 } 199 if (coeff.type==_SPOL1){ 200 gen remains=0; 201 return sparse_poly12gen(*coeff._SPOL1ptr,x,remains,true)+remains; 202 } 203 if (coeff.type!=_SYMB) 204 return coeff; 205 return symbolic(coeff._SYMBptr->sommet,spol12gen(coeff._SYMBptr->feuille,x)); 206 } 207 sparse_poly12gen(const sparse_poly1 & p,const gen & x,gen & remains,bool with_order_size)208 gen sparse_poly12gen(const sparse_poly1 & p,const gen & x,gen & remains,bool with_order_size){ 209 gen res; 210 remains=0; 211 sparse_poly1::const_iterator it=p.begin(),itend=p.end(); 212 for (;it!=itend;++it){ 213 gen coeff=it->coeff; 214 if (is_undef(coeff)){ 215 remains=pow(x,it->exponent,context0); // ok 216 if (with_order_size) 217 return res+remains*order_size(x,context0); 218 else 219 return res; 220 } 221 coeff=spol12gen(coeff,x); 222 res = res + coeff * pow(x,it->exponent,context0); // ok 223 } 224 return res; 225 } 226 227 ptruncate(sparse_poly1 & p,const gen & ordre,GIAC_CONTEXT)228 bool ptruncate(sparse_poly1 & p,const gen & ordre,GIAC_CONTEXT){ 229 if ( (series_flags(contextptr) & 0x2) || p.empty() ){ 230 sparse_poly1::iterator it=p.begin(),itend=p.end(); 231 gen first=it->exponent; 232 for (;it!=itend;++it){ 233 if (is_undef(it->coeff)) 234 return true; 235 if (ck_is_strictly_greater(it->exponent-first,ordre,contextptr)){ 236 it->coeff=undef; 237 p.erase(it+1,itend); 238 return true; 239 } 240 } 241 } 242 return true; 243 } 244 poly_truncate(sparse_poly1 & p,int ordre,GIAC_CONTEXT)245 void poly_truncate(sparse_poly1 & p,int ordre,GIAC_CONTEXT){ 246 if ( (series_flags(contextptr) & 0x2) || p.empty() ){ 247 sparse_poly1::iterator it=p.begin(),itend=p.end(); 248 for (;it!=itend;++it){ 249 if (is_undef(it->coeff)) 250 return ; 251 if (ck_is_strictly_greater(it->exponent,ordre,contextptr)){ 252 it->coeff=undef; 253 p.erase(it+1,itend); 254 return ; 255 } 256 } 257 } 258 return ; 259 } 260 261 static gen remove_lnexp(const gen & e,GIAC_CONTEXT); 262 padd(const sparse_poly1 & a,const sparse_poly1 & b,sparse_poly1 & res,GIAC_CONTEXT)263 bool padd(const sparse_poly1 & a,const sparse_poly1 &b, sparse_poly1 & res,GIAC_CONTEXT){ 264 // Series addition 265 if (a.empty()){ 266 if (&b!=&res) 267 res=b; 268 return true; 269 } 270 if (b.empty()){ 271 if (&a!=&res) 272 res=a; 273 return true; 274 } 275 sparse_poly1::const_iterator a_cur,a_end,b_cur,b_end; 276 sparse_poly1 temp_a,temp_b; 277 if (&res==&a){ // must make a copy of a 278 temp_a=a; 279 a_cur=temp_a.begin(); 280 a_end=temp_a.end(); 281 } 282 else { 283 a_cur=a.begin(); 284 a_end=a.end(); 285 } 286 if (&res==&b){ // must make a copy of b 287 temp_b=b; 288 b_cur=temp_b.begin(); 289 b_end=temp_b.end(); 290 } 291 else { 292 b_cur=b.begin(); 293 b_end=b.end(); 294 } 295 res.clear(); 296 res.reserve((a_end-a_cur)+(b_end-b_cur)); 297 for (;(a_cur!=a_end) && (b_cur!=b_end) ;) { 298 gen a_pow=a_cur->exponent; 299 gen b_pow=b_cur->exponent; 300 // a and b are non-empty, compare powers 301 if (ck_is_strictly_greater(b_pow,a_pow,contextptr)) { 302 // get coefficient from a 303 res.push_back(*a_cur); 304 if (is_undef(a_cur->coeff)){ 305 return true; 306 } 307 ++a_cur; 308 continue; 309 } 310 if (ck_is_strictly_greater(a_pow,b_pow,contextptr)) { 311 // get coefficient from b 312 res.push_back(*b_cur); 313 if (is_undef(b_cur->coeff)){ 314 return true; 315 } 316 ++b_cur; 317 continue; 318 } 319 // Add coefficient of a and b 320 gen sum=a_cur->coeff+b_cur->coeff; 321 if (sum.type>_POLY && sum.type!=_FRAC &&(res.empty() || (series_flags(contextptr) & 0x1) ) ){ 322 //cerr << sum << " "; 323 sum=recursive_normal(remove_lnexp(sum,contextptr),contextptr); 324 //cerr << sum << '\n'; 325 } 326 // gen sum=(a_cur->coeff+b_cur->coeff); 327 if (!is_zero(sum)) 328 res.push_back(monome(sum,a_pow)); 329 if (is_undef(sum)){ 330 return true; 331 } 332 ++a_cur; 333 ++b_cur; 334 } 335 for (;a_cur!=a_end;++a_cur) 336 res.push_back(*a_cur); 337 for (;b_cur!=b_end;++b_cur) 338 res.push_back(*b_cur); 339 return true; 340 } 341 spadd(const sparse_poly1 & a,const sparse_poly1 & b,GIAC_CONTEXT)342 sparse_poly1 spadd(const sparse_poly1 & a,const sparse_poly1 &b,GIAC_CONTEXT){ 343 sparse_poly1 res; 344 padd(a,b,res,contextptr); 345 return res; 346 } 347 spsub(const sparse_poly1 & a,const sparse_poly1 & b,GIAC_CONTEXT)348 sparse_poly1 spsub(const sparse_poly1 & a,const sparse_poly1 &b,GIAC_CONTEXT){ 349 sparse_poly1 res(b); 350 pneg(b,res,contextptr); 351 padd(a,res,res,contextptr); 352 return res; 353 } 354 355 pmul(const sparse_poly1 & a,const gen & b_orig,sparse_poly1 & res,GIAC_CONTEXT)356 bool pmul(const sparse_poly1 & a,const gen & b_orig, sparse_poly1 & res,GIAC_CONTEXT){ 357 gen b(b_orig); 358 if (&a==&res){ 359 if (is_one(b_orig)) return true; 360 sparse_poly1::iterator it=res.begin(),itend=res.end(); 361 for (;it!=itend;++it){ 362 gen g=it->coeff * b; 363 if (g.type>_POLY && g.type!=_FRAC) 364 g=ratnormal(g,contextptr) ; 365 it->coeff = g; 366 } 367 return true; 368 } 369 sparse_poly1::const_iterator it=a.begin(),itend=a.end(); 370 res.clear(); 371 res.reserve(itend-it); 372 for (;it!=itend;++it) 373 res.push_back(monome(ratnormal(it->coeff * b,contextptr), it->exponent)); 374 return true; 375 } 376 pmul(const gen & b,const sparse_poly1 & a,sparse_poly1 & res,GIAC_CONTEXT)377 bool pmul(const gen & b, const sparse_poly1 & a,sparse_poly1 & res,GIAC_CONTEXT){ 378 return pmul(a,b,res,contextptr); 379 } 380 spmul(const sparse_poly1 & a,const gen & b,GIAC_CONTEXT)381 sparse_poly1 spmul(const sparse_poly1 & a,const gen &b,GIAC_CONTEXT){ 382 sparse_poly1 res; 383 if (!pmul(a,b,res,contextptr)) 384 res=sparse_poly1(1,monome(1,undef)); 385 return res; 386 } 387 spmul(const gen & a,const sparse_poly1 & b,GIAC_CONTEXT)388 sparse_poly1 spmul(const gen & a,const sparse_poly1 &b,GIAC_CONTEXT){ 389 sparse_poly1 res; 390 if (!pmul(a,b,res,contextptr)) 391 res=sparse_poly1(1,monome(1,undef)); 392 return res; 393 } 394 395 struct monome_less { monome_lessgiac::monome_less396 monome_less() {} operator ()giac::monome_less397 bool operator () (const monome & a,const monome & b){ 398 return ck_is_strictly_greater(b.exponent,a.exponent,context0); 399 } 400 }; 401 402 struct symb_size_less_t { symb_size_less_tgiac::symb_size_less_t403 symb_size_less_t() {} operator ()giac::symb_size_less_t404 bool operator () (const gen &a,const gen &b){ 405 return symb_size_less(a,b); 406 } 407 }; 408 pmul(const sparse_poly1 & celuici,const sparse_poly1 & other,sparse_poly1 & final_seq,bool n_truncate,const gen & n_valuation,GIAC_CONTEXT)409 bool pmul(const sparse_poly1 & celuici,const sparse_poly1 &other, sparse_poly1 & final_seq,bool n_truncate,const gen & n_valuation,GIAC_CONTEXT){ 410 #ifdef TIMEOUT 411 control_c(); 412 #endif 413 if (ctrl_c || interrupted) { 414 interrupted=ctrl_c=true; 415 return false; 416 } 417 int asize=int(celuici.size()); 418 int bsize=int(other.size()); 419 if ( (!asize) || (!bsize) ) { 420 final_seq.clear(); 421 return true; 422 } 423 if (asize==1){ 424 gen temp(celuici.front().coeff); 425 pshift(other,celuici.front().exponent,final_seq,contextptr); 426 // COUT << other << "Shifted" << final_seq << '\n'; 427 return pmul(final_seq,temp,final_seq,contextptr); 428 // COUT << other << "Multiplied" << final_seq << '\n'; 429 } 430 if (bsize==1){ 431 gen temp(other.front().coeff); 432 pshift(celuici,other.front().exponent,final_seq,contextptr); 433 return pmul(final_seq,temp,final_seq,contextptr); 434 } 435 sparse_poly1 new_seq; 436 new_seq.reserve(asize*bsize); 437 // General sparse series multiplication: complexity is N*M*ln(N*M) 438 // Storage capacity 2*N*M expair 439 // That's much more than O(N+M) for dense poly *but* 440 // it works for non integer powers 441 // COUT << celuici << "pmul" << other << '\n'; 442 // First find the order product 443 gen a_max = porder(celuici); 444 gen b_max = porder(other); 445 gen a_min = celuici.front().exponent; 446 gen b_min = other.front().exponent; 447 gen c_min = normal(a_min + b_min,contextptr); 448 gen c_max = min(normal(a_min + b_max,contextptr),normal(b_min + a_max,contextptr),contextptr); 449 if (c_max.type==_SYMB && c_max._SYMBptr->sommet==at_max) 450 return false; // setsizeerr(gettext("series.cc/pmul")); 451 // compute all products term by term, with optimization for dense poly 452 // (coeff are sorted for dense poly) 453 sparse_poly1::const_iterator itb = other.begin(),itbend = other.end(); 454 sparse_poly1::const_iterator ita = celuici.begin(),ita_end=celuici.end(); 455 sparse_poly1::const_iterator itabegin = ita; 456 gen old_pow=normal(ita->exponent+itb->exponent,contextptr); 457 gen res(0); 458 for ( ; ita!=ita_end; ++ita ){ 459 sparse_poly1::const_iterator itacur=ita; 460 sparse_poly1::const_iterator itbcur=itb; 461 for (;;) { 462 #ifdef TIMEOUT 463 control_c(); 464 #endif 465 if (ctrl_c || interrupted) { 466 interrupted=ctrl_c=true; 467 return false; 468 } 469 gen cur_pow=normal(itacur->exponent+itbcur->exponent,contextptr); 470 if ((n_truncate && ck_is_strictly_greater(n_valuation,cur_pow,contextptr)) || ck_is_greater(c_max,cur_pow,contextptr)){ 471 if (cur_pow!=old_pow){ 472 new_seq.push_back( monome(res,old_pow )); 473 res=itacur->coeff * itbcur->coeff; 474 old_pow=cur_pow; 475 } 476 else 477 res=res+ itacur->coeff * itbcur->coeff; 478 } 479 if (itacur==itabegin) 480 break; 481 --itacur; 482 ++itbcur; 483 if (itbcur==itbend) 484 break; 485 } 486 } 487 --ita; 488 ++itb; 489 for ( ; itb!=itbend;++itb){ 490 sparse_poly1::const_iterator itacur=ita; 491 sparse_poly1::const_iterator itbcur=itb; 492 for (;;) { 493 #ifdef TIMEOUT 494 control_c(); 495 #endif 496 if (ctrl_c || interrupted) { 497 interrupted=ctrl_c=true; 498 return false; 499 } 500 gen cur_pow=normal(itacur->exponent + itbcur->exponent,contextptr); 501 if ((n_truncate && ck_is_strictly_greater(n_valuation,cur_pow,contextptr)) || ck_is_greater(c_max,cur_pow,contextptr)){ 502 if (cur_pow!=old_pow){ 503 new_seq.push_back( monome(res ,old_pow )); 504 res= itacur->coeff * itbcur->coeff ; 505 old_pow=cur_pow; 506 } 507 else 508 res=res+ itacur->coeff * itbcur->coeff ; 509 } 510 if (itacur==itabegin) 511 break; 512 --itacur; 513 ++itbcur; 514 if (itbcur==itbend) 515 break; 516 } 517 } 518 new_seq.push_back( monome(res ,old_pow )); 519 // COUT << new_seq << '\n'; 520 // sort by asc. power 521 sort( new_seq.begin(),new_seq.end(),monome_less()); 522 // COUT << "Sorted" << new_seq << '\n'; 523 // add terms with same power 524 sparse_poly1::const_iterator it=new_seq.begin(); 525 sparse_poly1::const_iterator itend=new_seq.end(); 526 final_seq.clear(); 527 final_seq.reserve(itend-it); 528 while (it!=itend){ 529 gen res=it->coeff; 530 gen pow=it->exponent; 531 if (is_undef(res)){ 532 final_seq.push_back(*it); 533 return true; 534 } 535 ++it; 536 while ( (it!=itend) && (it->exponent==pow)){ 537 #ifdef TIMEOUT 538 control_c(); 539 #endif 540 if (ctrl_c || interrupted) { 541 interrupted=ctrl_c=true; 542 return false; 543 } 544 if (is_undef(it->coeff)){ 545 final_seq.push_back(*it); 546 return true; 547 } 548 res=res+it->coeff; 549 ++it; 550 } 551 if (series_flags(contextptr) & 0x1) 552 res=recursive_normal(res,contextptr); 553 if (!is_zero(res)) 554 final_seq.push_back(monome(res, pow)); 555 } 556 if (c_max!=plus_inf) 557 final_seq.push_back(monome(undef, c_max)); 558 return true; 559 //COUT << final_seq.back().coeff << '\n'; 560 //return true; 561 } 562 spmul(const sparse_poly1 & a,const sparse_poly1 & b,GIAC_CONTEXT)563 sparse_poly1 spmul(const sparse_poly1 & a,const sparse_poly1 &b,GIAC_CONTEXT){ 564 sparse_poly1 res; 565 if (!pmul(a,b,res,false,0,contextptr)) 566 res=sparse_poly1(1,monome(1,undef)); 567 return res; 568 } 569 pneg(const sparse_poly1 & a,sparse_poly1 & res,GIAC_CONTEXT)570 bool pneg(const sparse_poly1 & a,sparse_poly1 & res,GIAC_CONTEXT){ 571 if (&a==&res){ 572 sparse_poly1::iterator it=res.begin(),itend=res.end(); 573 for (;it!=itend;++it) 574 it->coeff=-it->coeff; 575 return true; 576 } 577 sparse_poly1::const_iterator it=a.begin(),itend=a.end(); 578 res.clear(); 579 res.reserve(itend-it); 580 for (;it!=itend;++it) 581 res.push_back(monome(-it->coeff, it->exponent)); 582 return true; 583 } 584 spneg(const sparse_poly1 & a,GIAC_CONTEXT)585 sparse_poly1 spneg(const sparse_poly1 & a,GIAC_CONTEXT){ 586 sparse_poly1 res; 587 pneg(a,res,contextptr); 588 return res; 589 } 590 pshift(const sparse_poly1 & a,const gen & b_orig,sparse_poly1 & res,GIAC_CONTEXT)591 bool pshift(const sparse_poly1 & a,const gen & b_orig, sparse_poly1 & res,GIAC_CONTEXT){ 592 if (is_zero(b_orig)){ 593 if (&a!=&res) 594 res=a; 595 return true; 596 } 597 gen b(b_orig); 598 if (&a==&res){ 599 sparse_poly1::iterator it=res.begin(),itend=res.end(); 600 for (;it!=itend;++it) 601 it->exponent = normal(it->exponent + b,contextptr); 602 return true; 603 } 604 sparse_poly1::const_iterator it=a.begin(),itend=a.end(); 605 res.clear(); 606 res.reserve(itend-it); 607 for (;it!=itend;++it) 608 res.push_back(monome(it->coeff , normal(it->exponent +b,contextptr))); 609 return true; 610 } 611 612 // ascending order division pdiv(const sparse_poly1 & a,const sparse_poly1 & b_orig,sparse_poly1 & res,int ordre_orig,GIAC_CONTEXT)613 bool pdiv(const sparse_poly1 & a,const sparse_poly1 &b_orig, sparse_poly1 & res,int ordre_orig,GIAC_CONTEXT){ 614 #ifdef TIMEOUT 615 control_c(); 616 #endif 617 if (ctrl_c || interrupted) { 618 interrupted=ctrl_c=true; 619 return false; 620 } 621 //if (debug_infolevel) CERR << CLOCK()*1e-6 << " pdiv begin" <<'\n'; 622 sparse_poly1 b(b_orig); 623 ptruncate(b,ordre_orig,contextptr); 624 if (b.empty()){ 625 // divisionby0err(a); 626 return false; 627 } 628 gen b0=b.front().coeff; 629 if (is_undef(b0)){ 630 if (&b!=&res) 631 res=b; 632 return true; 633 } 634 if (b.size()==1){ 635 pshift(a,-b.front().exponent,res,contextptr); 636 return pdiv(res,b0,res,contextptr); 637 } 638 // COUT << a << "/" << b << '\n'; 639 if (&res==&b){ 640 // setsizeerr(gettext("series.cc/pdiv")); 641 return false; 642 } 643 gen e0=b.front().exponent; 644 gen ordre=min(min(porder(a),porder(b)-e0,contextptr),ordre_orig,contextptr); 645 if (ordre==plus_inf) 646 ordre=series_default_order(contextptr); 647 // COUT << ordre << '\n'; 648 if (ordre.type==_SYMB && ordre._SYMBptr->sommet==at_max) 649 return false; // setsizeerr(gettext("series.cc/pdiv")); 650 sparse_poly1 rem(a); 651 res.clear(); 652 sparse_poly1 bshift; 653 gen q_cur,e_cur; // current quotient, current exponent 654 for (;;){ 655 if (is_undef(rem.front().coeff)){ 656 res.push_back(monome(undef,rem.front().exponent-e0)); 657 // COUT << "=" << res << '\n'; 658 return true; 659 } 660 q_cur=rdiv(rem.front().coeff,b0,contextptr); 661 e_cur=rem.front().exponent-e0; 662 res.push_back(monome(q_cur,e_cur)); 663 pshift(b,e_cur,bshift,contextptr); 664 sparse_poly1::iterator it=bshift.begin(),itend=bshift.end(); 665 for (;it!=itend;++it){ 666 if (is_undef(it->coeff)) 667 break; 668 if (ck_is_strictly_greater(it->exponent,ordre,contextptr)){ 669 it->coeff=undef; 670 bshift.erase(it+1,itend); 671 break; 672 } 673 } 674 if (!pmul(-q_cur,bshift,bshift,contextptr)) 675 return false; 676 padd(rem,bshift,rem,contextptr); 677 // COUT << rem.front().exponent << " " << e0+ordre << '\n'; 678 if (ck_is_strictly_greater(rem.front().exponent,a.front().exponent+ordre,contextptr)){ 679 res.push_back(monome(undef,a.front().exponent+ordre+1-e0)); 680 return true; 681 } 682 } 683 return true; 684 } 685 spdiv(const sparse_poly1 & a,const sparse_poly1 & b,GIAC_CONTEXT)686 sparse_poly1 spdiv(const sparse_poly1 & a,const sparse_poly1 &b,GIAC_CONTEXT){ 687 sparse_poly1 res; 688 gen og=min(porder(a),porder(b),contextptr); 689 int o=series_default_order(contextptr); 690 if (og.type==_INT_) 691 o=og.val; 692 if (!pdiv(a,b,res,o,contextptr)) 693 res=sparse_poly1(1,monome(1,undef)); 694 return res; 695 } 696 pdiv(const sparse_poly1 & a,const gen & b_orig,sparse_poly1 & res,GIAC_CONTEXT)697 bool pdiv(const sparse_poly1 & a,const gen & b_orig, sparse_poly1 & res,GIAC_CONTEXT){ 698 #ifdef TIMEOUT 699 control_c(); 700 #endif 701 if (ctrl_c || interrupted) { 702 interrupted=ctrl_c=true; 703 return false; 704 } 705 if (is_zero(b_orig)) 706 return false; // divisionby0err(a); 707 if (is_one(b_orig)){ 708 if (&a!=&res) 709 res=a; 710 return true; 711 } 712 gen b(b_orig); 713 if (&a==&res){ 714 sparse_poly1::iterator it=res.begin(),itend=res.end(); 715 for (;it!=itend;++it){ 716 it->coeff=rdiv(it->coeff, b,contextptr); 717 if (series_flags(contextptr) & 0x1) 718 it->coeff=normal(it->coeff,contextptr); 719 } 720 // it->coeff=rdiv(it->coeff, b,contextptr); 721 return true; 722 } 723 sparse_poly1::const_iterator it=a.begin(),itend=a.end(); 724 res.clear(); 725 res.reserve(itend-it); 726 gen tmp; 727 for (;it!=itend;++it){ 728 tmp=rdiv(it->coeff,b,contextptr); 729 if (series_flags(contextptr) & 0x1) 730 tmp=normal(tmp,contextptr); 731 res.push_back(monome(tmp , it->exponent)); 732 } 733 // res.push_back(monome(rdiv(it->coeff,b,contextptr) , it->exponent)); 734 return true; 735 } 736 spdiv(const sparse_poly1 & a,const gen & b,GIAC_CONTEXT)737 sparse_poly1 spdiv(const sparse_poly1 & a,const gen &b,GIAC_CONTEXT){ 738 sparse_poly1 res; 739 if (!pdiv(a,b,res,contextptr)) 740 res=sparse_poly1(1,undef); 741 return res; 742 } 743 744 // v is replaced by e*v where e*v has no denominator lcmdeno(vecteur & v,gen & e,GIAC_CONTEXT)745 void lcmdeno(vecteur &v,gen & e,GIAC_CONTEXT){ 746 if (v.empty()){ 747 e=1; 748 return; 749 } 750 if (is_undef(v.front())){ 751 v.erase(v.begin()); 752 lcmdeno(v,e,contextptr); 753 v.insert(v.begin(),undef); 754 return; 755 } 756 vecteur l; 757 lvar(v,l); 758 //int l_size(l.size()); 759 vecteur w; 760 w.reserve(2*v.size()); 761 gen common=1,f,num,den; 762 // compute lcm of denominators in common 763 vecteur::iterator it=v.begin(),itend=v.end(); 764 for (;it!=itend;++it){ 765 if (is_integer(*it)){ 766 num=*it; den=1; 767 } 768 else { 769 if (it->type==_FRAC && is_integer(it->_FRACptr->num) && is_integer(it->_FRACptr->den)) 770 f=*it; 771 else 772 f=e2r(*it,l,contextptr); 773 fxnd(f,num,den); 774 } 775 w.push_back(num); 776 w.push_back(den); 777 // replace common by lcm of common and den 778 #ifndef USE_GMP_REPLACEMENTS 779 if (common.type==_ZINT && common.ref_count()==1 && is_integer(den)){ 780 if (den.type==_ZINT) 781 mpz_lcm(*common._ZINTptr,*common._ZINTptr,*den._ZINTptr); 782 else 783 mpz_lcm_ui(*common._ZINTptr,*common._ZINTptr,absint(den.val)); 784 } 785 else 786 common = lcm(common,den); 787 #else 788 common = lcm(common,den); 789 #endif 790 } 791 // compute e and recompute v 792 e=r2sym(common,l,contextptr); 793 it=v.begin(); 794 for (int i=0;it!=itend;++it,i=i+2){ 795 if (it->type==_FRAC && is_integer(it->_FRACptr->num) && is_integer(it->_FRACptr->den) && is_integer(common)) 796 *it=w[i]*rdiv(common,w[i+1],contextptr); 797 else 798 *it=r2sym(w[i]*rdiv(common,w[i+1],contextptr),l,contextptr); 799 } 800 } 801 lcmdeno_converted(vecteur & v,gen & e,GIAC_CONTEXT)802 void lcmdeno_converted(vecteur &v,gen & e,GIAC_CONTEXT){ 803 if (v.empty()){ 804 e=1; 805 return; 806 } 807 if (is_undef(v.front())){ 808 v.erase(v.begin()); 809 lcmdeno_converted(v,e,contextptr); 810 v.insert(v.begin(),undef); 811 return; 812 } 813 vecteur w; 814 w.reserve(2*v.size()); 815 gen common=1,f,num,den; 816 // compute lcm of denominators in common 817 vecteur::iterator it=v.begin(),itend=v.end(); 818 for (;it!=itend;++it){ 819 fxnd(*it,num,den); 820 w.push_back(num); 821 w.push_back(den); 822 // replace common by lcm of common and den 823 common = lcm(common,den); 824 } 825 // compute e and recompute v 826 e=common; 827 it=v.begin(); 828 for (int i=0;it!=itend;++it,i=i+2) 829 *it=w[i]*rdiv(common,w[i+1],contextptr); 830 } 831 lcmdeno(sparse_poly1 & v,gen & e,GIAC_CONTEXT)832 void lcmdeno(sparse_poly1 &v,gen & e,GIAC_CONTEXT){ 833 if (v.empty()){ 834 e=1; 835 return; 836 } 837 if (is_undef(v.back().coeff)){ 838 monome last=v.back(); 839 v.pop_back(); 840 lcmdeno(v,e,contextptr); 841 v.push_back(last); 842 return; 843 } 844 vecteur l; 845 lvar(v,l); 846 int l_size(int(l.size())); 847 vector<gen> w; 848 w.reserve(2*l_size); 849 gen common=1,num,den,f; 850 // compute lcm of denominators in common 851 sparse_poly1::iterator it=v.begin(),itend=v.end(); 852 for (;it!=itend;++it){ 853 f=e2r(it->coeff,l,contextptr); 854 fxnd(f,num,den); 855 w.push_back(num); 856 w.push_back(den); 857 common=lcm(common,den); 858 } 859 // compute e and recompute v 860 e=r2sym(common,l,contextptr); 861 it=v.begin(); 862 for (int i=0;it!=itend;++it,i=i+2){ 863 it->coeff=r2sym(w[i]*rdiv(common,w[i+1],contextptr),l,contextptr); 864 } 865 } 866 vreverse(iterateur a,iterateur aend)867 void vreverse(iterateur a,iterateur aend){ 868 iterateur b=aend-1; 869 for (;a<b;++a,--b){ 870 swapgen(*a,*b); 871 } 872 } 873 pcompose(const vecteur & v,const sparse_poly1 & p,sparse_poly1 & res,GIAC_CONTEXT)874 bool pcompose(const vecteur & v,const sparse_poly1 & p, sparse_poly1 & res,GIAC_CONTEXT){ 875 #ifdef TIMEOUT 876 control_c(); 877 #endif 878 if (ctrl_c || interrupted) { 879 interrupted=ctrl_c=true; 880 return false; 881 } 882 if (v.empty()){ 883 res.clear(); 884 return true; 885 } 886 if ( p.empty() ){ 887 res.clear(); 888 if (!is_zero(v.front())) 889 res.push_back(monome(v.front(),0)); 890 return true; 891 } 892 // Conversion of p and v to "internal" polynomial form 893 vecteur l; // will contain the list of variables common to v and p 894 alg_lvar(v,l); 895 alg_lvar(p,l); 896 // int l_size(l.size()); 897 gen plcm=plus_one,vlcm=plus_one,f,num,den; 898 // compute lcm of denominators of p in plcm 899 sparse_poly1::const_iterator its=p.begin(),itsend=p.end(); 900 vecteur ptemp; 901 ptemp.reserve(2*(itsend-its)); 902 for (;its!=itsend;++its){ 903 f=e2r(its->coeff,l,contextptr); 904 fxnd(f,num,den); 905 ptemp.push_back(num); 906 ptemp.push_back(den); 907 plcm=lcm(den,plcm); 908 } 909 // compute pcopy such that pcopy/plcm=p 910 its=p.begin(); 911 sparse_poly1 pcopy; 912 pcopy.reserve(itsend-its); 913 for (int i=0;its!=itsend;++its,i=i+2){ 914 num=ptemp[i]*rdiv(plcm,ptemp[i+1],contextptr); 915 pcopy.push_back(monome(num,its->exponent)); 916 } 917 // do the same thing on v 918 vecteur w; 919 // compute lcm of denominators in common 920 vecteur::const_iterator it=v.begin(),itend=v.end(); 921 w.reserve(2*(itend-it)); 922 for (;it!=itend;++it){ 923 f=e2r(*it,l,contextptr); 924 fxnd(f,num,den); 925 w.push_back(num); 926 w.push_back(den); 927 vlcm=lcm(vlcm,den); 928 } 929 // compute vcopy 930 it=v.begin(); 931 vecteur vcopy; 932 vcopy.reserve(itend-it); 933 for (int i=0;it!=itend;++it,i=i+2) 934 vcopy.push_back(w[i]*rdiv(vlcm,w[i+1],contextptr)); 935 vreverse(vcopy.begin(),vcopy.end()); 936 if (vcopy.empty() ){ 937 res=sparse_poly1(1,monome(undef,minus_inf)); 938 return true; 939 } 940 // COUT << "compose " << vcopy << " with " << pcopy << '\n'; 941 it=vcopy.begin(),itend=vcopy.end(); 942 int n=int(itend-it)-1; 943 bool n_truncate=false; 944 gen n_valuation; 945 if (is_undef(*it)){ 946 ++it; 947 n_truncate=true; 948 n_valuation=gen(n)*p.front().exponent; 949 // add undef order term 950 gen cur_ordre=porder(p); 951 // compare cur_ordre with n*valuation(pcopy) 952 if ( (cur_ordre==plus_inf) || (ck_is_strictly_greater(cur_ordre,n_valuation,contextptr)) ){ 953 // remove greater order terms from pcopy 954 for (;!pcopy.empty();){ 955 if (ck_is_strictly_greater(pcopy.back().exponent,n_valuation,contextptr)) 956 pcopy.pop_back(); 957 else 958 break; 959 } 960 // insert undef 961 if (pcopy.empty() || (!is_undef(pcopy.back().coeff)) ) 962 pcopy.push_back(monome(undef,n_valuation)); 963 } 964 } 965 // Skip 0 coeffs in the reverse list vcopy 966 for (;it!=itend;++it){ 967 if (!is_zero(*it)) 968 break; 969 } 970 if (it==itend){ 971 res=sparse_poly1(1,monome(undef)); 972 return true; 973 } 974 res=sparse_poly1(1,monome(*it)); 975 // COUT << res << '\n'; 976 ++it; 977 if (it==itend && is_undef(pcopy.back().coeff)) 978 res.push_back(monome(undef,min(n_valuation,pcopy.back().exponent,contextptr))); 979 gen plcmn=plus_one; 980 for (;it!=itend;++it){ 981 plcmn=plcmn*plcm; 982 // COUT << res << "*" << pcopy << '\n' ; 983 if (!pmul(res,pcopy,res,n_truncate,n_valuation,contextptr)) 984 return false; 985 if (n_truncate){ // Remove all terms of order > n_valuation 986 sparse_poly1::iterator sit=res.begin(),sitend=res.end(); 987 for (;sit!=sitend;++sit){ 988 if (ck_is_greater(sit->exponent,n_valuation,contextptr)){ 989 res.erase(sit,sitend); 990 res.push_back(monome(undef,n_valuation)); 991 break; 992 } 993 } 994 } 995 // COUT << res << '\n'; 996 if (!is_zero(*it)) 997 padd(res,sparse_poly1(1,monome(*it*plcmn)),res,contextptr); 998 // COUT << res << '\n'; 999 } 1000 den=vlcm*plcmn; 1001 // back conversion from res to symbolic form 1002 sparse_poly1::iterator sit=res.begin(),sitend=res.end(); 1003 for (;sit!=sitend;++sit){ 1004 num=den; 1005 sit->coeff=r2sym(fraction(sit->coeff,num).normal(),l,contextptr); 1006 } 1007 return true; 1008 } 1009 ppow(const sparse_poly1 & base,int m,int ordre,sparse_poly1 & res,GIAC_CONTEXT)1010 bool ppow(const sparse_poly1 & base,int m,int ordre,sparse_poly1 & res,GIAC_CONTEXT){ 1011 #ifdef TIMEOUT 1012 control_c(); 1013 #endif 1014 if (ctrl_c || interrupted) { 1015 interrupted=ctrl_c=true; 1016 return false; 1017 } 1018 if (m==0){ 1019 res.clear(); 1020 return true; 1021 } 1022 if (m==1){ 1023 if (&base!=&res) 1024 res=base; 1025 return true; 1026 } 1027 sparse_poly1 temp; 1028 if (!pmul(base,base,temp,true,ordre,contextptr)) 1029 return false; 1030 ptruncate(temp,ordre,contextptr); 1031 if (m%2){ 1032 if (!ppow(temp,m/2,ordre,temp,contextptr) || 1033 !pmul(temp,base,res,true,ordre,contextptr)) 1034 return false; 1035 } 1036 else { 1037 if (!ppow(temp,m/2,ordre,res,contextptr)) 1038 return false; 1039 } 1040 ptruncate(res,ordre,contextptr); 1041 return true; 1042 } 1043 1044 // constant power, otherwise use exp(ln) ppow(const sparse_poly1 & base,const gen & e,int ordre,int direction,sparse_poly1 & res,GIAC_CONTEXT)1045 bool ppow(const sparse_poly1 & base,const gen & e,int ordre,int direction,sparse_poly1 & res,GIAC_CONTEXT){ 1046 #ifdef TIMEOUT 1047 control_c(); 1048 #endif 1049 if (ctrl_c || interrupted) { 1050 interrupted=ctrl_c=true; 1051 return false; 1052 } 1053 if (base.size()==1){ 1054 gen basepow; 1055 if (e.type==_FRAC && e._FRACptr->den==2 && is_positive(-base.front().coeff,contextptr)) 1056 basepow=pow(cst_i,e._FRACptr->num,contextptr)*pow(-base.front().coeff,e,contextptr); 1057 else 1058 basepow=pow(base.front().coeff,e,contextptr); 1059 if (&base==&res){ 1060 res.front().coeff=basepow; 1061 res.front().exponent=res.front().exponent*e; 1062 } 1063 else 1064 res=sparse_poly1(1,monome(basepow,base.front().exponent*e)); 1065 return true; 1066 } 1067 gen n=porder(base); 1068 if ((n==plus_inf) && (e.type==_INT_) && (e.val>=0) ){ // exact power 1069 int m=e.val; 1070 return ppow(base,m,ordre,res,contextptr); 1071 } 1072 if (base.empty()){ 1073 if (ck_is_positive(e,contextptr)) 1074 res.clear(); 1075 else 1076 return false; // divisionby0err(base); 1077 return true; 1078 } 1079 // series expansion to a constant power 1080 monome first(base.front()); 1081 sparse_poly1 basecopy(base); 1082 basecopy.erase(basecopy.begin()); 1083 pshift(basecopy,-first.exponent,basecopy,contextptr); 1084 if (!pdiv(basecopy,first.coeff,basecopy,contextptr)) 1085 return false; 1086 if (n==plus_inf && !basecopy.empty()){ // add an O() error term 1087 monome last(undef,ordre+1); 1088 basecopy.push_back(last); 1089 } 1090 // If first.exponent!=0 and direction==0 we can not find 1091 // first.exponent^e consistently around 0 1092 if (!direction && !is_integer(e) && !is_zero(first.exponent) ){ 1093 *logptr(contextptr) << gettext("Warning: vanishing non integral power expansion") << '\n'; 1094 /* 1095 res.clear(); 1096 first.coeff=pow(first.coeff,e,contextptr); 1097 first.exponent = first.exponent*e; 1098 res.push_back(first); 1099 first.coeff=undef; 1100 first.exponent += basecopy[0].exponent; 1101 res.push_back(first); 1102 return; 1103 */ 1104 } 1105 // answer=first.coeff^e*x^(first.exponent*e)*(1+base)^e 1106 // first (1+base)^e -> compose( [1,e,e(e-1)/2,...], base) 1107 vecteur v(1,plus_one); 1108 gen produit(e),factorielle(1); 1109 for (int i=1;i<=ordre;++i){ 1110 v.push_back(rdiv(produit,factorielle,contextptr)); 1111 produit=produit*(e-gen(i)); 1112 factorielle=factorielle*gen(i+1); 1113 } 1114 if (e.type!=_INT_ || e.val>ordre) 1115 v.push_back(undef); 1116 // COUT << v << '\n'; 1117 if (!pcompose(v,basecopy,res,contextptr)) 1118 return false; 1119 // COUT << res << '\n'; 1120 // final multiplication ans shift 1121 pshift(res,first.exponent*e,res,contextptr); 1122 return pmul(res,normalize_sqrt(pow(first.coeff,e,contextptr),contextptr),res,contextptr); 1123 } 1124 sppow(const sparse_poly1 & a,const gen & b,GIAC_CONTEXT)1125 sparse_poly1 sppow(const sparse_poly1 & a,const gen &b,GIAC_CONTEXT){ 1126 sparse_poly1 res; 1127 if (!ppow(a,b,series_default_order(contextptr),0,res,contextptr)) 1128 res=sparse_poly1(1,monome(1,undef)); 1129 return res; 1130 } 1131 pintegrate(sparse_poly1 & p,const gen & t,GIAC_CONTEXT)1132 bool pintegrate(sparse_poly1 & p,const gen & t,GIAC_CONTEXT){ 1133 sparse_poly1::iterator it=p.begin(),itend=p.end(); 1134 identificateur idu("u"); gen u(idu); 1135 for (;it!=itend;++it){ 1136 #if 1 1137 it->coeff=integrate_gen(it->coeff,t,contextptr); 1138 #else 1139 gen tmp=subst(it->coeff,t,u,false,contextptr); 1140 it->coeff=_integrate(makesequence(tmp,u,0,t),contextptr); 1141 #endif 1142 } 1143 return true; 1144 } 1145 1146 // find q such that pcompose(p,q)=x 1147 // does not take care of cst coeff of p prevert(const sparse_poly1 & p_orig,sparse_poly1 & q,GIAC_CONTEXT)1148 bool prevert(const sparse_poly1 & p_orig,sparse_poly1 & q,GIAC_CONTEXT){ 1149 sparse_poly1 p(p_orig); 1150 if (p.empty()) 1151 return false; // setsizeerr(gettext("prevert")); 1152 if (p.front().exponent==0) 1153 p.erase(p.begin()); 1154 gen ak,k,invk,b1; 1155 if (p.empty() || is_undef( (ak=p.front().coeff) ) || ck_is_positive(- (k=p.front().exponent) ,contextptr) || k.type!=_INT_ ) 1156 return false; // setsizeerr(gettext("prevert")); 1157 invk=gen(1)/k; 1158 b1=pow(ak,invk,contextptr); 1159 vecteur pv(1); 1160 sparse_poly1::const_iterator it=p.begin(),itend=p.end(); 1161 int N=0; 1162 for (;it!=itend;++it){ 1163 gen Ng=it->exponent; 1164 if (Ng.type!=_INT_) 1165 return false; // setsizeerr(); 1166 N=Ng.val; 1167 if (is_undef(it->coeff)) 1168 break; 1169 for (int n=int(pv.size());n<N;++n){ 1170 pv.push_back(0); 1171 } 1172 pv.push_back(it->coeff); 1173 } 1174 if (it==itend) 1175 N++; 1176 N=k.val*N; 1177 q.clear(); 1178 q.push_back(monome(gen(1)/b1,invk)); 1179 for (int n=2;n<N;++n){ 1180 sparse_poly1 qtemp(q),res; 1181 qtemp.push_back(monome(undef,(n+1)*invk)); 1182 if (!pcompose(pv,qtemp,res,contextptr)) 1183 return false; 1184 // find coeff of order (n+k-1)/k 1185 sparse_poly1::const_iterator jt=res.begin(),jtend=res.end(); 1186 for (;jt!=jtend;++jt){ 1187 if (jt->exponent==(n+k-1)/k) 1188 break; 1189 } 1190 if (jt!=jtend) 1191 q.push_back(monome(-jt->coeff*invk/b1,gen(n)/k)); 1192 } 1193 q.push_back(monome(undef,gen(N)/k)); 1194 return true; 1195 } 1196 in_series__SPOL1(const gen & e,const identificateur & x,const vecteur & lvx,const vecteur & lvx_s,int ordre,int direction,sparse_poly1 & s,GIAC_CONTEXT)1197 static bool in_series__SPOL1(const gen & e,const identificateur & x,const vecteur & lvx, const vecteur & lvx_s,int ordre,int direction,sparse_poly1 & s,GIAC_CONTEXT){ 1198 #ifdef TIMEOUT 1199 control_c(); 1200 #endif 1201 if (ctrl_c || interrupted) { 1202 interrupted=ctrl_c=true; 1203 return false; 1204 } 1205 s.clear(); 1206 int pos=equalposcomp(lvx,e); 1207 if (pos){ 1208 gen f=lvx_s[pos-1]; // since vectors begin at position 0 1209 if (is_zero(f)){ 1210 return true; 1211 } 1212 if (f.type==_SPOL1){ 1213 s=*(f._SPOL1ptr); 1214 return true; 1215 } 1216 if (f.type!=_VECT) 1217 return false; // settypeerr(); 1218 vecteur2sparse_poly1(*f._VECTptr,s); 1219 return true; 1220 } 1221 if ( (e.type!=_SYMB) || !contains(e,x) ){ 1222 gen en=normal(e,contextptr); 1223 if (!is_zero(en)) 1224 s.push_back(monome(en,0)); 1225 return true; 1226 } 1227 // do rational operations 1228 if (e._SYMBptr->sommet==at_plus){ 1229 if (e._SYMBptr->feuille.type!=_VECT){ 1230 return in_series__SPOL1(e._SYMBptr->feuille,x,lvx,lvx_s,ordre,direction,s,contextptr); 1231 } 1232 const_iterateur it=e._SYMBptr->feuille._VECTptr->begin(),itend=e._SYMBptr->feuille._VECTptr->end(); 1233 sparse_poly1 temp; 1234 for (;it!=itend;++it){ 1235 if (!in_series__SPOL1(*it,x,lvx,lvx_s,ordre,direction,temp,contextptr)) 1236 return false; 1237 padd(s,temp,s,contextptr); 1238 } 1239 return true; 1240 } 1241 if (e._SYMBptr->sommet==at_neg){ 1242 if (e._SYMBptr->feuille.type!=_VECT){ 1243 if (!in_series__SPOL1(e._SYMBptr->feuille,x,lvx,lvx_s,ordre,direction,s,contextptr)) 1244 return false; 1245 pneg(s,s,contextptr); 1246 return true; 1247 } 1248 const_iterateur it=e._SYMBptr->feuille._VECTptr->begin(),itend=e._SYMBptr->feuille._VECTptr->end(); 1249 sparse_poly1 temp; 1250 for (;it!=itend;++it){ 1251 if (!in_series__SPOL1(*it,x,lvx,lvx_s,ordre,direction,temp,contextptr)) 1252 return false; 1253 pneg(temp,temp,contextptr); 1254 padd(s,temp,s,contextptr); 1255 } 1256 return true; 1257 } 1258 if (e._SYMBptr->sommet==at_prod){ 1259 if (e._SYMBptr->feuille.type!=_VECT){ 1260 if (!in_series__SPOL1(e._SYMBptr->feuille,x,lvx,lvx_s,ordre,direction,s,contextptr)) 1261 return false; 1262 return true; 1263 } 1264 const_iterateur it=e._SYMBptr->feuille._VECTptr->begin(),itend=e._SYMBptr->feuille._VECTptr->end(); 1265 sparse_poly1 temp; 1266 s=sparse_poly1(1,monome(1,0)); 1267 for (;it!=itend;++it){ 1268 if (!in_series__SPOL1(*it,x,lvx,lvx_s,ordre,direction,temp,contextptr) || 1269 !pmul(s,temp,s,true,ordre,contextptr)) 1270 return false; 1271 } 1272 return true; 1273 } 1274 if (e._SYMBptr->sommet==at_inv){ 1275 if (e._SYMBptr->feuille.type==_VECT) 1276 return false; // setsizeerr(gettext("series.cc/in_series__SPOL1")); 1277 sparse_poly1 temp; 1278 if (!in_series__SPOL1(e._SYMBptr->feuille,x,lvx,lvx_s,ordre,direction,temp,contextptr)) 1279 return false; 1280 return pdiv(sparse_poly1(1,monome(1,0)),temp,s,ordre,contextptr); 1281 } 1282 if (e._SYMBptr->sommet==at_pow){ 1283 // the power is independent on x 1284 gen base=(*(e._SYMBptr->feuille._VECTptr))[0]; 1285 gen exponent=(*(e._SYMBptr->feuille._VECTptr))[1]; 1286 if (!in_series__SPOL1(base,x,lvx,lvx_s,ordre,direction,s,contextptr)) 1287 return false; 1288 return ppow(s,exponent,ordre,direction,s,contextptr); 1289 } 1290 // unknown rational operator 1291 invalidserieserr(gettext("unknown rational operator")); 1292 return false; // 1293 } 1294 find_image(const symbolic & temp__SYMB,gen & image_of_lim_point,sparse_poly1 & s,int direction,GIAC_CONTEXT)1295 static void find_image(const symbolic & temp__SYMB,gen & image_of_lim_point,sparse_poly1 & s,int direction,GIAC_CONTEXT){ 1296 if (!s.empty()){ 1297 if (s.begin()->exponent==0){ 1298 image_of_lim_point=s.begin()->coeff; 1299 // ?? FIXME ??? 1300 if (temp__SYMB.sommet!=at_abs){ 1301 s.erase(s.begin()); // remove cst coeff from s 1302 for (;!s.empty();){ 1303 s.front().coeff=normal(s.front().coeff,contextptr); 1304 if (!is_zero(s.front().coeff)) 1305 break; 1306 s.erase(s.begin()); 1307 } 1308 } 1309 } 1310 else { 1311 if (ck_is_strictly_positive(s.begin()->exponent,contextptr)) 1312 image_of_lim_point=0; 1313 else { 1314 image_of_lim_point=unsigned_inf; 1315 if ( (s.begin()->exponent.type==_INT_) && !(s.begin()->exponent.val%2) ){ // odd negative exponent 1316 if (is_strictly_positive(s.begin()->coeff,contextptr)) 1317 image_of_lim_point=plus_inf; 1318 if (is_strictly_positive(-s.begin()->coeff,contextptr)) 1319 image_of_lim_point=minus_inf; 1320 } 1321 else { // other negative exponent 1322 if (direction){ 1323 if (is_strictly_positive(s.begin()->coeff,contextptr)) 1324 image_of_lim_point=plus_inf; 1325 if (is_strictly_positive(-s.begin()->coeff,contextptr)) 1326 image_of_lim_point=minus_inf; 1327 if (direction<0){ 1328 if (s.begin()->exponent.type==_INT_) 1329 image_of_lim_point=-image_of_lim_point; 1330 else 1331 image_of_lim_point=unsigned_inf; 1332 } 1333 } 1334 } 1335 } // end negative leading exponent 1336 } // end non-zero leading exponent 1337 } // end non-empty series 1338 } 1339 find_direction(const sparse_poly1 & s,int direction,GIAC_CONTEXT)1340 static int find_direction(const sparse_poly1 & s,int direction,GIAC_CONTEXT){ 1341 int image_of_direction=0; 1342 if (!s.empty() && fastsign(s.front().coeff,0)){ 1343 if (direction) 1344 image_of_direction=1; 1345 else { 1346 if (s.front().exponent.type==_INT_) { 1347 if (s.front().exponent.val %2) 1348 image_of_direction=direction; 1349 else 1350 image_of_direction=1; 1351 } 1352 } 1353 image_of_direction=image_of_direction*fastsign(s.front().coeff,0); 1354 return image_of_direction; 1355 } 1356 return 0; 1357 } 1358 ck_is_greater(const sparse_poly1 & s1,sparse_poly1 & s2,int direction,GIAC_CONTEXT)1359 static int ck_is_greater(const sparse_poly1 & s1, sparse_poly1 & s2,int direction,GIAC_CONTEXT){ 1360 sparse_poly1 s(s2); 1361 pneg(s,s,contextptr); 1362 padd(s1,s,s,contextptr); 1363 int image_of_direction=find_direction(s,direction,contextptr); 1364 if (!image_of_direction){ 1365 cksignerr(s); 1366 return -1; 1367 } 1368 return image_of_direction==1; 1369 } 1370 1371 static gen in_limit(const gen & e,const identificateur & x,const gen & lim_point,int direction,GIAC_CONTEXT); 1372 1373 static bool mrv_lead_term(const gen & e,const identificateur & x,gen & coeff, gen & mrv_var, gen & exponent,sparse_poly1 & q,int begin_ordre,GIAC_CONTEXT,bool series); 1374 1375 vecteur integrate(const vecteur & p,const gen & shift_coeff); 1376 series(const sparse_poly1 & s_,const unary_function_ptr & u,int direction,sparse_poly1 & res,GIAC_CONTEXT)1377 bool series(const sparse_poly1 & s_,const unary_function_ptr & u,int direction,sparse_poly1 & res,GIAC_CONTEXT){ 1378 sparse_poly1 s(s_); 1379 if (s.empty()) 1380 return false; 1381 gen shift_coeff=0; 1382 gen o=porder(s); 1383 if (o==plus_inf) 1384 o=series_default_order(contextptr); 1385 else 1386 o=_floor(o,contextptr); 1387 if (o.type!=_INT_) 1388 return false; 1389 gen exponent=s.front().exponent; 1390 gen c=s.front().coeff; 1391 if (is_undef(c) || is_strictly_positive(-exponent,contextptr)) 1392 return false; 1393 if (exponent==0) 1394 s.erase(s.begin()); 1395 else 1396 c=0; 1397 gen se=u.ptr()->series_expansion(c,o.val,u,direction,shift_coeff,contextptr); 1398 if (se.type==_SPOL1){ 1399 return false; 1400 } 1401 if (se.type!=_VECT || shift_coeff!=0) 1402 return false; 1403 if (!pcompose(*se._VECTptr,s,res,contextptr)) 1404 return false; 1405 return true; 1406 } 1407 series(const sparse_poly1 & s,const unary_function_ptr & u,int direction,GIAC_CONTEXT)1408 sparse_poly1 series(const sparse_poly1 & s,const unary_function_ptr & u,int direction,GIAC_CONTEXT){ 1409 sparse_poly1 res; 1410 if (!series(s,u,direction,res,contextptr)) 1411 return sparse_poly1(1,monome(undef,undef)); 1412 return res; 1413 } 1414 series__SPOL1(const gen & e_orig,const identificateur & x,const gen & lim_point,int ordre,int direction,sparse_poly1 & s,GIAC_CONTEXT)1415 bool series__SPOL1(const gen & e_orig,const identificateur & x,const gen & lim_point,int ordre,int direction,sparse_poly1 & s,GIAC_CONTEXT){ 1416 gen e(e_orig); 1417 // fast check first 1418 if (!contains(e,x)){ 1419 s.push_back(monome(e,0)); 1420 return true; 1421 } 1422 if (e.type==_IDNT){ 1423 if (!is_zero(lim_point)) 1424 s.push_back(monome(lim_point)); 1425 if (ordre) 1426 s.push_back(monome(1,1)); 1427 else 1428 s.push_back(monome(undef,1)); 1429 return true; 1430 } 1431 if (e.type!=_SYMB) 1432 return false; // settypeerr(); // comp not allowed 1433 // rewrite cos/sin/tan constants using rootof 1434 vecteur lv1(lvar(e)),lva,lvb; 1435 gen tmp; 1436 const_iterateur lv1_it=lv1.begin(),lv1_itend=lv1.end(); 1437 for (;lv1_it!=lv1_itend;++lv1_it){ 1438 if (lv1_it->type==_SYMB){ 1439 unary_function_ptr & u =lv1_it->_SYMBptr->sommet; 1440 if ( (u==at_cos || u==at_sin || u==at_tan) && has_evalf(*lv1_it,tmp,1,contextptr)){ 1441 tmp=normal(trig2exp(*lv1_it,contextptr),contextptr); 1442 if (lop(tmp,at_exp).empty()){ 1443 lva.push_back(*lv1_it); 1444 lvb.push_back(tmp); 1445 } 1446 } 1447 } 1448 } 1449 if (!lva.empty()) 1450 e=subst(e,lva,lvb,false,contextptr); 1451 // find list of vars depending on x 1452 vecteur lvx(rlvarx(e,x)); 1453 iterateur lvx_it=lvx.begin(),lvx_end=lvx.end(); 1454 // find asymptotic series expansion of vars in lvx 1455 vecteur lvx_s; 1456 lvx_s.reserve(lvx_end-lvx_it); 1457 for (;lvx_it!=lvx_end;++lvx_it){ 1458 if (lvx_it->type==_IDNT){ 1459 sparse_poly1 tmp; 1460 if (!is_zero(lim_point)) 1461 tmp.push_back(monome(lim_point)); 1462 if (ordre) 1463 tmp.push_back(monome(1,1)); 1464 else 1465 tmp.push_back(monome(undef,1)); 1466 lvx_s.push_back(tmp); 1467 continue; 1468 } 1469 if (lvx_it->type!=_SYMB) // just in case... 1470 return false; // settypeerr(); 1471 // test for a^b 1472 symbolic temp__SYMB=*lvx_it->_SYMBptr; 1473 if (temp__SYMB.sommet==at_order_size){ 1474 if (!is_zero(lim_point)) 1475 return false; 1476 sparse_poly1 tmp; 1477 tmp.push_back(monome(undef,0)); 1478 lvx_s.push_back(tmp); 1479 continue; 1480 } 1481 if ( (temp__SYMB.sommet==at_pow) && (!contains((*temp__SYMB.feuille._VECTptr)[1],x) ) ){ 1482 if (!in_series__SPOL1((*temp__SYMB.feuille._VECTptr)[0],x,lvx,lvx_s,ordre,direction,s,contextptr)|| 1483 !ppow(s,(*temp__SYMB.feuille._VECTptr)[1],ordre,direction,s,contextptr)) 1484 return false; 1485 lvx_s.push_back(s); 1486 continue; 1487 } 1488 if (temp__SYMB.sommet==at_pow) 1489 temp__SYMB=symbolic(at_exp,(*temp__SYMB.feuille._VECTptr)[1]*ln((*temp__SYMB.feuille._VECTptr)[0],contextptr)); 1490 // Check here for logarithms if image_of_lim_point=+/-inf 1491 // In such case we must factor x^s.begin().exponent and add 1492 // it to the exponent 0 term of the series expansion of the ln 1493 if (temp__SYMB.sommet==at_ln){ 1494 // recursive call, works since lvx is sorted by increasing size 1495 if (!in_series__SPOL1(temp__SYMB.feuille,x,lvx,lvx_s,ordre,direction,s,contextptr)) 1496 return false; 1497 gen exponent=s.front().exponent; 1498 gen c=s.front().coeff; 1499 s.erase(s.begin()); 1500 bool adjust=false; 1501 if (is_positive(-c,contextptr)){ 1502 // im(ln(c)) is i*pi, but im(ln(c*x^exposant+...)) might be -i*pi 1503 // check sign of imaginary part of expansion 1504 sparse_poly1::iterator it=s.begin(),itend=s.end(); 1505 for (;it!=itend;++it){ 1506 if (is_undef(it->coeff)) 1507 break; 1508 gen tmp=im(it->coeff,contextptr); 1509 if (!is_zero(tmp)){ 1510 if (is_positive(-tmp,contextptr)) 1511 adjust=true; 1512 break; 1513 } 1514 } 1515 } 1516 if (!s.empty()){ 1517 pshift(s,-exponent,s,contextptr); 1518 if (!pdiv(s,c,s,contextptr)) 1519 return false; 1520 vecteur expansion(1,zero); 1521 expansion.reserve(ordre); 1522 for (int i=1;i<=ordre;i++){ 1523 if (i%2) 1524 expansion.push_back(inv(gen(i),contextptr)); 1525 else 1526 expansion.push_back(-inv(gen(i),contextptr)); 1527 } 1528 expansion.push_back(undef); 1529 if (!pcompose(expansion,s,s,contextptr)) 1530 return false; 1531 } 1532 if (!is_zero(exponent) || !is_one(c)){ 1533 c=ln(c,contextptr); 1534 if (adjust) 1535 c-=cst_two_pi*cst_i; 1536 s.insert(s.begin(),monome(exponent*ln(x,contextptr)+c)); 1537 } 1538 lvx_s.push_back(s); 1539 continue; 1540 COUT << s.back() << '\n'; 1541 } 1542 // test for the special case var=f(x) 1543 if ((temp__SYMB.feuille.type==_IDNT) && (temp__SYMB.sommet!=at_abs)){ 1544 // Since e contains x feuille of e must be x 1545 if (!temp__SYMB.sommet.ptr()->series_expansion){ 1546 *logptr(contextptr) << gettext("no taylor method for ") << temp__SYMB.sommet.ptr()->print(contextptr) << '\n'; 1547 return false; 1548 } 1549 gen shift_coeff; 1550 gen res=temp__SYMB.sommet.ptr()->series_expansion(lim_point,ordre,temp__SYMB.sommet,direction,shift_coeff,contextptr); 1551 if (res.type==_SPOL1) 1552 lvx_s.push_back(*res._SPOL1ptr); 1553 else {// res must be a vecteur 1554 if (res.type!=_VECT) 1555 return false; // settypeerr(gettext("series.cc 1066")); 1556 sparse_poly1 temp(vecteur2sparse_poly1(*res._VECTptr)); 1557 if (!is_zero(shift_coeff)){ 1558 pshift(temp,shift_coeff,temp,contextptr); 1559 if (is_positive(shift_coeff,contextptr)) 1560 temp.insert(temp.begin(),monome(temp__SYMB.sommet(lim_point,contextptr))); 1561 } 1562 lvx_s.push_back(temp); 1563 } 1564 continue; 1565 } 1566 // Taylor not successful: find series_expansion of arg, 1567 // compose with sommet expansion 1568 // fixme: multiargs disabled, should return a vecteur of series_exp 1569 if (temp__SYMB.sommet != at_of && temp__SYMB.feuille.type==_VECT){ 1570 int nargs=int(temp__SYMB.feuille._VECTptr->size()); 1571 if (nargs==4 && temp__SYMB.sommet==at_sum){ 1572 vecteur & tempfv=*temp__SYMB.feuille._VECTptr; 1573 gen k=tempfv[1],lo=tempfv[2],up=tempfv[3]; 1574 if (contains(k,x)){ 1575 invalidserieserr(gettext("Summation variable must be != from series expansion variable")); 1576 return false; 1577 } 1578 if (k.type!=_IDNT || derive(lo,x,contextptr)!=0 || derive(up,x,contextptr)!=0) 1579 return false; // setsizeerr(); 1580 sparse_poly1 p; 1581 if (!in_series__SPOL1(tempfv[0],x,lvx,lvx_s,ordre,direction,p,contextptr)) 1582 return false; 1583 // sum lvx_s 1584 sparse_poly1::iterator it=p.begin(),itend=p.end(); 1585 for (;it!=itend;++it){ 1586 if (is_undef(it->coeff)) 1587 break; 1588 it->coeff=_sum(makesequence(it->coeff,k,lo,up),contextptr); 1589 } 1590 lvx_s.push_back(p); 1591 continue; 1592 //invalidserieserr(gettext("taylor of sum not implemented")); 1593 //return false; 1594 } 1595 if (temp__SYMB.sommet==at_euler_mac_laurin){ 1596 if (nargs!=5){ 1597 invalidserieserr(gettext("Integral must be definite")); 1598 return false; 1599 } 1600 vecteur & tempfv=*temp__SYMB.feuille._VECTptr; 1601 gen k=tempfv[2]; 1602 if (contains(k,x)){ 1603 invalidserieserr(gettext("Summation variable must be != from series expansion variable")); 1604 return false; 1605 } 1606 if (k.type!=_IDNT) 1607 return false; // setsizeerr(); 1608 // find upper and lower bound limit 1609 gen lower=tempfv[3],upper=tempfv[4],f=tempfv[0],F=tempfv[1]; 1610 gen fdiff=derive(f,k,contextptr); 1611 // first add integral part 1612 gen eff = preval(F,k,lower,upper,contextptr); 1613 eff += (limit(f,*k._IDNTptr,upper,-1,contextptr)+limit(f,*k._IDNTptr,lower,1,contextptr))/2; 1614 // then Bernoulli part 1615 if (is_undef(fdiff)) 1616 return false; 1617 for (int i=1;i<ordre;i++){ 1618 gen add= limit(fdiff,*k._IDNTptr,upper,-1,contextptr)-limit(fdiff,*k._IDNTptr,lower,1,contextptr); 1619 add=add*bernoulli(2*i)/factorial(2*i); 1620 eff += add; // fdiff flimdiff 2 fois 1621 fdiff=derive(fdiff,k,contextptr); 1622 fdiff=ratnormal(derive(fdiff,k,contextptr),contextptr); 1623 if (is_undef(fdiff)) 1624 return false; 1625 } 1626 // must do a recursive call since eff may contain new functions 1627 gen coeff,mrv_var,exponent; 1628 eff =subst(eff,x,inv(x,contextptr),true,contextptr); 1629 if (!mrv_lead_term(eff,x,coeff,mrv_var,exponent,s,ordre,contextptr,true)) 1630 return false; 1631 lvx_s.push_back(s); 1632 continue; 1633 // never reached setsizeerr(); 1634 } 1635 if (temp__SYMB.sommet==at_igamma_exp){ 1636 if (nargs!=2){ 1637 invalidserieserr(gettext("igamma: bad arg number")); 1638 return false; 1639 } 1640 vecteur & tempfv=*temp__SYMB.feuille._VECTptr; 1641 gen a=tempfv[0]; 1642 sparse_poly1 p; 1643 if (!in_series__SPOL1(tempfv[1],x,lvx,lvx_s,ordre,direction,p,contextptr)) 1644 return false; 1645 gen image_of_lim_point; 1646 find_image(temp__SYMB,image_of_lim_point,p,direction,contextptr); 1647 if (!is_inf(image_of_lim_point)) 1648 return false; 1649 // invert series expansion 1650 sparse_poly1 stmp; 1651 if (!pdiv(sparse_poly1(1,monome(1,0)),p,stmp,ordre,contextptr)) 1652 return false; 1653 p=stmp; 1654 // igamma(a,x)=Gamma(a)-int_x^inf exp(-t)*t^(a-1) dt 1655 // =Gamma(a)-exp(-x)*x^(a-1)-(a-1)*int_x^inf exp(-t)*t^(a-2) dt 1656 // ... 1657 // igamma_replace(a,x)=Gamma(a)-exp(-x)*igamma_exp(a,x) 1658 // therefore expansion of igamma_exp(a,x)=x^(a-1)[1+(a-1)/x+(a-1)*(a-2)/x^2...] 1659 vecteur v(ordre+1); gen facti(1); 1660 for (int i=0;i<=ordre;++i){ 1661 v[i]=facti; 1662 facti=(a-i-1)*facti; 1663 } 1664 v.push_back(undef); 1665 if (!pcompose(v,p,s,contextptr)) 1666 return false; 1667 pshift(s,1-a,s,contextptr); 1668 lvx_s.push_back(s); 1669 continue; 1670 } 1671 if (temp__SYMB.sommet==at_lower_incomplete_gamma){ 1672 if (nargs!=2){ 1673 invalidserieserr(gettext("igamma: bad arg number")); 1674 return false; 1675 } 1676 // series expansion of derivative 1677 vecteur & tempfv=*temp__SYMB.feuille._VECTptr; 1678 gen a=tempfv[0]; 1679 sparse_poly1 p; 1680 if (!in_series__SPOL1(tempfv[1],x,lvx,lvx_s,ordre,direction,p,contextptr)) 1681 return false; 1682 gen image_of_lim_point; 1683 find_image(temp__SYMB,image_of_lim_point,p,direction,contextptr); 1684 if (is_inf(image_of_lim_point)) 1685 return false; // inf is prevented by limit_symbolic_preprocessing 1686 if (is_zero(image_of_lim_point)){ 1687 int image_of_direction=0; 1688 image_of_direction = find_direction(p,direction,contextptr); 1689 if (image_of_direction==0) 1690 return false; 1691 vecteur v(ordre+1); gen facti(1); 1692 for (int i=0;i<=ordre;++i){ 1693 v[i]=inv(facti*(a+i),contextptr); 1694 facti=-(i+1)*facti; 1695 } 1696 if (!pcompose(v,p,s,contextptr)) 1697 return false; 1698 if (image_of_direction==-1) 1699 pneg(p,p,contextptr); 1700 if (!ppow(p,a,ordre,image_of_direction,p,contextptr)) 1701 return false; 1702 if (image_of_direction==-1) 1703 pneg(p,p,contextptr); 1704 if (!pmul(p,s,s,true,ordre,contextptr)) 1705 return false; 1706 } 1707 else { 1708 vecteur v; 1709 gen der; 1710 if (is_positive(image_of_lim_point,contextptr)) 1711 der=pow(x,a-1,contextptr)*exp(-x,contextptr); 1712 else 1713 der=-pow(-x,a-1,contextptr)*exp(-x,contextptr); 1714 if (!taylor(der,x,image_of_lim_point,ordre,v,contextptr)) 1715 return false; 1716 v=integrate(v,1); 1717 gen intcst=_lower_incomplete_gamma(makesequence(a,lim_point),contextptr); 1718 v.insert(v.begin(),intcst); 1719 if (!pcompose(v,p,s,contextptr)) 1720 return false; 1721 } 1722 lvx_s.push_back(s); 1723 continue; 1724 } 1725 if (temp__SYMB.sommet==at_integrate){ 1726 if (nargs!=4){ 1727 invalidserieserr(gettext("Integral must be definite")); 1728 return false; 1729 } 1730 vecteur & tempfv=*temp__SYMB.feuille._VECTptr; 1731 gen t=tempfv[1]; 1732 if (contains(t,x)){ 1733 invalidserieserr(gettext("Integration variable must be != from series expansion variable")); 1734 return false; 1735 } 1736 #if 1 1737 if (!contains(tempfv[0],x)){ 1738 vecteur v; 1739 // int(f(t),t,m(x),M(x))=F(M(x))-F(m(x)) 1740 // hence derivative is M'(x)*f(M(x))-m'(x)*f(m(x)) 1741 gen g=derive(tempfv[3],x,contextptr)*subst(tempfv[0],tempfv[1],tempfv[3],false,contextptr)-derive(tempfv[2],x,contextptr)*subst(tempfv[0],tempfv[1],tempfv[2],false,contextptr); 1742 g=exp2pow(g,contextptr); 1743 s.clear(); 1744 if (!series__SPOL1(g,x,lim_point,ordre,direction,s,contextptr)){ 1745 return false; 1746 } 1747 gen s0=eval(subst(temp__SYMB,x,lim_point,false,contextptr),1,contextptr); 1748 sparse_poly1 p; 1749 if (!is_zero(s0)) 1750 p.push_back(monome(s0,0)); 1751 sparse_poly1::const_iterator it=s.begin(),itend=s.end(); 1752 for (;it!=itend;++it){ 1753 gen i=it->exponent; 1754 if (i==-1) // should be i<=-1? 1755 return false; 1756 p.push_back(monome(it->coeff/(i+1),i+1)); 1757 } 1758 lvx_s.push_back(p); 1759 continue; 1760 } 1761 #else // old code 1762 if (!contains(tempfv[0],x)){ 1763 vecteur v; 1764 if (!taylor(temp__SYMB,x,lim_point,ordre,v,contextptr)) 1765 return false; 1766 s.clear(); 1767 for (int i=0;i<v.size();++i){ 1768 s.push_back(monome(v[i],i)); 1769 } 1770 lvx_s.push_back(s); 1771 continue; 1772 } 1773 #endif 1774 if (!in_series__SPOL1(tempfv[0],x,lvx,lvx_s,ordre,direction,s,contextptr)) 1775 return false; 1776 // FIXME if tempfv[3] and tempfv[2] tends to the same limit l 1777 // we may expand tempfv[0] w.r.t. t at l before integration 1778 1779 // integrate s term by term wrt t 1780 if (!pintegrate(s,t,contextptr)) 1781 return false; 1782 gen remains,primit=sparse_poly12gen(s,x,remains,false); 1783 // then compose primit at bounds and substract 1784 primit=subst(primit,t,tempfv[3],false,contextptr)-subst(primit,t,tempfv[2],false,contextptr); 1785 if (!series__SPOL1(primit,x,lim_point,ordre,direction,s,contextptr)) 1786 return false; 1787 // add remains to s: 1788 // int(remains,t,tempfv[2],tempfv[3]) = 1789 // (tempfv[3]-tempfv[2])*remains(theta) with theta in interval 1790 remains=remains*(tempfv[3]-tempfv[2]); 1791 sparse_poly1 p; 1792 if (!series__SPOL1(remains,x,lim_point,ordre,direction,p,contextptr)) 1793 return false; 1794 if (p.empty()){ 1795 invalidserieserr(gettext("Can not expand remainder of integrand")); 1796 return false; 1797 } 1798 p=sparse_poly1(p.begin(),p.begin()+1); 1799 p.front().coeff=undef; 1800 s=spadd(p,s,contextptr); 1801 lvx_s.push_back(s); 1802 continue; 1803 } 1804 const_iterateur fit=temp__SYMB.feuille._VECTptr->begin(),fitend=temp__SYMB.feuille._VECTptr->end(); 1805 if (temp__SYMB.sommet==at_Psi || temp__SYMB.sommet==at_Eta || temp__SYMB.sommet==at_Zeta){ 1806 if (!in_series__SPOL1(*fit,x,lvx,lvx_s,ordre,direction,s,contextptr)) 1807 return false; 1808 } 1809 else { 1810 bool ok=false; 1811 dbgprint_vector <sparse_poly1> vs; 1812 vs.reserve(fitend-fit); 1813 for (;fit!=fitend;++fit){ 1814 if (!in_series__SPOL1(*fit,x,lvx,lvx_s,ordre,direction,s,contextptr)) 1815 return false; 1816 vs.push_back(s); 1817 } 1818 if (temp__SYMB.sommet==at_max){ 1819 ok=true; 1820 int testck=ck_is_greater(vs.front(),vs.back(),direction,contextptr); 1821 if (testck==-1) 1822 return false; 1823 if (testck) 1824 s=vs.front(); 1825 else 1826 s=vs.back(); 1827 } 1828 if (temp__SYMB.sommet==at_min){ 1829 ok=true; 1830 int testck=ck_is_greater(vs.front(),vs.back(),direction,contextptr); 1831 if (testck==-1) 1832 return false; 1833 if (testck) 1834 s=vs.back(); 1835 else 1836 s=vs.front(); 1837 } 1838 if (!ok){ 1839 invalidserieserr(gettext(" multiargs not implemented")); 1840 return false; 1841 } 1842 lvx_s.push_back(s); 1843 continue; 1844 } 1845 } 1846 else { // 1-arg function 1847 if (temp__SYMB.sommet==at_of){ 1848 gen & tf=temp__SYMB.feuille; 1849 if (tf.type==_VECT && tf._VECTptr->size()==2){ 1850 gen tff=tf._VECTptr->front(); 1851 gen tfx=tf._VECTptr->back(); 1852 if (!in_series__SPOL1(tfx,x,lvx,lvx_s,ordre,direction,s,contextptr)) return false; // s<-arg 1853 } 1854 } 1855 else { 1856 if (!in_series__SPOL1(temp__SYMB.feuille,x,lvx,lvx_s,ordre,direction,s,contextptr)) return false; // s<-arg 1857 } 1858 } // end 1-arg function 1859 gen image_of_lim_point; 1860 find_image(temp__SYMB,image_of_lim_point,s,direction,contextptr); 1861 int image_of_direction=0; 1862 image_of_direction = find_direction(s,direction,contextptr); 1863 // Symbolic series expansion f(x), f is assumed to be analytic 1864 if (temp__SYMB.sommet==at_of){ 1865 gen & tf=temp__SYMB.feuille; 1866 if (tf.type==_VECT && tf._VECTptr->size()==2){ 1867 gen tff=tf._VECTptr->front(); 1868 gen tfx=tf._VECTptr->back(); 1869 // Symbolic Taylor expansion of tff 1870 vecteur expansion(1,symbolic(at_of,makesequence(tff,image_of_lim_point))); 1871 expansion.reserve(ordre); 1872 for (int i=1;i<=ordre;i++){ 1873 gen fn; 1874 if (i==1) 1875 fn=symbolic(at_of,makesequence(symbolic(at_function_diff,tff),image_of_lim_point)); 1876 else 1877 fn=symbolic(at_of,makesequence(symbolic(at_of,makesequence(symbolic(at_composepow,makesequence(at_function_diff,i)),tff)),image_of_lim_point)); 1878 expansion.push_back(fn/factorial(i)); 1879 } 1880 expansion.push_back(undef); 1881 if (!pcompose(expansion,s,s,contextptr)) 1882 return false; 1883 lvx_s.push_back(s); 1884 continue; 1885 } 1886 } 1887 if (temp__SYMB.sommet==at_abs){ 1888 if (!image_of_direction){ 1889 *logptr(contextptr) << gettext("Sign error ") << s << '\n'; 1890 return false; // cksignerr(s); 1891 } 1892 if (image_of_direction==-1) 1893 pneg(s,s,contextptr); 1894 lvx_s.push_back(s); 1895 continue; 1896 } 1897 if (is_inf(image_of_lim_point)){ 1898 // check for sin/cos 1899 if (temp__SYMB.sommet==at_cos || temp__SYMB.sommet==at_sin){ 1900 // split the series expansion in two parts, one tending -> 0 1901 sparse_poly1::iterator it=s.begin(),itend=s.end(); 1902 for (;it!=itend;++it){ 1903 if (ck_is_strictly_greater(it->exponent,zero,contextptr)) 1904 break; 1905 } 1906 sparse_poly1 s0(it,s.end()),s1(s.begin(),it); 1907 // expansion is done at s0 1908 image_of_lim_point=s1; 1909 s=s0; 1910 } 1911 else { 1912 // the function is assumed to have an expansion at infinity 1913 // invert series expansion 1914 sparse_poly1 stmp; 1915 if (!pdiv(sparse_poly1(1,monome(1,0)),s,stmp,ordre,contextptr)) 1916 return false; 1917 s=stmp; 1918 } 1919 } 1920 gen shift_coeff; 1921 if (!temp__SYMB.sommet.ptr()->series_expansion){ 1922 *logptr(contextptr) << string(gettext("Not expandable "))+temp__SYMB.sommet.ptr()->s << '\n'; 1923 return false; 1924 } 1925 int addorder=0; 1926 gen expansion; 1927 if ( (temp__SYMB.sommet==at_Psi ||temp__SYMB.sommet==at_Eta ||temp__SYMB.sommet==at_Zeta) && temp__SYMB.feuille.type==_VECT){ 1928 if (temp__SYMB.feuille._VECTptr->size()!=2 ) 1929 return false; // setsizeerr(); 1930 addorder=temp__SYMB.feuille._VECTptr->back().val; 1931 if (addorder<=0){ 1932 *logptr(contextptr) << gettext("Psi/Zeta/Eta: bad second argument") << '\n'; 1933 return false; 1934 } 1935 if (temp__SYMB.sommet==at_Psi) 1936 expansion=at_Psi_minus_ln->ptr()->series_expansion(image_of_lim_point,ordre+addorder,temp__SYMB.sommet,image_of_direction,shift_coeff,contextptr); 1937 } 1938 if (expansion==0) 1939 expansion=temp__SYMB.sommet.ptr()->series_expansion(image_of_lim_point,ordre+addorder,temp__SYMB.sommet,image_of_direction,shift_coeff,contextptr); 1940 if (expansion.type==_VECT){ 1941 if (addorder){ 1942 // derive expansion 1943 vecteur & v =*expansion._VECTptr; 1944 for (int i=0;i<addorder;++i){ 1945 int vs=int(v.size()); 1946 if (is_zero(shift_coeff)){ 1947 vecteur w(vs-1); 1948 for (int j=1;j<vs;++j){ 1949 w[j-1]=j*v[j]; 1950 } 1951 v=w; 1952 } 1953 else { 1954 for (int j=0;j<vs;++j){ 1955 v[j]=-v[j]*(j+shift_coeff); 1956 } 1957 shift_coeff += 1; 1958 } 1959 } 1960 // final correction for Psi 1961 if (temp__SYMB.sommet==at_Psi){ 1962 v.insert(v.begin(),gen((addorder%2)?1:-1)/factorial(addorder)); 1963 shift_coeff -= 1; 1964 } 1965 } 1966 if (is_zero(shift_coeff)){ 1967 if (!pcompose(*expansion._VECTptr,s,s,contextptr)) 1968 return false; 1969 } 1970 else { 1971 sparse_poly1 temp; 1972 if (!ppow(s,shift_coeff,ordre,direction,temp,contextptr)) 1973 return false; 1974 if (!pcompose(*expansion._VECTptr,s,s,contextptr) || 1975 !pmul(s,temp,s,true,ordre,contextptr)) 1976 return false; 1977 if (is_positive(shift_coeff,contextptr)){ 1978 gen imtemp; 1979 if (addorder>0) 1980 imtemp=temp__SYMB.sommet(makesequence(image_of_lim_point,addorder),contextptr); 1981 else 1982 imtemp=temp__SYMB.sommet(image_of_lim_point,contextptr); 1983 if (!is_zero(imtemp)) 1984 s.insert(s.begin(),monome(imtemp)); 1985 } 1986 } 1987 } 1988 else { 1989 s.clear(); 1990 s.push_back(monome(undef,minus_inf)); 1991 return true; 1992 } 1993 lvx_s.push_back(s); 1994 // fixme: add support for sparse_poly1 composition 1995 } // end loop lvx_it!=lvx_end 1996 return in_series__SPOL1(e,x,lvx,lvx_s,ordre,direction,s,contextptr); 1997 } 1998 series__SPOL1(const gen & e,const identificateur & x,const gen & lim_point,int ordre,int direction,GIAC_CONTEXT)1999 sparse_poly1 series__SPOL1(const gen & e,const identificateur & x,const gen & lim_point,int ordre,int direction,GIAC_CONTEXT){ 2000 sparse_poly1 s; 2001 if (!series__SPOL1(e,x,lim_point,ordre,direction,s,contextptr)) 2002 s=sparse_poly1(1,monome(1,undef)); 2003 return s; 2004 } 2005 ck_series__SPOL1(const gen & e,const identificateur & x,const gen & lim_point,int ordre,int direction,GIAC_CONTEXT)2006 static sparse_poly1 ck_series__SPOL1(const gen & e,const identificateur & x,const gen & lim_point,int ordre,int direction,GIAC_CONTEXT){ 2007 sparse_poly1 s; 2008 if (!series__SPOL1(e,x,lim_point,ordre,direction,s,contextptr)){ 2009 s=sparse_poly1(1,monome(1,undef)); 2010 return s; 2011 } 2012 // if s is not at order ordre, ask again with a modified ordre 2013 gen true_order=porder(s); 2014 if (true_order.type==_INT_ && true_order.val<=ordre){ 2015 if (series__SPOL1(e,x,lim_point,ordre+1+ordre-true_order.val,direction,s,contextptr)){ 2016 if (!s.empty()) 2017 ptruncate(s,ordre-s.front().exponent,contextptr); 2018 } 2019 else 2020 s=sparse_poly1(1,monome(1,undef)); 2021 } 2022 // truncate s 2023 for (unsigned i=0;i<s.size();++i){ 2024 if (is_strictly_greater(s[i].exponent,ordre,contextptr)){ 2025 s[i].coeff=undef; 2026 if (i<s.size()-1) 2027 s.erase(s.begin()+i+1,s.end()); 2028 } 2029 } 2030 return s; 2031 } 2032 2033 #ifdef DEBUG_SUPPORT inutile(sparse_poly1 & s)2034 static void inutile(sparse_poly1 & s){ 2035 s.dbgprint(); 2036 } 2037 #endif 2038 2039 // *********************** 2040 // LIMITS 2041 // *********************** 2042 contains(const vecteur & v,const gen & elem)2043 bool contains(const vecteur & v,const gen & elem){ 2044 vecteur::const_iterator it=v.begin(),itend=v.end(); 2045 for (;it!=itend;++it) 2046 if (contains(*it,elem)) 2047 return true; 2048 return false; 2049 } 2050 contains(const gen & e,const gen & elem)2051 bool contains(const gen & e,const gen & elem){ 2052 if (e==elem) 2053 return true; 2054 if (e.type==_VECT){ 2055 return contains(*e._VECTptr,elem); 2056 } 2057 if (e.type==_SYMB){ 2058 return contains(e._SYMBptr->feuille,elem); 2059 } 2060 if (e.type==_FRAC) 2061 return contains(e._FRACptr->num,elem) || contains(e._FRACptr->den,elem); 2062 #if defined HAVE_LIBMPFI && !defined NO_RTTI 2063 if (e.type==_REAL){ 2064 if (real_interval * ptr=dynamic_cast<real_interval *>(e._REALptr)){ 2065 mpfr_t tmp; mpfr_init2(tmp,mpfi_get_prec(ptr->infsup)); 2066 mpfi_get_left(tmp,ptr->infsup); 2067 gen einf=real_object(tmp); 2068 mpfi_get_right(tmp,ptr->infsup); 2069 gen esup=real_object(tmp); 2070 gen eleminf,elemsup; 2071 if (elem.type!=_REAL) 2072 ptr=0; 2073 else 2074 ptr=dynamic_cast<real_interval *>(elem._REALptr); 2075 if (ptr){ 2076 mpfi_get_left(tmp,ptr->infsup); 2077 eleminf=real_object(tmp); 2078 mpfi_get_right(tmp,ptr->infsup); 2079 elemsup=real_object(tmp); 2080 } 2081 else { 2082 eleminf=elem; elemsup=elem; 2083 } 2084 mpfr_clear(tmp); 2085 return is_greater(esup,elemsup,context0) && is_greater(eleminf,einf,context0); 2086 } 2087 } 2088 #endif 2089 return false; 2090 } 2091 lvarx(const gen & e,const gen & x,bool test)2092 vecteur lvarx(const gen &e,const gen & x,bool test){ 2093 vecteur v(lvar(e)); 2094 vecteur res; 2095 vecteur::const_iterator it=v.begin(),itend=v.end(); 2096 for (;it!=itend;++it){ 2097 // remove at_of if the function of of is x 2098 vecteur l=lop(*it,at_of); 2099 int i; 2100 for (i=0;i<l.size();++i){ 2101 if (contains(l[i]._SYMBptr->feuille[0],x)) 2102 break; 2103 } 2104 if (i<l.size()) 2105 continue; 2106 // remove ^ if exponent does not depend on x 2107 if ( (it->type==_SYMB) 2108 && ( (it->_SYMBptr->sommet==at_pow && !contains((*(it->_SYMBptr->feuille._VECTptr))[1],x)) || 2109 (it->_SYMBptr->sommet==at_NTHROOT && !contains((*(it->_SYMBptr->feuille._VECTptr))[0],x)) ) 2110 ){ 2111 vecteur tmp(lvarx((*(it->_SYMBptr->feuille._VECTptr))[(it->_SYMBptr->sommet==at_pow)?0:1],x)); 2112 const_iterateur it=tmp.begin(),itend=tmp.end(); 2113 for (;it!=itend;++it){ 2114 if (!equalposcomp(res,*it)) 2115 res.push_back(*it); 2116 } 2117 } 2118 else { 2119 if ( (!test || res.empty() || *it!=x ) && contains(*it,x) && !equalposcomp(res,*it)) 2120 res.push_back(*it); 2121 } 2122 } 2123 return res; 2124 } 2125 rlvarx(const gen & e,const gen & xgen,vecteur & res)2126 void rlvarx(const gen &e,const gen & xgen,vecteur & res){ 2127 const vecteur & v=lvar(e); 2128 vecteur::const_iterator it=v.begin(),itend=v.end(); 2129 for (;it!=itend;++it){ 2130 if (!contains(*it,xgen) || equalposcomp(res,*it)) 2131 continue; 2132 // recursive call 2133 res.push_back(*it); 2134 if (it->is_symb_of_sommet(at_derive) && it->_SYMBptr->feuille.type==_VECT && it->_SYMBptr->feuille._VECTptr->size()==3 && it->_SYMBptr->feuille._VECTptr->back().type==_INT_){ 2135 int n=it->_SYMBptr->feuille._VECTptr->back().val; 2136 for (--n;n>1;--n){ 2137 res.push_back(symbolic(at_derive,makesequence(it->_SYMBptr->feuille._VECTptr->front(),(*it->_SYMBptr->feuille._VECTptr)[1],n))); 2138 } 2139 res.push_back(symbolic(at_derive,makesequence(it->_SYMBptr->feuille._VECTptr->front(),(*it->_SYMBptr->feuille._VECTptr)[1]))); 2140 } 2141 if (it->type==_SYMB) { 2142 rlvarx(it->_SYMBptr->feuille,xgen,res); 2143 if ( (it->_SYMBptr->sommet==at_pow) 2144 && contains((*(it->_SYMBptr->feuille._VECTptr))[1],xgen) ) 2145 rlvarx(symbolic(at_ln,(*(it->_SYMBptr->feuille._VECTptr))[0]),xgen,res); 2146 } 2147 } 2148 } 2149 tri_rlvarx(vecteur & v)2150 void tri_rlvarx(vecteur & v){ 2151 int s=v.size(); 2152 for (;;){ 2153 bool ok=true; 2154 for (int i=0;i<s-1;++i){ 2155 if (symb_size_less(v[i+1],v[i])){ 2156 ok=false; 2157 swapgen(v[i],v[i+1]); 2158 } 2159 } 2160 if (ok) 2161 return; 2162 } 2163 } 2164 rlvarx(const gen & e,const gen & x)2165 vecteur rlvarx(const gen &e,const gen & x){ 2166 vecteur res; 2167 rlvarx(e,x,res); 2168 #ifdef FXCG 2169 tri_rlvarx(res); 2170 #else 2171 gen_sort_f(res.begin(),res.end(),symb_size_less); 2172 #endif 2173 return res; 2174 } 2175 upscale(gen & e,const identificateur & x,GIAC_CONTEXT)2176 static void upscale(gen & e,const identificateur & x,GIAC_CONTEXT){ 2177 vecteur a_remplacer,remplacer_par; 2178 a_remplacer.push_back(ln(x,contextptr)); 2179 remplacer_par.push_back(x); 2180 a_remplacer.push_back(x); 2181 remplacer_par.push_back(exp(x,contextptr)); 2182 e=subst(e,a_remplacer,remplacer_par,false,contextptr); 2183 } 2184 downscale(gen & e,const identificateur & x,GIAC_CONTEXT)2185 static void downscale(gen & e,const identificateur & x,GIAC_CONTEXT){ 2186 vecteur a_remplacer,remplacer_par; 2187 a_remplacer.push_back(exp(x,contextptr)); 2188 remplacer_par.push_back(x); 2189 a_remplacer.push_back(x); 2190 remplacer_par.push_back(ln(x,contextptr)); 2191 e=subst(e,a_remplacer,remplacer_par,false,contextptr); 2192 } 2193 2194 /* 2195 gen pow2exp(const gen & e,const identificateur & x){ 2196 if (e.type==_VECT){ 2197 const_iterateur it=e._VECTptr->begin(),itend=e._VECTptr->end(); 2198 vecteur v; 2199 v.reserve(itend-it); 2200 for (;it!=itend;++it) 2201 v.push_back(pow2exp(*it,x)); 2202 return v; 2203 } 2204 if (e.type!=_SYMB) 2205 return e; 2206 if ( e._SYMBptr->sommet==at_pow && contains((*(e._SYMBptr->feuille._VECTptr))[1],x)) 2207 return exp(pow2exp((*(e._SYMBptr->feuille._VECTptr))[1],x)*pow2exp(ln((*(e._SYMBptr->feuille._VECTptr))[0]),x)); 2208 if ( e._SYMBptr->sommet==at_tan && contains (e._SYMBptr->feuille,x)) 2209 return symbolic(at_sin,pow2exp(e._SYMBptr->feuille,x))/symbolic(at_cos,pow2exp(e._SYMBptr->feuille,x)); 2210 return e._SYMBptr->sommet(pow2exp(e._SYMBptr->feuille,x),contextptr); 2211 } 2212 */ 2213 check_bounded(const gen & g,GIAC_CONTEXT)2214 static int check_bounded(const gen & g,GIAC_CONTEXT){ 2215 vecteur v=loptab(g,sincostan_tab); 2216 int vs=int(v.size()); 2217 vecteur w; 2218 for (int i=0;i<vs;++i){ 2219 if (v[i].type==_SYMB && v[i]._SYMBptr->feuille.type==_SPOL1) 2220 w.push_back(v[i]); 2221 } 2222 if (w.empty()) 2223 return 0; 2224 for (unsigned i=0;i<w.size();++i){ 2225 vecteur lv=lvarxwithinv(g,w[i],contextptr); 2226 if (lv.empty()) 2227 continue; 2228 if (lv.size()>=2 || lv[0]!=w[i]){ 2229 //gensizeerr("Limit probably undefined, algorithm unable to handle "+g.print(contextptr)); 2230 return -1; 2231 } 2232 } 2233 return 1; 2234 } 2235 2236 // specialization equalposcomp(const std::vector<const unary_function_ptr * > & v,unary_function_ptr * w)2237 static int equalposcomp(const std::vector<const unary_function_ptr *> & v,unary_function_ptr * w){ 2238 int n=1; 2239 for (std::vector<const unary_function_ptr *>::const_iterator it=v.begin();it!=v.end();++it){ 2240 if (*(*it)==*w) 2241 return n; 2242 else 2243 n++; 2244 } 2245 return 0; 2246 } 2247 ln_expand0_(const gen & e,GIAC_CONTEXT)2248 static gen ln_expand0_(const gen & e,GIAC_CONTEXT){ 2249 if (e.type!=_SYMB) 2250 return ln(e,contextptr); 2251 if (e._SYMBptr->sommet==at_exp) 2252 return e._SYMBptr->feuille; 2253 if (e._SYMBptr->sommet==at_prod) 2254 return symbolic(at_plus,apply(e._SYMBptr->feuille,ln_expand0_,contextptr)); 2255 if (e._SYMBptr->sommet==at_inv) 2256 return -ln_expand0_(e._SYMBptr->feuille,contextptr); 2257 if (e._SYMBptr->sommet==at_pow){ 2258 gen & tmp=e._SYMBptr->feuille; 2259 if (tmp.type==_VECT && tmp._VECTptr->size()==2) 2260 return tmp._VECTptr->back()*ln_expand0_(tmp._VECTptr->front(),contextptr); 2261 } 2262 return ln(e,contextptr); 2263 } 2264 ln_expand_(const gen & e0,GIAC_CONTEXT)2265 static gen ln_expand_(const gen & e0,GIAC_CONTEXT){ 2266 gen e(factor(e0,false,contextptr)); 2267 return ln_expand0_(e,contextptr); 2268 } 2269 exp_series_(const gen & e0,GIAC_CONTEXT)2270 static gen exp_series_(const gen & e0,GIAC_CONTEXT){ 2271 vecteur v=lop(e0,at_ln); 2272 if (v.size()==1 && is_integer(v.front()._SYMBptr->feuille)){ 2273 gen a,b; 2274 if (is_linear_wrt(e0,v.front(),a,b,contextptr)) 2275 return exp(b,contextptr)*pow(v.front()._SYMBptr->feuille,a,contextptr); 2276 } 2277 return exp(e0,contextptr); 2278 } 2279 remove_lnexp(const gen & e,GIAC_CONTEXT)2280 static gen remove_lnexp(const gen & e,GIAC_CONTEXT){ 2281 vector<const unary_function_ptr *> v(1,at_ln); 2282 v.push_back(at_exp); 2283 vector< gen_op_context > w(1,&ln_expand_); 2284 w.push_back(&exp_series_); 2285 return subst(e,v,w,false,contextptr); 2286 } 2287 limit_symbolic_preprocess(const gen & e0,const identificateur & x,const gen & lim_point,int direction,GIAC_CONTEXT)2288 gen limit_symbolic_preprocess(const gen & e0,const identificateur & x,const gen & lim_point,int direction,GIAC_CONTEXT){ 2289 // FIXME: add support for int and sum 2290 gen e=factorial2gamma(e0,contextptr); 2291 int trigs=loptab(e,sincostan_tab).size(); 2292 if (trigs>=2){ 2293 gen tmp=trigcos(e,contextptr); 2294 int trigtmps=loptab(tmp,sincostan_tab).size(); 2295 if (trigtmps<trigs){ 2296 e=ratnormal(tmp); 2297 trigs=loptab(e,sincostan_tab).size(); 2298 } 2299 tmp=trigsin(e,contextptr); 2300 trigtmps=loptab(tmp,sincostan_tab).size(); 2301 if (trigtmps<trigs){ 2302 e=ratnormal(tmp); 2303 trigs=loptab(e,sincostan_tab).size(); 2304 } 2305 tmp=trigtan(e,contextptr); 2306 trigtmps=loptab(tmp,sincostan_tab).size(); 2307 if (trigtmps<trigs){ 2308 e=ratnormal(tmp); 2309 } 2310 } 2311 gen first_try=subst(e,x,lim_point,false,contextptr); 2312 first_try=simplifier(first_try,contextptr); 2313 if (!contains(lidnt(first_try),unsigned_inf)){ 2314 gen chknum; 2315 bool hasnum=has_evalf(first_try,chknum,1,contextptr); 2316 first_try=recursive_ratnormal(first_try,contextptr); 2317 gen chk=recursive_normal(first_try,contextptr); 2318 if (hasnum && !is_undef(chk) && abs(chk-chknum,contextptr)>1e-10 && abs(1-chk/chknum,contextptr)>1e-10) 2319 e=_simplify(e,contextptr); 2320 } 2321 // Find functions depending of x in e which are in the list 2322 // If their argument tends to +/-infinity, replace these functions 2323 vecteur v=rlvarx(e,x); 2324 int vs=int(v.size()),pos1,pos2=0; 2325 vecteur v1,v2; 2326 for (int i=0;i<vs;++i){ 2327 if (v[i].type==_SYMB){ 2328 if (v[i].is_symb_of_sommet(at_ln)){ 2329 gen g=limit(v[i]._SYMBptr->feuille,x,lim_point,direction,contextptr); 2330 if (is_inf(g) && g!=plus_inf) 2331 return gensizeerr(gettext("ln of unsigned or minus infinity")); 2332 } 2333 if (v[i].is_symb_of_sommet(at_sinh) || v[i].is_symb_of_sommet(at_cosh) || v[i].is_symb_of_sommet(at_tanh)){ 2334 gen g=limit(v[i]._SYMBptr->feuille,x,lim_point,direction,contextptr); 2335 if (is_inf(g)){ 2336 v1.push_back(v[i]); 2337 v2.push_back(hyp2exp(v[i],contextptr)); 2338 continue; 2339 } 2340 } 2341 if (v[i].is_symb_of_sommet(at_sum)){ 2342 gen tmp; 2343 if (v[i]._SYMBptr->feuille.type==_VECT && v[i]._SYMBptr->feuille._VECTptr->size()==4){ 2344 vecteur & vv=*v[i]._SYMBptr->feuille._VECTptr; 2345 if (derive(vv[2],x,contextptr)==0 && derive(vv[2],x,contextptr)==0) 2346 continue; 2347 } 2348 if (!convert_to_euler_mac_laurin(v[i],x,tmp,contextptr)) 2349 return gensizeerr(gettext("Unable to convert sum to Euler Mac-Laurin ")+v[i].print(contextptr)); 2350 v1.push_back(v[i]); 2351 v2.push_back(tmp); 2352 } 2353 if ( ( (pos1=equalposcomp(limit_tab,v[i]._SYMBptr->sommet)) || (pos2=equalposcomp(limit_tractable_functions(),&v[i]._SYMBptr->sommet)) ) ){ 2354 gen g=limit(v[i]._SYMBptr->feuille,x,lim_point,direction,contextptr); 2355 if ( is_inf(g) || (g.type==_VECT && !g._VECTptr->empty() && is_inf(g._VECTptr->back())) ){ 2356 v1.push_back(v[i]); 2357 v2.push_back(pos1?limit_replace[pos1-1](v[i]._SYMBptr->feuille,contextptr):limit_tractable_replace()[pos2-1](v[i]._SYMBptr->feuille,contextptr)); 2358 } 2359 if (is_zero(g)){ 2360 if (v[i]._SYMBptr->sommet==at_Ci){ 2361 v1.push_back(v[i]); 2362 v2.push_back(Ci_replace0(v[i]._SYMBptr->feuille,contextptr)); 2363 } 2364 if (v[i]._SYMBptr->sommet==at_Ei){ 2365 v1.push_back(v[i]); 2366 v2.push_back(Ei_replace0(v[i]._SYMBptr->feuille,contextptr)); 2367 } 2368 } 2369 } 2370 } 2371 } 2372 if (!v1.empty()){ 2373 v2=subst(v2,v1,v2,false,contextptr); 2374 e=remove_lnexp(subst(e,v1,v2,false,contextptr),contextptr); 2375 } 2376 v=rlvarx(e,x); 2377 vs=int(v.size()); 2378 v1.clear(); v2.clear(); 2379 for (int i=0;i<vs;++i){ 2380 #if 1 2381 if (v[i].is_symb_of_sommet(at_sign)){ 2382 gen g=v[i]._SYMBptr->feuille,tmp; 2383 unsigned j=0; 2384 for (;j<5;++j){ 2385 tmp=simplify(limit(g,x,lim_point,direction,contextptr),contextptr); 2386 if (!is_zero(tmp)){ 2387 break; 2388 } 2389 g=derive(g,x,contextptr); 2390 } 2391 if (direction==0 && j!=0) { 2392 continue; 2393 } 2394 if (j!=5){ 2395 // sign( x^j*tmp) 2396 tmp=sign(tmp,contextptr); 2397 if (direction==-1 && j%2) 2398 tmp=-tmp; 2399 v1.push_back(v[i]); 2400 v2.push_back(tmp); 2401 } 2402 } 2403 #else 2404 if (v[i].is_symb_of_sommet(at_sign) && v[i]!=e){ 2405 gen g; 2406 #ifndef NO_STDEXCEPT 2407 try { 2408 #endif 2409 g=limit(v[i],x,lim_point,direction,contextptr); 2410 v1.push_back(v[i]); 2411 v2.push_back(g); 2412 #ifndef NO_STDEXCEPT 2413 } catch (std::runtime_error & ) { 2414 } 2415 #endif 2416 } 2417 #endif 2418 } 2419 e=subst(e,v1,v2,false,contextptr); 2420 return e; 2421 } 2422 is_analytic(const gen & g)2423 bool is_analytic(const gen & g){ 2424 if (g.type==_VECT){ 2425 const_iterateur it=g._VECTptr->begin(),itend=g._VECTptr->end(); 2426 for (;it!=itend;++it){ 2427 if (!is_analytic(*it)) 2428 return false; 2429 } 2430 } 2431 if (g.type!=_SYMB) 2432 return true; 2433 if (equalposcomp(analytic_sommets,g._SYMBptr->sommet)) 2434 return is_analytic(g._SYMBptr->feuille); 2435 return false; 2436 } 2437 unidirectional_limit(const gen & e0,const identificateur & x,const gen & lim_point,int direction,GIAC_CONTEXT)2438 static gen unidirectional_limit(const gen & e0,const identificateur & x,const gen & lim_point,int direction,GIAC_CONTEXT){ 2439 gen e_copy=e0; 2440 // Unidirectional limit, rewrite first (if needed) 2441 if (is_inf(lim_point)){ 2442 if (lim_point==minus_inf) 2443 e_copy=subst(e_copy,x,-x,false,contextptr); 2444 } 2445 else { 2446 if (direction>0) 2447 e_copy=subst(e_copy,x,lim_point+inv(x,contextptr),false,contextptr); 2448 else 2449 e_copy=subst(e_copy,x,lim_point-inv(x,contextptr),false,contextptr); 2450 } 2451 gen coeff,mrv_var,exponent; 2452 sparse_poly1 p; 2453 if (!mrv_lead_term(e_copy,x,coeff,mrv_var,exponent,p,mrv_begin_order,contextptr,false) || is_undef(coeff)){ 2454 gensizeerr("Limit: Max order reached or unable to make series expansion"); 2455 return undef; 2456 } 2457 // check added for limit((tan(x)-x)/x^3,x=inf) 2458 for (unsigned i=0;i<p.size();++i){ 2459 if (check_bounded(p[i].coeff,contextptr)==-1) 2460 return gensizeerr("Limit probably undefined, algorithm unable to handle "+p[i].coeff.print(contextptr)); 2461 } 2462 if (ck_is_strictly_positive(exponent,contextptr)) 2463 return 0; 2464 if (is_zero(exponent)){ 2465 int l=check_bounded(coeff,contextptr); 2466 if (l==-1) 2467 return gensizeerr("Limit probably undefined, algorithm unable to handle "+coeff.print(contextptr)); 2468 return l==1?bounded_function(contextptr):coeff; 2469 } 2470 // check sign of coeff, if coeff depends on x first find equivalent 2471 gen essai=subst(coeff,x,plus_inf,false,contextptr); 2472 if (is_undef(essai) || is_zero(essai) || (essai==unsigned_inf)){ 2473 while (contains(coeff,x)){ 2474 e_copy=coeff; 2475 if (!mrv_lead_term(e_copy,x,coeff,mrv_var,exponent,p,mrv_begin_order,contextptr,false)) 2476 return gensizeerr(contextptr); 2477 } 2478 essai=coeff; 2479 } 2480 gen s=0; 2481 if (calc_mode(contextptr)!=1 || !has_i(p)) // should do it only up to order 0 terms 2482 s=sign(essai,contextptr); 2483 if (s==plus_one) 2484 return plus_inf; 2485 if (s==minus_one) 2486 return minus_inf; 2487 // ? FIXME: if essai=cos(<pi,-1>)+i*sin(<pi,-1>) unsigned_inf better 2488 // limit((-2)^n,n,inf) 2489 int l=check_bounded(essai,contextptr); 2490 if (l==-1) 2491 return gensizeerr("Limit probably undefined, algorithm unable to handle "+essai.print(contextptr)); 2492 return l==1?undef:unsigned_inf; 2493 /* 2494 essai=eval(subst(essai,sincosinf,vecteur(sincosinf.size(),undef))); 2495 if (is_undef(essai)) 2496 return undef; 2497 else 2498 return unsigned_inf; 2499 */ 2500 } 2501 in_limit(const gen & e0,const identificateur & x,const gen & lim_point,int direction,GIAC_CONTEXT)2502 static gen in_limit(const gen & e0,const identificateur & x,const gen & lim_point,int direction,GIAC_CONTEXT){ 2503 if (direction==-2) 2504 return gensizeerr(contextptr); 2505 vecteur vint=lop(rlvarx(e0,x),at_integrate); 2506 for (unsigned i=0;i<vint.size();++i){ 2507 gen f=vint[i]._SYMBptr->feuille; 2508 if (f.type!=_VECT || f._VECTptr->size()!=4) 2509 return gensizeerr(gettext("Undefined integral")); 2510 if ((*f._VECTptr)[1]==x) 2511 return gensizeerr(gettext("Integration variable and limit variable are the same")); 2512 if (!is_zero(derive(f._VECTptr->front(),x,contextptr))) 2513 return gensizeerr(gettext("Integral in limit not implemented yet")); 2514 } 2515 if (e0.type==_VECT){ 2516 const_iterateur it=e0._VECTptr->begin(),itend=e0._VECTptr->end(); 2517 vecteur res; 2518 res.reserve(itend-it); 2519 for (;it!=itend;++it){ 2520 res.push_back(in_limit(*it,x,lim_point,direction,contextptr)); 2521 } 2522 return gen(res,e0.subtype); 2523 } 2524 if (_about(x,contextptr)!=x){ 2525 identificateur xprime(" "+print_INT_(giac_rand(contextptr))); 2526 return in_limit(quotesubst(e0,x,xprime,contextptr),xprime,lim_point,direction,contextptr); 2527 } 2528 gen e=Heavisidetosign(when2sign(piecewise2when(e0,contextptr),contextptr),contextptr); 2529 // Adjust direction for +/- inf limits 2530 if (lim_point==plus_inf) 2531 direction=1; 2532 if (lim_point==minus_inf) 2533 direction=-1; 2534 // First try substitution 2535 if (has_i(lop(e,at_ln))) 2536 e=recursive_normal(expln2trig(e,contextptr),contextptr); 2537 vecteur vsign(loptab(e,sign_floor_ceil_round_tab)); 2538 if (0 && direction && !vsign.empty() && !equalposcomp(sign_floor_ceil_round_tab,e._SYMBptr->sommet)){ 2539 // evaluate vsign first 2540 vecteur res; 2541 for (int i=0;i<int(vsign.size());++i){ 2542 res[i]=in_limit(vsign[i],x,lim_point,direction,contextptr); 2543 } 2544 e=subst(e,vsign,res,false,contextptr); 2545 vsign.clear(); 2546 } 2547 if (vsign.empty()){ 2548 gen first_try=subst(e,x,lim_point,false,contextptr); 2549 first_try=eval(first_try,1,contextptr); 2550 first_try=simplifier(first_try,contextptr); 2551 // if (first_try==plus_inf || first_try==minus_inf) return first_try; 2552 if (!contains(lidnt(first_try),unsigned_inf)){ 2553 if (has_num_coeff(first_try)) 2554 return first_try; 2555 gen chknum; 2556 bool hasnum=has_evalf(first_try,chknum,1,contextptr); 2557 first_try=recursive_ratnormal(first_try,contextptr); 2558 gen chk=recursive_normal(first_try,contextptr); 2559 if (hasnum && !is_undef(chk) && abs(chk-chknum,contextptr)>1e-10 && abs(1-chk/chknum,contextptr)>1e-10){ 2560 chk=undef; 2561 e=_simplify(e,contextptr); 2562 } 2563 if (!is_undef(chk)){ 2564 /* 2565 if (!lop(chk,at_rootof).empty()) 2566 chk=ratnormal(first_try); 2567 */ 2568 if (!is_undef(chk) && !contains(lidnt(chk),unsigned_inf)){ 2569 chk=first_try; 2570 return taille(chk,100)<taille(first_try,100)?chk:first_try; 2571 } 2572 } 2573 } 2574 if (lim_point==unsigned_inf){ 2575 *logptr(contextptr) << gettext("Warning, infinity is unsigned, perhaps you meant +infinity")<< '\n'; 2576 first_try = subst(partfrac(e,false,contextptr),x,lim_point,false,contextptr); 2577 // first_try = subst(ratnormal(e,contextptr),x,lim_point,false,contextptr); 2578 } 2579 else { 2580 //bool b=assume_t_in_ab(x,direction==1?lim_point:lim_point-1,direction==-1?lim_point:lim_point+1,true,true,contextptr); 2581 first_try = quotesubst(partfrac(e,false,contextptr),x,lim_point,contextptr); 2582 // if (b) purgenoassume(x,contextptr); 2583 // first_try = quotesubst(ratnormal(e,contextptr),x,lim_point,contextptr); 2584 } 2585 bool absb=eval_abs(contextptr); 2586 first_try=eval(first_try,eval_level(contextptr),contextptr); // moved before eval_abs(false,contextptr), must be before simplifier below 2587 first_try=simplifier(first_try,contextptr); // for assume(a>0); limit((sqrt(2*a^3*x-x^4)-a*root(3,a^2*x))/(a-root(4,a*x^3)),x=a); 2588 eval_abs(false,contextptr); 2589 first_try = recursive_normal(first_try,contextptr); 2590 //first_try=eval(first_try,1,contextptr); 2591 eval_abs(absb,contextptr); 2592 if (is_undef(first_try) && first_try.type==_STRNG) 2593 return first_try; 2594 if (!is_undef(first_try)){ 2595 // if (!direction) return first_try; 2596 if (first_try!=unsigned_inf) 2597 return first_try; 2598 } 2599 } // end if vsign.empty() 2600 if (has_op(e,*at_surd) || has_op(e,*at_NTHROOT)){ 2601 // FIXME: adjust using limit/direction information 2602 vecteur subst1,subst2; 2603 surd2pow(e,subst1,subst2,contextptr); 2604 gen g=subst(e,subst1,subst2,false,contextptr); 2605 g=limit(g,x,lim_point,direction,contextptr); 2606 return subst(g,subst2,subst1,false,contextptr); 2607 } 2608 e=limit_symbolic_preprocess(e,x,lim_point,direction,contextptr); 2609 if (is_undef(e)) return e; 2610 gen errcode=checkanglemode(contextptr); 2611 if (is_undef(errcode)) 2612 return errcode; 2613 if (e.type!=_SYMB) // e might be an _IDNT equal to x (limit(x*sign(x),x,0,1)) 2614 return subst(e,x,lim_point,false,contextptr); 2615 if (e._SYMBptr->sommet==at_exp) 2616 return exp(in_limit(e._SYMBptr->feuille,x,lim_point,direction,contextptr),contextptr); 2617 if (e._SYMBptr->sommet==at_ln){ 2618 gen tmp=in_limit(e._SYMBptr->feuille,x,lim_point,direction,contextptr); 2619 if (is_undef(tmp)) return tmp; 2620 if (!is_positive(-tmp,contextptr)) 2621 return ln(tmp,contextptr); 2622 } 2623 gen e_copy; 2624 // Rewrite non rational ^ and tan 2625 e_copy=_pow2exp(tan2sincos(exact(e,contextptr),contextptr),contextptr); 2626 // FIXME: this translate exp(i*...) to sin/cos without bugging for 2627 // exp(exp(exp(x)/(1-1/x)))-exp(exp(exp(x)/(1-1/x-exp((-(ln(x)))*ln(ln(x)))))) 2628 if (has_i(e_copy)) { 2629 e_copy=subst(e_copy,tan_tab,tan2sincos_tab,true,contextptr); 2630 e_copy=subst(e_copy,exp_tab,exp2sincos_tab,true,contextptr); 2631 if (has_i(lop(e_copy,at_erfs))) 2632 return gensizeerr(gettext("erf/erfc/erfs with complex argument not yet implemented in limit")); 2633 } 2634 // Rewrite constants 2635 vecteur rv=rlvar(e_copy,false),cv; 2636 for (int i=0;i<rv.size();++i){ 2637 if (evalf(rv[i],1,contextptr).type<_CPLX) 2638 cv.push_back(rv[i]); 2639 } 2640 if (!cv.empty()){ 2641 gen cvg=tsimplify(cv,contextptr); 2642 if (cvg.type==_VECT && cvg._VECTptr->size()==cv.size()) 2643 e_copy=subst(e_copy,cv,*cvg._VECTptr,false,contextptr); 2644 } 2645 if (!direction) { 2646 if (!is_analytic(e_copy) || has_op(e_copy,*at_acos) || has_op(e_copy,*at_asin)){ 2647 gen g1=unidirectional_limit(e_copy,x,lim_point,1,contextptr); 2648 if (is_undef(g1)) 2649 return g1; 2650 gen g2=unidirectional_limit(e_copy,x,lim_point,-1,contextptr); 2651 if (is_undef(g2)) 2652 return g2; 2653 if (g1==g2 || is_zero(ratnormal(g1-g2,contextptr))) 2654 return g1; 2655 return gensizeerr("Unidirectional limits are distinct "+g2.print(contextptr)+","+g1.print(contextptr)); 2656 } 2657 // supposed to be analytic, try first series expansion 2658 sparse_poly1 p; 2659 p.push_back(monome(undef,0)); 2660 double ordre=mrv_begin_order; 2661 for ( ; !p.empty() && is_undef(p.front().coeff) && (ordre<mrv_begin_order*4);ordre=1.5*ordre+1) { 2662 p=series__SPOL1(e_copy,x,lim_point,int(ordre),0,contextptr); 2663 if (!p.empty() && !is_undef(p.front().coeff)){ 2664 if (p.front().coeff.type==_FRAC && is_strictly_positive(-p.front().coeff._FRACptr->den,contextptr)) 2665 p.front().coeff=fraction(-p.front().coeff._FRACptr->num,-p.front().coeff._FRACptr->den); 2666 break; 2667 } 2668 } 2669 // COUT << p << '\n'; 2670 if (ordre>=mrv_begin_order*4){ 2671 gen g1=unidirectional_limit(e_copy,x,lim_point,1,contextptr); 2672 if (is_undef(g1)) 2673 return g1; 2674 gen g2=unidirectional_limit(e_copy,x,lim_point,-1,contextptr); 2675 if (is_undef(g2)) 2676 return g2; 2677 if (is_zero(ratnormal(g1-g2,contextptr))) 2678 return g1; 2679 return gensizeerr("Unidirectional limits are distincts "+g2.print(contextptr)+","+g1.print(contextptr)); 2680 } 2681 if (p.empty() || ck_is_strictly_positive(p.front().exponent,contextptr) ){ 2682 if (check_bounded(p.front().coeff,contextptr)==-1) 2683 return gensizeerr("Unable to bound coefficient"); 2684 return 0; 2685 } 2686 if (ck_is_strictly_positive(-p.front().exponent,contextptr)){ 2687 if (p.front().exponent.type==_INT_ && p.front().exponent.val%2==0){ 2688 int s=fastsign(p.front().coeff,contextptr); 2689 if (s==1) 2690 return plus_inf; 2691 if (s==-1) 2692 return minus_inf; 2693 } 2694 return unsigned_inf; 2695 } 2696 if (is_zero(p.front().exponent)){ 2697 if (contains(p.front().coeff,x)){ 2698 return gensizeerr(gettext("Try unidirectional series")); 2699 } 2700 int l=check_bounded(p.front().coeff,contextptr); 2701 if (l==-1) 2702 return gensizeerr("Limit probably undefined, algorithm unable to handle "+p.front().coeff.print(contextptr)); 2703 return l==1?bounded_function(contextptr):p.front().coeff; 2704 } 2705 return gensizeerr(gettext("Series internal bug")); 2706 } 2707 return unidirectional_limit(e_copy,x,lim_point,direction,contextptr); 2708 } 2709 2710 // return plus_inf if a > b (at x=+infinity), !0 if a#b, 0 if a < b mrv_compare(const gen & a,const gen & b,const identificateur & x,GIAC_CONTEXT)2711 static gen mrv_compare(const gen & a,const gen & b,const identificateur & x,GIAC_CONTEXT){ 2712 if ((a.type!=_SYMB) && (b.type!=_SYMB)) 2713 return 1; 2714 gen lna,lnb,l; 2715 if ((a.type==_SYMB) && (a._SYMBptr->sommet==at_exp)) 2716 lna=a._SYMBptr->feuille; 2717 else 2718 lna=ln(a,contextptr); 2719 if ((b.type==_SYMB) && (b._SYMBptr->sommet==at_exp)) 2720 lnb=b._SYMBptr->feuille; 2721 else 2722 lnb=ln(b,contextptr); 2723 gen coeff,mrv_var,exponent; 2724 sparse_poly1 p; 2725 if (!mrv_lead_term(rdiv(lna,lnb,contextptr),x,coeff,mrv_var,exponent,p,mrv_begin_order,contextptr,false)) 2726 return gensizeerr(contextptr); 2727 if (ck_is_strictly_positive(exponent,contextptr)) 2728 return 0; 2729 if (is_zero(exponent)) 2730 return coeff; 2731 return plus_inf; 2732 } 2733 mrv_max(const vecteur & a_faster_var,const vecteur & a_coeff_ln,const vecteur & a_slower_var,const vecteur & b_faster_var,const vecteur & b_coeff_ln,const vecteur & b_slower_var,const identificateur & x,vecteur & faster_var,vecteur & coeff_ln,vecteur & slower_var,GIAC_CONTEXT)2734 static bool mrv_max(const vecteur & a_faster_var, const vecteur & a_coeff_ln, const vecteur & a_slower_var, const vecteur & b_faster_var,const vecteur & b_coeff_ln, const vecteur & b_slower_var,const identificateur & x, vecteur & faster_var, vecteur & coeff_ln, vecteur & slower_var,GIAC_CONTEXT){ 2735 int pos_a,pos_b; 2736 gen s; 2737 if (intersect(a_faster_var,b_faster_var,pos_a,pos_b)) 2738 s=normal(rdiv(a_faster_var[pos_a],b_faster_var[pos_b],contextptr),contextptr); 2739 else { 2740 if (a_faster_var.empty() || intersect(a_faster_var,b_slower_var,pos_a,pos_b)) 2741 s=0; 2742 else { 2743 if (b_faster_var.empty() || intersect(b_faster_var,a_slower_var,pos_a,pos_b) ) 2744 s=plus_inf; 2745 else 2746 s=mrv_compare(a_faster_var.front(),b_faster_var.front(),x,contextptr); 2747 } 2748 if (is_undef(s)) 2749 return false; 2750 if (s==plus_inf){ 2751 slower_var=mergevecteur(a_slower_var,b_slower_var); 2752 slower_var=mergevecteur(b_faster_var,slower_var); 2753 faster_var=a_faster_var; 2754 coeff_ln=a_coeff_ln; 2755 return true; 2756 } 2757 if (is_zero(s)){ 2758 slower_var=mergevecteur(a_slower_var,b_slower_var); 2759 slower_var=mergevecteur(a_faster_var,slower_var); 2760 faster_var=b_faster_var; 2761 coeff_ln=b_coeff_ln; 2762 return true; 2763 } 2764 } 2765 // s!=0 && s!=plus_inf 2766 // size test used to get w or inv(w) at the front() 2767 if (a_faster_var.front().symb_size()>b_faster_var.front().symb_size()){ 2768 coeff_ln=mergevecteur(b_coeff_ln,multvecteur(s,a_coeff_ln)); 2769 faster_var=mergevecteur(b_faster_var,a_faster_var); 2770 slower_var=mergevecteur(b_slower_var,a_slower_var); 2771 } 2772 else { 2773 coeff_ln=mergevecteur(a_coeff_ln,multvecteur(inv(s,contextptr),b_coeff_ln)); 2774 faster_var=mergevecteur(a_faster_var,b_faster_var); 2775 slower_var=mergevecteur(a_slower_var,b_slower_var); 2776 } 2777 return true; 2778 } 2779 2780 // Find most rapidly varying subexpression of e and res mrv(const gen & e,const identificateur & x,vecteur & faster_var,vecteur & coeff_ln,vecteur & slower_var,GIAC_CONTEXT)2781 static bool mrv(const gen & e,const identificateur & x,vecteur & faster_var,vecteur & coeff_ln, vecteur & slower_var,GIAC_CONTEXT){ 2782 // Find all var of e depending on x 2783 vecteur v0(lvarx(e,x)); 2784 // Find mrv of these vars 2785 vecteur::const_iterator it=v0.begin(),itend=v0.end(); 2786 for (;it!=itend;++it){ 2787 if (equalposcomp(faster_var,*it) || equalposcomp(slower_var,*it)) 2788 continue; 2789 if (it->type!=_SYMB){ 2790 if (!mrv_max(faster_var,coeff_ln,slower_var, 2791 vecteur(1,*it),vecteur(1,plus_one),vecteur(0), 2792 x,faster_var,coeff_ln,slower_var,contextptr)) 2793 return false; 2794 continue; 2795 } 2796 gen temp=*it; 2797 if (temp._SYMBptr->sommet==at_ln){ 2798 if (!mrv(temp._SYMBptr->feuille,x,faster_var,coeff_ln,slower_var,contextptr)) 2799 return false; 2800 continue; 2801 } 2802 if (temp._SYMBptr->sommet==at_Psi && temp._SYMBptr->feuille.type==_VECT){ 2803 if (!mrv(temp._SYMBptr->feuille[0],x,faster_var,coeff_ln,slower_var,contextptr)) 2804 return false; 2805 continue; 2806 } 2807 if (temp._SYMBptr->sommet==at_lower_incomplete_gamma || temp._SYMBptr->sommet==at_igamma_exp){ 2808 gen & f = temp._SYMBptr->feuille; 2809 if (depend(f[0],x)) 2810 return false; 2811 if (!mrv(f[1],x,faster_var,coeff_ln,slower_var,contextptr)) 2812 return false; 2813 continue; 2814 } 2815 if (temp._SYMBptr->sommet==at_pow) 2816 temp=symbolic(at_exp,(*(temp._SYMBptr->feuille._VECTptr))[1]*ln((*(temp._SYMBptr->feuille._VECTptr))[0],contextptr)); 2817 if (temp._SYMBptr->sommet==at_euler_mac_laurin){ 2818 gen & f = temp._SYMBptr->feuille; 2819 if (f.type==_VECT && f._VECTptr->size()==5){ 2820 if (!mrv(f[0],x,faster_var,coeff_ln,slower_var,contextptr) || 2821 !mrv(f[3],x,faster_var,coeff_ln,slower_var,contextptr) || 2822 !mrv(f[4],x,faster_var,coeff_ln,slower_var,contextptr)) 2823 return false; 2824 continue; 2825 } 2826 } 2827 if (temp._SYMBptr->sommet==at_integrate){ 2828 gen & f = temp._SYMBptr->feuille; 2829 if (f.type==_VECT && f._VECTptr->size()==4){ 2830 if (!mrv(f[0],x,faster_var,coeff_ln,slower_var,contextptr) || 2831 !mrv(f[2],x,faster_var,coeff_ln,slower_var,contextptr) || 2832 !mrv(f[3],x,faster_var,coeff_ln,slower_var,contextptr)) 2833 return false; 2834 continue; 2835 } 2836 } 2837 if (temp._SYMBptr->feuille.type==_VECT){ 2838 *logptr(contextptr) << gettext("Limit probably undefined, algorithm unable to handle ")+temp.print(contextptr) << '\n'; 2839 return false; 2840 } 2841 gen l=in_limit(temp._SYMBptr->feuille,x,plus_inf,0,contextptr); 2842 if (is_undef(l) || (l==unsigned_inf && temp._SYMBptr->sommet!=at_cos && temp._SYMBptr->sommet!=at_sin && temp._SYMBptr->sommet!=at_erfs)){ 2843 *logptr(contextptr) << gettext("Undef/Unsigned Inf encountered in limit") << '\n'; 2844 return false; 2845 } 2846 if (!is_inf(l)){ 2847 if (!mrv(temp._SYMBptr->feuille,x,faster_var,coeff_ln,slower_var,contextptr)) 2848 return false; 2849 continue; 2850 } 2851 if (temp._SYMBptr->sommet==at_exp){ 2852 if (!mrv(temp._SYMBptr->feuille,x,faster_var,coeff_ln,slower_var,contextptr) || 2853 !mrv_max(faster_var,coeff_ln,slower_var, 2854 vecteur(1,temp),vecteur(1,plus_one),vecteur(0), 2855 x,faster_var,coeff_ln,slower_var,contextptr)) 2856 return false; 2857 continue; 2858 } 2859 // (semi-)tractable functions? 2860 gen shift_coeff; 2861 if (!temp._SYMBptr->sommet.ptr()->series_expansion){ 2862 invalidserieserr(string(gettext("no taylor method for "))+temp._SYMBptr->sommet.ptr()->print(contextptr)); 2863 return false; 2864 } 2865 gen test=temp._SYMBptr->sommet.ptr()->series_expansion(l,0,temp._SYMBptr->sommet,0,shift_coeff,contextptr); // fixme: 0 should be image_of_direction 2866 if (!is_undef(test)){ 2867 if (!mrv(temp._SYMBptr->feuille,x,faster_var,coeff_ln,slower_var,contextptr)) 2868 return false; 2869 continue; 2870 } 2871 // fixme: add support for integrals and sums 2872 if (temp._SYMBptr->sommet!=at_exp) 2873 return false; // setsizeerr(gettext("series.cc/mrv")); 2874 } 2875 return true; 2876 } 2877 pnormal(sparse_poly1 & v,GIAC_CONTEXT)2878 bool pnormal(sparse_poly1 & v,GIAC_CONTEXT){ 2879 sparse_poly1::const_iterator it=v.begin(),itend=v.end(); 2880 sparse_poly1 p; 2881 gen e; 2882 for (;it!=itend;++it){ 2883 e=recursive_normal(it->coeff,contextptr); 2884 if (!is_zero(e)) 2885 p.push_back(monome(e,it->exponent)); 2886 } 2887 swap(p,v); 2888 return true; 2889 } 2890 2891 // find asymptotic equivalent of e in terms of the mrv var of e mrv_lead_term(const gen & e,const identificateur & x,gen & coeff,gen & mrv_var,gen & exponent,sparse_poly1 & q,int begin_ordre,GIAC_CONTEXT,bool series)2892 static bool mrv_lead_term(const gen & e,const identificateur & x,gen & coeff, gen & mrv_var, gen & exponent,sparse_poly1 & q,int begin_ordre,GIAC_CONTEXT,bool series){ 2893 if (!contains(e,x)){ 2894 coeff=ratnormal(e,contextptr); 2895 mrv_var=x; 2896 exponent=0; 2897 q.clear(); 2898 q.push_back(monome(coeff,0)); 2899 return true; 2900 } 2901 vecteur faster_var,coeff_ln,slower_var; 2902 gen ecopy(e); 2903 if (!mrv(ecopy,x,faster_var,coeff_ln,slower_var,contextptr)) 2904 return false; 2905 int rescaling=0; // number of ln(x) -> x substitution 2906 // if the mrv contains x, we must change scale ln(x) -> x and x -> e^x 2907 for (;equalposcomp(faster_var,x);++rescaling){ 2908 upscale(ecopy,x,contextptr); 2909 // next test needed to avoid infinite recursion 2910 // if there is an e^x inside the new mrv is m rescaled 2911 if (contains(ecopy,exp(x,contextptr))){ 2912 gen temp(faster_var); 2913 upscale(temp,x,contextptr); 2914 faster_var=*temp._VECTptr; 2915 } 2916 else { 2917 faster_var.clear(); 2918 slower_var.clear(); 2919 coeff_ln.clear(); 2920 if (!mrv(ecopy,x,faster_var,coeff_ln,slower_var,contextptr)) 2921 return false; 2922 } 2923 } 2924 if (faster_var.empty()){ 2925 coeff=ratnormal(ecopy,contextptr); 2926 mrv_var=x; 2927 exponent=0; 2928 q.clear(); 2929 q.push_back(monome(coeff,0)); 2930 return true; 2931 } 2932 // now find w, the mrv element -> 0 and express other elements of the mrv 2933 // algo: sort faster_var by symb_size 2934 // w is the shortest one, set g=ln(w) 2935 // replace the next shortest at position pos 2936 // by w^coeff_ln[pos]* exp(ln(faster_var[pos])-coeff_ln[pos]*g) 2937 // then replace faster_var[0] by w above 2938 // go on, the replace operation should be done with 2939 // previous ordered_faster by their expression in terms of w 2940 // At the end replace w by 1/w if w -> plus_inf 2941 bool dont_invert=is_zero(in_limit(faster_var.front(),x,plus_inf,0,contextptr)); 2942 vecteur faster_var_tmp(faster_var); 2943 stable_sort(faster_var.begin(),faster_var.end(),symb_size_less_t()); 2944 identificateur w(" w"); 2945 vecteur faster_var_subst(1,w); 2946 gen g=faster_var.front()._SYMBptr->feuille; 2947 iterateur it=faster_var.begin()+1,itend=faster_var.end(); 2948 gen c,f; 2949 for (;it!=itend;++it){ 2950 int p=equalposcomp(faster_var_tmp,*it); 2951 // assert(p); 2952 c=coeff_ln[p-1]; 2953 f=subst(it->_SYMBptr->feuille,faster_var,faster_var_subst,false,contextptr); 2954 faster_var_subst.push_back(pow(w,c,contextptr)*exp(normal(f-c*g,contextptr),contextptr)); 2955 } 2956 // subst in original expression and make the asymptotic expansion 2957 double ordre=begin_ordre; 2958 f=subst(ecopy,faster_var,faster_var_subst,false,contextptr); 2959 if (!dont_invert) 2960 f=subst(f,w,inv(w,contextptr),false,contextptr); 2961 if (faster_var.front().is_symb_of_sommet(at_exp)){ 2962 // replace ln(exp(g)^k*...) by k*g+ln(...) 2963 vecteur lf(lop(f,at_ln)),lf1,lf2; 2964 iterateur it=lf.begin(),itend=lf.end(); 2965 for (;it!=itend;++it){ 2966 gen argln=it->_SYMBptr->feuille; 2967 sparse_poly1 p=series__SPOL1(argln,w,0,int(ordre),1,contextptr); 2968 if (!p.empty() && !is_undef(p.front().coeff) ){ 2969 if (!is_zero(p.front().exponent)) 2970 argln=argln*symbolic(at_pow,gen(makevecteur(w,-p.front().exponent),_SEQ__VECT)); 2971 lf1.push_back(*it); 2972 lf2.push_back(p.front().exponent*(dont_invert?g:-g)+symbolic(at_ln,argln)); 2973 } 2974 } 2975 if (!lf1.empty()) 2976 f=subst(f,lf1,lf2,false,contextptr); 2977 } 2978 sparse_poly1 p; 2979 p.push_back(monome(undef,0)); 2980 // FIXME: if ordre>max_series it might return a wrong answer here 2981 for (unsigned count=0; (ordre<max_series_expansion_order) && !p.empty() && 2982 (p.size()>=1) && is_undef(p.front().coeff) ;++count,ordre=ordre*1.5+1){ 2983 bool inv=false; 2984 p=series__SPOL1(f,w,0,int(ordre),1,contextptr); 2985 // if (count==2 && p.size()==1 && is_undef(p.front().coeff)){ f=ratnormal(f,contextptr); p=series__SPOL1(f,w,0,int(ordre),1,contextptr); } 2986 #ifdef TIMEOUT 2987 control_c(); 2988 #endif 2989 if (ctrl_c || interrupted || (!p.empty() && is_undef(p.front().exponent))) 2990 return false; 2991 if (!p.empty() && !is_undef(p.front().coeff) ){ 2992 // substitution of ln(w) by +-g should not be useful anymore 2993 gen tmp=ratnormal(subst(p.front().coeff,ln(w,contextptr),(dont_invert?g:-g),false,contextptr),contextptr); 2994 if (is_undef(tmp) ){ 2995 inv=true; 2996 p=spdiv(sparse_poly1(1,monome(1,0)),p,contextptr); 2997 if (is_undef(p)) 2998 return false; 2999 pnormal(p,contextptr); 3000 } 3001 } 3002 // cerr << p << '\n'; 3003 // replace ln( w) in coeff by g or -g 3004 if (dont_invert) 3005 p=subst(p,ln(w,contextptr),g,false,contextptr); 3006 else 3007 p=subst(p,ln(w,contextptr),-g,false,contextptr); 3008 if (inv){ 3009 p=spdiv(sparse_poly1(1,monome(1,0)),p,contextptr); 3010 if (is_undef(p)) 3011 return false; 3012 pnormal(p,contextptr); 3013 } 3014 } 3015 if (!p.empty()) 3016 p.front().exponent=simplify(p.front().exponent,contextptr); 3017 q=p; 3018 if (!p.empty()){ 3019 bool done=false; 3020 // check for exponent 0 at front() 3021 if (is_zero(p.front().exponent) && contains(p.front().coeff,x) && !is_zero(derive(p.front().coeff,x,contextptr))){ 3022 // if p is uniquely composed of this coeff, expand 3023 if (!mrv_lead_term(p.front().coeff,x,coeff,mrv_var,exponent,p,mrv_begin_order,contextptr,false)) 3024 return false; 3025 if (p.empty() 3026 // test below allow expanding series(ln(x+1),x=inf) 3027 // the boolean series was added 3028 // otherwise testlimit would not work correctly 3029 || (series && p.size()==1 && q.size()>1 &&!is_undef(p.front().coeff)) 3030 ){ 3031 done=false; 3032 p=q; 3033 } 3034 else { 3035 if (q.size()>1 && !is_undef(p.back().coeff)) 3036 p.push_back(monome(undef,p.back().exponent+begin_ordre)); 3037 q=p; 3038 done=true; 3039 } 3040 } 3041 if (!done) { 3042 if (dont_invert) 3043 mrv_var=faster_var.front(); 3044 else 3045 mrv_var=inv(faster_var.front(),contextptr); 3046 coeff = p.front().coeff; 3047 exponent = p.front().exponent; 3048 } 3049 } 3050 // mrv_var is w as function of x, rescaled using the int variable rescaling 3051 // coeff must be rescaled as well 3052 sparse_poly1::iterator it0=q.begin(),it1=q.end(); 3053 for (;rescaling;--rescaling){ 3054 downscale(mrv_var,x,contextptr); 3055 downscale(coeff,x,contextptr); 3056 for (sparse_poly1::iterator it=it0;it!=it1;++it) 3057 downscale(it->coeff,x,contextptr); 3058 } 3059 return true; 3060 } 3061 intersect(const vecteur & a,const vecteur & b,int & pos_a,int & pos_b)3062 bool intersect(const vecteur & a,const vecteur &b,int & pos_a,int & pos_b){ 3063 vecteur res; 3064 if (a.empty() || b.empty()) 3065 return false; 3066 vecteur::const_iterator it=a.begin(),itend=a.end(); 3067 for (;it!=itend;++it) 3068 pos_b=equalposcomp(b,*it); 3069 if (pos_b){ 3070 --pos_b; 3071 pos_a=int(it-a.begin()); 3072 return true; 3073 } 3074 return false; 3075 } 3076 convert_to_direction(const gen & l)3077 static int convert_to_direction(const gen & l){ 3078 if (is_one(l) || l==at_plus) 3079 return 1; 3080 if (is_minus_one(l) || l==at_binary_minus || l==at_neg) 3081 return -1; 3082 if (is_zero(l)) 3083 return 0; 3084 return -2; 3085 } 3086 3087 // Main limit entry point limit(const gen & e,const identificateur & x,const gen & lim_point_,int direction,GIAC_CONTEXT)3088 gen limit(const gen & e,const identificateur & x,const gen & lim_point_,int direction,GIAC_CONTEXT){ 3089 gen lim_point(eval(lim_point_,1,contextptr)); 3090 if (is_undef(lim_point)) 3091 return lim_point; 3092 if (lim_point.type==_DOUBLE_){ 3093 double d=lim_point._DOUBLE_val; 3094 if (d==int(d)) 3095 return limit(e,x,int(d),direction,contextptr); 3096 } 3097 if (has_num_coeff(lim_point)) // otherwise A:=conic(x^2/4+y^2/3=1); B:=element(A,1)+nop(-1.666-1.243*i) runs forever 3098 return subst(e,x,lim_point,false,contextptr); 3099 // Insert here code for cleaning limit remember 3100 int save_series_flags=series_flags(contextptr); 3101 series_flags(save_series_flags | 8,contextptr); 3102 // sincosinf.clear(); 3103 gen e_exact=exact(e,contextptr); 3104 gen lim_point_exact=exact(lim_point,contextptr); 3105 gen l=in_limit(e_exact,x,lim_point_exact,direction,contextptr); 3106 if (e.is_approx() || lim_point.is_approx()) 3107 l=evalf(l,1,contextptr); 3108 series_flags(save_series_flags,contextptr); 3109 // vecteur sincosinfsub(sincosinf.size(),undef); 3110 // l=eval(subst(l,sincosinf,sincosinfsub)); 3111 return l; 3112 } 3113 quotedlimit(const gen & e,const identificateur & x,const gen & lim_point,int direction,GIAC_CONTEXT)3114 gen quotedlimit(const gen & e,const identificateur & x,const gen & lim_point,int direction,GIAC_CONTEXT){ 3115 vecteur v(1,exact(e,contextptr)); 3116 v=quote_eval(v,vecteur(1,x),contextptr); 3117 return limit(v[0],x,lim_point,direction,contextptr); 3118 } 3119 3120 // "unary" version 3121 static const char _limit_s []="limit"; _limit(const gen & args,GIAC_CONTEXT)3122 gen _limit(const gen & args,GIAC_CONTEXT){ 3123 if ( args.type==_STRNG && args.subtype==-1) return args; 3124 if (args.type!=_VECT) 3125 return quotedlimit(args,*vx_var._IDNTptr,0,0,contextptr); 3126 vecteur v =*args._VECTptr; 3127 int s=int(v.size()); 3128 if (!s) 3129 toofewargs(_limit_s); 3130 gen G=v[0]; 3131 if (s==1) 3132 return quotedlimit(G,*vx_var._IDNTptr,0,0,contextptr); 3133 gen e=v[1]; 3134 if (s==2){ 3135 if (calc_mode(contextptr)==1 && !v[1].is_symb_of_sommet(at_equal)) 3136 return _limit(makesequence(G,ggb_var(G),e),contextptr); 3137 if (e.type==_IDNT) 3138 return quotedlimit(G,*e._IDNTptr,0,0,contextptr); 3139 if (e.type!=_SYMB) 3140 return gentypeerr(contextptr); 3141 if (!is_equal(e)) 3142 return gensizeerr(contextptr); 3143 gen x=(*(e._SYMBptr->feuille._VECTptr))[0]; 3144 if (x.type!=_IDNT) 3145 return gensizeerr(contextptr); 3146 return quotedlimit(G,*x._IDNTptr,(*(e._SYMBptr->feuille._VECTptr))[1],0,contextptr); 3147 } 3148 if (s>2) 3149 v[2]=eval(v[2],1,contextptr); 3150 if (s>3) 3151 v[3]=eval(v[3],1,contextptr); 3152 if (s==3){ 3153 gen arg3=v[2]; 3154 if (e.type==_IDNT) 3155 return quotedlimit(G,*e._IDNTptr,arg3,0,contextptr); 3156 if (e.type!=_SYMB){ 3157 if (is_one(arg3)||is_minus_one(arg3)) 3158 return quotedlimit(G,*ggb_var(G)._IDNTptr,e,int(evalf_double(arg3,1,contextptr)._DOUBLE_val),contextptr); 3159 return gentypeerr(contextptr); 3160 } 3161 if (!is_equal(e)){ 3162 if (is_one(arg3)||is_minus_one(arg3)) 3163 return quotedlimit(G,*ggb_var(G)._IDNTptr,e,int(evalf_double(arg3,1,contextptr)._DOUBLE_val),contextptr); 3164 return gensizeerr(contextptr); 3165 } 3166 gen x=(*(e._SYMBptr->feuille._VECTptr))[0]; 3167 if (x.type!=_IDNT) 3168 return gensizeerr(contextptr); 3169 return quotedlimit(G,*x._IDNTptr,(*(e._SYMBptr->feuille._VECTptr))[1],convert_to_direction(v[2]),contextptr); 3170 } 3171 if (s>4) 3172 return gentoomanyargs(_limit_s); 3173 if (e.type!=_IDNT) 3174 return gentypeerr(contextptr); 3175 return quotedlimit(G,*e._IDNTptr,v[2],convert_to_direction(v[3]),contextptr); 3176 } texprintaslimit(const gen & g,const char * orig_s,GIAC_CONTEXT)3177 static string texprintaslimit(const gen & g,const char * orig_s,GIAC_CONTEXT){ 3178 string s("\\lim "); 3179 if (g.type!=_VECT) 3180 return s+gen2tex(g,contextptr); 3181 vecteur v(*g._VECTptr); 3182 int l(int(v.size())); 3183 if (!l) 3184 return s; 3185 if (l==1) 3186 return s+gen2tex(v[0],contextptr); 3187 if (l==2) 3188 return s+"_{"+gen2tex(v[1],contextptr)+"}"+gen2tex(v[0],contextptr); 3189 // directional limit 3190 if (l==3){ 3191 if (is_one(v[2])) 3192 return s+"_{"+gen2tex(v[1],contextptr)+"^+}"+gen2tex(v[0],contextptr); 3193 if (is_minus_one(v[2])) 3194 return s+"_{"+gen2tex(v[1],contextptr)+"^-}"+gen2tex(v[0],contextptr); 3195 else return s+"_{"+gen2tex(v[1],contextptr)+"}"+gen2tex(v[0],contextptr); 3196 } 3197 return s; 3198 } 3199 static define_unary_function_eval4 (__limit,&_limit,_limit_s,0,&texprintaslimit); 3200 define_unary_function_ptr5( at_limit ,alias_at_limit,&__limit,_QUOTE_ARGUMENTS,true); 3201 3202 #if 0 3203 static const char _lim_s []="lim"; 3204 static define_unary_function_eval4 (__lim,&_limit,_lim_s,0,&texprintaslimit); 3205 define_unary_function_ptr5( at_lim ,alias_at_lim,&__lim,_QUOTE_ARGUMENTS,true); 3206 #endif 3207 3208 // like sparse_poly12gen, but if there is only 1 term and no remainder 3209 // expand it, l.1976 sparse_poly12gen_expand(const sparse_poly1 & s,const identificateur & x,const gen & mrv_var,int ordre,gen & remains,bool with_order_size,GIAC_CONTEXT)3210 static gen sparse_poly12gen_expand(const sparse_poly1 & s,const identificateur & x,const gen & mrv_var,int ordre,gen & remains,bool with_order_size,GIAC_CONTEXT){ 3211 if (s.size()!=1 || !contains(s.front().coeff,x) ) 3212 return sparse_poly12gen(s,mrv_var,remains,with_order_size); 3213 gen afaire=s.front().coeff,mrv_fait(mrv_var),a,b,c,exponent(s.front().exponent); 3214 if (mrv_fait.is_symb_of_sommet(at_inv) && mrv_fait._SYMBptr->feuille.is_symb_of_sommet(at_exp)) 3215 mrv_fait=symbolic(at_exp,symbolic(at_neg,mrv_fait._SYMBptr->feuille._SYMBptr->feuille)); 3216 if (mrv_fait.is_symb_of_sommet(at_exp)){ 3217 // search embedded ln if mrv_var is an exp var 3218 mrv_fait=mrv_fait._SYMBptr->feuille; 3219 vecteur l(lop(mrv_fait,at_ln)); 3220 int ls=int(l.size()); 3221 for (int i=0;i<ls;++i){ 3222 identificateur tx(" x"); 3223 gen tmpx(tx); 3224 gen mrv_temp=subst(mrv_fait,l[i],tmpx,true,contextptr); 3225 if (is_linear_wrt(mrv_temp,tmpx,a,b,contextptr)){ 3226 // extract constant part of a in c using decompose_plus 3227 if (a.is_symb_of_sommet(at_plus) && a._SYMBptr->feuille.type==_VECT){ 3228 c=0; 3229 vecteur non_constant; 3230 decompose_plus(*a._SYMBptr->feuille._VECTptr,x,non_constant,c,contextptr); 3231 a=_plus(non_constant,contextptr); 3232 afaire=afaire*pow(l[i]._SYMBptr->feuille,c*exponent,contextptr); 3233 mrv_fait=a*l[i]+b; 3234 } 3235 } 3236 } 3237 mrv_fait=symbolic(at_exp,mrv_fait); 3238 } 3239 gen coeff2,mrv_var2,exponent2; 3240 sparse_poly1 s2; 3241 if (!mrv_lead_term(afaire,x,coeff2,mrv_var2,exponent2,s2,ordre,contextptr,true)) 3242 return false; 3243 return sparse_poly12gen_expand(s2,x,mrv_var2,ordre,remains,with_order_size,contextptr)*pow(mrv_fait,exponent,contextptr); 3244 } 3245 in_series(const gen & e0,const identificateur & x,const gen & lim_point,int ordre,int direction,GIAC_CONTEXT)3246 static gen in_series(const gen & e0,const identificateur & x,const gen & lim_point,int ordre,int direction,GIAC_CONTEXT){ 3247 gen e=limit_symbolic_preprocess(e0,x,lim_point,direction,contextptr); 3248 if (is_undef(e)) return e; 3249 gen errcode=checkanglemode(contextptr); 3250 if (is_undef(errcode)) return errcode; 3251 if (lim_point==plus_inf){ 3252 gen coeff,mrv_var,exponent,remains; 3253 sparse_poly1 s; 3254 if (!mrv_lead_term(e,x,coeff,mrv_var,exponent,s,ordre,contextptr,true)){ 3255 #ifdef EMCC 3256 return undef; 3257 #else 3258 return gensizeerr(contextptr); 3259 #endif 3260 } 3261 if (series_flags(contextptr) & (1<<4) ) 3262 return s; // no back conversion if bit 4 is set 3263 return sparse_poly12gen_expand(s,x,mrv_var,ordre,remains,true,contextptr); 3264 } 3265 if (lim_point==minus_inf){ 3266 gen coeff,mrv_var,exponent,remains; 3267 sparse_poly1 s; 3268 if (!mrv_lead_term(subst(e,x,-x,false,contextptr),x,coeff,mrv_var,exponent,s,ordre,contextptr,true)) 3269 return gensizeerr(contextptr); 3270 if (series_flags(contextptr) & (1<<4) ) 3271 return s; // no back conversion if bit 4 is set 3272 return subst(sparse_poly12gen_expand(s,x,mrv_var,ordre,remains,true,contextptr),x,-x,false,contextptr); 3273 } 3274 if (direction==1){ 3275 gen ecopy=subst(e,x,lim_point+inv(x,contextptr),false,contextptr); 3276 gen coeff,mrv_var,exponent,remains; 3277 sparse_poly1 s; 3278 if (!mrv_lead_term(ecopy,x,coeff,mrv_var,exponent,s,ordre,contextptr,true)) 3279 return gensizeerr(contextptr); 3280 if (series_flags(contextptr) & (1<<4) ) 3281 return s; // no back conversion if bit 4 is set 3282 return subst(sparse_poly12gen_expand(s,x,mrv_var,ordre,remains,true,contextptr),x,inv(x-lim_point,contextptr),false,contextptr); 3283 } 3284 if (direction==-1){ 3285 gen ecopy=subst(e,x,lim_point-inv(x,contextptr),false,contextptr); 3286 gen coeff,mrv_var,exponent,remains; 3287 sparse_poly1 s; 3288 if (!mrv_lead_term(ecopy,x,coeff,mrv_var,exponent,s,ordre,contextptr,true)) 3289 return gensizeerr(contextptr); 3290 if (series_flags(contextptr) & (1<<4) ) 3291 return s; // no back conversion if bit 4 is set 3292 return subst(sparse_poly12gen_expand(s,x,mrv_var,ordre,remains,true,contextptr),x,inv(lim_point-x,contextptr),false,contextptr); 3293 } 3294 gen remains; 3295 switch (e.type){ 3296 case _INT_: case _ZINT: case _DOUBLE_: case _CPLX: 3297 return e; 3298 case _IDNT: 3299 if ( !ordre && (*e._IDNTptr==x) ) 3300 return lim_point+(x-lim_point)*order_size(x-lim_point,contextptr); 3301 else 3302 return e; 3303 case _SYMB: { 3304 sparse_poly1 s; 3305 s=ck_series__SPOL1(e,x,lim_point,ordre,direction,contextptr); 3306 if (series_flags(contextptr) & (1<<4) ) 3307 return s; // no back conversion if bit 4 is set 3308 return sparse_poly12gen(s,x-lim_point,remains,true); 3309 } 3310 default: 3311 return symbolic(at_series,makesequence(e,x,lim_point,ordre)); 3312 } 3313 } 3314 3315 // Main series entry point series(const gen & e,const identificateur & x,const gen & lim_point,int ordre,int direction,GIAC_CONTEXT)3316 gen series(const gen & e,const identificateur & x,const gen & lim_point,int ordre,int direction,GIAC_CONTEXT){ 3317 int save_series_flags=series_flags(contextptr); 3318 series_flags(save_series_flags | 8,contextptr); 3319 if (has_op(e,*at_surd) || has_op(e,*at_NTHROOT)){ 3320 vecteur subst1,subst2; 3321 surd2pow(e,subst1,subst2,contextptr); 3322 gen g=subst(e,subst1,subst2,false,contextptr); 3323 g=series(g,x,lim_point,ordre,direction,contextptr); 3324 series_flags(save_series_flags,contextptr); 3325 return subst(g,subst2,subst1,false,contextptr); 3326 } 3327 if (e.type==_VECT){ 3328 vecteur res = *e._VECTptr; 3329 int l=int(res.size()); 3330 for (int i=0;i<l;++i){ 3331 res[i]=in_series(_pow2exp(tan2sincos(res[i],contextptr),contextptr),x,lim_point,ordre,direction,contextptr); 3332 } 3333 series_flags(save_series_flags,contextptr); 3334 return res; 3335 } 3336 gen res=in_series(_pow2exp(tan2sincos(e,contextptr),contextptr),x,lim_point,ordre,direction,contextptr); 3337 series_flags(save_series_flags,contextptr); 3338 return res; 3339 } 3340 series(const gen & e,const gen & vars,const gen & lim_point,int ordre,int direction,GIAC_CONTEXT)3341 gen series(const gen & e,const gen & vars,const gen & lim_point,int ordre,int direction,GIAC_CONTEXT){ 3342 gen x,l; 3343 if (is_equal(vars)){ 3344 // vars= x==lim_point, overwrites lim_point definition 3345 // direction is given by lim_point (for interactive input) 3346 x = (*(vars._SYMBptr->feuille._VECTptr)) [0]; 3347 l = (*(vars._SYMBptr->feuille._VECTptr)) [1]; 3348 if (lim_point.type!=_INT_) 3349 return gensizeerr(contextptr); 3350 if (absint(lim_point.val)>0){ 3351 if (!direction && absint(ordre)<2) 3352 direction=ordre; 3353 ordre=absint(lim_point.val); 3354 } 3355 else 3356 direction = lim_point.val; 3357 } 3358 else { 3359 x=vars; 3360 l=lim_point; 3361 } 3362 if (x.type==_VECT && l.type==_VECT){ 3363 vecteur &v=*x._VECTptr; 3364 gen h(identificateur(" h")); 3365 vecteur w=addvecteur(*l._VECTptr,multvecteur(h,v)); 3366 gen newe=subst(e,v,w,false,contextptr); 3367 sparse_poly1 res=series__SPOL1(newe,*h._IDNTptr,zero,ordre,direction,contextptr); 3368 poly_truncate(res,ordre,contextptr); 3369 if (!res.empty() && is_undef(res.back().coeff)) 3370 res.pop_back(); 3371 // order term has been removed 3372 gen remains; 3373 gen r=sparse_poly12gen(res,1,remains,false); 3374 return subst(r,v,subvecteur(v,*l._VECTptr),false,contextptr); 3375 } 3376 if (x.type!=_IDNT){ 3377 identificateur xx("x"); 3378 gen res=series(subst(e,x,xx,false,contextptr),xx,l,ordre,direction,contextptr); 3379 return subst(res,xx,x,false,contextptr); 3380 } 3381 return series(e,*x._IDNTptr,l,ordre,direction,contextptr); 3382 } 3383 series(const gen & e,const gen & vars,const gen & lim_point,const gen & ordre0,GIAC_CONTEXT)3384 gen series(const gen & e,const gen & vars,const gen & lim_point,const gen & ordre0,GIAC_CONTEXT){ 3385 gen ordre(ordre0); 3386 if (!is_integral(ordre)) 3387 return gensizeerr(contextptr); 3388 return series(e,vars,lim_point,ordre.val,0,contextptr); // it's the direction 3389 } 3390 3391 // "unary" version 3392 static const char _series_s []="series"; _series(const gen & args,GIAC_CONTEXT)3393 gen _series(const gen & args,GIAC_CONTEXT){ 3394 if ( args.type==_STRNG && args.subtype==-1) return args; 3395 if (args.type==_SPOL1) 3396 return args; 3397 if (args.type==_VECT && args._VECTptr->size()==2 && args._VECTptr->front().type==_STRNG && args._VECTptr->back().type==_INT_){ 3398 int n=args._VECTptr->back().val; 3399 if (n<=0 || n>1024) 3400 return gensizeerr("Default series order must be >0 and <=1024"); 3401 series_default_order(n,contextptr); 3402 return _series(args._VECTptr->front(),contextptr); 3403 } 3404 if (args.type==_STRNG && args._STRNGptr->size()==1){ 3405 char ch=(*args._STRNGptr)[0]; 3406 gen h(*args._STRNGptr,contextptr); 3407 series_variable_name(ch,contextptr); 3408 sparse_poly1 s(1,monome(1,1)); 3409 sto(s,h,contextptr); 3410 series_flags(contextptr) = series_flags(contextptr) | (1<<5); 3411 *logptr(contextptr) << "Setting " << ch << " as series variable name" << '\n'; 3412 string Os=abs_calc_mode(contextptr)==38?"b":"O"; 3413 gen O(Os,contextptr); 3414 if (eval(O,1,contextptr)!=O) 3415 *logptr(contextptr) << "Purge "<<Os<<" if you want to use "<<Os<<"("<< h <<"^...) notation"<< '\n'; 3416 else { 3417 gen prog=symb_program(vx_var,0,vx_var*symbolic(at_order_size,h),contextptr); 3418 *logptr(contextptr) << "Assigning "<<Os<<" so that you can use use "<<Os<<"("<< h<<"^...) notation"<< '\n'; 3419 sto(prog,O,contextptr); 3420 series_flags(contextptr)=series_flags(contextptr) | (1<<6) ; 3421 } 3422 return s; 3423 } 3424 if (args.type!=_VECT) 3425 return series(args,vx_var,0,series_default_order(contextptr),0,contextptr); 3426 vecteur v=*args._VECTptr; 3427 if (v.empty()) 3428 return gensizeerr(contextptr); 3429 if (v.back().type==_INT_ && v.back().subtype==_INT_MAPLECONVERSION && v.back().val==_POLY1__VECT){ 3430 gen p=v.back(); 3431 v.pop_back(); 3432 gen res=_series(gen(v,_SEQ__VECT),contextptr); 3433 res=_convert(makesequence(res,p),contextptr); 3434 return res; 3435 } 3436 v[0]=Heavisidetosign(when2sign(piecewise2when(v[0],contextptr),contextptr),contextptr); 3437 int s=int(v.size()); 3438 if (!s) 3439 toofewargs(_series_s); 3440 if (s==1) 3441 return series( v[0],vx_var,0,series_default_order(contextptr),0,contextptr); 3442 if (s==2){ 3443 if (v[1].type==_INT_) 3444 return series( v[0],vx_var,0,v[1],contextptr); 3445 return series( v[0],v[1],0,series_default_order(contextptr),contextptr); 3446 } 3447 if (s==3){ 3448 if ( (v[1].type==_VECT && v[2].type==_VECT) || 3449 ( v[1].type==_IDNT || ( v[1].type==_SYMB && (v[1]._SYMBptr->sommet==at_equal || v[1]._SYMBptr->sommet==at_equal2 || v[1]._SYMBptr->sommet==at_at ) ) ) 3450 ) 3451 return series( v[0],v[1],v[2],series_default_order(contextptr),contextptr); 3452 if (v[1].type==_VECT){ 3453 vecteur vars=*v[1]._VECTptr,vals(vars.size()); 3454 for (int i=0;i<vars.size();++i){ 3455 if (vars[i].is_symb_of_sommet(at_equal)){ 3456 vals[i]=vars[i]._SYMBptr->feuille[1]; 3457 vars[i]=vars[i]._SYMBptr->feuille[0]; 3458 } 3459 } 3460 return series(v[0],vars,vals,v[2],contextptr); 3461 } 3462 return series( v[0],symbolic(at_equal,makesequence(vx_var,v[1])),v[2],series_default_order(contextptr),contextptr); 3463 } 3464 if (s==4) 3465 return series( v[0],v[1],v[2],v[3],contextptr); 3466 if (s>5 || v[3].type!=_INT_ || v[4].type!=_INT_) 3467 return gentoomanyargs(_series_s); 3468 return series(v[0],v[1],v[2],v[3].val,v[4].val,contextptr); 3469 } 3470 static define_unary_function_eval (__series,&_series,_series_s); 3471 define_unary_function_ptr5( at_series ,alias_at_series,&__series,0,true); 3472 3473 // "unary" version 3474 static const char _revert_s []="revert"; _revert(const gen & args,GIAC_CONTEXT)3475 gen _revert(const gen & args,GIAC_CONTEXT){ 3476 if ( args.type==_STRNG && args.subtype==-1) return args; 3477 if (args.type==_SPOL1){ 3478 const sparse_poly1 & s=*args._SPOL1ptr; 3479 if (s.empty() || s.front().exponent==0) 3480 return gensizeerr("Not invertible for composition"); 3481 sparse_poly1 res; 3482 if (!prevert(s,res,contextptr)) 3483 return gensizeerr("Not invertible for composition"); 3484 return res; 3485 } 3486 vecteur v=gen2vecteur(args); 3487 if (v.empty()) 3488 return gensizeerr(contextptr); 3489 gen g=v[0],x; 3490 if (v.size()==1) 3491 x=vx_var; 3492 else 3493 x=v[1]; 3494 if (x.type!=_IDNT){ 3495 identificateur idx(" trevert"); 3496 return _revert(subst(args,x,idx,false,contextptr),contextptr); 3497 } 3498 // find ordre 3499 int ordre=series_default_order(contextptr); 3500 if (v.size()>2 && v[2].type==_INT_) 3501 ordre=v[2].val; 3502 vecteur w=lop(g,at_order_size); 3503 if (w.size()==1){ 3504 gen xn=derive(g,w.front(),contextptr); 3505 if (is_undef(xn)) return xn; 3506 if (xn.is_symb_of_sommet(at_pow)){ 3507 gen & f=xn._SYMBptr->feuille; 3508 if (f.type==_VECT && f._VECTptr->size()==2 && f._VECTptr->back().type==_INT_){ 3509 ordre=f._VECTptr->back().val; 3510 w.clear(); 3511 g=subst(g,w.front(),0,false,contextptr); 3512 } 3513 } 3514 } 3515 if (!w.empty()) 3516 return gensizeerr(contextptr); 3517 sparse_poly1 p=series__SPOL1(g,*x._IDNTptr,zero,ordre,0,contextptr),res; 3518 if (!prevert(p,res,contextptr)) 3519 return gensizeerr(contextptr); 3520 gen remains; 3521 return sparse_poly12gen(res,x,remains,false); 3522 } 3523 static define_unary_function_eval (__revert,&_revert,_revert_s); 3524 define_unary_function_ptr5( at_revert ,alias_at_revert,&__revert,0,true); 3525 3526 static const char _bounded_function_s []="bounded_function"; _bounded_function(const gen & args,GIAC_CONTEXT)3527 gen _bounded_function(const gen & args,GIAC_CONTEXT){ 3528 if ( args.type==_STRNG && args.subtype==-1) return args; 3529 return symbolic(at_bounded_function,args); 3530 } bounded_function(GIAC_CONTEXT)3531 gen bounded_function(GIAC_CONTEXT){ 3532 int i=bounded_function_no(contextptr); 3533 ++i; 3534 bounded_function_no(i,contextptr); 3535 return symbolic(at_bounded_function,i); 3536 } 3537 static define_unary_function_eval (__bounded_function,&_bounded_function,_bounded_function_s); 3538 define_unary_function_ptr5( at_bounded_function ,alias_at_bounded_function,&__bounded_function,0,true); 3539 3540 // internal function, used to replace sum for limit/series 3541 // args = expression, antiderivative, variable, lower_bound, upper_bound 3542 static const char _euler_mac_laurin_s []="euler_mac_laurin"; _euler_mac_laurin(const gen & args,GIAC_CONTEXT)3543 gen _euler_mac_laurin(const gen & args,GIAC_CONTEXT){ 3544 if ( args.type==_STRNG && args.subtype==-1) return args; 3545 return symbolic(at_euler_mac_laurin,args); 3546 } 3547 static define_unary_function_eval (__euler_mac_laurin,&_euler_mac_laurin,_euler_mac_laurin_s); 3548 define_unary_function_ptr5( at_euler_mac_laurin ,alias_at_euler_mac_laurin,&__euler_mac_laurin,0,true); 3549 convert_to_euler_mac_laurin(const gen & g,const identificateur & n,gen & res,GIAC_CONTEXT)3550 bool convert_to_euler_mac_laurin(const gen & g,const identificateur & n,gen & res,GIAC_CONTEXT){ 3551 if (g.is_symb_of_sommet(at_sum)){ 3552 gen & f = g._SYMBptr->feuille; 3553 if (f.type!=_VECT || f._VECTptr->size()!=4) 3554 return false; 3555 gen l=in_limit((f[3]-f[2])/n,n,plus_inf,1,contextptr); 3556 if (is_zero(l) || is_undef(l) || is_inf(l)) 3557 return false; 3558 // check that the expression to be summed has a derivative 3559 // which is a o(expression) as n -> inf 3560 gen f0=f._VECTptr->front(); 3561 gen x =f[1]; 3562 if (x.type!=_IDNT){ 3563 *logptr(contextptr) << gettext("Unable to convert to euler mac laurin"); 3564 return false; 3565 } 3566 gen f0prime=derive(f0,x,contextptr), f03=derive(f0prime,x,contextptr); 3567 f03=derive(f03,x,contextptr); 3568 if (is_undef(f03)) return false; 3569 l=in_limit(f03/f0prime,n,plus_inf,1,contextptr); 3570 if (!is_zero(l)) 3571 return false; 3572 gen remains; 3573 gen F0=integrate_gen_rem(f0,x,remains,contextptr); 3574 if (!is_zero(remains) || is_undef(F0)) 3575 return false; 3576 res=symbolic(at_euler_mac_laurin,gen(makevecteur(f0,F0,x,f[2],f[3]),_SEQ__VECT)); 3577 return true; 3578 } 3579 vecteur v=lop(g,at_sum); 3580 vecteur w=v; 3581 int s=int(v.size()); 3582 for (int i=0;i<s;++i){ 3583 if (!convert_to_euler_mac_laurin(v[i],n,w[i],contextptr)) 3584 return false; 3585 } 3586 res=subst(g,v,w,false,contextptr); 3587 return true; 3588 } 3589 3590 #ifndef NO_NAMESPACE_GIAC 3591 } // namespace giac 3592 #endif // ndef NO_NAMESPACE_GIAC 3593 3594