1 /* -*- mode:C++ ; compile-command: "g++-3.4 -DHAVE_CONFIG_H -I. -I.. -DIN_GIAC -g -c derive.cc" -*- */ 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 "usual.h" 24 #include "symbolic.h" 25 #include "unary.h" 26 #include "poly.h" 27 #include "sym2poly.h" // for equalposcomp 28 #include "tex.h" 29 #include "prog.h" 30 #include "intg.h" 31 #include "subst.h" 32 #include "plot.h" 33 #include "modpoly.h" 34 #include "moyal.h" 35 #include "alg_ext.h" 36 #include "giacintl.h" 37 38 #ifndef NO_NAMESPACE_GIAC 39 namespace giac { 40 #endif // ndef NO_NAMESPACE_GIAC 41 eval_before_diff(const gen & expr,const gen & variable,GIAC_CONTEXT)42 gen eval_before_diff(const gen & expr,const gen & variable,GIAC_CONTEXT){ 43 identificateur tmp_x_id(" eval_before_diff_x"); 44 gen tmp_x(tmp_x_id); 45 gen res=subst(expr,variable,tmp_x,false,contextptr); // replace variable by a non affected identifier 46 gen save_vx_var=vx_var; 47 if (variable==vx_var) vx_var=tmp_x; 48 int m=calc_mode(contextptr); 49 calc_mode(0,contextptr); 50 res=eval(res,1,contextptr); // eval res (all identifiers except X will be replaced by their values) 51 res=eval(res,1,contextptr); // eval res (all identifiers except X will be replaced by their values) 52 calc_mode(m,contextptr); 53 vx_var=save_vx_var; 54 res=subst(res,tmp_x,variable,false,contextptr); 55 return res; 56 } 57 depend(const gen & g,const identificateur & i)58 bool depend(const gen & g,const identificateur & i){ 59 if (g.type==_IDNT) 60 return *g._IDNTptr==i; 61 if (g.type==_SYMB) 62 return depend(g._SYMBptr->feuille,i); 63 if (g.type!=_VECT) 64 return false; 65 const_iterateur it=g._VECTptr->begin(),itend=g._VECTptr->end(); 66 for (;it!=itend;++it){ 67 if (depend(*it,i)) 68 return true; 69 } 70 return false; 71 } 72 count_noncst(const gen & g,const identificateur & i)73 static int count_noncst(const gen & g,const identificateur & i){ 74 if (g.type!=_VECT) 75 return depend(g,i)?1:0; 76 int res=0; 77 for (unsigned j=0;j<g._VECTptr->size();++j){ 78 if (depend((*g._VECTptr)[j],i)) 79 ++res; 80 } 81 return res; 82 } 83 derive_SYMB(const gen & g_orig,const identificateur & i,GIAC_CONTEXT)84 static gen derive_SYMB(const gen &g_orig,const identificateur & i,GIAC_CONTEXT){ 85 const symbolic & s = *g_orig._SYMBptr; 86 if (s.sommet==at_pnt){ 87 gen f=g_orig._SYMBptr->feuille; 88 if (f.type==_VECT && !f._VECTptr->empty()){ 89 vecteur v=*f._VECTptr; 90 v[0]=derive(v[0],i,contextptr); 91 f=gen(v,f.subtype); 92 return symbolic(at_pnt,f); 93 } 94 } 95 // if s does not depend on i return 0 96 if (!depend(g_orig,i)) 97 return zero; 98 // rational operators are treated first for efficiency 99 if (s.sommet==at_plus){ 100 bool do_step=step_infolevel(contextptr)>1 && count_noncst(s.feuille,i)>1; 101 if (do_step) 102 gprintf(gettext("Derivative of %gen apply linearity: (u+v+...)'=u'+v'+..."),makevecteur(s),contextptr); 103 if (s.feuille.type!=_VECT) 104 return derive(s.feuille,i,contextptr); 105 vecteur::const_iterator iti=s.feuille._VECTptr->begin(),itend=s.feuille._VECTptr->end(); 106 int taille=int(itend-iti); 107 if (taille==2) 108 return derive(*iti,i,contextptr)+derive(*(iti+1),i,contextptr); 109 vecteur v; 110 v.reserve(taille); 111 gen e; 112 for (;iti!=itend;++iti){ 113 e=derive(*iti,i,contextptr); 114 if (is_undef(e)) 115 return e; 116 if (!is_zero(e)) 117 v.push_back(e); 118 } 119 if (v.size()==1) 120 return v.front(); 121 if (v.empty()) 122 return zero; 123 gen res=_plus(gen(v,_SEQ__VECT),contextptr); // symbolic(at_plus,v); 124 if (do_step) 125 gprintf(gettext("Hence derivative of %gen by linearity is %gen"),makevecteur(g_orig,res),contextptr); 126 return res; 127 } 128 if (s.sommet==at_prod){ 129 bool do_step=step_infolevel(contextptr)>1 && count_noncst(s.feuille,i)>1; 130 if (s.feuille.type==_VECT && s.feuille._VECTptr->size()==2 && s.feuille._VECTptr->back().is_symb_of_sommet(at_inv) && !is_constant_wrt(s.feuille._VECTptr->back()._SYMBptr->feuille,i,contextptr)){ 131 gen u=s.feuille._VECTptr->front(),v=s.feuille._VECTptr->back()._SYMBptr->feuille; 132 if (do_step) 133 gprintf(gettext("Derivative of %gen/%gen, a quotient: (u/v)'=(u'*v-u*v')/v^2"),makevecteur(u,v),contextptr); 134 return (derive(u,i,contextptr)*v-u*derive(v,i,contextptr))/pow(v,2,contextptr); 135 } 136 if (do_step) 137 gprintf(gettext("Derivative of %gen, apply product rule: (u*v*...)'=u'*v*...+u*v'*...+..."),makevecteur(s),contextptr); 138 if (s.feuille.type!=_VECT) 139 return derive(s.feuille,i,contextptr); 140 vecteur::const_iterator itbegin=s.feuille._VECTptr->begin(),itj,iti,itend=s.feuille._VECTptr->end(); 141 int taille=int(itend-itbegin); 142 // does not work because of is_linear_wrt e.g. for cos(3*pi/4) 143 // if (taille==2) return derive(*itbegin,i,contextptr)*(*(itbegin+1))+(*itbegin)*derive(*(itbegin+1),i,contextptr); 144 vecteur v,w; 145 v.reserve(taille); 146 w.reserve(taille); 147 gen e; 148 for (iti=itbegin;iti!=itend;++iti){ 149 w.clear(); 150 e=derive(*iti,i,contextptr); 151 if (is_undef(e)) 152 return e; 153 if (!is_zero(e)){ 154 for (itj=itbegin;itj!=iti;++itj) 155 w.push_back(*itj); 156 w.push_back(e); 157 ++itj; 158 for (;itj!=itend;++itj) 159 w.push_back(*itj); 160 v.push_back(_prod(w,contextptr)); 161 } 162 } 163 if (v.size()==1) 164 return v.front(); 165 if (v.empty()) 166 return zero; 167 gen res=symbolic(at_plus,gen(v,_SEQ__VECT)); 168 if (do_step) 169 gprintf(gettext("Hence derivative of %gen by product rule is %gen"),makevecteur(g_orig,res),contextptr); 170 return res; 171 } 172 if (s.sommet==at_neg) 173 return -derive(s.feuille,i,contextptr); 174 if (s.sommet==at_pow){ 175 if (s.feuille.type!=_VECT || s.feuille._VECTptr->size()!=2) 176 return gensizeerr(contextptr); 177 gen base = s.feuille._VECTptr->front(),exponent=s.feuille._VECTptr->back(); 178 if (step_infolevel(contextptr)>1){ 179 if (is_constant_wrt(exponent,i,contextptr)) 180 gprintf(gettext("Derivative of a power: (%gen)'=(%gen)*(%gen)'*%gen"),makevecteur(symb_pow(base,exponent),exponent,base,symb_pow(base,exponent-1)),contextptr); 181 else 182 gprintf(gettext("Derivative of a power: (%gen)'=%gen*(%gen)'*ln(%gen)+(%gen)*(%gen)'*%gen"),makevecteur(symb_pow(base,exponent),symb_pow(base,exponent),exponent,base,exponent,base,symb_pow(base,exponent-1)),contextptr); 183 } 184 gen dbase=derive(base,i,contextptr),dexponent=derive(exponent,i,contextptr); 185 // diff(base^exponent)=diff(exp(exponent*ln(base))) 186 // =base^exponent*diff(exponent)*ln(base)+base^(exponent-1)*exponent*diff(base) 187 gen expm1=exponent+gen(-1); 188 if (is_zero(dexponent)) 189 return exponent*dbase*pow(base,expm1,contextptr); 190 return dexponent*ln(base,contextptr)*s+exponent*dbase*pow(base,expm1,contextptr); 191 } 192 if (s.sommet==at_inv){ 193 if (step_infolevel(contextptr)>1) 194 gprintf(gettext("Derivative of inv(u)=-u'/u^2 with u=%gen"),makevecteur(s.feuille),contextptr); 195 if (s.feuille.is_symb_of_sommet(at_pow)){ 196 gen & f = s.feuille._SYMBptr->feuille; 197 if (f.type==_VECT && f._VECTptr->size()==2) 198 return derive(symb_pow(f._VECTptr->front(),-f._VECTptr->back()),i,contextptr); 199 } 200 return rdiv(-derive(s.feuille,i,contextptr),pow(s.feuille,2),contextptr); 201 } 202 if (equalposcomp(inequality_tab,s.sommet)) 203 return 0; 204 if (s.sommet==at_fsolve && s.feuille.type==_VECT && s.feuille._VECTptr->size()>=2){ 205 vecteur v=*s.feuille._VECTptr; 206 if (v[1].is_symb_of_sommet(at_equal) && v[1]._SYMBptr->feuille.type==_VECT && !v[1]._SYMBptr->feuille._VECTptr->empty()) 207 v[1]=v[1]._SYMBptr->feuille._VECTptr->front(); 208 gen eq=remove_equal(v[0]),y=v[1],x=i; // fsolve(eq(x,y),y) -> y(x), dy/dx=-(deq/dx)/(deq/dy) 209 gen res=-derive(eq,x,contextptr)/derive(eq,y,contextptr); 210 res=subst(res,y,s,false,contextptr); 211 return res; 212 } 213 if (s.sommet==at_rootof){ 214 gen f=s.feuille; 215 if (f.type==_VECT && f._VECTptr->size()==2 && f._VECTptr->front().type==_VECT && f._VECTptr->back().type==_VECT){ 216 vecteur P=*f._VECTptr->front()._VECTptr; 217 vecteur Q=*f._VECTptr->back()._VECTptr; 218 // d/dx(P(y))=(dP/dx)(y) + (dP/dy) (dy/dx) 219 // =(dP/dx)(y) - (dP/dy) (dQ/dx)/(dQ/dy) 220 // where y dependency is in the list polynomial P/Q 221 gen Px=trim(*derive(P,i,contextptr)._VECTptr,0); 222 Px=symb_rootof(Px,Q,contextptr); 223 gen Qx=trim(*derive(Q,i,contextptr)._VECTptr,0); 224 Qx=symb_rootof(Qx,Q,contextptr); 225 gen Py=derivative(P); 226 Py=symb_rootof(Py,Q,contextptr); 227 gen Qy=derivative(Q); 228 Qy=symb_rootof(Qy,Q,contextptr); 229 gen res=Px-(Py*Qx)/Qy; 230 res=normal(res,contextptr); 231 return res; 232 } 233 return gensizeerr(gettext("Derivative of rootof currently not handled")); 234 } 235 if (step_infolevel(contextptr)>1 && s.feuille.type!=_VECT){ 236 if (s.feuille==i){ 237 int save_step=step_infolevel(contextptr); 238 step_infolevel(contextptr)=0; 239 gen der=derive_SYMB(g_orig,i,contextptr); 240 step_infolevel(contextptr)=save_step; 241 gprintf(gettext("Derivative of elementary function %gen is %gen"),makevecteur(g_orig,der),contextptr); 242 } 243 else 244 gprintf(gettext("Derivative of a composition: (%gen)'=(%gen)'*f'(%gen) where f=%gen"),makevecteur(g_orig,s.feuille,s.feuille,s.sommet),contextptr); 245 } 246 if (s.sommet==at_UTPT){ 247 if (s.feuille.type!=_VECT || s.feuille._VECTptr->size()!=2) 248 return gensizeerr(contextptr); 249 gen & arg=s.feuille._VECTptr->back(); 250 return -derive(arg,i,contextptr)*_student(s.feuille,contextptr); 251 } 252 if (s.sommet==at_UTPC){ 253 if (s.feuille.type!=_VECT || s.feuille._VECTptr->size()!=2) 254 return gensizeerr(contextptr); 255 gen & arg=s.feuille._VECTptr->back(); 256 return -derive(arg,i,contextptr)*_chisquare(s.feuille,contextptr); 257 } 258 if (s.sommet==at_UTPF){ 259 if (s.feuille.type!=_VECT || s.feuille._VECTptr->size()!=3) 260 return gensizeerr(contextptr); 261 gen & arg=s.feuille._VECTptr->back(); 262 return -derive(arg,i,contextptr)*_snedecor(s.feuille,contextptr); 263 } 264 if (s.sommet==at_program){ 265 return gensizeerr(gettext("Expecting an expression, not a function")); 266 } 267 if (s.sommet==at_ln){ 268 if (s.feuille.is_symb_of_sommet(at_abs) ) 269 return rdiv(derive(s.feuille._SYMBptr->feuille,i,contextptr),s.feuille._SYMBptr->feuille,contextptr); 270 if (s.feuille.is_symb_of_sommet(at_inv)) 271 return -derive(symbolic(at_ln,s.feuille._SYMBptr->feuille),i,contextptr); 272 if (s.feuille.is_symb_of_sommet(at_prod)){ 273 gen res; 274 const gen &f=s.feuille._SYMBptr->feuille; 275 if (f.type==_VECT){ 276 const_iterateur it=f._VECTptr->begin(),itend=f._VECTptr->end(); 277 for (;it!=itend;++it) 278 res=res+derive(symbolic(at_ln,*it),i,contextptr); 279 return res; 280 } 281 } 282 } 283 if (s.feuille.type==_VECT){ 284 vecteur v=*s.feuille._VECTptr; 285 int vs=int(v.size()); 286 if (vs>=3 && (s.sommet==at_ifte || s.sommet==at_when) ){ 287 for (int j=1;j<vs;++j){ 288 gen & tmp=v[j]; 289 tmp=derive(tmp,i,contextptr); // v[j]=derive(v[j],i,contextptr); 290 // if (is_undef(tmp)) return tmp; 291 // commented otherwise diff(when(x<0,x^2+3,undef)) returns undef 292 } 293 return symbolic(s.sommet,gen(v,s.feuille.subtype)); 294 } 295 if (s.sommet==at_piecewise){ 296 for (int j=0;j<vs/2;++j){ 297 gen & tmp=v[2*j+1]; 298 tmp=derive(tmp,i,contextptr); // v[2*j+1]=derive(v[2*j+1],i,contextptr); 299 } 300 if (vs%2){ 301 gen & tmp=v[vs-1]; 302 tmp=derive(tmp,i,contextptr); // v[vs-1]=derive(v[vs-1],i,contextptr); 303 } 304 return symbolic(s.sommet,gen(v,s.feuille.subtype)); 305 } 306 if (vs==2 && s.sommet==at_NTHROOT){ 307 gen base = v[1],exponent=inv(v[0],contextptr); 308 gen dbase=derive(base,i,contextptr),dexponent=derive(exponent,i,contextptr); 309 // diff(base^exponent)=diff(exp(exponent*ln(base))) 310 // =base^exponent*diff(exponent)*ln(base)+base^(exponent-1)*exponent*diff(base) 311 if (is_zero(dexponent)) 312 return exponent*dbase*s/v[1]; 313 return dexponent*ln(base,contextptr)*s+exponent*dbase*s/v[1]; 314 } 315 if (vs>=3 && s.sommet==at_Beta){ 316 gen v0=v[0],v1=v[1],v2=v[2]; 317 if (!is_zero(derive(v0,i,contextptr)) || !is_zero(derive(v1,i,contextptr)) ) 318 return gensizeerr("diff of incomplete beta with respect to non constant 1st or 2nd arg not implemented"); 319 // diff/v2 of int_0^v2 t^(v0-1)*(1-t)^(v1-1) dt 320 gen tmp=pow(v2,v0-1,contextptr)*pow(1-v2,v1-1,contextptr)*derive(v2,i,contextptr); 321 if (vs==4){ 322 gen v3=v[3]; 323 if (is_one(v3)) 324 return tmp/Beta(v0,v1,contextptr); 325 return gensizeerr(contextptr); 326 gen v3p=derive(v3,i,contextptr); 327 if (!is_zero(v3p)) 328 return tmp-pow(v3,v0-1,contextptr)*pow(1-v3,v1-1,contextptr)*v3p; 329 } 330 return tmp; 331 } 332 if (vs==4 && s.sommet==at_sum){ 333 gen v0=v[0],v1=v[1],v2=v[2],v3=v[3]; 334 if (!is_zero(derive(v1,i,contextptr)) || !is_zero(derive(v2,i,contextptr)) || ! is_zero(derive(v3,i,contextptr)) ) 335 return gensizeerr(gettext("diff of sum with boundaries or mute variable depending on differentiation variable")); 336 if (is_inf(v2) || is_inf(v3)) 337 *logptr(contextptr) << gettext("Warning, assuming derivative commutes with infinite sum") << '\n'; 338 return _sum(makesequence(derive(v0,i,contextptr),v1,v2,v3),contextptr); 339 } 340 if ( (vs==2 || (vs==3 && is_zero(v[2]))) && (s.sommet==at_upper_incomplete_gamma || s.sommet==at_lower_incomplete_gamma || s.sommet==at_Gamma)){ 341 gen v0=v[0],v1=v[1]; 342 if (!is_zero(derive(v0,i,contextptr))) 343 return gensizeerr(gettext("diff of incomplete gamma with respect to non constant 1st arg not implemented")); 344 // diff(int_v1^inf exp(-t)*t^(v0-1) dt) 345 gen tmp1=exp(-v1,contextptr)*pow(v1,v0-1,contextptr)*derive(v1,i,contextptr); 346 return (s.sommet==at_lower_incomplete_gamma)?tmp1:-tmp1; 347 } 348 if (vs==3 && (s.sommet==at_upper_incomplete_gamma || s.sommet==at_lower_incomplete_gamma || s.sommet==at_Gamma)){ 349 return derive(symbolic(s.sommet,makesequence(v[0],v[1]))/symbolic(at_Gamma,v[0]),i,contextptr); 350 } 351 } 352 // now look at other operators, first onearg operator 353 if (s.sommet.ptr()->D){ 354 if (s.feuille.type!=_VECT) 355 return derive(s.feuille,i,contextptr)*(*s.sommet.ptr()->D)(1)(s.feuille,contextptr); 356 // multiargs operators 357 int taille=int(s.feuille._VECTptr->size()); 358 vecteur v; 359 v.reserve(taille); 360 vecteur::const_iterator iti=s.feuille._VECTptr->begin(),itend=s.feuille._VECTptr->end(); 361 gen e; 362 for (int j=1;iti!=itend;++iti,++j){ 363 e=derive(*iti,i,contextptr); 364 if (is_undef(e)) 365 return e; 366 if (!is_zero(e)) 367 v.push_back(e*(*s.sommet.ptr()->D)(j)(s.feuille,contextptr)); 368 } 369 if (v.size()==1) 370 return v.front(); 371 if (v.empty()) 372 return zero; 373 return symbolic(at_plus,gen(v,_SEQ__VECT)); 374 } 375 // integrate 376 if (s.sommet==at_integrate || s.sommet==at_HPINT){ 377 if (s.feuille.type!=_VECT) 378 return s.feuille; 379 vecteur v=*s.feuille._VECTptr; 380 int nargs=int(v.size()); 381 if (nargs<=1) 382 return s.feuille; 383 if (nargs==2 && is_equal(v[1])){ 384 gen v1f=v[1]._SYMBptr->feuille; 385 if (v1f.type==_VECT && v1f._VECTptr->size()==2){ 386 gen v1f1=v1f._VECTptr->front(); 387 gen v1f2=v1f._VECTptr->back(); 388 v[1]=v1f1; 389 if (v1f2.is_symb_of_sommet(at_interval)){ 390 v.push_back(v1f2._SYMBptr->feuille._VECTptr->front()); 391 v.push_back(v1f2._SYMBptr->feuille._VECTptr->back()); 392 nargs=4; 393 } 394 } 395 } 396 gen res,newint; 397 if (v[1]==i) 398 res=v[0]; 399 else { 400 res=subst(v[0],v[1],i,false,contextptr); 401 newint=derive(v[0],i,contextptr); 402 if (nargs<4) 403 newint=integrate_gen(newint,v[1],contextptr); 404 } 405 if (nargs==2) 406 return res+newint; 407 if (nargs==3) 408 return derive(v[2],i,contextptr)*subst(res,i,v[2],false,contextptr); 409 if (nargs==4){ 410 gen a3=derive(v[3],i,contextptr); 411 gen b3=is_zero(a3)?zero:limit(res,i,v[3],-1,contextptr); 412 gen a2=derive(v[2],i,contextptr); 413 gen b2=is_zero(a2)?zero:limit(res,i,v[2],1,contextptr); 414 return a3*b3-a2*b2+_integrate(gen(makevecteur(newint,v[1],v[2],v[3]),_SEQ__VECT),contextptr); 415 } 416 return gensizeerr(contextptr); 417 } 418 if (s.sommet==at_of && s.feuille.type==_VECT && s.feuille._VECTptr->size()==2){ 419 // assuming we do not have an index in a list or matrix! 420 gen f=s.feuille._VECTptr->front(); 421 gen arg=s.feuille._VECTptr->back(); 422 gen darg=derive(arg,i,contextptr); 423 if (!is_one(darg)){ 424 if (darg.type==_VECT){ 425 gen res=0; 426 for (int i=0;i<int(darg._VECTptr->size());++i){ 427 gen fprime=symbolic(at_derive,makesequence(f,i)); 428 res += darg[i]*symbolic(at_of,makesequence(fprime,arg)); 429 } 430 return res; 431 } 432 // f(arg)'=arg'*f'(arg) 433 gen fprime=symbolic(at_derive,f); 434 return darg*symbolic(at_of,makesequence(fprime,arg)); 435 } 436 } 437 // multi derivative and multi-indice derivatives 438 if (s.sommet==at_derive){ 439 if (s.feuille.type!=_VECT) 440 return symbolic(at_derive,gen(makevecteur(s.feuille,vx_var,2),_SEQ__VECT)); 441 if (s.feuille._VECTptr->size()==2){ // derive(f,x) 442 gen othervar=(*s.feuille._VECTptr)[1]; 443 if (othervar.type!=_IDNT) return gensizeerr(gettext("derive.cc/derive_SYMB")); 444 if (*othervar._IDNTptr==i){ // _FUNCnd derivative 445 vecteur res(*s.feuille._VECTptr); 446 symbolic sprime(s); 447 res.push_back(2); 448 return symbolic(at_derive,gen(res,_SEQ__VECT)); 449 } 450 else { 451 vecteur var; 452 var.push_back(othervar); 453 var.push_back(i); 454 vecteur nderiv; 455 nderiv.push_back(1); 456 nderiv.push_back(1); 457 return symbolic(at_derive,gen(makevecteur((*s.feuille._VECTptr)[0],var,nderiv),_SEQ__VECT)); 458 } 459 } 460 else { // derive(f,x,n) 461 if (s.feuille._VECTptr->size()!=3) return gensizeerr(gettext("derive.cc/derive_SYMB")); 462 gen othervar=(*s.feuille._VECTptr)[1]; 463 if (othervar.type==_IDNT){ 464 if (*othervar._IDNTptr==i){ // n+1 derivative 465 vecteur vprime=(*s.feuille._VECTptr); 466 vprime[2] += 1; 467 return symbolic(s.sommet,gen(vprime,_SEQ__VECT)); 468 } 469 else { 470 vecteur var; 471 var.push_back(othervar); 472 var.push_back(i); 473 vecteur nderiv; 474 nderiv.push_back((*s.feuille._VECTptr)[2]); 475 nderiv.push_back(1); 476 return symbolic(at_derive,gen(makevecteur((*s.feuille._VECTptr)[0],var,nderiv),_SEQ__VECT)); 477 } 478 } // end if othervar.type==_IDNT 479 else { // othervar.type must be _VECT 480 if (othervar.type!=_VECT) return gensizeerr(gettext("derive.cc/derive_SYMB")); 481 gen nder((*s.feuille._VECTptr)[2]); 482 if (nder.type!=_VECT || 483 nder._VECTptr->size()!=othervar._VECTptr->size()) return gensizeerr(gettext("derive.cc/derive_SYMB")); 484 vecteur nderiv(*nder._VECTptr); 485 int pos=equalposcomp(*othervar._VECTptr,i); 486 if (pos){ 487 nderiv[pos-1]=nderiv[pos-1]+1; 488 } 489 else { 490 othervar._VECTptr->push_back(i); 491 nderiv.push_back(1); 492 } 493 return symbolic(at_derive,gen(makevecteur((*s.feuille._VECTptr)[0],othervar,nderiv),_SEQ__VECT)); 494 } 495 } 496 } 497 if (s.sommet==at_re || s.sommet==at_im || s.sommet==at_conj){ 498 return s.sommet(derive(s.feuille,i,contextptr),contextptr); 499 } 500 // no info about derivative 501 return symbolic(at_derive,gen(makevecteur(s,i),_SEQ__VECT)); 502 //i.dbgprint(); 503 //s.dbgprint(); 504 } 505 derive_VECT(const vecteur & v,const identificateur & i,GIAC_CONTEXT)506 static gen derive_VECT(const vecteur & v,const identificateur & i,GIAC_CONTEXT){ 507 vecteur w; 508 w.reserve(v.size()); 509 vecteur::const_iterator it=v.begin(),itend=v.end(); 510 for (;it!=itend;++it){ 511 gen tmp=derive(*it,i,contextptr); 512 if (is_undef(tmp)) 513 return tmp; 514 w.push_back(tmp); 515 } 516 return w; 517 } 518 derive(const gen & e,const identificateur & i,GIAC_CONTEXT)519 gen derive(const gen & e,const identificateur & i,GIAC_CONTEXT){ 520 if (abs_calc_mode(contextptr)==38 && i.id_name[0]>='A' && i.id_name[0]<='Z'){ 521 identificateur tmp("xdiff"); 522 gen ee=subst(e,i,tmp,true,contextptr); 523 ee=eval(ee,1,contextptr); 524 ee=subst(ee,i,tmp,true,contextptr); 525 ee=derive(ee,tmp,contextptr); 526 ee=subst(ee,tmp,i,true,contextptr); 527 return ee; 528 } 529 switch (e.type){ 530 case _INT_: case _DOUBLE_: case _ZINT: case _CPLX: case _MOD: case _REAL: case _USER: case _FLOAT_: 531 return 0; 532 case _IDNT: 533 if (is_undef(e)) 534 return e; 535 if (*e._IDNTptr==i) 536 return 1; 537 else 538 return 0; 539 case _SYMB: 540 return derive_SYMB(e,i,contextptr); 541 case _VECT: { 542 gen res=derive_VECT(*e._VECTptr,i,contextptr); 543 if (res.type==_VECT) res.subtype=e.subtype; 544 return res; 545 } 546 case _FRAC: 547 return fraction(derive(e._FRACptr->num,i,contextptr)*e._FRACptr->den-(e._FRACptr->num)*derive(e._FRACptr->den,i,contextptr),e._FRACptr->den); 548 case _EXT: 549 if (is_zero(derive(*(e._EXTptr+1),i,contextptr))) 550 return algebraic_EXTension(derive(*e._EXTptr,i,contextptr),*(e._EXTptr+1)); 551 default: 552 return gentypeerr(contextptr); 553 } 554 return 0; 555 } 556 _VECTderive(const gen & e,const vecteur & v,GIAC_CONTEXT)557 static gen _VECTderive(const gen & e,const vecteur & v,GIAC_CONTEXT){ 558 vecteur w; 559 w.reserve(v.size()); 560 vecteur::const_iterator it=v.begin(),itend=v.end(); 561 for (;it!=itend;++it){ 562 gen tmp=derive(e,*it,contextptr); 563 if (is_undef(tmp)) 564 return tmp; 565 w.push_back(tmp); 566 } 567 return w; 568 } 569 derivesymb(const gen & e,const gen & var,GIAC_CONTEXT)570 static gen derivesymb(const gen& e,const gen & var,GIAC_CONTEXT){ 571 identificateur x(" x"); 572 gen xx(x); 573 gen f=subst(e,var,xx,false,contextptr); 574 f=derive(f,x,contextptr); 575 f=subst(f,xx,var,false,contextptr); 576 return f; 577 } derive(const gen & e,const gen & vars,GIAC_CONTEXT)578 gen derive(const gen & e,const gen & vars,GIAC_CONTEXT){ 579 // cout << e << " " << vars << '\n'; 580 if (is_equal(e)) 581 return symb_equal(derive(e._SYMBptr->feuille[0],vars,contextptr), 582 derive(e._SYMBptr->feuille[1],vars,contextptr)); 583 switch (vars.type){ 584 case _INT_: 585 return symbolic(at_derive,makesequence(e,vars)); 586 case _IDNT: 587 return derive(e,*vars._IDNTptr,contextptr); 588 case _VECT: 589 return _VECTderive(e,*vars._VECTptr,contextptr); 590 case _SYMB: 591 return derivesymb(e,vars,contextptr); 592 default: 593 return gensizeerr(contextptr); 594 } 595 return 0; 596 } 597 derive(const gen & e,const gen & vars,const gen & nderiv,GIAC_CONTEXT)598 gen derive(const gen & e,const gen & vars,const gen & nderiv,GIAC_CONTEXT){ 599 if (is_equal(e)) 600 return symb_equal(derive(e._SYMBptr->feuille[0],vars,nderiv,contextptr), 601 derive(e._SYMBptr->feuille[1],vars,nderiv,contextptr)); 602 if (nderiv.type==_INT_){ 603 int n=nderiv.val; 604 gen ecopie(e),eprime(e); 605 int j=1; 606 for (;j<=n;++j){ 607 eprime=ratnormal(derive(ecopie,vars,contextptr),contextptr); 608 if (is_undef(eprime)) 609 return eprime; 610 if ( (eprime.type==_SYMB) && (eprime._SYMBptr->sommet==at_derive)) 611 break; 612 ecopie=eprime; 613 } 614 if (j==n+1) 615 return eprime; 616 return symbolic(at_derive,gen(makevecteur(ecopie,vars,n+1-j),_SEQ__VECT)); 617 } 618 // multi-index derivation 619 if (nderiv.type!=_VECT || 620 vars.type!=_VECT) return gensizeerr(gettext("derive.cc/derive")); 621 int s=int(nderiv._VECTptr->size()); 622 if (s!=signed(vars._VECTptr->size())) return gensizeerr(gettext("derive.cc/derive")); 623 int j=0; 624 gen ecopie(e); 625 for (;j<s;++j){ 626 ecopie=derive(ecopie,(*vars._VECTptr)[j],(*nderiv._VECTptr)[j],contextptr); 627 } 628 return ecopie; 629 } 630 symb_derive(const gen & a,const gen & b)631 symbolic symb_derive(const gen & a,const gen & b){ 632 return symbolic(at_derive,gen(makevecteur(a,b),_SEQ__VECT)); 633 } 634 symb_derive(const gen & a,const gen & b,const gen & c)635 gen symb_derive(const gen & a,const gen & b,const gen &c){ 636 if (is_zero(c)) 637 return a; 638 if (is_one(c)) 639 return symb_derive(a,b); 640 return symbolic(at_derive,gen(makevecteur(a,b,c),_SEQ__VECT)); 641 } 642 _derive(const gen & args,GIAC_CONTEXT)643 gen _derive(const gen & args,GIAC_CONTEXT){ 644 if (args.type==_STRNG && args.subtype==-1) return args; 645 if (is_equal(args)) 646 return apply_to_equal(args,_derive,contextptr); 647 #ifndef NSPIRE 648 if (calc_mode(contextptr)==1 && args.type!=_VECT) 649 return _derive(makesequence(args,ggb_var(args)),contextptr); 650 #endif 651 vecteur v; 652 if (args.type==_VECT && args.subtype==_POLY1__VECT) 653 return gen(derivative(*args._VECTptr),_POLY1__VECT); 654 if (args.type==_VECT) 655 v=plotpreprocess(gen(*args._VECTptr,_SEQ__VECT),contextptr); 656 else { 657 if (args==vx_var){ // special handling for x(t):=t^2; x' 658 gen tmp=eval(args,1,contextptr); 659 if (tmp!=vx_var) 660 return _derive(tmp,contextptr); 661 } 662 v=plotpreprocess(makesequence(args,vx_var),contextptr); 663 } 664 if (v.size()>1 && v[1].is_symb_of_sommet(at_unquote)) 665 v[1]=eval(v[1],1,contextptr); 666 if (is_undef(v)) 667 return v; 668 if (step_infolevel(contextptr) && v.size()==2 && v[0].type==_SYMB) 669 gprintf(step_derive_header,gettext("===== Derive %gen with respect to %gen ====="),makevecteur(v[0],v[1]),contextptr); 670 gen var,res; 671 if (args.type!=_VECT && is_algebraic_program(v[0],var,res)){ 672 if (var.type==_VECT && var.subtype==_SEQ__VECT && var._VECTptr->size()==1) 673 var=var._VECTptr->front(); 674 res=derive(res,var,contextptr); 675 return symbolic(at_program,makesequence(var,0,res)); 676 } 677 int s=int(v.size()); 678 if (s==2){ 679 if (v[1].type==_VECT && v[1].subtype==_SEQ__VECT){ 680 vecteur & w=*v[1]._VECTptr; 681 int ss=int(w.size()); 682 gen res=v[0]; 683 for (int i=0;i<ss;++i) 684 res=ratnormal(derive(res,w[i],contextptr),contextptr); 685 return res; 686 } 687 if (v[0].type==_SPOL1){ 688 sparse_poly1 res=*v[0]._SPOL1ptr; 689 sparse_poly1::iterator it=res.begin(),itend=res.end(); 690 for (;it!=itend;++it){ 691 gen e=it->exponent; 692 it->coeff=it->coeff*e; 693 it->exponent=e-1; 694 } 695 return res; 696 } 697 if (args.type!=_VECT && v[0].type==_VECT && v[0].subtype==_POLY1__VECT) 698 return gen(derivative(*v[0]._VECTptr),_POLY1__VECT); 699 return derive(v[0],v[1],contextptr); 700 } 701 if (s==3 && (v[2].type==_INT_ || (v[2].type==_VECT && v[2].subtype!=_SEQ__VECT)) ) 702 return derive( v[0],v[1],v[2],contextptr); 703 if (s<3) 704 return gensizeerr(contextptr); 705 if (s>=3 && v.back().is_symb_of_sommet(at_equal)){ 706 gen v_=gen(vecteur(v.begin(),v.end()-1),_SEQ__VECT); 707 v_=_derive(v_,contextptr); 708 return _subst(makesequence(v_,v.back()),contextptr); 709 } 710 const_iterateur it=v.begin()+1,itend=v.end(); 711 res=v[0]; 712 for (;it!=itend;++it) 713 res=ratnormal(_derive(gen(makevecteur(res,*it),_SEQ__VECT),contextptr),contextptr); 714 return res; 715 } 716 // "unary" version step_derive(const gen & args,GIAC_CONTEXT)717 gen step_derive(const gen & args,GIAC_CONTEXT){ 718 if (step_infolevel(contextptr)) 719 ++step_infolevel(contextptr); 720 gen res; 721 #ifndef NO_STDEXCEPT 722 try { 723 res=_derive(args,contextptr); 724 } catch (std::runtime_error & e){ 725 last_evaled_argptr(contextptr)=NULL; 726 res=string2gen(e.what(),false); 727 res.subtype=-1; 728 } 729 #else 730 res=_derive(args,contextptr); 731 #endif 732 if (step_infolevel(contextptr)) 733 --step_infolevel(contextptr); 734 return res; 735 } _diff(const gen & g,GIAC_CONTEXT)736 gen _diff(const gen & g,GIAC_CONTEXT){ 737 return _derive(g,contextptr); 738 } 739 static const char _derive_s []="diff"; printasderive(const gen & feuille,const char * sommetstr,GIAC_CONTEXT)740 static string printasderive(const gen & feuille,const char * sommetstr,GIAC_CONTEXT){ 741 if (feuille.type!=_VECT){ 742 if (feuille.type>=_POLY && feuille.type!=_IDNT) 743 return "("+feuille.print()+")'"; 744 return feuille.print()+"'"; 745 } 746 return sommetstr+("("+feuille.print(contextptr)+")"); 747 } texprintasderive(const gen & feuille,const char * sommetstr,GIAC_CONTEXT)748 static string texprintasderive(const gen & feuille,const char * sommetstr,GIAC_CONTEXT){ 749 if (feuille.type!=_VECT) 750 return gen2tex(feuille,contextptr)+"'"; 751 return "\\frac{\\partial \\left("+gen2tex(feuille._VECTptr->front(),contextptr)+"\\right)}{\\partial "+gen2tex(feuille._VECTptr->back(),contextptr)+"}"; 752 } 753 static define_unary_function_eval4_quoted (__derive,&step_derive,_derive_s,printasderive,texprintasderive); 754 define_unary_function_ptr5( at_derive ,alias_at_derive,&__derive,_QUOTE_ARGUMENTS,true); 755 _grad(const gen & args,GIAC_CONTEXT)756 gen _grad(const gen & args,GIAC_CONTEXT){ 757 if ( args.type==_STRNG && args.subtype==-1) return args; 758 if (args.type!=_VECT) 759 return gensizeerr(contextptr); 760 if (args._VECTptr->size()==3){ 761 gen opt=args._VECTptr->back(); 762 if (opt.is_symb_of_sommet(at_equal) && (opt._SYMBptr->feuille[0]==at_coordonnees || (opt._SYMBptr->feuille[0].type==_INT_ && opt._SYMBptr->feuille[0].val==_COORDS))){ 763 gen coord=(*args._VECTptr)[1]; 764 gen res=_derive(makesequence(args._VECTptr->front(),coord),contextptr); 765 if (res.type==_VECT){ 766 vecteur resv=*res._VECTptr; 767 if (opt._SYMBptr->feuille[1]==at_sphere && resv.size()==3){ 768 resv[1]=resv[1]/coord[0]; 769 resv[2]=resv[2]/(coord[0]*sin(coord[1],contextptr)); 770 return resv; 771 } 772 if (opt._SYMBptr->feuille[1]==at_cylindre && resv.size()>=2){ 773 resv[1]=resv[1]/coord[0]; 774 return resv; 775 } 776 } 777 } 778 } 779 if (args._VECTptr->size()!=2) 780 return gensizeerr(contextptr); 781 return _derive(args,contextptr); 782 } 783 static const char _grad_s []="grad"; 784 static define_unary_function_eval_quoted (__grad,&_grad,_grad_s); 785 define_unary_function_ptr5( at_grad ,alias_at_grad,&__grad,_QUOTE_ARGUMENTS,true); 786 critical(const gen & g,bool extrema_only,GIAC_CONTEXT)787 gen critical(const gen & g,bool extrema_only,GIAC_CONTEXT){ 788 gen arg,var; 789 if (g.type!=_VECT){ 790 arg=g; 791 var=ggb_var(arg); 792 } 793 else { 794 if (g.subtype!=_SEQ__VECT || g._VECTptr->size()<2) 795 return gensizeerr(contextptr); 796 arg=g._VECTptr->front(); 797 var=(*g._VECTptr)[1]; 798 } 799 int savestep=step_infolevel(contextptr); 800 gprintf(gettext("===== Critical points for %gen ====="),makevecteur(arg),contextptr); 801 step_infolevel(contextptr)=0; 802 gen d=_derive(makesequence(arg,var),contextptr); 803 gen deq=_equal(makesequence(d,0*var),contextptr); 804 // *logptr(contextptr) << "Critical points for "<< arg <<": solving " << deq << " with respect to " << var << '\n'; 805 int c=calc_mode(contextptr); 806 calc_mode(0,contextptr); 807 gen s=_solve(makesequence(deq,var),contextptr); 808 step_infolevel(contextptr)=savestep; 809 gprintf(step_extrema1,gettext("Derivative of %gen with respect to %gen is %gen\nSolving %gen with respect to %gen answer %gen"),makevecteur(arg,var,d,deq,var,s.type==_VECT?change_subtype(s,_SEQ__VECT):s),contextptr); 810 calc_mode(c,contextptr); 811 vecteur ls=lidnt(s); 812 for (int i=0;i<int(ls.size());++i){ 813 if (ls[i]==var || (var.type==_VECT && equalposcomp(*var._VECTptr,ls[i]))) 814 return gensizeerr(gettext("solve error while finding critical points")); 815 } 816 if (s.type==_VECT){ 817 vecteur res; 818 step_infolevel(contextptr)=0; 819 gen d2=_derive(makesequence(d,var),contextptr); 820 step_infolevel(contextptr)=savestep; 821 gprintf(step_extrema2,gettext("Hessian at %gen : %gen"),makevecteur(var,d2),contextptr); 822 // *logptr(contextptr) << "Hessian " << d2 << '\n'; 823 vecteur v=*s._VECTptr; 824 int vs=int(v.size()); 825 for (int i=0;i<vs;++i){ 826 gen g=simplify(subst(d2,var,v[i],false,contextptr),contextptr); 827 gprintf(step_extrema3,gettext("Hessian at %gen : %gen"),makevecteur(v[i],g),contextptr); 828 // *logptr(contextptr) << "Hessian at " << v[i] << ": " << g << '\n'; 829 if (ckmatrix(g)){ 830 g=evalf(g,1,contextptr); 831 if (g.type!=_VECT || !is_numericm(*g._VECTptr,true)){ 832 gprintf(step_extrema4,gettext("%gen critical point (unknown type)"),makevecteur(v[i]),contextptr); 833 // *logptr(contextptr) << v[i] << " critical point (unknown type)" << '\n'; 834 continue; 835 } 836 vecteur w=megvl(*g._VECTptr,contextptr); 837 if (!ckmatrix(w)){ 838 gprintf(step_extrema4,gettext("%gen critical point (unknown type)"),makevecteur(v[i]),contextptr); 839 // *logptr(contextptr) << v[i] << " critical point (unknown type)" << '\n'; 840 continue; 841 } 842 int j=0,ws=int(w.size()); 843 for (;j<ws;++j){ 844 if (is_zero(w[0][0])){ 845 gprintf(step_extrema5,gettext("%gen critical point (0 as eigenvalue) %gen"),makevecteur(v[i],_diag(w,contextptr)),contextptr); 846 // *logptr(contextptr) << v[i] << " critical point (0 as eigenvalue) " << _diag(w,contextptr) << '\n'; 847 break; 848 } 849 if (is_positive(-w[0][0]*w[j][j],contextptr)){ 850 gprintf(step_extrema5,gettext("%gen is a saddle point (2 eigenvalues with opposite sign) %gen"),makevecteur(v[i],_diag(w,contextptr)),contextptr); 851 // *logptr(contextptr) << v[i] << " saddle point (2 eigenvalues with opposite sign) " << _diag(w,contextptr) << '\n'; 852 break; 853 } 854 } 855 if (j==ws){ 856 res.push_back(v[i]); 857 if (is_positive(w[0][0],contextptr)){ 858 gprintf(step_extrema6,gettext("%gen is a local minimum %gen"),makevecteur(v[i],_diag(w,contextptr)),contextptr); 859 // *logptr(contextptr) << v[i] << " local minimum " << _diag(w,contextptr) << '\n'; 860 } 861 else { 862 gprintf(step_extrema6,gettext("%gen is a local maximum %gen"),makevecteur(v[i],_diag(w,contextptr)),contextptr); 863 // *logptr(contextptr) << v[i] << " local maximum " << _diag(w,contextptr) << '\n'; 864 } 865 } 866 continue; 867 } // end multi-dimension 868 int d=2; 869 gen curd=d2; 870 if (is_zero(g)){ 871 for (++d;d< NEWTON_DEFAULT_ITERATION;++d){ 872 curd=ratnormal(_derive(makesequence(curd,var),contextptr),contextptr); 873 g=simplify(subst(curd,var,v[i],false,contextptr),contextptr); 874 if (!is_zero(g)) 875 break; 876 } 877 } 878 if (d%2==0 && is_strictly_positive(g,contextptr)){ 879 gprintf(step_extrema7,gettext("%gen is a local minimum"),makevecteur(v[i]),contextptr); 880 // *logptr(contextptr) << v[i] << " local minimum" << '\n'; 881 res.push_back(v[i]); 882 continue; 883 } 884 if (d%2==0 && is_strictly_positive(-g,contextptr)){ 885 gprintf(step_extrema7,gettext("%gen is a local maximum"),makevecteur(v[i]),contextptr); 886 // *logptr(contextptr) << v[i] << " local maximum" << '\n'; 887 res.push_back(v[i]); 888 continue; 889 } 890 if (d==NEWTON_DEFAULT_ITERATION){ 891 gprintf(step_extrema4,gettext("%gen is a critical point (unknown type)"),makevecteur(v[i]),contextptr); 892 //*logptr(contextptr) << v[i] << " critical point (unknown type)" << '\n'; 893 } 894 else { 895 gprintf(step_extrema8,gettext("%gen is an inflection point"),makevecteur(v[i]),contextptr); 896 // *logptr(contextptr) << v[i] << " inflection point" << '\n'; 897 } 898 } 899 if (extrema_only) 900 return res; 901 else 902 return s; 903 } 904 else 905 return s; 906 } 907 908 #if defined USE_GMP_REPLACEMENTS || defined GIAC_GGB _extrema(const gen & g,GIAC_CONTEXT)909 gen _extrema(const gen & g,GIAC_CONTEXT){ 910 return critical(g,true,contextptr); 911 } 912 static const char _extrema_s []="extrema"; 913 static define_unary_function_eval_quoted (__extrema,&_extrema,_extrema_s); 914 define_unary_function_ptr5( at_extrema ,alias_at_extrema,&__extrema,_QUOTE_ARGUMENTS,true); 915 #endif 916 _critical(const gen & g,GIAC_CONTEXT)917 gen _critical(const gen & g,GIAC_CONTEXT){ 918 return critical(g,false,contextptr); 919 } 920 static const char _critical_s []="critical"; 921 static define_unary_function_eval_quoted (__critical,&_critical,_critical_s); 922 define_unary_function_ptr5( at_critical ,alias_at_critical,&__critical,_QUOTE_ARGUMENTS,true); 923 924 // FIXME: This should not use any temporary identifier 925 // Should define the identity operator and write again all rules here 926 // NB: It requires all D operators for functions to be functions! _function_diff(const gen & g,GIAC_CONTEXT)927 gen _function_diff(const gen & g,GIAC_CONTEXT){ 928 if ( g.type==_STRNG && g.subtype==-1) return g; 929 if (g.is_symb_of_sommet(at_function_diff)){ 930 gen & f = g._SYMBptr->feuille; 931 return symbolic(at_of,makesequence(gen(symbolic(at_composepow,makesequence(at_function_diff,2))),f)); 932 } 933 if (g.is_symb_of_sommet(at_of)){ 934 gen & f = g._SYMBptr->feuille; 935 if (f.type==_VECT && f._VECTptr->size()==2){ 936 gen & f1=f._VECTptr->front(); 937 gen & f2=f._VECTptr->back(); 938 if (f1.is_symb_of_sommet(at_composepow)){ 939 gen & f1f=f1._SYMBptr->feuille; 940 if (f1f.type==_VECT && f1f._VECTptr->size()==2 && f1f._VECTptr->front()==at_function_diff){ 941 return symbolic(at_of,makesequence(gen(symbolic(at_composepow,makesequence(at_function_diff,f1f._VECTptr->back()+1))),f2)); 942 } 943 } 944 } 945 } 946 identificateur _tmpi(" _x"); 947 gen _tmp(_tmpi); 948 gen dg(derive(g(_tmp,contextptr),_tmp,contextptr)); 949 if (lop(dg,at_derive).empty()){ 950 identificateur tmpi(" x"); 951 gen tmp(tmpi); 952 gen res=symb_program(tmp,zero,quotesubst(dg,_tmp,tmp,contextptr),contextptr); 953 return res; 954 } 955 return symbolic(at_function_diff,g); 956 } 957 static const char _function_diff_s []="function_diff"; 958 static define_unary_function_eval (__function_diff,&_function_diff,_function_diff_s); 959 define_unary_function_ptr5( at_function_diff ,alias_at_function_diff,&__function_diff,0,true); 960 961 static const char _fonction_derivee_s []="fonction_derivee"; 962 static define_unary_function_eval (__fonction_derivee,&_function_diff,_fonction_derivee_s); 963 define_unary_function_ptr5( at_fonction_derivee ,alias_at_fonction_derivee,&__fonction_derivee,0,true); 964 _implicit_diff(const gen & args,GIAC_CONTEXT)965 gen _implicit_diff(const gen & args,GIAC_CONTEXT){ 966 if (is_undef(args)) return args; 967 if (args.type!=_VECT || (args._VECTptr->size()!=3 && args._VECTptr->size()!=4)) 968 return gensizeerr(contextptr); 969 int ndiff=1; 970 if (args._VECTptr->size()==4){ 971 gen g=args._VECTptr->back(); 972 if (!is_integral(g) || g.type!=_INT_ || g.val<1) 973 return gensizeerr(contextptr); 974 ndiff=g.val; 975 } 976 gen eq(remove_equal(args._VECTptr->front())),x((*args._VECTptr)[1]),y((*args._VECTptr)[2]); 977 gen yprime=-derive(eq,x,contextptr)/derive(eq,y,contextptr); 978 if (ndiff==1) 979 return yprime; 980 gen yn=yprime; 981 for (int n=2;n<=ndiff;++n){ 982 yn=ratnormal(derive(yn,x,contextptr)+derive(yn,y,contextptr)*yprime,contextptr); 983 } 984 return yn; 985 } 986 static const char _implicit_diff_s []="implicit_diff"; 987 static define_unary_function_eval (__implicit_diff,&_implicit_diff,_implicit_diff_s); 988 define_unary_function_ptr5( at_implicit_diff ,alias_at_implicit_diff,&__implicit_diff,0,true); 989 990 // mode==0 for domain, ==1 for singular values domain(const gen & f,const gen & x,vecteur & eqs,vecteur & excluded,int mode,GIAC_CONTEXT)991 void domain(const gen & f,const gen & x,vecteur & eqs,vecteur &excluded,int mode,GIAC_CONTEXT){ 992 vecteur v=lvarxwithinv(f,x,contextptr); 993 lvar(f,v); 994 for (int i=0;i<int(v.size());++i){ 995 gen g=v[i]; 996 if (g.is_symb_of_sommet(at_nop)) 997 g=g._SYMBptr->feuille; 998 if (is_constant_wrt(g,x,contextptr)) 999 continue; 1000 if (g.type!=_SYMB) 1001 continue; 1002 gen gf=g._SYMBptr->feuille; 1003 domain(gf,x,eqs,excluded,mode,contextptr); 1004 unary_function_ptr & u=g._SYMBptr->sommet; 1005 if (u==at_inv || u==at_Ei || (mode==1 && (u==at_ln || u==at_log10 || u==at_Ci))){ 1006 excluded=mergevecteur(excluded,gen2vecteur(_solve(makesequence(symb_equal(gf,0),x),contextptr))); 1007 continue; 1008 } 1009 if (u==at_pow){ 1010 if (mode==1){ 1011 excluded=mergevecteur(excluded,gen2vecteur(_solve(makesequence(symb_equal(gf[0],0),x),contextptr))); 1012 continue; 1013 } 1014 if (is_constant_wrt(gf[1],x,contextptr) && is_greater(gf[1],0,contextptr)) 1015 eqs.push_back(symb_superieur_egal(gf[0],0)); 1016 else 1017 eqs.push_back(symb_superieur_strict(gf[0],0)); 1018 continue; 1019 } 1020 if (u==at_ln || u==at_log10 || u==at_Ci){ 1021 eqs.push_back(symb_superieur_strict(gf,0)); 1022 continue; 1023 } 1024 if (u==at_acosh){ 1025 if (mode==1) 1026 excluded=mergevecteur(excluded,gen2vecteur(_solve(makesequence(symb_equal(gf,0),x),contextptr))); 1027 else 1028 eqs.push_back(symb_superieur_egal(gf,1)); 1029 continue; 1030 } 1031 if (u==at_asin || u==at_acos || u==at_atanh){ 1032 if (mode==1) 1033 excluded=mergevecteur(excluded,gen2vecteur(_solve(makesequence(symb_equal(pow(gf,2,contextptr),0),x),contextptr))); 1034 else 1035 eqs.push_back(symb_inferieur_egal(pow(gf,2,contextptr),1)); 1036 continue; 1037 } 1038 if (u==at_tan){ 1039 excluded=mergevecteur(excluded,gen2vecteur(_solve(makesequence(symb_equal(symb_cos(gf),0),x),contextptr))); 1040 continue; 1041 } 1042 if (u==at_sin || u==at_cos || u==at_exp || u==at_atan) 1043 continue; 1044 if (u==at_sinh || u==at_cosh || u==at_tanh) 1045 continue; 1046 if (u==at_floor || u==at_ceil || u==at_round || u==at_abs || u==at_sign || u==at_max || u==at_min) 1047 continue; 1048 *logptr(contextptr) << g << gettext(" function not supported, doing like if it was defined") << '\n'; 1049 } 1050 } domain(const gen & f,const gen & x,int mode,GIAC_CONTEXT)1051 gen domain(const gen & f,const gen & x,int mode,GIAC_CONTEXT){ 1052 // domain of expression f with respect to variable x 1053 if (x.type!=_IDNT){ 1054 gen domainx(identificateur("domainx")); 1055 return domain(subst(f,x,domainx,false,contextptr),domainx,mode,contextptr); 1056 } 1057 vecteur eqs,excluded,res; 1058 bool b=complex_mode(contextptr); 1059 complex_mode(false,contextptr); 1060 #ifndef NO_STDEXCEPT 1061 try { 1062 #endif 1063 domain(f,x,eqs,excluded,mode,contextptr); 1064 res=gen2vecteur(_solve(makesequence(eqs,x),contextptr)); 1065 #ifndef NO_STDEXCEPT 1066 } catch (std::runtime_error & e ) { 1067 last_evaled_argptr(contextptr)=NULL; 1068 *logptr(contextptr) << e.what() << '\n'; 1069 } 1070 #endif 1071 complex_mode(b,contextptr); 1072 comprim(excluded); 1073 if (mode==1) 1074 return excluded; 1075 if (excluded.empty()) 1076 return res.size()==1?res.front():res; 1077 vecteur tmp; 1078 for (int i=0;i<int(excluded.size());++i){ 1079 tmp.push_back(symbolic(at_different,makesequence(x,excluded[i]))); 1080 } 1081 if (res.size()==1 && res.front()==x) 1082 return tmp.size()==1?tmp.front():symbolic(at_and,gen(tmp,_SEQ__VECT)); 1083 else { 1084 // check if excluded values are solutions inside res 1085 for (int i=0;i<int(res.size());++i){ 1086 for (int j=0;j<int(excluded.size());++j){ 1087 gen resi=subst(res[i],x,excluded[j],false,contextptr); 1088 resi=eval(resi,1,contextptr); 1089 if (is_zero(resi)) 1090 continue; 1091 if (!res[i].is_symb_of_sommet(at_and)){ 1092 res[i]=symbolic(at_and,makesequence(res[i],tmp[j])); 1093 continue; 1094 } 1095 vecteur v=gen2vecteur(res[i]._SYMBptr->feuille); 1096 v.push_back(tmp[j]); 1097 res[i]=symbolic(at_and,gen(v,_SEQ__VECT)); 1098 } 1099 } 1100 return res; 1101 } 1102 // not reached 1103 if (res.size()==1){ 1104 tmp.insert(tmp.begin(),res.front()); 1105 return symbolic(at_and,gen(tmp,_SEQ__VECT)); 1106 } 1107 tmp.insert(tmp.begin(),symbolic(at_ou,gen(res,_SEQ__VECT))); 1108 return symbolic(at_and,gen(tmp,_SEQ__VECT)); 1109 } _domain(const gen & args,GIAC_CONTEXT)1110 gen _domain(const gen & args,GIAC_CONTEXT){ 1111 if (is_undef(args)) return args; 1112 if (args.type!=_VECT || args.subtype!=_SEQ__VECT) 1113 return domain(args,vx_var,0,contextptr); 1114 vecteur v=*args._VECTptr; 1115 if (v.size()<2) 1116 return gensizeerr(contextptr); 1117 if (is_integral(v[1])) 1118 v.insert(v.begin()+1,vx_var); 1119 if (v.size()==2) 1120 v.push_back(0); 1121 if (v[2].type!=_INT_) 1122 return gensizeerr(contextptr); 1123 return domain(v[0],v[1],v[2].val,contextptr); 1124 } 1125 static const char _domain_s []="domain"; 1126 static define_unary_function_eval (__domain,&_domain,_domain_s); 1127 define_unary_function_ptr5( at_domain ,alias_at_domain,&__domain,0,true); 1128 1129 #ifndef NO_NAMESPACE_GIAC 1130 } // namespace giac 1131 #endif // ndef NO_NAMESPACE_GIAC 1132