1 // -*- mode:C++ ; compile-command: "g++-3.4 -I.. -g -c intgab.cc -DHAVE_CONFIG_H -DIN_GIAC" -*- 2 #include "giacPCH.h" 3 /* 4 * Copyright (C) 2000,2014 B. Parisse, Institut Fourier, 38402 St Martin d'Heres 5 * 6 * This program is free software; you can redistribute it and/or modify 7 * it under the terms of the GNU General Public License as published by 8 * the Free Software Foundation; either version 3 of the License, or 9 * (at your option) any later version. 10 * 11 * This program is distributed in the hope that it will be useful, 12 * but WITHOUT ANY WARRANTY; without even the implied warranty of 13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 14 * GNU General Public License for more details. 15 * 16 * You should have received a copy of the GNU General Public License 17 * along with this program. If not, see <http://www.gnu.org/licenses/>. 18 */ 19 using namespace std; 20 #include <stdexcept> 21 #include "vector.h" 22 #include <cmath> 23 #include <cstdlib> 24 #include "sym2poly.h" 25 #include "usual.h" 26 #include "intgab.h" 27 #include "subst.h" 28 #include "derive.h" 29 #include "lin.h" 30 #include "vecteur.h" 31 #include "gausspol.h" 32 #include "plot.h" 33 #include "prog.h" 34 #include "modpoly.h" 35 #include "series.h" 36 #include "tex.h" 37 #include "ifactor.h" 38 #include "risch.h" 39 #include "solve.h" 40 #include "intg.h" 41 #include "desolve.h" 42 #include "alg_ext.h" 43 #include "misc.h" 44 #include "maple.h" 45 #include "rpn.h" 46 #include "giacintl.h" 47 #ifdef HAVE_LIBGSL 48 #include <gsl/gsl_math.h> 49 #include <gsl/gsl_sf_gamma.h> 50 #include <gsl/gsl_sf_psi.h> 51 #include <gsl/gsl_sf_zeta.h> 52 #include <gsl/gsl_odeiv.h> 53 #include <gsl/gsl_errno.h> 54 #endif 55 56 #ifndef NO_NAMESPACE_GIAC 57 namespace giac { 58 #endif // ndef NO_NAMESPACE_GIAC 59 60 // check whether an expression is meromorphic 61 // returns -1 if x is not an IDNT 62 // return 2 if rational 63 // return 3 if rational fraction of x, sin(a*x+b), cos(a*x+b) 64 // return 4 if rational fraction of x, exp(a*x+b) where re(a)!=0 65 // return 5 if ln(P)*A==rational fraction + B==rational fraction 66 // TODO: 6 if exp(a[0]*x^2+a[1]*x+a[2])*b+P, P,b =rational, a[0]<0 is_meromorphic(const gen & g,const gen & x,gen & a,gen & b,gen & P,GIAC_CONTEXT)67 int is_meromorphic(const gen & g,const gen & x,gen & a,gen & b,gen & P,GIAC_CONTEXT){ 68 if (x.type!=_IDNT) 69 return -1; 70 if (g.type<=_IDNT) 71 return 2; 72 if (g.type==_VECT){ 73 const_iterateur it=g._VECTptr->begin(),itend=g._VECTptr->end(); 74 for (;it!=itend;++it){ 75 if (!is_meromorphic(*it,x,a,b,P,contextptr)) 76 return 0; 77 } 78 return 1; 79 } 80 if (g.type==_SYMB){ 81 vecteur v; 82 rlvarx(g,x,v); 83 islesscomplexthanf_sort(v.begin(),v.end()); 84 if (v==vecteur(1,x)) 85 return 2; 86 if (v.size()<2 || v[1].type!=_SYMB) 87 return 0; 88 P=v[1]._SYMBptr->feuille; 89 if (v.size()==2) { 90 if (v[1].is_symb_of_sommet(at_exp)){ 91 if (is_linear_wrt(P,x,a,b,contextptr)){ 92 if (is_zero(re(a,contextptr))){ 93 a=a/cst_i; 94 b=b/cst_i; 95 return 3; 96 } 97 return 4; 98 } 99 identificateur t(" t"); 100 gen tt(t); 101 gen g2=subst(g,v[1],tt,false,contextptr),bcst; 102 if (is_linear_wrt(g2,t,b,bcst,contextptr)){ 103 gen A,B,C; 104 if (is_quadratic_wrt(P,x,A,B,C,contextptr) && is_strictly_positive(-A,contextptr)){ 105 a=makevecteur(A,B,C); 106 P=bcst; 107 return 6; 108 } 109 } 110 } 111 if ( (v[1].is_symb_of_sommet(at_cos) || v[1].is_symb_of_sommet(at_sin)) && is_linear_wrt(P,x,a,b,contextptr) && is_zero(im(a,contextptr)) ) 112 return 3; 113 if (v[1].is_symb_of_sommet(at_ln)){ 114 identificateur t(" t"); 115 gen tt(t); 116 gen g2=subst(g,v[1],tt,false,contextptr); 117 if (is_linear_wrt(g2,t,a,b,contextptr)) 118 return 5; 119 } 120 } 121 if (v.size()==3){ 122 if (v[1].is_symb_of_sommet(at_cos) && v[2].is_symb_of_sommet(at_sin)){ 123 gen v2f=v[2]._SYMBptr->feuille; 124 if (P==v2f && is_linear_wrt(P,x,a,b,contextptr) && is_zero(im(a,contextptr))) 125 return 3; 126 } 127 } 128 const_iterateur it=v.begin(),itend=v.end(); 129 for (;it!=itend;++it){ 130 if (it->type==_IDNT) 131 continue; 132 if (it->type!=_SYMB) 133 return 0; 134 unary_function_ptr * u=&it->_SYMBptr->sommet; 135 if (u==at_sin || u==at_cos || u==at_tan || u==at_exp || u==at_sinh || u==at_cosh || u==at_tanh){ 136 // the argument must not have singularities 137 if (find_singularities(it->_SYMBptr->feuille,*x._IDNTptr,1,contextptr).empty()) 138 continue; 139 } 140 return 0; 141 } 142 return 1; 143 } 144 return 0; 145 } 146 fast_is_even_odd(const gen & f,const gen & x,GIAC_CONTEXT)147 int fast_is_even_odd(const gen & f,const gen & x,GIAC_CONTEXT){ 148 if (f==x) // x is odd 149 return 2; 150 if (f.type==_VECT){ 151 vecteur & v =*f._VECTptr; 152 const_iterateur it=v.begin(),itend=v.end(); 153 int res=0,current; 154 for (;it!=itend;++it){ 155 if (! (current=fast_is_even_odd(*it,x,contextptr)) ) 156 return 0; 157 if (!res) 158 res=current; 159 if (res!=current) 160 return 0; 161 } 162 return res; 163 } 164 if (f.type!=_SYMB) // constant is even 165 return 1; 166 gen & ff=f._SYMBptr->feuille; 167 int res; 168 const unary_function_ptr & u = f._SYMBptr->sommet; 169 if (u==at_pow && ff.type==_VECT && ff._VECTptr->size()==2){ 170 gen ff2=ff._VECTptr->back(); 171 if (ff2.type==_INT_){ 172 gen ff1=ff._VECTptr->front(); 173 res=fast_is_even_odd(ff1,x,contextptr); 174 if (res<2) 175 return res; 176 return (ff2.val%2)?2:1; 177 } 178 return 0; 179 } 180 res=fast_is_even_odd(ff,x,contextptr); 181 if (res<2) 182 return res; 183 // ff is odd 184 if (u==at_plus || u==at_neg || u==at_inv || u==at_sin || u==at_tan || u==at_sinh || u==at_tanh || u==at_atan || u==at_atanh) 185 return res; 186 if (u==at_prod){ 187 if (ff.type!=_VECT || (ff._VECTptr->size()%2)) 188 return res; 189 return 1; 190 } 191 if (u==at_cos || u==at_cosh || u==at_abs) 192 return 1; 193 return 0; 194 } 195 196 // 0 none, 1 even, 2 odd is_even_odd(const gen & f,const gen & x,GIAC_CONTEXT)197 int is_even_odd(const gen & f,const gen & x,GIAC_CONTEXT){ 198 int res=fast_is_even_odd(f,x,contextptr); 199 if (res) 200 return res; 201 gen f1=f,f2=subst(f,x,-x,false,contextptr); 202 if (lvar(f)==vecteur(1,x)){ // rational case 203 if (is_zero(normal(f1-f2,contextptr))) 204 return 1; 205 if (is_zero(normal(f1+f2,contextptr))) 206 return 2; 207 return 0; 208 } 209 f1=normal(_texpand(f1,contextptr),contextptr); 210 f2=normal(_texpand(f2,contextptr),contextptr); 211 if (f1==f2) 212 return 1; 213 if (is_zero(ratnormal(f1+f2,contextptr))) 214 return 2; 215 return 0; 216 } 217 has_sparse_poly1(const gen & g)218 bool has_sparse_poly1(const gen & g){ 219 if (g.type==_SPOL1) 220 return true; 221 if (g.type==_VECT){ 222 vecteur & v=*g._VECTptr; 223 for (size_t i=0;i<v.size();++i){ 224 if (has_sparse_poly1(v[i])) 225 return true; 226 } 227 return false; 228 } 229 if (g.type==_SYMB) 230 return has_sparse_poly1(g._SYMBptr->feuille); 231 return false; 232 } 233 234 // residue of g at x=a residue(const gen & g_,const gen & x,const gen & a,GIAC_CONTEXT)235 gen residue(const gen & g_,const gen & x,const gen & a,GIAC_CONTEXT){ 236 if (x.type!=_IDNT) 237 return gensizeerr(contextptr); 238 gen xval=x._IDNTptr->eval(1,x,contextptr); 239 if (xval!=x){ 240 #if 1 241 identificateur tmpid(x._IDNTptr->id_name+string("_")); 242 gen tmp(tmpid); 243 gen g(subst(g_,x,tmp,false,contextptr)); 244 return residue(g,tmp,a,contextptr); 245 #else 246 _purge(x,contextptr); 247 xval=x._IDNTptr->eval(1,x,contextptr); 248 if (xval!=x){ 249 string s="Unable to purge "+x.print(contextptr)+ ", choose another free variable name"; 250 *logptr(contextptr) << s << '\n'; 251 return gensizeerr(s); 252 } 253 gen res=residue(g_,x,a,contextptr); 254 sto(xval,x,contextptr); 255 return res; 256 #endif 257 } 258 gen g1=fxnd(g_); 259 if (g1.type==_VECT && g1._VECTptr->size()==2){ 260 gen n=g1._VECTptr->front(),d=derive(g1._VECTptr->back(),x,contextptr); 261 gen da=subst(d,*x._IDNTptr,a,false,contextptr); 262 if (!is_zero(da)){ 263 gen na=subst(n,*x._IDNTptr,a,false,contextptr); 264 n=recursive_normal(na/da,contextptr); 265 if (!is_zero(n) && !is_undef(n) && !is_inf(n)) 266 return n; 267 } 268 } 269 gen g=_pow2exp(tan2sincos(g_,contextptr),contextptr); 270 for (int ordre=2;ordre<max_series_expansion_order;ordre=2*ordre){ 271 #ifdef TIMEOUT 272 control_c(); 273 #endif 274 if (ctrl_c || interrupted) { 275 interrupted = true; ctrl_c=false; 276 gen res; 277 gensizeerr(gettext("Stopped by user interruption."),res); 278 return res; 279 } 280 sparse_poly1 s=series__SPOL1(g,*x._IDNTptr,a,ordre,0,contextptr); 281 int n=int(s.size()); 282 if (n && (is_undef(s[0].coeff) || is_undef(s[0].exponent))) 283 continue; // stop the loop, try a larger order 284 for (int i=0;i<n;++i){ 285 gen e1=s[i].exponent+1; 286 if (is_strictly_positive(e1,contextptr)){ 287 gen res=s[i].coeff; 288 if (!has_sparse_poly1(res) && is_zero(derive(res,x,contextptr))) 289 return 0; 290 return gensizeerr(gettext("Non holomorphic function ")+g.print(contextptr)+" at "+x.print(contextptr)+"="+a.print(contextptr)); 291 } 292 if (e1.type==_FRAC) 293 return undef; 294 if (is_undef(s[i].coeff)) 295 break; // stop the loop, try a larger order 296 if (is_zero(e1)){ 297 gen res=s[i].coeff; 298 if (is_zero(derive(res,x,contextptr))) 299 return res; 300 return gensizeerr(gettext("Non holomorphic function ")+g.print(contextptr)+" at "+x.print(contextptr)+"="+a.print(contextptr)); 301 } 302 if (i==n-1) // exact expansion 303 return 0; 304 } 305 } 306 return genmaxordererr(contextptr); 307 } _residue(const gen & args,GIAC_CONTEXT)308 gen _residue(const gen & args,GIAC_CONTEXT){ 309 if ( args.type==_STRNG && args.subtype==-1) return args; 310 if (args.type!=_VECT) return gensizeerr(contextptr); 311 vecteur v(*args._VECTptr); 312 int s=int(v.size()); 313 if (s==2){ 314 if (is_equal(v[1])){ 315 vecteur & w=*v[1]._SYMBptr->feuille._VECTptr; 316 v.push_back(w[1]); 317 v[1]=w[0]; 318 } 319 else 320 return gensizeerr(gettext("Syntax residue(expr,x=a)")); 321 ++s; 322 } 323 if (s<3) 324 return gensizeerr(contextptr); 325 return residue(v[0],v[1],v[2],contextptr); 326 } 327 static const char _residue_s []="residue"; 328 static define_unary_function_eval (__residue,&_residue,_residue_s); 329 define_unary_function_ptr5( at_residue ,alias_at_residue,&__residue,0,true); 330 singular(const gen & g,const gen & x,GIAC_CONTEXT)331 vecteur singular(const gen & g,const gen & x,GIAC_CONTEXT){ 332 // FIXME handle set of variables 333 if (x.type!=_IDNT) 334 return vecteur(1,gentypeerr(contextptr)); 335 vecteur res=find_singularities(g,*x._IDNTptr,1,contextptr); 336 vecteur res2; 337 const_iterateur it=res.begin(),itend=res.end(); 338 for (;it!=itend;++it){ 339 if (!equalposcomp(res2,*it)) 340 res2.push_back(*it); 341 } 342 return res2; 343 } _singular(const gen & args,GIAC_CONTEXT)344 gen _singular(const gen & args,GIAC_CONTEXT){ 345 if ( args.type==_STRNG && args.subtype==-1) return args; 346 vecteur v(gen2vecteur(args)); 347 int s=int(v.size()); 348 if (s==1){ 349 v.push_back(vx_var); 350 ++s; 351 } 352 if (s<2) 353 return gensizeerr(contextptr); 354 if (s>2 && v[1].type==_IDNT){ 355 vecteur res=find_singularities(v[0],*v[1]._IDNTptr,(is_zero(v[2])?1:9),contextptr); 356 comprim(res); 357 return res; 358 } 359 return singular(v[0],v[1],contextptr); 360 } 361 static const char _singular_s []="singular"; 362 static define_unary_function_eval (__singular,&_singular,_singular_s); 363 define_unary_function_ptr5( at_singular ,alias_at_singular,&__singular,0,true); 364 assume_t_in_ab(const gen & t,const gen & a,const gen & b,bool exclude_a,bool exclude_b,GIAC_CONTEXT)365 bool assume_t_in_ab(const gen & t,const gen & a,const gen & b,bool exclude_a,bool exclude_b,GIAC_CONTEXT){ 366 vecteur v_interval(1,gen(makevecteur(a,b),_LINE__VECT)); 367 vecteur v_excluded; 368 if (exclude_a) 369 v_excluded.push_back(a); 370 if (exclude_b) 371 v_excluded.push_back(b); 372 return !is_undef(sto(gen(makevecteur(gen(_DOUBLE_).change_subtype(1),v_interval,v_excluded),_ASSUME__VECT),t,contextptr)); 373 } 374 375 // reduce g, a rational fraction wrt to x, to a sqff lnpart 376 // and adds the non sqff integrated part to ratpart intreduce(const gen & e,const gen & x,gen & lnpart,gen & ratpart,GIAC_CONTEXT)377 static bool intreduce(const gen & e,const gen & x,gen & lnpart,gen & ratpart,GIAC_CONTEXT){ 378 vecteur l; 379 l.push_back(x); // insure x is the main var 380 l=vecteur(1,l); 381 alg_lvar(e,l); 382 int s=int(l.front()._VECTptr->size()); 383 if (!s){ 384 l.erase(l.begin()); 385 s=int(l.front()._VECTptr->size()); 386 } 387 if (!s) 388 return false; 389 gen r=e2r(e,l,contextptr); 390 gen r_num,r_den; 391 fxnd(r,r_num,r_den); 392 if (r_num.type==_EXT){ 393 return false; 394 } 395 if (r_den.type!=_POLY){ 396 l.front()._VECTptr->front()=x; 397 lnpart=0; 398 if (r_num.type==_POLY) 399 ratpart=rdiv(r2e(r_num._POLYptr->integrate(),l,contextptr),r2sym(r_den,l,contextptr),contextptr); 400 else 401 ratpart=e*x; 402 return true; 403 } 404 polynome den(*r_den._POLYptr),num(s); 405 if (r_num.type==_POLY) 406 num=*r_num._POLYptr; 407 else 408 num=polynome(r_num,s); 409 l.front()._VECTptr->front()=x; 410 polynome p_content(lgcd(den)); 411 factorization vden(sqff(den/p_content)); // first square-free factorization 412 vector< pf<gen> > pfde_VECT; 413 polynome ipnum(s),ipden(s),temp(s),tmp(s); 414 partfrac(num,den,vden,pfde_VECT,ipnum,ipden); 415 vector< pf<gen> >::iterator it=pfde_VECT.begin(); 416 vector< pf<gen> >::const_iterator itend=pfde_VECT.end(); 417 vector< pf<gen> > ratpartv; 418 for (;it!=itend;++it){ 419 pf<gen> single(intreduce_pf(*it,ratpartv)); 420 lnpart += r2e(single.num,l,contextptr)/r2e(single.den,l,contextptr); 421 // FIXME: add ratpartv to ratpart 422 } 423 return true; 424 } 425 intgab_sincos(const gen & g,const gen & x,const gen & a,const gen & b,gen & A,gen & B,gen & res,GIAC_CONTEXT)426 static bool intgab_sincos(const gen & g,const gen & x,const gen & a,const gen & b,gen & A, gen & B,gen & res,GIAC_CONTEXT){ 427 gen g1=trig2exp(g,contextptr); 428 // write it as a rational fraction of x,X=exp(i*(A*x+b)) 429 identificateur Xid(" X"); 430 gen X(Xid),expx(exp(cst_i*(A*x+B),contextptr)); 431 g1=subst(g1,expx,X,false,contextptr); 432 if (is_positive(-A,contextptr)) 433 g1=subst(g1,X,inv(X,contextptr),false,contextptr); 434 // Separable variables? 435 vecteur f=factors(g1,x,contextptr); // Factor then split factors 436 gen xfact(plus_one),Xfact(plus_one); 437 if (separate_variables(f,Xid,x,Xfact,xfact,contextptr)){ 438 // xfact must be a proper fraction 439 if (!is_zero(limit(xfact,*x._IDNTptr,plus_inf,1,contextptr))){ 440 res=undef; 441 return true; 442 } 443 // Xfact must be rewritten as a generalized polynomial part 444 // + a rational part N/D, the roots of D are in pairs r, 1/conj(r) 445 // if there are roots with norm = 1, D vanishes inf. times on R 446 // hence the integral is undef 447 // we select all roots with norm > 1 and take the corresp. part of 448 // the partial fraction expansion, we have 449 // N/D=2*re(true_poly_part)+2*re(n/d)-re(n(0)/d(0)) 450 // where d has no poles in C^+ and X tends to 0 at inf in C^+ 451 vecteur vX(1,X); 452 lvar(Xfact,vX); 453 gen ND=sym2r(Xfact,vX,contextptr); 454 gen N,D; 455 fxnd(ND,N,D); 456 polynome Np,Dp,Q,R; 457 if (D.type!=_POLY) 458 return false; 459 Dp=*D._POLYptr; 460 if (N.type==_POLY) 461 Np=*N._POLYptr; 462 else 463 Np=polynome(N,1); 464 Np.TDivRem(Dp,Q,R); 465 int Qd=Q.degree(0); 466 // Q is the true poly part 467 if (Qd){ 468 // Dp must be divisible by Qd 469 int Dval=Dp.valuation(0); 470 if (Dval!=Qd) 471 return false; // setsizeerr(); 472 index_t decal(vX.size()); 473 decal[0]=-Dval; 474 Dp=Dp.shift(decal); 475 polynome XQd(gen(1),int(vX.size())),U(int(vX.size())),V(int(vX.size())),C(int(vX.size())); 476 decal[0]=Dval; 477 XQd=XQd.shift(decal); 478 Tabcuv(Dp,XQd,R,U,V,C); // C*Np=Dp*U+X^Dval*V 479 // Np/(Dp*X^Dval)=V/Dp/C+... 480 R=V; 481 Dp=Dp*C; 482 } 483 if (!Q.coord.empty() && Q.coord.back().index.is_zero()) 484 Q.coord.back().value=Q.coord.back().value/2; 485 // R/Dp is the true fractional part, Q is the true poly. part 486 // now check roots of norm=1 487 gen dp=r2sym(Dp,vX,contextptr); 488 identificateur XXi(" x"),XYi(" y"); 489 gen XX(XXi),XY(XYi); 490 dp=subst(dp,X,XX+cst_i*XY,false,contextptr); 491 dp=_resultant(gen(makevecteur(dp,XX*XX+XY*XY-1,XY),_SEQ__VECT),contextptr); 492 if (is_undef(dp)) return false; 493 dp=gcd(re(dp,contextptr),im(dp,contextptr),contextptr); 494 vecteur vdp=factors(dp,XX,contextptr); 495 int vdps=int(vdp.size()); 496 for (int i=0;i<vdps;i+=2){ 497 if (sturmab(vdp[i],XX,-1,1,contextptr)>0){ 498 res=undef; 499 return true; 500 } 501 } 502 // ok, now find roots of norm>1 503 factorization fd; 504 polynome Dp_content; 505 gen extra_div=1; 506 if (!factor(Dp,Dp_content,fd,false,true,true,1,extra_div) || extra_div!=1){ 507 *logptr(contextptr) << gettext("Unable to factor ") << r2sym(Dp,vX,contextptr) << '\n'; 508 res=undef; 509 return true; 510 } 511 // check that each factor has degree 1 512 polynome D1(gen(1),int(vX.size())); 513 factorization::const_iterator f_it=fd.begin(),f_itend=fd.end(); 514 for (;f_it!=f_itend;++f_it){ 515 if (f_it->fact.degree(0)>1){ 516 *logptr(contextptr) << gettext("Unable to factor ") << r2sym(f_it->fact,vX,contextptr) << '\n'; 517 res=undef; 518 return true; 519 } 520 if (f_it->fact.coord.size()==2){ 521 gen f1=f_it->fact.coord.front().value; 522 gen f2=f_it->fact.coord.back().value; 523 gen f21=r2sym(-f2/f1,vecteur(0),contextptr); 524 if (is_positive(abs(f21,contextptr)-1,contextptr)) 525 D1=D1*pow(f_it->fact,f_it->mult); 526 } 527 } // end f_it 528 // keep only those roots 529 polynome D2,tmp,U,V,C; 530 Dp.TDivRem(D1,D2,tmp); 531 Tabcuv(D1,D2,R,U,V,C); // R/D=(D1*U+D2*V)/(C*D1*D2) -> V/(C*D1) 532 Dp=C*D1; 533 R=V; // back to integrating 2*Re(R/Dp) 534 // find R/Dp at 0, real part should be substracted 535 vecteur Rv(polynome2poly1(R,1)),Dv(polynome2poly1(Dp,1)),Qv(polynome2poly1(Q,1)); 536 vecteur vX1(vX.begin()+1,vX.end()); 537 Rv=*r2sym(Rv,vX1,contextptr)._VECTptr; 538 Dv=*r2sym(Dv,vX1,contextptr)._VECTptr; 539 Qv=*r2sym(Qv,vX1,contextptr)._VECTptr; 540 gen correc=re(Rv.back()/Dv.back(),contextptr); 541 xfact=xfact*(2*(horner(Rv,expx)/horner(Dv,expx)+horner(Qv,expx))-correc); 542 res=0; 543 vecteur v=singular(xfact,x,contextptr); 544 if (!v.empty() && is_undef(v.front())) 545 return false; 546 int s=int(v.size()); 547 for (int i=0;i<s;++i){ 548 gen coeff=0,vi=v[i]; 549 if (is_real(vi,contextptr)){ 550 // check that xfact has a finite limit 551 gen test=limit(g,*x._IDNTptr,vi,0,contextptr); 552 if (is_inf(test)){ 553 res=undef; 554 return true; 555 } 556 coeff=1; 557 } 558 else { 559 if (ck_is_positive(im(v[i],contextptr),contextptr)) 560 coeff=2; 561 } 562 if (!is_zero(coeff)){ 563 coeff=coeff*residue(xfact,x,vi,contextptr); 564 if (is_undef(coeff)) 565 return false; 566 res=res+re(coeff*cst_i*cst_pi,contextptr); 567 } 568 } 569 res = normal(res,contextptr); 570 return true; 571 } // end if separate_variables 572 return false; 573 } 574 nvars_depend_x(const vecteur & v,const gen & x)575 static int nvars_depend_x(const vecteur & v,const gen & x){ 576 int res=0; 577 for (unsigned int i=0;i<v.size();i++){ 578 if (contains(v[i],x)) 579 ++res; 580 } 581 return res; 582 } 583 intgab_r(const gen & g0,const gen & x,const gen & a,const gen & b,bool rational,gen & res,GIAC_CONTEXT)584 bool intgab_r(const gen & g0,const gen & x,const gen & a,const gen & b,bool rational,gen & res,GIAC_CONTEXT){ 585 gen g(g0); 586 // check if g may be integrated using the residue formula 587 gen A,B,P; 588 int typeint=rational?2:is_meromorphic(g,x,A,B,P,contextptr); 589 if (typeint==6 && a==minus_inf && b==plus_inf && is_zero(P) && is_constant_wrt(A,x,contextptr) && lvarxwithinv(B,x,contextptr).size()<2 && is_strictly_positive(-A[0],contextptr)){ 590 gen A_(A[0]),B_(A[1]),C_(A[2]); 591 gen der(2*A_*x+B_),tmp(B); 592 B=0; 593 while (!is_zero(tmp,contextptr)){ 594 tmp=_quorem(makesequence(tmp,der,x),contextptr); 595 if (tmp.type!=_VECT || tmp._VECTptr->size()!=2) 596 return false; 597 B += tmp._VECTptr->back(); 598 tmp=-derive(tmp._VECTptr->front(),x,contextptr); 599 } 600 // int(exp(A_*x^2+B_*x+C_)*B,x,-inf,inf) 601 // =int(exp(A_*(x+B_/2/A_)^2+(-B_^2/4/A_+C_))*B,x,-inf,inf) 602 // =sqrt(pi/A_)*B*exp(B_^2/4/A_-C_) 603 res=sqrt(-cst_pi/A_,contextptr)*B*exp(ratnormal(C_-B_*B_/4/A_,contextptr),contextptr); 604 return true; 605 } 606 if (typeint==5){ 607 // A*ln(P(x))+B, A/B/P rational fractions 608 bool estreel=is_zero(im(A,contextptr))&&is_zero(im(P,contextptr)); 609 vecteur lv(1,x); 610 lvar(P,lv); 611 int lvs=int(lv.size()); 612 for (int i=0;i<lvs;++i){ 613 gen tmp=derive(lv[i],x,contextptr); 614 if (!is_zero(tmp) && !is_one(tmp)) 615 return false; 616 } 617 gen resB; 618 if (!intgab(B,x,a,b,resB,contextptr)) 619 return false; 620 gen Af=sym2r(A,lv,contextptr),Anum,Aden; 621 fxnd(Af,Anum,Aden); 622 int numdeg=0; 623 if (Anum.type==_POLY) 624 numdeg=Anum._POLYptr->lexsorted_degree(); 625 int dendeg=0; 626 if (Aden.type==_POLY) 627 dendeg=Aden._POLYptr->lexsorted_degree(); 628 if (numdeg>=dendeg-1) 629 return false; 630 // A must have non real roots 631 vecteur rA=singular(A,x,contextptr); 632 if (!rA.empty() && is_undef(rA)) 633 return false; 634 vecteur Pv=factors(P,x,contextptr); 635 int Pvs=int(Pv.size()/2); 636 for (int Pi=0;Pi<Pvs;++Pi){ 637 gen somme_residus=0; 638 P=Pv[2*Pi]; 639 int cur_mult=Pv[2*Pi+1].val; 640 if (!equalposcomp(lidnt(P),x)){ 641 gen resA; 642 if (!intgab(A,x,a,b,resA,contextptr)) 643 return false; 644 res += cur_mult*resA*ln(P,contextptr); 645 continue; 646 } 647 vecteur rP=singular(symbolic(at_ln,P),x,contextptr); 648 if (!rP.empty() && is_undef(rP)) 649 return false; 650 // compute constant coeff of P 651 gen P2=sym2r(P,lv,contextptr),N,D; 652 fxnd(P2,N,D); 653 gen lncst=ln(_lcoeff(gen(makevecteur(r2sym(N,lv,contextptr),x),_SEQ__VECT),contextptr)/_lcoeff(gen(makevecteur(r2sym(D,lv,contextptr),x),_SEQ__VECT),contextptr),contextptr); 654 // roots and poles of P are sorted: for im>0, contour is C- 655 // for im<=0, contour is C+, 656 // for roots of P take +residue(ln(x-r)*A) 657 // for poles of P take -residue(ln(x-r)*A) 658 int rAs=int(rA.size()),rPs=int(rP.size()); 659 for (int i=0;i<rAs;++i){ 660 gen racineA=rA[i]; 661 bool residucplus=ck_is_positive(im(racineA,contextptr),contextptr); 662 if (residucplus && !is_zero(lncst)) 663 somme_residus += lncst*residue(A,*x._IDNTptr,racineA,contextptr); 664 if (is_undef(somme_residus)) 665 return false; 666 for (int j=0;j<rPs;++j){ 667 gen racineP=rP[j]; 668 gen imracineP=im(racineP,contextptr); 669 bool racinePcplus=ck_is_strictly_positive(imracineP,contextptr); 670 if (racinePcplus){ // contour is C- for this pole/root of ln 671 if (!estreel && !residucplus){ 672 gen tmp=residue(A*ln(x-racineP,contextptr),*x._IDNTptr,racineA,contextptr); 673 if (is_undef(tmp)) 674 return false; 675 somme_residus -= tmp; 676 } 677 } 678 else { // contour is C+ for this pole/root of ln 679 if (residucplus){ 680 gen tmp=residue(A*ln(x-racineP,contextptr),*x._IDNTptr,racineA,contextptr); 681 if (is_undef(tmp)) 682 return false; 683 if (estreel && !is_zero(imracineP)) 684 tmp=2*cst_i*im(tmp,contextptr); 685 somme_residus += tmp; 686 } 687 } 688 } 689 } 690 somme_residus=normal(2*cur_mult*cst_pi*cst_i*somme_residus,contextptr); 691 res += somme_residus+resB; 692 } 693 return true; 694 } 695 if (typeint==4){ 696 // rational fraction of x and exp(A*x+B) 697 // the exp part is periodic if x -> x+2*i*pi/A 698 // if the rat frac of x and rat frac of exp are separate 699 // the problem is to find a function such that 700 // f(.+2*i*pi/A)-f(.)=ratfrac(x) 701 identificateur Xid(" X"); 702 gen X(Xid),expx(symb_exp(A*x+B)); 703 gen g1=subst(g,expx,X,false,contextptr); 704 // Separable variables? 705 vecteur f=factors(g1,x,contextptr); // Factor then split factors 706 gen xfact(plus_one),Xfact(plus_one),T(2*cst_i*cst_pi/A); 707 gen imT(im(T,contextptr)); 708 if (separate_variables(f,x,Xid,xfact,Xfact,contextptr)){ 709 // rescale xfact, in order to find a discrete antiderivative 710 gen xfactscaled=subst(xfact,x,x*T,false,contextptr),remains_to_sum,xfactint; 711 if (!rational_sum(xfactscaled,x,xfactint,remains_to_sum, 712 /* psi allowed */ true,contextptr)) 713 return false; 714 // psi function has poles at 0,-1,-2,... 715 // all have Laurent series -1/(x-pole) 716 xfactint=ratnormal(subst(xfactint,x,x/T,false,contextptr),contextptr); 717 // now int()==contour_integral of xfactint*Xfact over rectangle 718 // -inf .. + inf -> inf+T ..-inf+T -> 719 // just compute all residues at poles where im is in [0,im(T)] 720 // for xfactint, poles are the same as the poles of xfact translated 721 // for Xfact, find pole in X then take A*x+B=ln(pole in X) 722 vecteur rA=singular(xfact,x,contextptr); 723 vecteur rP=singular(Xfact,X,contextptr); 724 if (is_undef(rA) || is_undef(rP)) 725 return false; 726 gen tmp=xfactint*subst(Xfact,X,expx,false,contextptr); 727 gen somme_residus; 728 int rAs=int(rA.size()),rPs=int(rP.size()); 729 for (int i=0;i<rAs;++i){ 730 gen rac=rA[i]; 731 // adjust imaginary part 732 gen imrac=im(rac,contextptr); 733 gen k=_floor(imrac/imT,contextptr); 734 rac -= k*T; 735 if (is_positive(k,contextptr)) 736 somme_residus += residue(tmp,*x._IDNTptr,rac,contextptr); 737 if (is_undef(somme_residus)) 738 return false; 739 } 740 for (int i=0;i<rPs;++i){ 741 gen rac=(ln(rP[i],contextptr)-B)/A; 742 // adjust imaginary part between 0 and im(T) 743 gen imrac=im(rac,contextptr); 744 gen k=_floor(imrac/imT,contextptr); 745 rac -= k*T; 746 somme_residus += residue(tmp,*x._IDNTptr,rac,contextptr); 747 if (is_undef(somme_residus)) return false; 748 } 749 res=normal(-2*cst_pi*cst_i*somme_residus,contextptr); 750 // - because the integration on inf+T..-inf+T is done backward 751 return true; 752 } 753 return false; 754 } 755 if (typeint==3){ 756 // rational fraction of x and sin(A*x+B)|cos(A*x+B) or exp(i*(A*x+B)) 757 gen img=im(g,contextptr); 758 img=simplify(img,contextptr); 759 if (!is_zero(img)){ 760 gen resre,resim; 761 gen reg=simplify(re(g,contextptr),contextptr); 762 if (!intgab_sincos(reg,x,a,b,A,B,resre,contextptr)) 763 return false; 764 if (!intgab_sincos(img,x,a,b,A,B,resim,contextptr)) 765 return false; 766 res=resre+cst_i*resim; 767 return true; 768 } 769 return intgab_sincos(g,x,a,b,A,B,res,contextptr); 770 } 771 if (typeint==2){ 772 // FIXME replace by if (typeint>0) but should handle transc. func. 773 // correctly... 774 #ifndef NO_STDEXCEPT 775 try { 776 #endif 777 gen gl,glim; 778 identificateur t(" t"),r(" r"); 779 gen gt(t),gr(r),geff(g); 780 if (typeint==2){ 781 // replace g by the log part of g 782 if (intgab_ratfrac(g,x,res,contextptr)){ 783 return true; 784 } 785 gen ratpart,lnpart; 786 if (!intreduce(g,x,lnpart,ratpart,contextptr)) 787 return false; 788 gl=limit(g,*x._IDNTptr,plus_inf,1,contextptr); 789 geff=lnpart; 790 } 791 else { 792 // has limit 0 at infinity in the upper or lower half plane 793 // replace x by r*exp(i.t), assume(t in ]0,pi[ or ]-pi,0[) 794 // and look for limit(r*g) 795 glim=gr*subst(g,x,gr*symbolic(at_exp,cst_i*t),false,contextptr); 796 if (!assume_t_in_ab(gt,0,cst_pi,true,true,contextptr)) 797 return false; 798 gl=limit(glim,r,plus_inf,1,contextptr); 799 } 800 if (is_zero(gl)){ // use upper half plan 801 res=0; 802 vecteur v=singular(geff,x,contextptr); 803 if (is_undef(v)) 804 return false; 805 int s=int(v.size()),nresidue=0; 806 for (int i=0;i<s;++i){ 807 if (is_real(v[i],contextptr)){ 808 res=undef; // singularity on the real axis 809 *logptr(contextptr) << gettext("Warning: pole at ") << v[i] << '\n'; 810 purgenoassume(gt,contextptr); 811 return false; 812 } 813 if (is_positive(im(v[i],contextptr),contextptr)){ 814 nresidue++; 815 res=res+residue(geff,x,v[i],contextptr); 816 if (is_undef(res)) return false; 817 } 818 } 819 if (s && !nresidue) // e.g. int(1/(x^2+a^2),x,0,inf) 820 return false; 821 res = normal(2*cst_i*cst_pi*res,contextptr); 822 purgenoassume(gt,contextptr); 823 return true; 824 } 825 if (typeint==2){ 826 gen gl2=limit(g,*x._IDNTptr,minus_inf,1,contextptr); 827 if (is_strictly_positive(gl,contextptr) && is_strictly_positive(gl2,contextptr)) 828 res=plus_inf; 829 else { 830 if (is_strictly_positive(-gl,contextptr) && is_strictly_positive(-gl2,contextptr)) 831 res=minus_inf; 832 else 833 res=undef; 834 } 835 purgenoassume(gt,contextptr); 836 return true; 837 } 838 if (!assume_t_in_ab(gt,-cst_pi,0,true,true,contextptr)) 839 return false; 840 gl=limit(glim,r,plus_inf,1,contextptr); 841 if (is_zero(gl)){ // use lower half plan 842 res=0; 843 vecteur v=singular(g,x,contextptr); 844 if (is_undef(v)) 845 return false; 846 int s=int(v.size()); 847 for (int i=0;i<s;++i){ 848 if (is_real(v[i],contextptr)){ 849 res=undef; // singularity on the real axis 850 purgenoassume(gt,contextptr); 851 return true; 852 } 853 if (is_positive(-im(v[i],contextptr),contextptr)){ 854 res=res+residue(g,x,v[i],contextptr); 855 if (is_undef(res)) return false; 856 } 857 } 858 res = normal(2*cst_i*cst_pi*res,contextptr); 859 purgenoassume(gt,contextptr); 860 return true; 861 } 862 #ifndef NO_STDEXCEPT 863 } catch (...){ 864 last_evaled_argptr(contextptr)=NULL; 865 return false; 866 } 867 #endif 868 } 869 return false; 870 } 871 doapplyinv(const gen & g,GIAC_CONTEXT)872 gen doapplyinv(const gen &g,GIAC_CONTEXT){ 873 if (g.type!=_SYMB) 874 return inv(g,contextptr); 875 if (g._SYMBptr->sommet==at_pow && g._SYMBptr->feuille.type==_VECT && g._SYMBptr->feuille._VECTptr->size()==2) 876 return symb_pow(g._SYMBptr->feuille._VECTptr->front(),-g._SYMBptr->feuille._VECTptr->back()); 877 if (g._SYMBptr->sommet==at_prod && g._SYMBptr->feuille.type==_VECT){ 878 vecteur v=*g._SYMBptr->feuille._VECTptr; 879 for (unsigned i=0;i<v.size();++i){ 880 v[i]=doapplyinv(v[i],contextptr); 881 } 882 return _prod(v,contextptr); 883 } 884 return inv(g,contextptr); 885 } 886 applyinv(const gen & g,GIAC_CONTEXT)887 gen applyinv(const gen & g,GIAC_CONTEXT){ 888 vector<const unary_function_ptr *> inv_v(1,at_inv); 889 vector< gen_op_context > applyinv_v(1,doapplyinv); 890 return subst(g,inv_v,applyinv_v,false,contextptr); 891 } 892 intgab(const gen & g0,const gen & x,const gen & a,const gen & b,gen & res,bool nonrecursive,GIAC_CONTEXT)893 static bool intgab(const gen & g0,const gen & x,const gen & a,const gen & b,gen & res,bool nonrecursive,GIAC_CONTEXT){ 894 if (x.type!=_IDNT) 895 return false; 896 if (is_zero(g0)){ 897 res=zero; 898 return true; 899 } 900 if (a==unsigned_inf || b==unsigned_inf){ 901 *logptr(contextptr) << gettext("Please use +infinity or -infinity since infinity is unsigned") << '\n'; 902 return false; 903 } 904 if (is_strictly_greater(a,b,contextptr)){ 905 bool bo=intgab(g0,x,b,a,res,contextptr); 906 res=-res; 907 return bo; 908 } 909 if (!is_strictly_greater(b,a,contextptr)) 910 return false; 911 if (equalposcomp(lidnt(a),x) || equalposcomp(lidnt(b),x)) 912 return false; 913 vecteur subst1,subst2; 914 surd2pow(g0,subst1,subst2,contextptr); 915 gen g0_(subst(g0,subst1,subst2,false,contextptr)),g0mult(1); 916 // apply inv to pow and * 917 g0_=applyinv(g0_,contextptr); if (is_undef(g0_)) return false; 918 if (x.type==_IDNT && g0_.is_symb_of_sommet(at_prod) && g0_._SYMBptr->feuille.type==_VECT){ 919 // extract csts 920 vecteur v=*g0_._SYMBptr->feuille._VECTptr,v1,v2; 921 for (unsigned i=0;i<v.size();++i){ 922 if (depend(v[i],*x._IDNTptr)) 923 v2.push_back(v[i]); 924 else 925 v1.push_back(v[i]); 926 } 927 if (!v1.empty()){ 928 if (!intgab(_prod(v2,contextptr),x,a,b,res,contextptr)) 929 return false; 930 res=_prod(v1,contextptr)*res; 931 return true; 932 } 933 } 934 if (!is_inf(a) && !is_inf(b) && g0_.is_symb_of_sommet(at_pow)){ 935 g0_=g0_._SYMBptr->feuille; 936 if (g0_.type==_VECT && g0_._VECTptr->size()==2){ 937 gen expo=g0_._VECTptr->back(); 938 gen base=g0_._VECTptr->front(); 939 vecteur lv=lvarxwithinv(base,x,contextptr);//rlvarx(base,x); 940 if (lv.size()==1 && lv.front()==x){ 941 int na=0,nb=0; 942 for (;;){ 943 gen tmp=_quorem(makesequence(base,x-a,x),contextptr); 944 if (tmp.type==_VECT && tmp._VECTptr->size()==2 && tmp._VECTptr->back()==0){ 945 ++na; 946 base=tmp._VECTptr->front(); 947 continue; 948 } 949 tmp=_quorem(makesequence(base,b-x,x),contextptr); 950 if (tmp.type!=_VECT || tmp._VECTptr->size()!=2 || tmp._VECTptr->back()!=0) 951 break; 952 ++nb; 953 base=tmp._VECTptr->front(); 954 } 955 if (derive(base,x,contextptr)==0){ 956 g0mult=pow(base,expo,contextptr); 957 g0_=symbolic(at_pow,makesequence(x-a,na*expo))*symbolic(at_pow,makesequence(b-x,nb*expo)); 958 nb=0; // insure next test is not true 959 } 960 if (nb==1 && !na){ 961 base=g0_._VECTptr->front(); 962 gen tmp=_horner(makesequence(base,a,x),contextptr); 963 base=base-tmp; 964 for (;;){ 965 gen tmp=_quorem(makesequence(base,x-a,x),contextptr); 966 if (tmp.type==_VECT && tmp._VECTptr->size()==2 && tmp._VECTptr->back()==0){ 967 ++na; 968 base=tmp._VECTptr->front(); 969 continue; 970 } 971 break; 972 } 973 if (derive(base,x,contextptr)==0){ 974 // pow(-base,expo)*int(((b-a)^na-(x-a)^na)^expo,x,a,b) 975 // let x=a+(b-a)*t^(1/na) 976 // -> pow(-base,expo)*(b-a)^(1+na*expo)/na*int((1-t)^expo*t^(1/na-1),t,0,1) 977 res= pow(-base,expo,contextptr)*pow(b-a,1+na*expo,contextptr)*Gamma(inv(na,contextptr),contextptr)*Gamma(expo+1,contextptr)/Gamma(expo+1+inv(na,contextptr),contextptr)/na; 978 return true; 979 } 980 } // nb==1 && !na 981 } 982 } 983 } 984 if (!is_inf(a) && !is_inf(b) && g0_.is_symb_of_sommet(at_prod) && g0_._SYMBptr->feuille.type==_VECT && g0_._SYMBptr->feuille._VECTptr->size()==2){ // Beta? 985 // rewrite ^ of powers 986 vecteur v=*g0_._SYMBptr->feuille._VECTptr,v1; 987 for (unsigned i=0;i<v.size();++i){ 988 gen tmp=v[i]; 989 if (!is_zero(derive(tmp,x,contextptr))){ 990 v1.push_back(tmp); 991 } 992 } 993 if (v1.size()==2 && v1.front().is_symb_of_sommet(at_pow)&& v1.back().is_symb_of_sommet(at_pow)){ 994 gen va=v1.front()._SYMBptr->feuille,vb=v1.back()._SYMBptr->feuille; 995 if (va.type==_VECT && va._VECTptr->size()==2 && vb.type==_VECT && vb._VECTptr->size()==2){ 996 gen va1=va._VECTptr->front(),va1x,va1c,va2=va._VECTptr->back(),vb1=vb._VECTptr->front(),vb2=vb._VECTptr->back(),vb1x,vb1c; 997 if (va1.is_symb_of_sommet(at_pow) && va1._SYMBptr->feuille.type==_VECT && va1._SYMBptr->feuille._VECTptr->size()==2){ 998 va2=va1._SYMBptr->feuille._VECTptr->back()*va2; 999 va1=va1._SYMBptr->feuille._VECTptr->front(); 1000 } 1001 if (vb1.is_symb_of_sommet(at_pow) && vb1._SYMBptr->feuille.type==_VECT && vb1._SYMBptr->feuille._VECTptr->size()==2){ 1002 vb2=vb1._SYMBptr->feuille._VECTptr->back()*vb2; 1003 vb1=vb1._SYMBptr->feuille._VECTptr->front(); 1004 } 1005 if (is_linear_wrt(va1,x,va1x,va1c,contextptr) && is_linear_wrt(vb1,x,vb1x,vb1c,contextptr)){ 1006 if (is_zero(recursive_normal(va1c/va1x+a,contextptr),contextptr) && 1007 is_zero(recursive_normal(vb1c/vb1x+b,contextptr),contextptr)){ 1008 // int( (va1x*(x-a))^va2*(vb1x*(x-b))^vb2,a,b) 1009 res=g0mult*pow(va1x,va2,contextptr)*pow(-vb1x,vb2,contextptr)*pow(b-a,va2+vb2+1,contextptr)*Beta(va2+1,vb2+1,contextptr); 1010 return true; 1011 } 1012 if (is_zero(recursive_normal(va1c/va1x+b,contextptr),contextptr) && 1013 is_zero(recursive_normal(vb1c/vb1x+a,contextptr),contextptr)){ 1014 // int( (va1x*(x-b))^va2*(vb1x*(x-a))^vb2,a,b) 1015 res=g0mult*pow(-va1x,va2,contextptr)*pow(vb1x,vb2,contextptr)*pow(b-a,va2+vb2+1,contextptr)*Beta(va2+1,vb2+1,contextptr); 1016 return true; 1017 } 1018 } 1019 } 1020 } 1021 } 1022 // detect Dirac 1023 vecteur v=lop(g0,at_Dirac); 1024 if (!v.empty()){ 1025 gen A,B,a0,b0; 1026 #ifdef GIAC_HAS_STO_38 1027 identificateur t("tsumab_"); 1028 #else 1029 identificateur t(" tsumab"); 1030 #endif 1031 gen h=quotesubst(g0,v.front(),t,contextptr); 1032 if (!is_linear_wrt(h,t,A,B,contextptr)) 1033 return false; 1034 gen heav=v.front()._SYMBptr->feuille; 1035 if (heav.type==_VECT && heav._VECTptr->size()==2 && heav._VECTptr->back().type==_INT_ ){ 1036 int diracorder=heav._VECTptr->back().val; 1037 if (diracorder<0){ 1038 *logptr(contextptr) << gettext("Negative second Dirac argument") << '\n'; 1039 return false; 1040 } 1041 A=derive(A,x,diracorder,contextptr); 1042 if (is_undef(A)) 1043 return false; 1044 if (diracorder%2) 1045 A=-A; 1046 heav=heav._VECTptr->front(); 1047 } 1048 if (!is_linear_wrt(heav,x,a0,b0,contextptr) || is_zero(a0)) 1049 return false; 1050 if (!intgab(B,x,a,b,res,contextptr)) 1051 return false; 1052 gen c=-b0/a0; 1053 if (ck_is_greater(c,a,contextptr) && ck_is_greater(b,c,contextptr)) 1054 res += quotesubst(A,x,c,contextptr); 1055 else 1056 *logptr(contextptr) << gettext("Warning, Dirac function outside summation interval") << '\n'; 1057 return true; 1058 } 1059 if (a==b){ 1060 res=0; 1061 return true; 1062 } 1063 gen g=hyp2exp(g0,contextptr); 1064 vecteur lvarg = lvar(g); 1065 bool rational = lvarg==vecteur(1,x); 1066 if (!rational){ 1067 int s1=nvars_depend_x(loptab(g,sincostan_tab),x); 1068 // rewrite cos/sin/tan if more than 1 available, 1069 // do not rewrite atan/asin/acos 1070 if (s1) // check added otherwise int(1/(x-a)^999,x,a-1,a+1) takes forever 1071 g=tsimplify_noexpln(g,s1,0,contextptr); 1072 } 1073 // FIXME should check integrability at -/+inf 1074 if (a==minus_inf){ 1075 gen A=limit(g,*x._IDNTptr,a,1,contextptr); 1076 if (!is_zero(A) && b!=plus_inf){ 1077 res=-a*A; 1078 return !is_undef(res); // true; 1079 } 1080 if (b==plus_inf){ 1081 vecteur singu=find_singularities(g,*x._IDNTptr,0 /* real singularities*/,contextptr); 1082 if (!singu.empty()){ 1083 *logptr(contextptr) << "Warning, singularities at " << singu << '\n'; 1084 if (calc_mode(contextptr)==1 || abs_calc_mode(contextptr)==38){ 1085 res=undef; 1086 return true; 1087 } 1088 } 1089 gen B=limit(g,*x._IDNTptr,b,-1,contextptr); 1090 if (!is_zero(B)){ 1091 if (is_zero(A)) 1092 res=b*B; 1093 else 1094 res=b*B-a*A; 1095 return !is_undef(res); // true; 1096 } 1097 int ieo=is_even_odd(g,x,contextptr); 1098 if (ieo==2){ 1099 res=0; 1100 return true; 1101 } 1102 if (is_zero(A) && intgab_r(g,x,a,b,rational,res,contextptr)) 1103 return true; 1104 if (ieo==1){ 1105 // simplify g on 0..inf 1106 assumesymbolic(symb_superieur_egal(x,0),0,contextptr); 1107 gen g1=eval(g,1,contextptr); 1108 purgenoassume(x,contextptr); 1109 if (g1!=eval(g,1,contextptr)){ 1110 res=2*_integrate(makesequence(g1,x,0,plus_inf),contextptr); 1111 return true; 1112 } 1113 } 1114 return false; 1115 } // end b==plus_inf (a is still minus_inf) 1116 // subst x by x+b, check parity: even -> 1/2 int(-inf,+inf) 1117 gen gb=subst(g,x,x+b,false,contextptr); 1118 int eo=is_even_odd(gb,x,contextptr); 1119 if (eo==1){ 1120 if ( (rational && intgab_ratfrac(gb,x,res,contextptr)) || 1121 intgab(gb,x,a,plus_inf,res,contextptr) ){ 1122 if (!is_inf(res)) 1123 res=ratnormal(res/2,contextptr); 1124 return true; 1125 } 1126 } 1127 vecteur v; 1128 rlvarx(gb,x,v); 1129 int vs=int(v.size()); 1130 for (int i=0;i<vs;++i){ 1131 if (v[i].is_symb_of_sommet(at_ln)){ 1132 gen f=v[i]._SYMBptr->feuille,a,b; 1133 // if f is a*x make the change of var x=-exp(t) 1134 if (is_linear_wrt(f,x,a,b,contextptr) && is_zero(b)){ 1135 vecteur vin=makevecteur(v[i],x); 1136 vecteur vout=makevecteur(ln(-a,contextptr)+x,-exp(x,contextptr)); 1137 gb=quotesubst(gb,vin,vout,contextptr)*exp(x,contextptr); 1138 return intgab(gb,x,minus_inf,plus_inf,res,contextptr); 1139 } 1140 } 1141 } 1142 return false; 1143 } // end a==minus_inf 1144 if (b==plus_inf){ 1145 gen ga_orig=subst(g,x,x+a,false,contextptr),ga(ga_orig); 1146 // additional check for int(t^n/(exp(alpha*t)-1),t,0,inf)=n!/alpha^(n+1)*Zeta(n+1) 1147 vecteur vax=rlvarx(ga,x); 1148 if (vax.size()==2 && vax.front()==x && vax.back().is_symb_of_sommet(at_exp)){ 1149 gen expo=vax.back(),expo_a,expo_b; 1150 if (is_linear_wrt(expo._SYMBptr->feuille,x,expo_a,expo_b,contextptr) && is_strictly_positive(expo_a,contextptr)){ 1151 gen gand=_fxnd(ga,contextptr); 1152 if (gand.type==_VECT && gand._VECTptr->size()==2){ 1153 gen gan=gand._VECTptr->front(),gad=gand._VECTptr->back(),gad_a,gad_b; 1154 // gad must be a power of expo-exp(expo_b), starting with linear 1155 if (rlvarx(gan,x).size()==1 && is_linear_wrt(gad,expo,gad_a,gad_b,contextptr)){ 1156 // gad=gad_a*(expo+gad_b/gad_a) 1157 gen test=ratnormal(gad_a*exp(expo_b,contextptr)+gad_b,contextptr); 1158 if (is_zero(test)){ 1159 // -1/gad_b*int(gan/(exp(expo_a*t)-1),t,0,inf) 1160 gen ganv=_coeff(makesequence(gan,x),contextptr); 1161 if (ganv.type==_VECT && !ganv._VECTptr->empty() && is_zero(ganv._VECTptr->back())){ 1162 res=0; 1163 vecteur v=*ganv._VECTptr; 1164 gen facti=pow(expo_a,-2,contextptr); 1165 for (int i=1;i<int(v.size());++i){ 1166 gen coeff=v[v.size()-i-1]; 1167 if (!is_zero(coeff)) 1168 res += coeff*facti*Zeta(i+1,contextptr); 1169 facti=facti*gen(i+1)/expo_a; 1170 } 1171 res=ratnormal(-res/gad_b,contextptr); 1172 return true; 1173 } 1174 } // end if is_zero(test) 1175 } // end if (rlvarx(gan,x).size()==1 1176 } // end if gand.type==_VECT 1177 identificateur t(" tintgab"); 1178 gen y(t); 1179 ga=subst(ga,expo,y,false,contextptr); 1180 vecteur f=factors(ga,x,contextptr); // Factor then split factors 1181 gen xfact(plus_one),yfact(plus_one); 1182 if (separate_variables(f,x,y,xfact,yfact,contextptr)){ 1183 // yfact must be a fraction with denominator a power of expo-exp(expo_b) 1184 gen y0=exp(expo_b,contextptr); 1185 gen expofact=_fxnd(yfact,contextptr); 1186 if (expofact.type==_VECT && expofact._VECTptr->size()==2){ 1187 gen exponum=expofact._VECTptr->front(),expoden=expofact._VECTptr->back(); 1188 int n=_degree(makesequence(expoden,y),contextptr).val; 1189 gen coeffden=ratnormal(expoden/pow(y-y0,n,contextptr),contextptr); 1190 if (n>1 && is_zero(derive(coeffden,y,contextptr))){ 1191 // 1/expoden*xfact*exponum(y)/(y-1)^n, y'=expo_a*y 1192 gen additional=_quorem(makesequence(exponum,pow(y-1,n-1),y),contextptr); 1193 exponum=additional[1]; 1194 gen xfactc=_coeff(makesequence(xfact,x),contextptr); 1195 if (xfactc.type==_VECT && !xfactc._VECTptr->empty()){ 1196 vecteur & xfactv=*xfactc._VECTptr; 1197 // check cancellation of xfact at 0 at least order n 1198 bool check=true; 1199 for (int i=1;i<=n;++i){ 1200 if (xfactv[xfactv.size()-i]!=0){ 1201 check=false; 1202 break; 1203 } 1204 } 1205 if (check){ 1206 // reduce degree by integration by part 1207 // int(P(t)*Q(e^at)/(e^at-1)^n= 1208 // int( (1/(n-1)*P'/a-P)*Q(e^at)+1/(n-1)P*Q'(e^at))/(e^at-1)^(n-1) 1209 gen int1=(derive(xfact,x,contextptr)/((n-1)*expo_a)-xfact)*subst(exponum/pow(y-1,n-1,contextptr),y,exp(expo_a*x,contextptr),false,contextptr); 1210 int1=_integrate(makesequence(int1,x,0,plus_inf),contextptr); 1211 gen int2=xfact/gen(n-1)*subst(derive(exponum,y,contextptr)/pow(y-1,n-1,contextptr),y,exp(expo_a*x,contextptr),false,contextptr); 1212 int2=_integrate(makesequence(int2,x,0,plus_inf),contextptr); 1213 gen int3=xfact*subst(additional[0]/(y-1),y,exp(expo_a*x,contextptr),false,contextptr); 1214 int3=_integrate(makesequence(int3,x,0,plus_inf),contextptr); 1215 res=(int1+int2+int3)/coeffden; 1216 return true; 1217 } 1218 } 1219 } 1220 } 1221 } 1222 } // end if (is_linear_wrt(expo...)) 1223 } // end varx.size()==2 1224 ga=ga_orig; 1225 int eo=is_even_odd(ga,x,contextptr); 1226 if (eo==1){ 1227 vecteur singu=find_singularities(g,*x._IDNTptr,0 /* real singularities*/,contextptr); 1228 if (singu.empty()){ 1229 if ( (rational && intgab_ratfrac(ga,x,res,contextptr)) || 1230 intgab(ga,x,minus_inf,plus_inf,res,contextptr) ){ 1231 if (!is_inf(res)) 1232 res=ratnormal(res/2,contextptr); 1233 return !is_undef(res); 1234 } 1235 } 1236 } 1237 vecteur v; 1238 rlvarx(ga,x,v); 1239 int vs=int(v.size()); 1240 for (int i=0;i<vs;++i){ 1241 if (v[i].is_symb_of_sommet(at_ln)){ 1242 gen f=v[i]._SYMBptr->feuille,a,b; 1243 // if f is a*x make the change of var x=exp(t) 1244 if (is_linear_wrt(f,x,a,b,contextptr) && is_zero(b)){ 1245 vecteur vin=makevecteur(v[i],x); 1246 vecteur vout=makevecteur(ln(a,contextptr)+x,exp(x,contextptr)); 1247 ga=quotesubst(ga,vin,vout,contextptr)*exp(x,contextptr); 1248 return intgab(ga,x,minus_inf,plus_inf,res,contextptr); 1249 } 1250 } 1251 } 1252 return false; 1253 } 1254 gen gab=subst(g0,x,b,false,contextptr)-subst(g0,x,a,false,contextptr); 1255 gen gabd; 1256 if (!has_evalf(gab,gabd,1,contextptr) || is_zero(gabd)) 1257 gab=simplify(gab,contextptr); 1258 gen gm=subst(g0,x,b,false,contextptr)+subst(g0,x,a,false,contextptr); 1259 if (!has_evalf(gm,gabd,1,contextptr) || is_zero(gabd)) 1260 gm=simplify(gm,contextptr); 1261 if (is_constant_wrt(g,x,contextptr)){ 1262 if (contains(g,x)) 1263 g=ratnormal(g,contextptr); 1264 res=g*(b-a); 1265 return true; 1266 } 1267 int eo=0; 1268 if (is_zero(gab) || is_zero(gm) ){ 1269 identificateur t("tintgab_"); 1270 gen tt(t); 1271 gm=subst(g0,x,tt+(a+b)/2,false,contextptr); 1272 eo=is_even_odd(gm,tt,contextptr); 1273 } 1274 if (!nonrecursive && eo==1){ 1275 if (!intgab(g0,x,a,(a+b)/2,res,true,contextptr)) 1276 return false; 1277 res=2*res; 1278 return true; 1279 } 1280 if (eo==2){ 1281 #if 0 // set to 1 if you want to check for singularities before returning 0 1282 vecteur sp=find_singularities(g,*x._IDNTptr,false,contextptr); 1283 for (int i=0;i<sp.size();++i){ 1284 if (is_greater(sp[i],a,contextptr) && is_greater(b,sp[i],contextptr)) 1285 return false; 1286 } 1287 #endif 1288 res=0; 1289 return true; 1290 } 1291 // check periodicity 1292 gm=gab; 1293 if (is_zero(gm)){ 1294 // do it only if g0 depends on potentially periodical functions exp/sin/cos/tan 1295 vecteur vx=lvarx(g0,x); 1296 for (unsigned i=0;i<vx.size();++i){ 1297 if (vx[i].type!=_SYMB || (vx[i]._SYMBptr->sommet!=at_exp && vx[i]._SYMBptr->sommet!=at_sin && vx[i]._SYMBptr->sommet!=at_cos && vx[i]._SYMBptr->sommet!=at_tan)) 1298 gm=1; 1299 } 1300 if (is_zero(gm)){ 1301 gm=subst(g0,x,x+(b-a),false,contextptr); 1302 gm=simplify(gm-g0,contextptr); 1303 } 1304 } 1305 #ifndef NO_STDEXCEPT 1306 try { 1307 #endif 1308 if (is_zero(gm)){ 1309 // try to rewrite g as a function of exp(2*i*pi*x/(b-a)) 1310 g=_lin(trig2exp(g0,contextptr),contextptr); 1311 vecteur v; 1312 rlvarx(g,x,v); 1313 islesscomplexthanf_sort(v.begin(),v.end()); 1314 int i,s=int(v.size()); 1315 if (s>=2){ 1316 gen v0,alpha,beta,alphacur,betacur,gof,periode,periodecur; 1317 for (i=0;i<s;++i){ 1318 if (v[i].is_symb_of_sommet(at_exp)){ 1319 v0=v[i]; 1320 gen v0arg=v0._SYMBptr->feuille; 1321 if (is_linear_wrt(v0arg,x,alphacur,betacur,contextptr) && is_integer( (periodecur=normal(alphacur*(b-a)/cst_two_pi/cst_i,contextptr)) )){ 1322 periode=gcd(periode,periodecur,contextptr); 1323 } 1324 } 1325 } 1326 if (!is_zero(periode)){ 1327 alpha=normal(periode*cst_two_pi/(b-a)*cst_i,contextptr); 1328 if (is_zero(re(alpha,contextptr))){ 1329 beta=normal(betacur*alpha/alphacur,contextptr); 1330 gen radius=exp(re(beta,contextptr),contextptr); 1331 // vO=exp(alpha*x+beta) -> x=(ln(v0)-beta)/alpha 1332 vecteur vin=makevecteur(x); 1333 vecteur vout=makevecteur((ln(x,contextptr)-beta)/alpha); 1334 // check for essential singularities 1335 vecteur v2=lop(recursive_normal(rlvarx(subst(v,vin,vout,false,contextptr),x),contextptr),at_exp); 1336 vecteur w2=singular(v2,x,contextptr); 1337 unsigned w2i=0; 1338 for (;w2i<w2.size();++w2i){ 1339 if (is_greater(radius,w2[w2i],contextptr)) 1340 break; 1341 } 1342 bool changesign=false; 1343 if (w2i!=w2.size()){ 1344 changesign=true; 1345 // try again after doing alpha->-alpha 1346 alpha=-alpha; 1347 beta=normal(betacur*alpha/alphacur,contextptr); 1348 radius=exp(re(beta,contextptr),contextptr); 1349 // vO=exp(alpha*x+beta) -> x=(ln(v0)-beta)/alpha 1350 vin=makevecteur(x); 1351 vout=makevecteur((ln(x,contextptr)-beta)/alpha); 1352 // check for essential singularities 1353 v2=lop(recursive_normal(rlvarx(subst(v,vin,vout,false,contextptr),x),contextptr),at_exp); 1354 w2=singular(v2,x,contextptr); 1355 w2i=0; 1356 for (;w2i<w2.size();++w2i){ 1357 if (is_greater(radius,w2[w2i],contextptr)) 1358 break; 1359 } 1360 if (w2i!=w2.size()) 1361 return false; 1362 } 1363 gof=quotesubst(g,vin,vout,contextptr); 1364 gof=_exp2pow(gof/x,contextptr); 1365 // bool b=complex_mode(contextptr); 1366 // complex_mode(true,contextptr); 1367 gof=recursive_normal(gof,contextptr); 1368 // complex_mode(b,contextptr); 1369 // Sucess! Now integration of gof on the unit circle using residues 1370 if (debug_infolevel) 1371 *logptr(contextptr) << gettext("Searching int of ") << gof << gettext(" where ") << x << gettext(" is on the unit circle, using residues") << '\n'; 1372 // replace x by another variable because we might have assumptions on x 1373 identificateur tmpid("_intgab38"); 1374 vecteur w=singular(subst(gof,x,tmpid,false,contextptr),tmpid,contextptr); 1375 if (is_undef(w)) 1376 return false; 1377 int s=int(w.size()); 1378 gen somme_residues=0; 1379 for (int i=0;i<s;++i){ 1380 #ifdef TIMEOUT 1381 control_c(); 1382 #endif 1383 if (ctrl_c || interrupted) 1384 return false; 1385 gen wabs=normal(w[i]*conj(w[i],contextptr),contextptr); 1386 if ((wabs==radius) 1387 // && !is_zero(tmpresidue) // commented otherwise int(1/cos(x)^2,x,0,pi) returns 0 1388 ){ 1389 res = unsigned_inf; 1390 return true; 1391 } 1392 if (is_strictly_greater(radius,wabs,contextptr)){ 1393 gen tmpresidue=residue(gof,x,w[i],contextptr); 1394 if (is_undef(tmpresidue)) 1395 return false; 1396 somme_residues = somme_residues + tmpresidue; // was residue(gof,x,w[i],contextptr) but no reason to compute it twice... 1397 if (is_undef(somme_residues)) 1398 return false; 1399 } 1400 } 1401 res = ratnormal(normal(periode*cst_two_pi/alpha*cst_i,contextptr)*somme_residues,contextptr); 1402 if (changesign) 1403 res=-res; 1404 return true; 1405 } // end if(is_zero(re(alpha))) 1406 } 1407 } 1408 } 1409 #ifndef NO_STDEXCEPT 1410 } catch (std::runtime_error & ){ 1411 last_evaled_argptr(contextptr)=NULL; 1412 } 1413 #endif 1414 return false; 1415 } 1416 1417 // if true put int(g,x=a..b) into res intgab(const gen & g0,const gen & x,const gen & a,const gen & b,gen & res,GIAC_CONTEXT)1418 bool intgab(const gen & g0,const gen & x,const gen & a,const gen & b,gen & res,GIAC_CONTEXT){ 1419 return intgab(g0,x,a,b,res,false,contextptr); 1420 } 1421 quotesubstcheck(const gen & g,const gen & x,const gen & i,const vecteur & v,GIAC_CONTEXT)1422 static gen quotesubstcheck(const gen & g,const gen & x,const gen & i,const vecteur & v,GIAC_CONTEXT){ 1423 if (!equalposcomp(v,i)) 1424 return quotesubst(g,x,i,contextptr); 1425 return 0; 1426 } 1427 1428 // from csturm.cc 1429 vecteur crationalroot(polynome & p,bool complexe); 1430 1431 // check whether P has only integer roots is_admissible_poly(const polynome & P,int & deg,polynome & lcoeff,vecteur & roots,GIAC_CONTEXT)1432 bool is_admissible_poly(const polynome & P,int & deg,polynome & lcoeff,vecteur & roots,GIAC_CONTEXT){ 1433 lcoeff=Tfirstcoeff(P); 1434 index_t degs=P.degree(); 1435 deg=degs[0]; 1436 for (unsigned int i=1;i<degs.size();++i){ 1437 if (degs[i]) 1438 return false; 1439 } 1440 polynome PP=poly12polynome(polynome2poly1(P)); 1441 polynome P1=PP.derivative(); 1442 if (gcd(PP,P1).degree(0)>0) 1443 return false; 1444 roots.clear(); 1445 if (deg<1) 1446 return true; 1447 roots=crationalroot(PP,false); 1448 roots=*_sort(roots,contextptr)._VECTptr; 1449 if (int(roots.size())!=deg) 1450 return false; 1451 return true; 1452 } 1453 1454 // tmp1=sum_k a_k x^k, find sum_k a_k/s_k x^k 1455 // where s_x=product(k-decals[i],i) and decals[i] is rationnal 1456 // if decals[i] is an integer, multiply by x^(-1-decals[i]), int 1457 // and mult by x^decals[i] in_sumab_int(gen & tmp1,const gen & gx,const vecteur & decals,const gen & lcoeff,GIAC_CONTEXT)1458 static bool in_sumab_int(gen & tmp1,const gen & gx,const vecteur & decals,const gen & lcoeff,GIAC_CONTEXT){ 1459 int nstep=int(decals.size()); 1460 gen coeff=lcoeff; 1461 gen remains; 1462 for (int i=0;i<nstep;++i){ 1463 if (decals[i].type==_FRAC){ 1464 gen n=decals[i]._FRACptr->num; 1465 gen d=decals[i]._FRACptr->den; 1466 // coeff*(k-n/d)=coeff/d*(d*k-n) 1467 coeff = coeff/d; 1468 // sum a_k/(d*k-n)*gx^k 1469 // set gx=X^d : sum_ a_k/(d*k-n)*X^(d*k) 1470 tmp1=subst(tmp1,gx,pow(gx,d,contextptr),false,contextptr); 1471 tmp1=tmp1*pow(gx,-1-n,contextptr); 1472 tmp1=integrate_id_rem(tmp1,gx,remains,contextptr,0); 1473 if (is_undef(tmp1)) return false; 1474 tmp1=tmp1-limit(tmp1,*gx._IDNTptr,0,1,contextptr); 1475 if (is_inf(tmp1)) return false; // for sum(1/((n+1)*(2*n-1)),n,0,inf); 1476 tmp1=ratnormal(tmp1*pow(gx,n,contextptr),contextptr); 1477 tmp1=ratnormal(subst(tmp1,gx,pow(gx,inv(d,contextptr),contextptr),false,contextptr),contextptr); 1478 } 1479 else { 1480 tmp1=tmp1*pow(gx,-1-decals[i],contextptr); 1481 tmp1=integrate_id_rem(tmp1,gx,remains,contextptr,0); 1482 if (is_undef(tmp1)) return false; 1483 tmp1=tmp1-limit(tmp1,*gx._IDNTptr,0,1,contextptr); 1484 tmp1=tmp1*pow(gx,decals[i],contextptr); 1485 } 1486 if (!is_zero(remains)) 1487 return false; 1488 } 1489 tmp1=tmp1/coeff; 1490 return true; 1491 // do_lnabs(b,contextptr); 1492 } 1493 sumab_int(gen & tmp1,const gen & gx,const vecteur & decals,const gen & lcoeff,GIAC_CONTEXT)1494 static bool sumab_int(gen & tmp1,const gen & gx,const vecteur & decals,const gen & lcoeff,GIAC_CONTEXT){ 1495 bool b=do_lnabs(contextptr); 1496 do_lnabs(false,contextptr); 1497 // bool c=complex_mode(contextptr); 1498 // complex_mode(true,contextptr); 1499 bool bres=in_sumab_int(tmp1,gx,decals,lcoeff,contextptr); 1500 do_lnabs(b,contextptr); 1501 // complex_mode(c,contextptr); 1502 return bres; 1503 } 1504 sumab_ps(const polynome & Q,const polynome & R,const vecteur & v,const gen & a,const gen & x,const gen & g,bool est_reel,const polynome & p,const polynome & s,gen & res,GIAC_CONTEXT)1505 static bool sumab_ps(const polynome & Q,const polynome & R,const vecteur & v,const gen & a,const gen & x,const gen & g,bool est_reel,const polynome & p,const polynome & s,gen & res,GIAC_CONTEXT){ 1506 // p corresponds to derivation, s to integration 1507 // cerr << "p=" << p << " s=" << s << " Q=" << Q << " R=" << R << '\n'; 1508 // Q must be independant of x 1509 // If R is independant of x we use the geometric series 1510 // If R=x-integer the exponential (must change bounds by integer) 1511 // If R=2x(2x+1) sinh/cosh etc. 1512 if (Q.degree(0)==0){ 1513 // count "integrations" step in s 1514 int intstep; 1515 vecteur decals; 1516 polynome lcoeffs; 1517 if (is_admissible_poly(s,intstep,lcoeffs,decals,contextptr)){ 1518 gen lcoeff=r2e(lcoeffs,v,contextptr); 1519 #ifdef GIAC_HAS_STO_38 1520 identificateur idx("sumw_"); // identificateur idx(" x"); // 1521 #else 1522 identificateur idx(" sumw"); // identificateur idx(" x"); // 1523 #endif 1524 gen gx(idx); // ("` sumw`",contextptr); 1525 // parser instead of temporary otherwise bug with a:=1; ZT(f,z):=sum(f(n)/z^n,n,0,inf); ZT(k->c^k,z); ZT; 1526 // otherwise while purge(gx) happens, the string `sumw` is destroyed 1527 // and the global map is not sorted correctly anymore 1528 if (!assume_t_in_ab(gx,0,1,true,true,contextptr)) 1529 return false; 1530 // R must be the product of degree(R) consecutive terms 1531 int r=R.degree(0); 1532 if (r==1){ 1533 vecteur Rv=iroots(R); 1534 if (Rv.size()!=1){ 1535 purgenoassume(gx,contextptr); 1536 return false; 1537 } 1538 gen R0=Rv[0]; 1539 if (is_strictly_greater(R0,a,contextptr)){ 1540 res=undef; 1541 purgenoassume(gx,contextptr); 1542 return true; 1543 } 1544 // Q=Q/R.coord.front().value; 1545 index_t ind=R.coord.front().index.iref(); 1546 ind[0]=0; 1547 gen Qg=r2e(Q,v,contextptr)/r2e(polynome(monomial<gen>(R.coord.front().value,ind)),v,contextptr); 1548 // (g|x=a)*s(a)/p(a)/Q^(a-R0)/(a-R0)!*sum(p(n)/s(n)*Q^(n-R0)/(n-R0)!,n=a..inf); 1549 // first compute sum(p(n)*Q^(n-R0)/(n-R0)!,n=R0..inf) 1550 // = sum(p(n+R0)*Q^(n)/n!,n=0..inf) 1551 int d=p.degree(0); 1552 gen Pg=r2e(p,v,contextptr); 1553 vecteur vx(d+1),vy(d+1); 1554 for (int i=0;i<=d;++i){ 1555 vx[i]=i; 1556 vy[i]=ratnormal(subst(Pg,x,R0+i,false,contextptr),contextptr); 1557 } 1558 vecteur w=divided_differences(vx,vy); 1559 // p(n+R0)=w[0]+w[1]*n+w[2]*n*(n-1)+... 1560 // hence the sum is exp(Q)*(w[0]+w[1]*Q+...) 1561 reverse(w.begin(),w.end()); 1562 gen tmp1=symb_horner(w,gx)*exp(gx,contextptr),remains; 1563 // substract sum(p(n)*Q^(n-R0)/(n-R0)!,n=R0..a-1) 1564 for (int n=R0.val;n<a.val;++n){ 1565 tmp1 -= subst(Pg,x,n,false,contextptr)*pow(Qg,n-R0.val)/factorial(n-R0.val); 1566 } 1567 // Integration step 1568 decals = subvecteur(decals,vecteur(decals.size(),R0)); 1569 if (sumab_int(tmp1,gx,decals,lcoeff,contextptr)){ 1570 gen coeffa=limit(g*r2e(s,v,contextptr)/r2e(p,v,contextptr)*factorial(a.val-R0.val),*x._IDNTptr,a,1,contextptr); 1571 res=limit(coeffa*tmp1,*gx._IDNTptr,Qg,1,contextptr)/pow(Qg,a-R0,contextptr); 1572 if (est_reel) 1573 res=re(res,contextptr); 1574 res=ratnormal(res,contextptr); 1575 purgenoassume(gx,contextptr); 1576 return true; 1577 } 1578 } 1579 if (r>=2){ 1580 polynome Rc=lgcd(R); 1581 vecteur Rv=polynome2poly1(R/Rc,1); 1582 // Rv should be a multiple of (r*x-R0)*(r*x-(R0+1))*... 1583 // -Rv[1]/Rv[0]= sum of roots = r*R0 + sum(j,j=0..r-1) 1584 // R0 = -Rv[1]/Rv[0] - (r-1)/2 1585 gen R0=-Rv[1]/Rv[0]-gen(r-1)/gen(2); 1586 if (R0.type!=_INT_){ 1587 purgenoassume(gx,contextptr); 1588 return false; 1589 } 1590 // check that Rv = cst*product(r*x-(R0+j),j=0..r-1) 1591 vecteur test(1,1); 1592 for (int j=0;j<r;++j){ 1593 test=operator_times(test,makevecteur(r,-(R0+j)),0); 1594 } 1595 if (Rv[0]*test!=test[0]*Rv){ 1596 purgenoassume(gx,contextptr); 1597 return false; 1598 } 1599 int r0=R0.val; 1600 if (r0>r*a.val){ 1601 res=undef; 1602 purgenoassume(gx,contextptr); 1603 return true; 1604 } 1605 Rc=Rv[0]/pow(gen(r),r)*Rc; 1606 gen Qg=r2e(Q,v,contextptr)/r2e(Rc,v,contextptr); 1607 if (r==2) 1608 Qg=sqrt(Qg,contextptr); 1609 else 1610 Qg=pow(Qg,inv(gen(r),contextptr),contextptr); 1611 // (g|x=a)*s(a)/p(a)/Qg^(r*a-R0)*(r*a-R0)!*sum(p(n)/s(n)*Qg^(r*n-R0)/(r*n-R0)!,n=a..inf); 1612 gen coeffa=g*r2e(s,v,contextptr)/r2e(p,v,contextptr)/pow(Qg,r*a-R0,contextptr)*factorial(r*a.val-r0); 1613 coeffa=limit(coeffa,*x._IDNTptr,a,1,contextptr); 1614 // Set k=r*n-R0 and compute sum((p/s)((k+R0)/r)*X^k/k!,k=r*a-R0..inf) 1615 int d=p.degree(0); 1616 gen Pg=r2e(p,v,contextptr); 1617 vecteur vx(d+1),vy(d+1); 1618 for (int i=0;i<=d;++i){ 1619 vx[i]=i; 1620 vy[i]=ratnormal(subst(Pg,x,(i+R0)/r,false,contextptr),contextptr); 1621 } 1622 vecteur w=divided_differences(vx,vy); 1623 reverse(w.begin(),w.end()); 1624 gen tmp=symb_horner(w,gx)*exp(gx,contextptr),remains; 1625 // substract sum(...,k=0..r*a-R0-1) 1626 for (int k=0;k<r*a.val-r0;++k){ 1627 // pow(Qg,k) replaced by pow(gx,k) for assume(x>0);somme(x^(4n+1)/(4n+1)!,n,1,inf); 1628 tmp -= subst(Pg,x,(k+R0)/r,false,contextptr)*pow(gx,k)/factorial(k); 1629 } 1630 // keep terms which are = -R0 mod r = N 1631 // for example if r=2 and R0 even, keep even terms 1632 // that is (f(X)+f(-X))/2 1633 // more generally take 1634 // 1/r*sum(f(X*exp(2i pi*k/r))*exp(-2i pi*k*N/r),k=0..r-1) 1635 int N= -r0 % r; 1636 gen tmp1=0,tmpadd; 1637 for (int k=0;k<r;++k){ 1638 tmpadd = subst(tmp,gx,gx*exp(2*k*cst_i*cst_pi/r,contextptr),false,contextptr)*exp(-2*k*N*cst_i*cst_pi/r,contextptr); 1639 tmp1 += tmpadd; 1640 } 1641 tmp1=tmp1/r; 1642 // Integration step 1643 if (sumab_int(tmp1,gx,decals,lcoeff,contextptr)){ 1644 res=limit(coeffa*tmp1,*gx._IDNTptr,Qg,1,contextptr); 1645 if (est_reel) 1646 res=re(res,contextptr); 1647 res=ratnormal(res,contextptr); 1648 purgenoassume(gx,contextptr); 1649 return true; 1650 } 1651 } 1652 if (!r){ 1653 // (g|x=a)*s(a)/p(a)/Q^a*sum(p(n)/s(n)*X^n)|X=Q, will work 1654 // first compute sum(p(n)*X^n) 1655 // then multiply by X^sdecal and integrate intstep times 1656 // then subst X by Q 1657 gen tmp=r2e(p,v,contextptr)*pow(gx,x,contextptr); 1658 gen remains,tmp1=sum(tmp,x,remains,contextptr); 1659 if (!is_zero(remains) || is_undef(tmp1)){ 1660 *logptr(contextptr) << gettext("Unable to sum ")+remains.print(contextptr) << '\n'; 1661 return false; 1662 } 1663 tmp1=-subst(tmp1,x,a,false,contextptr); 1664 // bool b=do_lnabs(contextptr); 1665 // do_lnabs(false,contextptr); 1666 if (sumab_int(tmp1,gx,decals,lcoeff,contextptr)){ 1667 gen Qg=r2e(Q,v,contextptr)/r2e(R,v,contextptr); 1668 gen coeffa=limit(g*r2e(s,v,contextptr)/r2e(p,v,contextptr),*x._IDNTptr,a,1,contextptr); 1669 res=limit(coeffa*tmp1,*gx._IDNTptr,Qg,1,contextptr)/pow(Qg,a,contextptr); 1670 if (est_reel) 1671 res=re(res,contextptr); 1672 res=ratnormal(res,contextptr); 1673 purgenoassume(gx,contextptr); 1674 return true; 1675 } 1676 } 1677 } 1678 } 1679 return false; 1680 } 1681 in_sumab(const gen & g,const gen & x,const gen & a_orig,const gen & b_orig,gen & res,bool testi,bool dopartfrac,GIAC_CONTEXT)1682 static bool in_sumab(const gen & g,const gen & x,const gen & a_orig,const gen & b_orig,gen & res,bool testi,bool dopartfrac,GIAC_CONTEXT){ 1683 //grad 1684 if (x.type!=_IDNT || !angle_radian(contextptr)) //if not in radians, exit 1685 return false; 1686 if (is_zero(g)){ 1687 res=zero; 1688 return true; 1689 } 1690 if (dopartfrac){ 1691 gen gp=_partfrac(gen(makevecteur(g,*x._IDNTptr),_SEQ__VECT),contextptr); 1692 if (gp.is_symb_of_sommet(at_plus) && gp._SYMBptr->feuille.type==_VECT){ 1693 vecteur vp=*gp._SYMBptr->feuille._VECTptr; 1694 res=0; 1695 int i=0; 1696 for (;i<int(vp.size());++i){ 1697 gen resi; 1698 if (!in_sumab(vp[i],x,a_orig,b_orig,resi,testi,false,contextptr)) 1699 break; 1700 res += resi; 1701 if (is_undef(res)){ 1702 res=0; 1703 return in_sumab(g,x,a_orig,b_orig,res,testi,false,contextptr); 1704 } 1705 } 1706 if (i==int(vp.size())) 1707 return true; 1708 } 1709 } 1710 vecteur v=lvarx(g,x); 1711 vecteur vD=lop(v,at_Kronecker); 1712 if (!vD.empty()){ 1713 gen vD0=vD.front(),A,B,a,b; 1714 if (is_linear_wrt(g,vD0,A,B,contextptr) && is_linear_wrt(vD0._SYMBptr->feuille,x,a,b,contextptr) && in_sumab(B,x,a_orig,b_orig,res,testi,false,contextptr)){ 1715 // sum(A*Kronecker(a*x+b),x,a_orig,b_orig)+res 1716 gen xval=-b/a; 1717 if (is_integer(xval) && is_greater(xval,a_orig,contextptr) && is_greater(b_orig,xval,contextptr)) 1718 res += subst(A,x,xval,false,contextptr); 1719 return true; 1720 } 1721 } 1722 vD=lop(v,at_Heaviside); 1723 if (!vD.empty()){ 1724 gen vD0=vD.front(),A,B,a,b; 1725 if (is_linear_wrt(g,vD0,A,B,contextptr) && is_linear_wrt(vD0._SYMBptr->feuille,x,a,b,contextptr) && in_sumab(B,x,a_orig,b_orig,res,testi,false,contextptr)){ 1726 // sum(A*Heaviside(a*x+b),x,a_orig,b_orig)+res 1727 gen xval=-b/a,resadd; 1728 if (is_integer(xval) && is_strictly_positive(a,contextptr) && in_sumab(A,x,max(a_orig,xval,contextptr),b_orig,resadd,testi,false,contextptr)){ 1729 res += resadd; 1730 return true; 1731 } 1732 } 1733 } 1734 v=loptab(v,sincostan_tab); 1735 bool est_reel=testi?!has_i(g):true; 1736 if (!v.empty()){ 1737 gen w=trig2exp(v,contextptr); 1738 vecteur vexp; 1739 lin(subst(g,v,*w._VECTptr,true,contextptr),vexp,contextptr); 1740 const_iterateur it=vexp.begin(),itend=vexp.end(); 1741 for (;it!=itend;){ 1742 gen coeff=*it,tmp; 1743 ++it; // it -> on the arg of the exp that must be linear 1744 gen axb=coeff*symbolic(at_exp,*it); 1745 ++it; 1746 if (!sumab(axb,*x._IDNTptr,a_orig,b_orig,tmp,!est_reel,contextptr)){ 1747 return false; 1748 } 1749 res += tmp; 1750 } 1751 return true; 1752 } 1753 v.clear(); 1754 polynome p,q,r; 1755 gen a(a_orig),b(b_orig); 1756 if (!est_reel || complex_mode(contextptr)){ 1757 bool b=complex_mode(contextptr); 1758 complex_mode(contextptr)=false; 1759 gen reg(g),img(0),reres,imres; 1760 if (!est_reel){ 1761 reg=re(g,contextptr),img=im(g,contextptr); 1762 } 1763 if (!in_sumab(reg,x,a_orig,b_orig,reres,false,false,contextptr) || !in_sumab(img,x,a_orig,b_orig,imres,false,false,contextptr)){ 1764 complex_mode(contextptr)=b; 1765 return false; 1766 } 1767 complex_mode(contextptr)=b; 1768 res=reres+cst_i*imres; 1769 return true; 1770 } 1771 bool Hyper=is_hypergeometric(g,*x._IDNTptr,v,p,q,r,contextptr); 1772 // g(x+1)/g(x) as p(x+1)/p(x)*q(x)/r(x+1) 1773 if (Hyper){ 1774 // Newton binomial: sum_{x=a}^{b} comb(b-a,x-a)*p^x = (p+1)^(b-a)*p^a 1775 // n=b-a 1776 // comb(n,x+1-a)*p^(x+1)/comb(n,x-a)/p^x 1777 // = p*(x-a)!*(n-x+a)!/(x+1-a)!/(n-x-1+a)!=p*(n-x+a)/(x+1-a) 1778 // q=(-qa)*(n-x+a)=(-qa)*(b-x), r=ra*(x-a) -> -q/qa+r/ra=n 1779 // can be generallized with j-unitroots to 1780 // sum_{x=a}^{b} comb(b-a,j*x-j*a)*p^x 1781 // 1782 gen Q=r2sym(q,v,contextptr),R=r2sym(r,v,contextptr),Qa,Qb,Ra,Rb; 1783 if (a.type==_INT_ && b==plus_inf && p.lexsorted_degree()==0 && r.coord.size()==1 && q+r==0 ){ 1784 // gen coeff=inv(r.coord.front().value,contextptr); 1785 int pui=r.lexsorted_degree(); 1786 // coeff*sum((-1)^k/k^pui,k,a,inf) 1787 // if a is even set res=0, if a is odd set res=-1/a^pui and a++ 1788 res=0; R=inv(R,contextptr); 1789 if (pui==1){ 1790 if (a.val<=0){ 1791 res=R*plus_inf; 1792 return true; 1793 } 1794 res=symbolic(at_ln,2); 1795 if (a.val>1) 1796 res = res-_sum(makesequence(R,x,1,a.val-1),contextptr); 1797 res=res*subst(g/R,x,1,false,contextptr); 1798 return true; 1799 } 1800 // add sum(1/k^pui,k,a,inf) 1801 // -> 2*sum(1/(2*k)^pui,k,a/2,inf)=2^(1-pui)*sum(1/k^pui,k,a/2,inf) 1802 res = - _sum(makesequence(R,x,a,plus_inf),contextptr) + pow(2,1-pui,contextptr)*_sum(makesequence(R,x,(a.val+1)/2,plus_inf),contextptr); 1803 res = -res*subst(g/R,x,a,false,contextptr); 1804 return true; 1805 } 1806 gen n=b-a; 1807 if (is_linear_wrt(Q,x,Qa,Qb,contextptr) && is_linear_wrt(R,x,Ra,Rb,contextptr)){ 1808 // Q/R=(Qa*x+Qb)/(Ra*x+Rb)=(-Qa/Ra)*(-x-Qb/Qa)/(x+Rb/Ra) 1809 gen trueb=normal(-Qb/Qa,contextptr); 1810 gen truea=normal(-Rb/Ra,contextptr); 1811 gen truen=normal(trueb-truea,contextptr); 1812 gen diffa=normal(a-truea,contextptr); 1813 gen diffb=normal(b-trueb,contextptr); 1814 if (diffa.type==_INT_ && diffb.type==_INT_){ 1815 gen P=r2sym(p,v,contextptr); 1816 if (p.lexsorted_degree()==0){ 1817 res = P*(-Qa/Ra)+1; 1818 if (!is_zero(res)) 1819 res=simplify(pow(res,truen,contextptr)*limit(g,*x._IDNTptr,truea,0,contextptr),contextptr); 1820 if (absint(diffb.val)>100 || absint(diffa.val)>100) 1821 return false; 1822 if (diffb.val>0){ // b>trueb: add sum(g,x,trueb+1,b-1) 1823 for (int i=0;i<diffb.val;++i) 1824 res += simplify(limit(g,*x._IDNTptr,trueb+1+i,0,contextptr),contextptr); 1825 } 1826 else { // b<=trueb substract sum(g,x,b+1,trueb) 1827 for (int i=0;i<-diffb.val;++i) 1828 res -= simplify(limit(g,*x._IDNTptr,b+1+i,0,contextptr),contextptr); 1829 } 1830 if (diffa.val>0){ // a>truea : substract sum(g,x,truea,a-1) 1831 for (int i=0;i<diffa.val;++i) 1832 res -= simplify(limit(g,*x._IDNTptr,truea+i,0,contextptr),contextptr); 1833 } 1834 else { // a<=truea: add sum(g,x,a,truea-1) 1835 for (int i=0;i<-diffa.val;++i) 1836 res += simplify(limit(g,*x._IDNTptr,a+i,0,contextptr),contextptr); 1837 } 1838 return true; 1839 } 1840 // recompute the sum by repeted division, first divide g by P 1841 // then divide P by R, by R-1, by R-2, etc. 1842 gen gsurP=ratnormal(g/P,contextptr),prod=1; 1843 while (!is_zero(P)){ 1844 gen tmp=_quorem(gen(makevecteur(P,R,x),_SEQ__VECT),contextptr),tmpres; 1845 if (tmp.type!=_VECT || tmp._VECTptr->size()!=2) 1846 return false; // setsizeerr(); 1847 if (!in_sumab(gsurP*prod,x,a_orig,b_orig,tmpres,false,false,contextptr)) 1848 return false; 1849 res += tmp._VECTptr->back()*tmpres; 1850 prod = prod*R; 1851 R=subst(R,x,x-1,false,contextptr); // R=R-1; 1852 P=tmp._VECTptr->front(); 1853 } 1854 return true; 1855 } 1856 } 1857 } 1858 if (!is_inf(a) && !is_inf(b)) 1859 return false; 1860 if (b==plus_inf && a.type==_INT_ && Hyper){ 1861 polynome s,Q,R; 1862 // limit of q/r at infinity must be 0 1863 gen test=r2e(q,v,contextptr)/r2e(r,v,contextptr); 1864 test=ratnormal(test*test-1,contextptr); 1865 if (is_zero(test)){ 1866 res=_limit(gen(makevecteur(g,x,plus_inf),_SEQ__VECT),contextptr); 1867 if (!is_zero(res)){ 1868 return is_inf(res); 1869 } 1870 } 1871 if (is_strictly_positive(test,contextptr)){ 1872 res=_limit(gen(makevecteur(g,x,plus_inf),_SEQ__VECT),contextptr); 1873 return true; 1874 } 1875 r=taylor(r,1); 1876 // r(x)/q(x)=s(x+1)/s(x)*R(x)/Q(x+1) 1877 AB2PQR(r,q,s,R,Q); 1878 R=taylor(R,-1); 1879 simplify(p,s); 1880 // IMPROVE: make a partial fraction decomposition of p(n)/s(n) 1881 // [could also make ln return ln(1-x) instead of ln(x-1)] 1882 if (sumab_ps(Q,R,v,a,x,g,est_reel,p,s,res,contextptr)) 1883 return true; 1884 } 1885 gen A,B,P; 1886 int type=is_meromorphic(g,x,A,B,P,contextptr); 1887 int even=is_even_odd(g,x,contextptr); 1888 bool complete=(a==minus_inf && b==plus_inf); 1889 if (type==2){ 1890 if (!is_zero(limit(g,*x._IDNTptr,plus_inf,0,contextptr))){ 1891 res=undef; 1892 return true; 1893 } 1894 if ( complete || even==1){ 1895 if (!complete){ 1896 if (a==minus_inf){ 1897 a=-b; 1898 b=plus_inf; 1899 } 1900 if (a.type!=_INT_) 1901 return false; 1902 } 1903 int ai=a.val; 1904 // find the value of int(cos(pi*x)/sin(pi*x)*g,circle(0,R)) 1905 vecteur v=singular(g,x,contextptr); 1906 if (is_undef(v)) 1907 return false; 1908 gen g2=g*cst_pi*symb_cos(cst_pi*x)/symb_sin(cst_pi*x); 1909 // if root is integer it must not be inside a..b 1910 // the sum of residues of v + sum(g,n=-inf..inf, n not in v) will be 0 1911 gen somme_residus; 1912 int vs=int(v.size()); 1913 gen correc=0; 1914 for (int i=0;i<vs;++i){ 1915 gen vi=v[i]; 1916 if (is_integer(vi)){ 1917 if (is_greater(vi,a,contextptr)){ 1918 res=undef; 1919 return true; 1920 } 1921 } 1922 somme_residus += residue(g2,x,vi,contextptr); 1923 if (is_undef(somme_residus)) 1924 return false; 1925 } 1926 if (complete) 1927 res=-somme_residus; 1928 else { // even fraction 1929 if (ai>0){ 1930 for (int i=1;i<ai;++i) 1931 correc -= quotesubstcheck(g,x,i,v,contextptr); 1932 } 1933 if (ai<=0){ 1934 for (int i=0;i>=ai;--i) 1935 correc += quotesubstcheck(g,x,i,v,contextptr); 1936 } 1937 res=(-quotesubstcheck(g,x,0,v,contextptr)-somme_residus)/2; 1938 } 1939 res=correc+res; 1940 res=recursive_normal(trig2exp(res,contextptr),contextptr); 1941 return true; 1942 } 1943 } 1944 return false; 1945 } 1946 1947 // if true put int(g,x=a..b) into res 1948 // a or b must be +/-infinity and a<b sumab(const gen & g,const gen & x,const gen & a_orig,const gen & b_orig,gen & res,bool testi,GIAC_CONTEXT)1949 bool sumab(const gen & g,const gen & x,const gen & a_orig,const gen & b_orig,gen & res,bool testi,GIAC_CONTEXT){ 1950 if (x.type!=_IDNT) 1951 return false; 1952 if (g.is_symb_of_sommet(at_plus)){ 1953 vecteur argv=gen2vecteur(g._SYMBptr->feuille); 1954 int args=int(argv.size()),i; 1955 res=0; 1956 gen tmp; 1957 for (i=0;i<args;++i){ 1958 if (!sumab(argv[i],x,a_orig,b_orig,tmp,testi,contextptr)) 1959 break; 1960 res += tmp; 1961 } 1962 if (i==args) 1963 return true; 1964 } 1965 // detect when 1966 vecteur v=lop(g,at_when); 1967 if (!v.empty()){ 1968 gen A,B,a,b; 1969 identificateur t(" tsumab"); 1970 gen h=quotesubst(g,v.front(),t,contextptr); 1971 if (!is_linear_wrt(h,t,A,B,contextptr)) 1972 return false; 1973 gen heav=v.front()._SYMBptr->feuille; 1974 if (heav.type==_VECT && heav._VECTptr->size()==3){ 1975 B=B+A*heav._VECTptr->back(); 1976 A=A*((*heav._VECTptr)[1]-heav._VECTptr->back()); 1977 heav=heav._VECTptr->front(); 1978 } 1979 else return false; // setsizeerr(); 1980 if (!is_linear_wrt(heav,x,a,b,contextptr) || is_zero(a)) 1981 return false; 1982 if (!sumab(B,x,a_orig,b_orig,res,testi,contextptr)) 1983 return false; 1984 gen c=-b/a; 1985 if (ck_is_greater(c,a_orig,contextptr) && ck_is_greater(b_orig,c,contextptr)) 1986 res += quotesubst(A,x,c,contextptr); 1987 else 1988 *logptr(contextptr) << gettext("Warning, Dirac function outside summation interval") << '\n'; 1989 return true; 1990 } 1991 // detect Heaviside 1992 v=lop(g,at_Heaviside); 1993 if (v.empty()) 1994 return in_sumab(g,x,a_orig,b_orig,res,testi,true /* do partfrac */,contextptr); 1995 gen A,B,a,b; 1996 identificateur t(" tsumab"); 1997 gen h=quotesubst(g,v.front(),t,contextptr); 1998 if (!is_linear_wrt(h,t,A,B,contextptr)) 1999 return false; 2000 // A*Heaviside()+B 2001 gen heav=v.front()._SYMBptr->feuille; 2002 if (!is_linear_wrt(heav,x,a,b,contextptr) || is_zero(a)) 2003 return false; 2004 if (!sumab(B,x,a_orig,b_orig,res,testi,contextptr)) 2005 return false; 2006 // for A additional condition a*x+b>=0 : if a>0 x>=-b/a else x<=-b/a 2007 gen c=-b/a; 2008 gen newa,newb; 2009 if (ck_is_positive(a,contextptr)){ 2010 if (ck_is_greater(a_orig,c,contextptr)){ 2011 newa=a_orig; newb=b_orig; 2012 } 2013 else { 2014 if (ck_is_greater(b_orig,c,contextptr)){ 2015 newa=_ceil(c,contextptr); 2016 newb=b_orig; 2017 } 2018 else 2019 return true; 2020 } 2021 } 2022 else { 2023 if (ck_is_greater(a_orig,c,contextptr)) 2024 return true; 2025 if (ck_is_greater(b_orig,c,contextptr)){ 2026 newa=a_orig; 2027 newb=_floor(c,contextptr); 2028 } 2029 else { 2030 newa=a_orig; 2031 newb=b_orig; 2032 } 2033 } 2034 gen sumA; 2035 if (!sumab(A,x,newa,newb,sumA,testi,contextptr)) 2036 return false; 2037 res += sumA; 2038 return true; 2039 } 2040 2041 #ifndef NO_NAMESPACE_GIAC 2042 } // namespace giac 2043 #endif // ndef NO_NAMESPACE_GIAC 2044 2045