1 // -*- mode:C++ ; compile-command: "g++-3.4 -I.. -g -c risch.cc -DHAVE_CONFIG_H -DIN_GIAC" -*- 2 #include "giacPCH.h" 3 /* 4 * Copyright (C) 2003,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 "vector.h" 22 #include <cmath> 23 #include <cstdlib> 24 #include "sym2poly.h" 25 #include "usual.h" 26 #include "intg.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 "misc.h" 40 #include "giacintl.h" 41 42 #ifndef NO_NAMESPACE_GIAC 43 namespace giac { 44 #endif // ndef NO_NAMESPACE_GIAC 45 46 static bool in_risch(const gen & e,const identificateur & x,const vecteur & v,const gen & allowed_lnarg,gen & prim,gen & lncoeff,gen & remains_to_integrate,GIAC_CONTEXT); 47 static bool risch_poly_part(const vecteur & e,int shift,const identificateur & x,const vecteur & v,const gen & allowed_lnarg,gen & prim,gen & lncoeff,gen & remains_to_integrate,GIAC_CONTEXT); 48 static bool risch_desolve(const gen & f,const gen & g,const identificateur & x,const vecteur & v,gen & y,bool f_is_derivative,GIAC_CONTEXT); 49 ln_exp(const gen & g,GIAC_CONTEXT)50 static gen ln_exp(const gen & g,GIAC_CONTEXT){ 51 if (g.is_symb_of_sommet(at_exp)) 52 return g._SYMBptr->feuille; 53 return symbolic(at_ln,g); 54 } 55 56 // returns true & the tower of extension if g is elementary 57 // false otherwise risch_tower(const identificateur & x,gen & g,vecteur & v,GIAC_CONTEXT)58 static bool risch_tower(const identificateur & x,gen &g, vecteur & v,GIAC_CONTEXT){ 59 g=tsimplify(pow2expln(g,x,contextptr),contextptr); 60 // ln(exp(...))-> ... 61 vector<const unary_function_ptr *> vu; 62 vu.push_back(at_ln); 63 vector <gen_op_context> vv; 64 vv.push_back(ln_exp); 65 g=ratnormal(subst(g,vu,vv,true,contextptr)); 66 v=rlvarx(g,x); 67 const_iterateur it=v.begin(),itend=v.end(); 68 for (;it!=itend;++it){ 69 if (*it==x) 70 continue; 71 if (!it->is_symb_of_sommet(at_exp) && (!it->is_symb_of_sommet(at_ln)) ) 72 return false; 73 } 74 reverse(v.begin(),v.end()); // most complex var at the beginning 75 return true; 76 } 77 78 // Compute the derivative of a poly or fraction 79 // The derivative wrt to x_i (i-th index) is the i-th element of v 80 // The last element of v should normally be 1 (derivative of x) diff(const polynome & f,const vecteur & v)81 static fraction diff(const polynome & f,const vecteur & v){ 82 int s=int(v.size()); 83 if (f.dim<s) 84 return fraction(gensizeerr(gettext("Risch diff dimension"))); 85 fraction res(zero); 86 std::vector< monomial<gen> >::const_iterator it=f.coord.begin(),itend=f.coord.end(); 87 for (;it!=itend;++it){ 88 index_t i= it->index.iref(); 89 fraction tmp(zero); 90 for (int n=0;n<s;++n){ 91 if (i[n]){ 92 --i[n]; 93 if (v[n].type==_FRAC) 94 tmp=tmp+gen(polynome(monomial<gen>(i[n]+1,i)))*(*v[n]._FRACptr); 95 else 96 tmp=tmp+fraction(gen(polynome(monomial<gen>(i[n]+1,i)))*v[n]); 97 ++i[n]; 98 } 99 } 100 res=res+it->value*tmp; 101 } 102 return res; 103 } 104 105 rothstein_trager_resultant(const polynome & num,const polynome & den,const vecteur & vl,polynome & p1,GIAC_CONTEXT)106 static polynome rothstein_trager_resultant(const polynome & num,const polynome & den,const vecteur & vl,polynome & p1,GIAC_CONTEXT){ 107 int s=num.dim; 108 fraction denprime=diff(den,vl); 109 if (is_undef(denprime.num)){ 110 return gen2poly(denprime.num,s); 111 } 112 p1=num*gen2poly(denprime.den,s); 113 p1=p1.untrunc1(); 114 polynome p2(monomial<gen>(plus_one,1,s+1)); 115 p2=p2*gen2poly(denprime.num,s).untrunc1(); 116 p1=p1-p2; 117 p2=den.untrunc1(); 118 // exchange var 1 (parameter t) and 2 (top tower variable) 119 vector<int> i=transposition(0,1,s+1); 120 p1.reorder(i); 121 // Change sign of p1 if first coeff is negative 122 if (is_positive(-p1.coord.front().value,context0)) // ok 123 p1=-p1; 124 p2.reorder(i); 125 polynome pres=Tresultant<gen>(p1,p2),pcontent(s); 126 p1.reorder(i); 127 return pres.trunc1(); // pres top var is the parameter t 128 } 129 130 /* 131 fraction diff(const fraction & f,const vecteur & v){ 132 if (f.num.type!=_POLY){ 133 if (f.den.type!=_POLY) 134 return zero; 135 return -diff(*f.den._POLYptr,v)*fraction(f.num,f.den*f.den); 136 } 137 polynome & num = *f.num._POLYptr; 138 if (f.den.type!=_POLY) 139 return diff(*f.num._POLYptr,v)/f.num; 140 polynome & den = *f.den._POLYptr; 141 return diff(num,v)*f.den-diff(den,v)*fraction(f.num,f.den*f.den); 142 } 143 */ 144 145 // diff(n/d*exp(a*x)) diff(const polynome & n,const polynome & d,const gen & a,const vecteur & v,polynome & resn,polynome & resd)146 static bool diff(const polynome & n,const polynome & d,const gen & a,const vecteur & v,polynome & resn,polynome & resd){ 147 int s=n.dim; 148 fraction nprime(diff(n,v)),dprime(diff(d,v)); 149 if (is_undef(nprime.num) || is_undef(dprime.num)) 150 return false; 151 polynome nn(gen2polynome(nprime.num,s)),nd(gen2polynome(nprime.den,s)); 152 polynome dn(gen2polynome(dprime.num,s)),dd(gen2polynome(dprime.den,s)); 153 // ((n/d)*exp(a*x))'=(n'*d-d'*n)/d^2+a*n/d = (nn/nd*d-dn/dd*n)/d^2+a*n/d 154 resn = nn*d*dd-dn*n*nd; 155 resd = d*d*dd*nd; 156 if (!is_zero(a)){ 157 // resn/resd + a*n/d = (resn*aden*d+resd*anum*n)/(d*resd) 158 gen num,den; 159 fxnd(a,num,den); 160 polynome anum=gen2poly(num,resd.dim),aden=gen2poly(den,resd.dim); 161 resn = resn*aden*d + resd*anum*n; 162 resd = resd*d*aden; 163 } 164 simplify(resn,resd); 165 return true; 166 } 167 168 rischde_simplify(polynome & R,polynome & S,polynome & T)169 static bool rischde_simplify(polynome & R,polynome & S, polynome & T){ 170 polynome lcmdeno=gcd(gcd(R,S),T); 171 R=R/lcmdeno; 172 S=S/lcmdeno; 173 T=T/lcmdeno; 174 // if gcd(R,S) does not divide T then there is no solution 175 lcmdeno=gcd(R,S); 176 return lcmdeno.lexsorted_degree()==0; 177 } 178 spde_x(const polynome & S0,const polynome & R0,const polynome & T0,const vecteur & vdiff,polynome & N1,polynome & N2)179 static bool spde_x(const polynome & S0,const polynome & R0,const polynome & T0,const vecteur & vdiff,polynome & N1,polynome & N2){ // Solve S*N+R*N'=T, R constant poly 180 gen con=gcd(Tcontent(S0),gcd(Tcontent(R0),Tcontent(T0))); 181 polynome S(S0/con),R(R0/con),T(T0/con); 182 int s=T.dim; 183 if (T.coord.empty()){ 184 N1=T; 185 N2=gen2poly(plus_one,s); 186 return true; 187 } 188 if (T.lexsorted_degree()<S.lexsorted_degree()) 189 return false; 190 polynome quo(s),rem(s),a(s); 191 T.TPseudoDivRem(S,quo,rem,a); 192 // a*T=quo*S+rem -> quo/a is the quotient T/S, 193 // it has the same high degree polynomial part as N 194 // We set N=quo/a+ M, where M satisfies 195 // S*M+R*M'=T-S*(quo/a)-R*(quo/a)' 196 // or a*dNden*(S*M+R*M')= a*dNden*(T-S*(quo/a)-R*(quo/a)') 197 // Let R*(quo/a)'=dNnum/dNden 198 // newS=a*dNden*S, newR=a*dNden*R, newT=a*dNden*T-S*quo*dNden-a*dNnum 199 fraction quo_a=fraction(quo,a).normal(); 200 // fraction dN(diff(quo,vdiff)/a-(diff(a,vdiff)*quo)/(a*a)); 201 fraction tmp1(diff(quo,vdiff)); 202 tmp1.den=a*tmp1.den; 203 fraction tmp2(diff(a,vdiff)); 204 if (is_undef(tmp1.num) || is_undef(tmp2.num)) 205 return false; 206 tmp2.num=tmp2.num*quo; 207 tmp2.den=tmp2.den*a*a; 208 fraction dN((tmp1-tmp2).normal()); 209 polynome dNnum(R*gen2poly(dN.num,s)),dNden(gen2poly(dN.den,s)); 210 polynome newT(T*dNden*a-quo*S*dNden-a*dNnum); 211 if (!spde_x(a*dNden*S,a*dNden*R,newT,vdiff,N1,N2)) 212 return false; 213 // M=N1/N2 hence N=quo/a+M=quo/a+N1/N2 214 fraction Nres=quo_a+fraction(N1,N2); 215 N1=gen2poly(Nres.num,s); 216 N2=gen2poly(Nres.den,s); 217 return true; 218 } 219 SPDE(const polynome & R0,const polynome & S0,const polynome & T0,const identificateur & x,const vecteur & v,const vecteur & vdiff,const vecteur & lv,int ydeg,gen & prim,GIAC_CONTEXT)220 static bool SPDE(const polynome &R0,const polynome & S0,const polynome & T0,const identificateur & x,const vecteur & v,const vecteur & vdiff,const vecteur & lv,int ydeg, gen & prim,GIAC_CONTEXT){ 221 // SPDE algorithm to reduce R to a constant polynomial wrt to Z 222 // this will also reduce ydeg, initial equation is Ry'+Sy=T 223 // Principe: if degree(R)>0, find U and V s.t. RU+SV=T and deg(V)<deg(R) 224 // Then y = V + R*z and y'=U - S*z 225 // hence V'+Rz'+R'z=U-S*z -> R z' + (R' - S) z = U - V' 226 if (T0.coord.empty()){ 227 prim=zero; 228 return true; 229 } 230 polynome R(R0),S(S0),T(T0); 231 int s=int(lv.size()); 232 if (ydeg<0 || !rischde_simplify(R,S,T)) 233 return false; 234 int r=R.lexsorted_degree(); 235 if (R.lexsorted_degree()){ 236 polynome U(s),V(s),c(s); 237 // U,V / R*U+S*V=c*T, c cst wrt main var, degV<degR 238 Tabcuv<gen>(S,R,T,V,U,c); 239 fraction dR(diff(R,vdiff)),dV(diff(V,vdiff)); 240 if (is_undef(dR.num)||is_undef(dV.num)) 241 return false; 242 dV.den=dV.den*c; 243 fraction tmpfrac=fraction(V,c*c)*diff(c,vdiff); 244 if (is_undef(tmpfrac.num)) 245 return false; 246 dV=dV-tmpfrac; // now dV=(V/c)' 247 polynome dRnum(gen2poly(dR.num,s)),dRden(gen2poly(dR.den,s)),dVnum(gen2poly(dV.num,s)),dVden(gen2poly(dV.den,s)); 248 polynome newR(R*dRden*dVden),newS((S*dRden+dRnum)*dVden),newT((((U*dVden)/c)-dVnum)*dRden); 249 ydeg=ydeg-r; 250 if(!SPDE(newR,newS,newT,x,v,vdiff,lv,ydeg,prim,contextptr)) 251 return false; 252 prim=r2sym(R,lv,contextptr)*prim+r2sym(V,lv,contextptr)/r2sym(c,lv,contextptr); 253 return true; 254 } 255 // degree(R)=0, Final resolution 256 gen Z=v.front(); 257 if (Z==x || S.lexsorted_degree()){ // S*N+N'=T 258 polynome N1,N2; 259 if (!spde_x(S,R,T,vdiff,N1,N2)) 260 return false; 261 prim=r2sym(N1,lv,contextptr)/r2sym(N2,lv,contextptr); 262 return true; 263 } 264 vecteur v1(v.begin()+1,v.end()); 265 vecteur lv1(lv.begin()+1,lv.end()); 266 vecteur t=polynome2poly1(T,1); 267 t=*r2sym(t,lv1,contextptr)._VECTptr; 268 gen rr=r2sym(R,lv,contextptr); 269 t=divvecteur(t,rr); 270 if (S.coord.empty()){ // solve y'=T -> polynomial part with 1 less var 271 gen lncoeff,remains; 272 return risch_poly_part(t,0,x,v1,plus_one,prim,lncoeff,remains,contextptr); 273 } 274 // for each power solve a risch de with 1 less var 275 if (ydeg<signed(t.size())-1 || Z.type!=_SYMB) 276 return false; 277 // t=mergevecteur(vecteur(ydeg-t.size()-1),t); 278 gen z=Z._SYMBptr->feuille; 279 gen dz=derive(z,x,contextptr); 280 if (is_undef(dz)) 281 return false; 282 gen b=r2sym(S,lv,contextptr)/rr; 283 int tdeg=int(t.size())-1; 284 gen previous,sol; 285 bool ok; 286 for (int k=tdeg;k>=0;--k){ 287 if (Z.is_symb_of_sommet(at_exp)) 288 ok=risch_desolve(k*dz+b,t[tdeg-k],x,v1,sol,false,contextptr); 289 else 290 ok=risch_desolve(b,t[tdeg-k]-(k+1)*previous*dz/z,x,v1,sol,false,contextptr); 291 if (!ok) 292 return false; 293 prim=prim+sol*pow(Z,k); 294 previous=sol; 295 } 296 return true; 297 } 298 299 // solve y'+f*y=g in rational fraction of v 300 // returns true and y or false risch_desolve(const gen & f,const gen & g,const identificateur & x,const vecteur & v,gen & y,bool f_is_derivative,GIAC_CONTEXT)301 static bool risch_desolve(const gen & f,const gen & g,const identificateur & x,const vecteur & v,gen & y,bool f_is_derivative,GIAC_CONTEXT){ 302 gen Z=v.front(); 303 int s=int(v.size()); 304 vecteur lv(lvar(v)); 305 vecteur v1(v.begin()+1,v.end()); 306 vecteur vdiff(s); 307 for (int i=0;i<s;++i){ 308 vdiff[i]=derive(v[i],x,contextptr); 309 if (is_undef(vdiff[i])) 310 return false; 311 } 312 lvar(vdiff,lv); 313 lvar(f,lv); 314 lvar(g,lv); 315 vecteur lv1(lv.begin()+1,lv.end()); 316 int ss=int(lv.size()); 317 for (int i=0;i<s;++i) 318 vdiff[i]=sym2r(vdiff[i],lv,contextptr); 319 // identificateur Zi(" Z"); 320 // gen Ze(Zi); 321 // gen fZ=quotesubst(f,Z,Ze); 322 fraction ff(sym2r(f,lv,contextptr)); 323 polynome fden(gen2poly(ff.den,ss)),fnum(gen2poly(ff.num,ss)),D(gen2poly(plus_one,ss)); 324 polynome fdenred(fden); 325 // compute denominator of y 326 fraction gg(sym2r(g,lv,contextptr)); 327 polynome gden(gen2poly(gg.den,ss)),gnum(gen2poly(gg.num,ss)); 328 polynome gdenred(gden); 329 if (Z.is_symb_of_sommet(at_exp)){ 330 // if Z is an exp, eliminate it in fden/gden and compute Z in D 331 int fdenv=fden.valuation(0),fnumv=fnum.valuation(0),gdenv=gden.valuation(0),gnumv=gnum.valuation(0); 332 index_t ii(ss); 333 ii[0]=-fdenv; 334 fdenred=fdenred.shift(ii); 335 ii[0]=-gdenv; 336 gdenred=gdenred.shift(ii); 337 int alpha=0,beta=fnumv-fdenv,gamma=gnumv-gdenv; 338 if (gamma<0){ 339 alpha=gamma; 340 if (beta<=0){ 341 if (beta!=0) 342 alpha=giacmin(0,gamma-beta); 343 else { // possible cancellation case depend of cst coeff of f 344 vecteur vtmp(polynome2poly1(fnum,1)); 345 gen f0=r2sym(vtmp[0],lv1,contextptr); 346 vtmp=polynome2poly1(fdenred,1); 347 f0=f0/r2sym(vtmp[0],lv1,contextptr); 348 gen lnc,prim,remains; 349 if (in_risch(f0,x,v1,Z._SYMBptr->feuille,prim,lnc,remains,contextptr)&&lnc.type==_INT_) 350 alpha=giacmin(lnc.val,alpha); 351 } 352 } 353 } 354 if (gamma>=0 && beta==0){ 355 // possible cancellation case depend of cst coeff of f 356 vecteur vtmp(polynome2poly1(fnum,1)); 357 gen f0=r2sym(vtmp[0],lv1,contextptr); 358 if (f0.type==_INT_ && f0.val>0) 359 alpha=-f0.val; 360 } 361 D=polynome(monomial<gen>(plus_one,-alpha,1,ss)); // Z^(-alpha) 362 } 363 if (!f_is_derivative){ 364 // Fixme: eliminate residues in fden -> new fdenred 365 // Find degree 1 factors of fdenred 366 polynome tmpy(fdenred.derivative()),tmpw(fdenred); 367 polynome tmpc(simplify(tmpy,tmpw)); 368 tmpy=tmpy-tmpw.derivative(); 369 polynome f1=simplify(tmpw,tmpy); 370 if (f1.lexsorted_degree()){ /// FIXME 371 polynome f1cofact(fdenred/f1); 372 polynome N1(ss),N2(ss),c(ss); 373 Tabcuv<gen>(f1,f1cofact,fnum,N1,N2,c); // fnum/fdenred=N2/f1+... 374 // find resultant_Z(N-t f1' , f1) 375 polynome p1(ss); 376 polynome pres=rothstein_trager_resultant(N2,f1,vdiff,p1,contextptr); 377 // for each negative integer root of pres, multiply D 378 // find linear factors of pres -> FIXME does not work 379 factorization vden; 380 gen extra_div=1; 381 factor(pres,N1,vden,false,false,false,1,extra_div); 382 factorization::const_iterator f_it=vden.begin(),f_itend=vden.end(); 383 // bool ok=true; 384 for (;f_it!=f_itend;++f_it){ 385 int deg=f_it->fact.lexsorted_degree(); 386 if (deg!=1) 387 continue; 388 // extract the root 389 vecteur vtmp=polynome2poly1(f_it->fact,1); 390 gen root=-r2sym(vtmp.back()/vtmp.front(),lv1,contextptr); 391 if (root.type==_INT_ && root.val<0){ 392 identificateur t(" t"); 393 gen tmp1=r2sym(p1,mergevecteur(vecteur(1,t),lv),contextptr); 394 tmp1=subst(tmp1,t,root,false,contextptr); 395 polynome p1subst(gen2poly(sym2r(tmp1,lv,contextptr),ss)); 396 p1subst=gcd(p1subst,f1); 397 D=D*pow(p1subst,-root.val); 398 } 399 } 400 } 401 } 402 polynome c(gcd(fdenred,gdenred)); 403 D=D*gcd(gdenred,gdenred.derivative())/gcd(c,c.derivative()); 404 // y'+f*y=g -> new equation is Ry'+Sy=T, compute R=D,S=fD-D',T=gD^2 405 fraction dD(diff(D,vdiff)); 406 if (is_undef(dD.num)) 407 return false; 408 polynome dDnum(gen2poly(dD.num,ss)),dDden(gen2poly(dD.den,ss)); 409 // then multiply by the lcm of denominators of S and T 410 // simplify by gcd of R, S, T 411 polynome lcmdeno(gden*fden/gcd(gden,fden)*dDden); 412 polynome R(D*lcmdeno),S(fnum*(lcmdeno/fden)*D-dDnum*(lcmdeno/dDden)),T(D*D*gnum*(lcmdeno/gden)); 413 if (!rischde_simplify(R,S,T)) 414 return false; 415 int rd=R.lexsorted_degree(); 416 int sd=S.lexsorted_degree(); 417 int td=T.lexsorted_degree(); 418 polynome Rr(Tfirstcoeff<gen>(R)),Ss(Tfirstcoeff<gen>(S)); 419 // compute max possible degree of y: it depends on Z type 420 int ydeg=td-sd; 421 gen expshift=plus_one; // multiplicative change of variable 422 if (Z==x){ 423 ydeg=td-giacmax(rd-1,sd); 424 if (rd-1==sd){ // test whether S_s/R_r is an integer 425 gen n=Ss.coord.front().value/Rr.coord.front().value; 426 if (n.type==_INT_ && n.val>ydeg && (Ss-n*Rr).coord.empty()) 427 ydeg=giacmax(ydeg,n.val); 428 } 429 } 430 else { 431 if (Z.type!=_SYMB) 432 return false; 433 gen z=Z._SYMBptr->feuille; 434 ydeg=td-giacmax(rd,sd); 435 if (Z.is_symb_of_sommet(at_exp)){ 436 if (rd==sd){ // test whether int S_s/R_r is elementary with n*z coeff 437 gen lnc,prim,remains,tmp=r2sym(Rr,lv,contextptr)/r2sym(Ss,lv,contextptr); 438 if (in_risch(tmp,x,v1,z,prim,lnc,remains,contextptr)&&lnc.type==_INT_) 439 ydeg=giacmax(lnc.val,ydeg); 440 } 441 } 442 else { 443 if (rd==sd+1){ 444 gen lnc,prim,remains,tmp=r2sym(Rr,lv,contextptr)/r2sym(Ss,lv,contextptr); 445 // test whether int S_s/R_r is elementary with exp also elementary 446 if (in_risch(tmp,x,v1,plus_one,prim,lnc,remains,contextptr)){ 447 prim=tsimplify(exp(prim,contextptr),contextptr); 448 vecteur lv2(lv); 449 lvar(prim,lv2); 450 if (lv2==lv){ // the exp is elementary, change variables 451 expshift=prim; 452 fraction expshiftf(sym2r(prim,lv,contextptr)); 453 polynome expnum(gen2poly(expshiftf.num,ss)),expden(gen2poly(expshiftf.den,ss)); 454 R=R*Rr*expden; 455 S=(S*Rr-R*Ss)*expden; 456 // sd=S.lexsorted_degree(); 457 T=(T*Rr)*expnum; 458 } 459 } 460 } 461 } 462 } 463 bool ok=SPDE(R,S,T,x,v,vdiff,lv,ydeg,y,contextptr); 464 y=y*expshift/r2sym(D,lv,contextptr); 465 return ok; 466 } 467 hermite_reduce(const pf<gen> & p_cst,const gen & a,const vecteur & v_derivatives,const vecteur & lv,gen & prim,GIAC_CONTEXT)468 static pf<gen> hermite_reduce(const pf<gen> & p_cst,const gen & a,const vecteur & v_derivatives,const vecteur & lv,gen & prim,GIAC_CONTEXT){ 469 pf<gen> p(p_cst); 470 if (p.mult<=0){ 471 prim=gensizeerr(gettext("risch.cc/hermite_reduce")); 472 return p; 473 } 474 if (p.mult==1) 475 return p_cst; 476 gen expax=exp(r2sym(a,lv,contextptr)*lv.front(),contextptr); 477 fraction pprime=diff(p.fact,v_derivatives); 478 if (is_undef(pprime.num)){ 479 prim=pprime.num; 480 return p; 481 } 482 int s=int(lv.size()); 483 polynome fprime=gen2poly(pprime.num,s),fprimeden=gen2poly(pprime.den,s); 484 polynome d(s),u(s),v(s),C(s); 485 polynome resnum(s),resden(plus_one,s),numtemp(s),dentemp(s); 486 Tegcdpsr(fprime,p.fact,v,u,d); // f*u+f'.num*v=d 487 polynome usave(u),vsave(v),pdensave(s); 488 // reduce p.den to the cofactor 489 p.den=p.den/pow(p.fact,p.mult); 490 // now we are integrating p.num/(p.den*f^p.mult) 491 while (p.mult>1){ 492 pdensave=p.den; 493 Tegcdtoabcuv(fprime,p.fact,p.num,v,u,d,C); // f*u+f'.num*v=C*p.num 494 v=v*fprimeden; // f*u+f'*v=C*p.num 495 // p.num/(p.den*f^p.mult)=(f*u+f'*v)/(C*p.den*f^p.mult) 496 p.mult--; 497 // int(f'/f^(p.mult+1) * v/(C*p.den)*exp(a*x) ) 498 // = 1/p.mult*[-1/f^(p.mult)*v/(C*p.den)*exp(a*x) 499 // + int(1/f^(p.mult)*(v/C*p.den*exp(a*x))')] 500 // update non integrated term 501 if (!diff(v,C*p.den,a,v_derivatives,numtemp,dentemp)){ 502 prim=gensizeerr(gettext("risch.cc/hermite_reduce")); 503 return p; 504 } 505 dentemp=dentemp*p.mult; 506 Tfracadd<gen>(numtemp,dentemp,u,C*p.den,p.num,p.den); 507 simplify(p.num,p.den); 508 // update integrated term 509 pdensave=-C*pdensave; 510 TsimplifybyTlgcd(pdensave,v); 511 pdensave=pdensave*pow(p.fact,p.mult)*p.mult; 512 Tfracadd<gen>(resnum,resden,v,pdensave,numtemp,dentemp); 513 resnum=numtemp; 514 resden=dentemp; 515 // finished? 516 if (p.mult==1) 517 break; 518 // restore Bezout coeffs 519 u=usave; 520 v=vsave; 521 } 522 prim=prim+r2sym(resnum,lv,contextptr)/r2sym(resden,lv,contextptr)*expax; 523 // restore the factor in p.den 524 p.den=p.den*p.fact; 525 return pf<gen>(p); 526 } 527 528 // int(exp(a*x)*coeff,x) where coeff is a rational fraction wrt x 529 // polynomial part: recursive call of int 530 // remaining part: reduce to square free denom by int by part 531 // partial fraction decomp -> Ei 532 // Ei with imaginary arguments -> Ei(i*t)=i*(-pi/2)+i*Si(x)+Ci(x) integrate_ei(const gen & a,const gen & coeff,const identificateur & x,gen & prim,gen & remains_to_integrate,GIAC_CONTEXT)533 static bool integrate_ei(const gen & a,const gen & coeff,const identificateur & x,gen & prim,gen & remains_to_integrate,GIAC_CONTEXT){ 534 gen expax=exp(a*x,contextptr),ima,rea=re(a,contextptr); 535 bool imaneg=false; 536 if (is_zero(rea)){ 537 ima=im(a,contextptr); 538 imaneg=is_positive(-ima,contextptr); 539 } 540 vecteur l(1,x); 541 lvar(coeff,l); 542 lvar(a,l); 543 int s=int(l.size()); 544 vecteur l1(l.begin()+1,l.end()); 545 gen r=e2r(coeff,l,contextptr),ar=e2r(a,l,contextptr); 546 // cout << "Int " << r << '\n'; 547 gen r_num,r_den; 548 fxnd(r,r_num,r_den); 549 if (r_num.type==_EXT) 550 return false; 551 polynome den(gen2poly(r_den,s)),num(gen2poly(r_num,s)); 552 polynome p_content(lgcd(den)); 553 // Square-free factorization 554 factorization vden(sqff(den/p_content)); 555 vector< pf<gen> > pfdecomp; 556 polynome ipnum(s),ipden(s); 557 gen ipshift; 558 partfrac(num,den,vden,pfdecomp,ipnum,ipden); 559 // int( ipnum/ipden*exp(a*x),x) 560 int save=calc_mode(contextptr); 561 calc_mode(0,contextptr); 562 prim += _integrate(gen(makevecteur(expax*r2sym(ipnum/ipden,l,contextptr),x),_SEQ__VECT),contextptr); 563 calc_mode(save,contextptr); 564 if (is_undef(prim)) return false; 565 // Hermite reduction 566 vector< pf<gen> >::iterator it=pfdecomp.begin(); 567 vector< pf<gen> >::const_iterator itend=pfdecomp.end(); 568 for (;it!=itend;++it){ 569 pf<gen> tmp(*it); 570 if (it->mult>1){ 571 vecteur vl(1,1); 572 tmp=hermite_reduce(*it,ar,vl,l,prim,contextptr); 573 if (is_undef(prim)) 574 return false; 575 } 576 if (is_zero(tmp.num)) 577 continue; 578 // ei part for it->num/it->den 579 vecteur itnum=polynome2poly1(tmp.num,1); 580 vecteur itden=derivative(polynome2poly1(tmp.den,1)); 581 factorization vden; 582 gen extra_div=1; 583 factor(tmp.fact,p_content,vden,false,true,true,1,extra_div); // complex+sqrt ok 584 factorization::const_iterator f_it=vden.begin(),f_itend=vden.end(); 585 gen add_prim; 586 // bool ok=true; 587 for (;f_it!=f_itend;++f_it){ 588 int deg=f_it->fact.lexsorted_degree(); 589 if (!deg) 590 continue; 591 if (deg!=1){ 592 remains_to_integrate=remains_to_integrate+expax*r2sym(it->num,l,contextptr)/r2sym(it->den,l,contextptr); 593 break; 594 } 595 // extract the root 596 vecteur vtmp=polynome2poly1(f_it->fact,1); 597 gen root=-vtmp.back()/vtmp.front(); 598 gen rootn=horner(itnum,root); 599 gen rootd=horner(itden,root); 600 root = r2sym(root,l1,contextptr); 601 if (is_zero(im(root,contextptr)) && is_zero(rea)){ 602 prim += (cos(ima*root,contextptr)+cst_i*sin(ima*root,contextptr))*(r2sym(rootn/rootd,l1,contextptr)*(symbolic(at_Ci,(imaneg?(-ima):ima)*(x-root))+(imaneg?(-cst_i):cst_i)*symbolic(at_Si,(imaneg?(-ima):ima)*(x-root)))); 603 } 604 else 605 prim += r2sym(rootn/rootd,l1,contextptr)*exp(a*root,contextptr)*symbolic(at_Ei,a*(x-root)); 606 } 607 } 608 return true; 609 } 610 611 // e is assumed to be a (generallized) poly wrt the top var of v 612 // shift is 0 for a poly or the max Laurent exponent for a generallized poly risch_poly_part(const vecteur & e,int shift,const identificateur & x,const vecteur & v,const gen & allowed_lnarg,gen & prim,gen & lncoeff,gen & remains_to_integrate,GIAC_CONTEXT)613 static bool risch_poly_part(const vecteur & e,int shift,const identificateur & x,const vecteur & v,const gen & allowed_lnarg,gen & prim,gen & lncoeff,gen & remains_to_integrate,GIAC_CONTEXT){ 614 if (v.size()==1){ // the shift must be 1 for integration of a polynomial 615 vecteur tmp=e; 616 reverse(tmp.begin(),tmp.end()); 617 tmp=integrate(tmp,1); 618 if (is_undef(tmp)) 619 return false; 620 reverse(tmp.begin(),tmp.end()); 621 tmp.push_back(zero); 622 prim=prim+symb_horner(tmp,x); 623 return true; 624 } 625 int s=int(e.size()); // degree is s-1 626 vecteur v1(v.begin()+1,v.end()); 627 gen X=v.front(); 628 if (X.is_symb_of_sommet(at_ln)){ 629 gen dX=ratnormal(derive(X,x,contextptr),contextptr); 630 if (is_undef(dX)) 631 return false; 632 // log extension 633 vecteur eprim(s+1); 634 gen lnc,remains; 635 for (int j=s-1;j>0;--j){ 636 // eprim[s-j] ' = e[s-1-j] - (j+1) eprim[s-j-1] * v.front()' 637 if (!in_risch(e[s-1-j]-(j+1)*eprim[s-1-j]*dX,x,v1,X._SYMBptr->feuille,eprim[s-j],lnc,remains,contextptr)){ 638 remains_to_integrate=remains_to_integrate+symb_horner(e,X); 639 return false; 640 } 641 eprim[s-1-j]=eprim[s-1-j]+rdiv(lnc,j+1,contextptr); 642 } 643 gen prim_add; 644 bool ok=in_risch(e[s-1]-eprim[s-1]*dX,x,v1,zero,prim_add,lnc,remains,contextptr); 645 prim=prim+prim_add+symb_horner(eprim,X); 646 remains_to_integrate=remains_to_integrate+remains; 647 return ok; 648 } 649 // exp extension: we have to solve a Risch diff equation for each power 650 gen prim_add,remains; 651 // exp extension 652 vecteur eprim(s); 653 if (!X.is_symb_of_sommet(at_exp)) 654 return false; 655 gen dY=derive(X._SYMBptr->feuille,x,contextptr); 656 if (is_undef(dY)) 657 return false; 658 for (int j=s-1;j>=0;--j){ 659 if (j+shift==0){ 660 bool ok=in_risch(e[s-1-j],x,v1,allowed_lnarg,prim_add,lncoeff,remains,contextptr); 661 prim=prim+prim_add; 662 remains_to_integrate=remains_to_integrate+remains; 663 if (!ok &&!is_zero(allowed_lnarg)) 664 return ok; 665 } 666 else { 667 if (!risch_desolve((j+shift)*dY,e[s-1-j],x,v1,eprim[s-1-j],true,contextptr)) { 668 gen coeff=e[s-1-j],pui=j+shift; 669 gen a,b,c,d; 670 if (is_zero(allowed_lnarg)&&is_linear_wrt(X._SYMBptr->feuille,x,a,b,contextptr) && lvarx(coeff,x)==vecteur(1,x) && integrate_ei(pui*a,coeff,x,c,d,contextptr)){ 671 // add to prim int(exp(pui*(a*x+b))*coeff) 672 // expressed with Ei 673 prim += exp(pui*b,contextptr)*c; 674 remains_to_integrate += exp(pui*b,contextptr)*d; 675 continue; 676 } 677 remains_to_integrate=remains_to_integrate+coeff*pow(X,pui,contextptr); 678 eprim[s-1-j]=0; 679 if (!is_zero(allowed_lnarg)) 680 return false; 681 } 682 } 683 } 684 prim=prim+symb_horner(eprim,X)*pow(X,shift); 685 return true; 686 } 687 688 // Inner Risch algorithm call, v is the tower extension 689 // allowed_lnarg is zero if all ln are allowed or the arg 690 // of the allowed ln if only 1 ln is allowed 691 // lncoeff will contain this coeff if it's the case 692 // prim is the antiderivative 693 // returns false if only 1 ln allowed and no elementary integral found in_risch(const gen & e,const identificateur & x,const vecteur & v,const gen & allowed_lnarg,gen & prim,gen & lncoeff,gen & remains_to_integrate,GIAC_CONTEXT)694 static bool in_risch(const gen & e,const identificateur & x,const vecteur & v,const gen & allowed_lnarg,gen & prim,gen & lncoeff,gen & remains_to_integrate,GIAC_CONTEXT){ 695 prim=zero; 696 lncoeff=zero; 697 remains_to_integrate= zero; 698 int vs=int(v.size()); 699 vecteur l(v); 700 // add top non-x vars 701 lvar(e,l); 702 gen diffv=derive(v,x,contextptr); 703 if (is_undef(diffv) || diffv.type!=_VECT) 704 return false; 705 vecteur vx=*diffv._VECTptr; 706 lvar(vx,l); 707 vecteur vl; 708 int s=int(l.size()); 709 for (int i=0;i<vs;++i) 710 vl.push_back(e2r(vx[i],l,contextptr)); 711 vecteur l1(l.begin()+1,l.end()); 712 gen r=e2r(e,l,contextptr); 713 // cout << "Int " << r << '\n'; 714 gen r_num,r_den; 715 fxnd(r,r_num,r_den); 716 if (r_num.type==_EXT){ 717 remains_to_integrate= e; 718 return false; 719 } 720 polynome den(gen2poly(r_den,s)),num(gen2poly(r_num,s)); 721 polynome p_content(lgcd(den)); 722 // Square-free factorization 723 // FIXME: if ex[ extension, should always treat pole 0 apart 724 factorization vden; 725 int zeromult=den.coord.back().index.front(); 726 if (zeromult){ 727 index_t sh(den.dim); 728 sh[0]=-zeromult; 729 polynome dens=den.shift(sh); 730 vden=sqff(dens/p_content); 731 vden.push_back(facteur<polynome>(polynome(monomial<gen>(1,1,den.dim)),zeromult)); 732 } 733 else 734 vden=sqff(den/p_content); 735 vector< pf<gen> > pfdecomp; 736 polynome ipnum(s),ipden(s); 737 gen ipshift; 738 partfrac(num,den,vden,pfdecomp,ipnum,ipden); 739 // Hermite/ln reduction and 0 isolation for exp extensions 740 vector< pf<gen> >::iterator it=pfdecomp.begin(); 741 vector< pf<gen> >::const_iterator itend=pfdecomp.end(); 742 for (;it!=itend;++it){ 743 if (v.front().is_symb_of_sommet(at_exp) && it->fact.coord.size()==1){ 744 // generalized polynomial part, fact must be 1 [1,0,0...] 745 index_t i(s); 746 i.front()=-it->mult; 747 vecteur tmp=polynome2poly1(it->num,1); 748 tmp=*r2sym(tmp,l1,contextptr)._VECTptr; 749 tmp=divvecteur(tmp,r2sym(it->den.shift(i),l,contextptr)); 750 if (!risch_poly_part(tmp,-it->mult,x,v,allowed_lnarg,prim,lncoeff,remains_to_integrate,contextptr) && !is_zero(allowed_lnarg)) 751 return false; 752 } 753 else { // Hermite reduce it->num/it->den 754 pf<gen> tmp(*it); 755 if (it->mult>1){ 756 tmp=hermite_reduce(*it,0,vl,l,prim,contextptr); 757 if (is_undef(prim)) 758 return false; 759 } 760 if (is_zero(tmp.num)) 761 continue; 762 if (!is_zero(allowed_lnarg)){ // if only 1 log allowed 763 if (pfdecomp.size()>1) 764 return false; 765 // compute u'/u, must be a multiple of it->num/it->den 766 lncoeff=normal(allowed_lnarg*r2sym(it->num,l,contextptr)/(derive(allowed_lnarg,x,contextptr)*r2sym(it->den,l,contextptr)),contextptr); 767 if (!is_zero(derive(lncoeff,x,contextptr))) 768 return false; 769 continue; 770 } 771 // Logarithmic part: compute resultant of num - t * den 772 polynome p1(s); 773 polynome pres=rothstein_trager_resultant(tmp.num,tmp.den,vl,p1,contextptr); 774 // Factorization, should return 1st order factor independant of 775 // the tower variables 776 factorization vden; 777 gen extra_div=1; 778 factor(pres,p_content,vden,false,withsqrt(contextptr),true,1,extra_div); 779 factorization::const_iterator f_it=vden.begin(),f_itend=vden.end(); 780 gen add_prim; 781 bool ok=true; 782 for (;f_it!=f_itend;++f_it){ 783 int deg=f_it->fact.lexsorted_degree(); 784 if (!deg) 785 continue; 786 if (deg!=1){ 787 ok=false; 788 remains_to_integrate=remains_to_integrate+r2sym(it->num,l,contextptr)/r2sym(it->den,l,contextptr); 789 break; 790 } 791 // extract the root 792 vecteur vtmp=polynome2poly1(f_it->fact,1); 793 gen root=-r2sym(vtmp.back()/vtmp.front(),l1,contextptr); 794 if (!is_zero(derive(root,x,contextptr))){ 795 ok=false; 796 // remains_to_integrate=remains_to_integrate+r2sym(it->num,l,contextptr)/r2sym(it->den,l,contextptr); 797 remains_to_integrate=remains_to_integrate+r2sym(tmp.num,l,contextptr)/r2sym(tmp.den,l,contextptr); 798 break; 799 } 800 identificateur t(" t"); 801 gen tmp1=r2sym(p1,mergevecteur(vecteur(1,t),l),contextptr); 802 tmp1=subst(tmp1,t,root,false,contextptr); 803 gen tmparg=_gcd(makevecteur(recursive_normal(tmp1,contextptr),recursive_normal(r2sym(it->fact,l,contextptr),contextptr)),contextptr); 804 add_prim=add_prim+root*ln(tmparg,contextptr); 805 // If tmparg is of maximal degree, this will change the integral 806 // part of the fraction by root*_quo(tmparg',tmparg,X) 807 ipshift=ipshift-root*_quo(makevecteur(derive(tmparg,x,contextptr),tmparg,recursive_normal(v.front(),contextptr)),contextptr); 808 if (is_undef(ipshift)) 809 return false; 810 } 811 if (ok) 812 prim=prim+add_prim; 813 } // end else case (non 0 pole, ie Hermite reduction) 814 } // end for (all denominateurs) 815 // ipnum/ipden = the polynomial part -> integrate it 816 simplify(ipnum,ipden); 817 vecteur tmp=polynome2poly1(ipnum,1); 818 tmp=*r2sym(tmp,l1,contextptr)._VECTptr; 819 tmp=divvecteur(tmp,r2sym(ipden,l,contextptr)); 820 if (tmp.empty()) 821 tmp=vecteur(1,ipshift); 822 else 823 tmp.back() += ipshift; 824 if (!risch_poly_part(tmp,0,x,v,allowed_lnarg,prim,lncoeff,remains_to_integrate,contextptr) && !is_zero(allowed_lnarg)) 825 return false; 826 return true; 827 } 828 risch_lin(const gen & e_orig,const identificateur & x,gen & remains_to_integrate,GIAC_CONTEXT)829 static gen risch_lin(const gen & e_orig,const identificateur & x,gen & remains_to_integrate,GIAC_CONTEXT){ 830 vecteur v; 831 gen e=e_orig; 832 vecteur vatan=lop(e,at_atan); 833 vecteur vln; 834 for (int i=0;i<int(vatan.size());++i){ 835 gen ga=vatan[i]; 836 gen g=ga._SYMBptr->feuille; 837 vln.push_back(ln(ratnormal(rdiv(cst_i+g,cst_i-g),contextptr),contextptr)); 838 vatan[i]=-2*ga*cst_i; 839 } 840 #if 1 841 // extract cst part of exponentials 842 vecteur v1=lop(e,at_exp),v2=v1; 843 for (unsigned i=0;i<v1.size();++i){ 844 gen tmp=v1[i]._SYMBptr->feuille; 845 vecteur tmpv=rlvarx(tmp,x); 846 vecteur tmpw(tmpv.size()); 847 gen tmpcst=subst(tmp,tmpv,tmpw,false,contextptr); 848 tmpcst=exp(tmpcst,contextptr)*exp(ratnormal(tmp-tmpcst,contextptr),contextptr); 849 if (!is_inf(tmpcst) && !is_undef(tmpcst)) 850 v2[i]=tmpcst; 851 } 852 e=subst(e,v1,v2,false,contextptr); 853 #else 854 // texpand added for integrate(x *(x - (exp(x) - exp(-x)) / 2 / ((exp(1) - exp(-1)) / 2))); 855 gen e2=_texpand(e,contextptr); 856 if (rlvarx(e2,x).size()<rlvarx(e,x).size()) 857 e=e2; 858 #endif 859 if (!risch_tower(x,e,v,contextptr)){ 860 remains_to_integrate=e_orig; 861 return zero; 862 } 863 if (v.empty()) 864 return e*x; 865 gen prim,lncoeff; 866 in_risch(e,x,v,zero,prim,lncoeff,remains_to_integrate,contextptr); 867 vector<const unary_function_ptr *> SiCi(1,at_Si); 868 SiCi.push_back(at_Ci); 869 if (!lop(prim,at_Si).empty()){ 870 prim=recursive_normal(prim,contextptr); 871 if (has_i(prim)){ 872 prim=_exp2trig(prim,contextptr); 873 prim=recursive_normal(prim,contextptr); 874 } 875 } 876 if (!vatan.empty()) 877 prim=ratnormal(subst(recursive_ratnormal(prim,contextptr),vln,vatan,false,contextptr),contextptr); 878 if (is_zero(prim)) 879 remains_to_integrate=e_orig; 880 return prim; 881 } 882 risch(const gen & e_orig,const identificateur & x,gen & remains_to_integrate,GIAC_CONTEXT)883 gen risch(const gen & e_orig,const identificateur & x,gen & remains_to_integrate,GIAC_CONTEXT){ 884 #if 0 // def GIAC_HAS_STO_38 885 remains_to_integrate=e_orig; 886 return 0; 887 #endif 888 vecteur vexp; 889 lin(trig2exp(e_orig,contextptr),vexp,contextptr); 890 const_iterateur it=vexp.begin(),itend=vexp.end(); 891 gen rem,remsum,res; 892 for (;it!=itend;){ 893 gen coeff=*it; 894 ++it; 895 gen expo=*it; 896 ++it; rem=0; 897 res = res+risch_lin(coeff*exp(expo,contextptr),x,rem,contextptr); 898 remsum += rem; 899 } 900 res = res+risch_lin(remsum,x,remains_to_integrate,contextptr); 901 if (is_zero(res)){ // perhaps we should derive res and substract it from e_orig 902 remains_to_integrate=e_orig; 903 } 904 else { 905 if (!has_i(e_orig) && has_i(remains_to_integrate)){ 906 // ?fails for x/(2*i)/(exp(i*x)+exp(-i*x))*2*exp(i*x)+(-i)*x/(exp((-i)*x)^2+1) 907 gen tmp=_exp2trig(remains_to_integrate,contextptr),r,i; 908 reim(tmp,r,i,contextptr); 909 if (is_zero(ratnormal(i,contextptr))) 910 remains_to_integrate=ratnormal(r,contextptr); 911 else { 912 res=0; 913 remains_to_integrate=e_orig; 914 } 915 } 916 } 917 vector<const unary_function_ptr *> SiCiexp(1,at_Si); 918 SiCiexp.push_back(at_Ci); 919 SiCiexp.push_back(at_exp); 920 if (!lop(res,SiCiexp).empty()){ 921 res=recursive_normal(res,contextptr); 922 if (!has_i(e_orig) && has_i(res)){ 923 res=_exp2trig(res,contextptr); 924 res=recursive_normal(res,contextptr); 925 if (has_i(res)) 926 res=recursive_normal(re(halftan(res,contextptr),contextptr),contextptr); 927 } 928 } 929 return res; 930 } 931 _risch(const gen & g,GIAC_CONTEXT)932 gen _risch(const gen & g,GIAC_CONTEXT){ 933 if ( g.type==_STRNG && g.subtype==-1) return g; 934 if (g.type!=_VECT) 935 return _risch(vecteur(1,g),contextptr); 936 vecteur & v=*g._VECTptr; 937 int s=int(v.size()); 938 if (s>2) 939 return _integrate(g,contextptr); 940 gen tmp; 941 gen var=x__IDNT_e; 942 if (s==2 && v.back().type==_IDNT) 943 var=v.back(); 944 gen res=risch(v.front(),*var._IDNTptr,tmp,contextptr); 945 if (is_zero(tmp)) 946 return res; 947 return res+symbolic(at_integrate,makesequence(tmp,var)); 948 } 949 static const char _risch_s []="risch"; 950 static define_unary_function_eval (__risch,&_risch,_risch_s); 951 define_unary_function_ptr5( at_risch ,alias_at_risch ,&__risch,0,true); 952 953 // integer roots of a polynomial iroots(const polynome & p)954 vecteur iroots(const polynome & p){ 955 int s=p.dim; 956 vecteur zerozero(s-1); 957 vecteur P0(polynome2poly1(p/lgcd(p),1)); 958 vecteur P(P0); 959 // eval every coeff at (0,...,0) 960 int d=int(P.size()); 961 for (int i=0;i<d;++i){ 962 if (P[i].type==_POLY) 963 P[i]=peval(*P[i]._POLYptr,zerozero,0); 964 } 965 // now search the integer roots of this polynomial 966 polynome p0(poly12polynome(P,1,1)); 967 polynome p1=p0.derivative(); 968 p1=gcd(p0,p1); 969 p0=p0/p1; // p0 is now squarefree with the same roots as initial p0 970 // check that all coeffs are integer, if not call normal factorizatio 971 vector< monomial<gen> >::const_iterator it=p0.coord.begin(),itend=p0.coord.end(); 972 vecteur res; 973 for (;it!=itend;++it){ 974 #ifndef HAVE_LIBNTL // with LIBNTL, linearfind does nothing! 975 if (!is_integer(it->value)){ 976 #endif 977 factorization vden; 978 gen extra_div=1; 979 factor(p0,p1,vden,false,false,false,1,extra_div); 980 factorization::const_iterator f_it=vden.begin(),f_itend=vden.end(); 981 // bool ok=true; 982 for (;f_it!=f_itend;++f_it){ 983 int deg=f_it->fact.lexsorted_degree(); 984 if (deg!=1) 985 continue; 986 // extract the root 987 vecteur vtmp=polynome2poly1(f_it->fact,1); 988 gen root=-vtmp.back()/vtmp.front(); 989 if (root.type==_INT_) 990 res.push_back(root); 991 } 992 return res; 993 #ifndef HAVE_LIBNTL 994 } 995 #endif 996 } 997 environment * env=new environment; 998 polynome temp(1); 999 vectpoly v; 1000 int ithprime=1; 1001 if (!linearfind(p0,env,temp,v,ithprime)) // FIXME?? 1002 res.clear();// int bound= 1003 delete env; 1004 d=int(v.size()); 1005 for (int i=0;i<d;++i){ 1006 vecteur tmpv=polynome2poly1(v[i]); 1007 if (tmpv.size()!=2) 1008 continue; 1009 gen g=-tmpv[1]/tmpv[0]; 1010 if (g.type!=_INT_) 1011 continue; 1012 gen tmp=horner(P0,g); 1013 if (is_zero(tmp)) 1014 res.push_back(g); 1015 } 1016 return res; 1017 } 1018 #ifndef NO_NAMESPACE_GIAC 1019 } // namespace giac 1020 #endif // ndef NO_NAMESPACE_GIAC 1021 1022