1 // -*- mode:C++ ; compile-command: "g++ -I.. -I../include -g -c -Wall -DHAVE_CONFIG_H -DIN_GIAC csturm.cc" -*- 2 #include "giacPCH.h" 3 4 /* 5 * Copyright (C) 2000,14 B. Parisse, Institut Fourier, 38402 St Martin d'Heres 6 * 7 * This program is free software; you can redistribute it and/or modify 8 * it under the terms of the GNU General Public License as published by 9 * the Free Software Foundation; either version 3 of the License, or 10 * (at your option) any later version. 11 * 12 * This program is distributed in the hope that it will be useful, 13 * but WITHOUT ANY WARRANTY; without even the implied warranty of 14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 15 * GNU General Public License for more details. 16 * 17 * You should have received a copy of the GNU General Public License 18 * along with this program. If not, see <http://www.gnu.org/licenses/>. 19 */ 20 using namespace std; 21 #include <cmath> 22 #include <stdexcept> 23 #include <map> 24 #include "gen.h" 25 #include "csturm.h" 26 #include "vecteur.h" 27 #include "modpoly.h" 28 #include "unary.h" 29 #include "symbolic.h" 30 #include "usual.h" 31 #include "sym2poly.h" 32 #include "solve.h" 33 #include "prog.h" 34 #include "subst.h" 35 #include "permu.h" 36 #include "series.h" 37 #include "alg_ext.h" 38 #include "ti89.h" 39 #include "plot.h" 40 #include "modfactor.h" 41 #include"giacintl.h" 42 43 #ifndef NO_NAMESPACE_GIAC 44 namespace giac { 45 #endif // ndef NO_NAMESPACE_GIAC 46 47 // compute Sturm sequence of r0 and r1, 48 // returns gcd (without content) 49 // and compute list of quotients, coeffP, coeffR 50 // such that coeffR*r_(k+2) = Q_k*r_(k+1) - coeffP_k*r_k csturm_seq(modpoly & r0,modpoly & r1,vecteur & listquo,vecteur & coeffP,vecteur & coeffR,GIAC_CONTEXT)51 gen csturm_seq(modpoly & r0,modpoly & r1,vecteur & listquo,vecteur & coeffP, vecteur & coeffR,GIAC_CONTEXT){ 52 listquo.clear(); 53 coeffP.clear(); 54 coeffR.clear(); 55 if (r0.empty()) 56 return r1; 57 if (r1.empty()) 58 return r0; 59 gen tmp; 60 lcmdeno(r0,tmp,contextptr); 61 if (ck_is_positive(-tmp,contextptr)) 62 r0=-r0; 63 r0=r0/abs(lgcd(r0),contextptr); 64 lcmdeno(r1,tmp,contextptr); 65 if (ck_is_positive(-tmp,contextptr)) 66 r1=-r1; 67 r1=r1/abs(lgcd(r0),contextptr); 68 // set auxiliary constants g and h to 1 69 gen g(1),h(1); 70 modpoly a(r0),b(r1),quo,r; 71 gen b0(1); 72 for (int loop_counter=0;;++loop_counter){ 73 int m=int(a.size())-1; 74 int n=int(b.size())-1; 75 int ddeg=m-n; // should be 1 generically 76 if (!n) { // if b is constant, gcd=1 77 return 1; 78 } 79 b0=b.front(); 80 if (b.front().type==_VECT) { 81 // ddeg should be even if b0 is a _POLY1 82 if (ddeg%2==0) 83 *logptr(contextptr) << gettext("Singular parametric Sturm sequence ") << a << "/" << b << '\n'; 84 } 85 else 86 b0=abs(b.front(),contextptr); 87 coeffP.push_back(pow(b0,ddeg+1)); 88 DivRem(coeffP.back()*a,b,0,quo,r); 89 listquo.push_back(quo); 90 coeffR.push_back(g*pow(h,ddeg)); 91 if (r.empty()){ 92 return b/abs(lgcd(b),contextptr); 93 } 94 // remainder is non 0, loop continue: a <- b 95 a=b; 96 // now divides r by g*h^(m-n) and change sign, result is the new b 97 b= -r/coeffR.back(); 98 g=b0; 99 h=pow(b0,ddeg)/pow(h,ddeg-1); 100 } // end while loop 101 } 102 csturm_horner(const modpoly & p,const gen & a)103 static gen csturm_horner(const modpoly & p,const gen & a){ 104 if (p.size()==1 && p.front().type==_POLY && p.front()._POLYptr->dim==1){ 105 // patch for "sparse modpoly" 106 vector< monomial<gen> >::const_iterator it=p.front()._POLYptr->coord.begin(),itend=p.front()._POLYptr->coord.end(); 107 gen res=0,anum=a,aden=1,den=1; 108 int oldpui=0,pui; 109 if (a.type==_FRAC){ 110 anum=a._FRACptr->num; 111 aden=a._FRACptr->den; 112 } 113 for (;it!=itend;++it){ 114 pui=it->index.front(); 115 if (oldpui){ 116 res = res * pow(anum,oldpui-pui,context0); 117 den = den * pow(aden,oldpui-pui,context0); 118 } 119 res += it->value*den; 120 oldpui=pui; 121 } 122 if (oldpui){ 123 res = res *pow(a,oldpui,context0); 124 den = den *pow(a,oldpui,context0); 125 } 126 return res/den; 127 } 128 else 129 return horner(p,a); 130 } 131 csturm_vertex_ab(const modpoly & r0,const modpoly & r1,const vecteur & listquo,const vecteur & coeffP,const vecteur & coeffR,const gen & a,int start,GIAC_CONTEXT)132 static int csturm_vertex_ab(const modpoly & r0,const modpoly & r1,const vecteur & listquo,const vecteur & coeffP, const vecteur & coeffR,const gen & a,int start,GIAC_CONTEXT){ 133 int n=int(listquo.size()),j,k,res=0; 134 vecteur R(n+2); 135 R[0]=csturm_horner(r0,a); 136 R[1]=csturm_horner(r1,a); 137 for (j=0;j<n;j++) 138 R[j+2]=(-coeffP[j]*R[j]+R[j+1]*horner(listquo[j],a))/coeffR[j]; 139 // signes 140 for (j=start;j<n+2;j++){ 141 if (R[j]!=0) break; 142 } 143 for (k=j+1;k<n+2;k++){ 144 if (is_zero(R[k])) continue; 145 if (fastsign(R[j],context0)*fastsign(R[k],context0)<0 146 // is_positive(-R[j]*R[k],contextptr) 147 ){ 148 res++; 149 j=k; 150 } 151 } 152 return res; 153 } 154 155 #if 0 156 // compute vertex index at a (==0 unless s(a)==0) 157 static int csturm_vertex_a(const modpoly & s,const modpoly & r,const gen & a,int direction,GIAC_CONTEXT){ 158 int j; 159 modpoly s1,s2; 160 gen sa=horner(s,a,0,s1); 161 if (!is_zero(sa)) return 0; 162 for (j=1;;j++){ 163 sa=horner(s1,a,0,s2); 164 if (!is_zero(sa)) 165 break; 166 s1=s2; 167 } 168 if (direction==1) j=0; 169 gen tmp=sign(sa,contextptr)*sign(horner(r,a),contextptr); 170 return tmp.val*((j%2)?-1:1); 171 } 172 #endif 173 change_scale(modpoly & p,const gen & l)174 void change_scale(modpoly & p,const gen & l){ 175 int n=int(p.size()); 176 gen lton(l); 177 for (int i=n-2;i>=0;--i){ 178 p[i] = p[i] * lton; 179 lton = lton * l; 180 } 181 } 182 back_change_scale(modpoly & p,const gen & l)183 void back_change_scale(modpoly & p,const gen & l){ 184 int n=int(p.size()); 185 gen lton(l); 186 for (int i=n-2;i>=0;--i){ 187 p[i] = p[i] / lton; 188 lton = lton * l; 189 } 190 } 191 192 // p(x)->p(a*x+b) linear_changevar(const modpoly & p,const gen & a,const gen & b)193 modpoly linear_changevar(const modpoly & p,const gen & a,const gen & b){ 194 modpoly res(taylor(p,b)); 195 change_scale(res,a); 196 return res; 197 } 198 199 // p(a*x+b)->p(x) 200 // t=a*x+b -> pgcd(t)=g((t-b)/a) inv_linear_changevar(const modpoly & p,const gen & a,const gen & b)201 modpoly inv_linear_changevar(const modpoly & p,const gen & a,const gen & b){ 202 gen A=inv(a,context0); 203 gen B=-b/a; 204 modpoly res(taylor(p,B)); 205 change_scale(res,A); 206 return res; 207 } 208 209 // Find roots of R, S=R' at precision eps, returns number of roots 210 // if eps==0 does not compute intervals for roots csturm_realroots(const modpoly & S,const modpoly & R,const vecteur & listquo,const vecteur & coeffP,const vecteur & coeffR,const gen & a,const gen & b,const gen & t0,const gen & t1,vecteur & realroots,double eps,GIAC_CONTEXT)211 static int csturm_realroots(const modpoly & S,const modpoly & R,const vecteur & listquo,const vecteur & coeffP, const vecteur & coeffR,const gen & a,const gen & b,const gen & t0, const gen & t1,vecteur & realroots,double eps,GIAC_CONTEXT){ 212 if (is_inf(t0)) // replace with max(R) 213 return csturm_realroots(S,R,listquo,coeffP,coeffR,a,b,-linfnorm(R,contextptr),t1,realroots,eps,contextptr); 214 if (is_inf(t1)) // replace with max(R) 215 return csturm_realroots(S,R,listquo,coeffP,coeffR,a,b,t0,linfnorm(R,contextptr),realroots,eps,contextptr); 216 int n1=csturm_vertex_ab(S,R,listquo,coeffP,coeffR,t0,1,contextptr); 217 int n2=csturm_vertex_ab(S,R,listquo,coeffP,coeffR,t1,1,contextptr); 218 int n=(n2-n1); 219 if (!eps || !n) 220 return n; 221 /* disabled localization of roots, do isolation of roots instead 222 if (is_strictly_greater(eps,(t1-t0)*abs(b,contextptr),contextptr)){ 223 realroots.push_back(makevecteur(makevecteur(a+t0*b,a+t1*b),n)); 224 return n; 225 } 226 */ 227 if (n==1){ 228 gen T0=t0,T1=t1,T2; 229 int s0=fastsign(csturm_horner(R,T0),contextptr); 230 // int s1=fastsign(csturm_horner(R,T1),contextptr); 231 int s2; 232 gen delta=evalf_double(log((T1-T0)*abs(b,contextptr)/eps,contextptr)/log(2.,contextptr),1,contextptr); 233 if (delta.type!=_DOUBLE_){ 234 realroots=vecteur(1,gentypeerr(contextptr)); 235 return -2; 236 } 237 int nstep=int(delta._DOUBLE_val+1); 238 for (int step=0;step<nstep;++step){ 239 T2=(T0+T1)/2; 240 s2=fastsign(csturm_horner(R,T2),contextptr); 241 if (!s2){ 242 realroots.push_back(makevecteur(a+T2*b,n)); 243 return n; 244 } 245 if (s2==s0) 246 T0=T2; 247 else 248 T1=T2; 249 } 250 realroots.push_back(makevecteur(makevecteur(a+T0*b,a+T1*b),n)); 251 return n; 252 } 253 gen T0=t0,T1=t1,t01; 254 for (;;){ 255 t01=(T0+T1)/2; 256 int n01=csturm_vertex_ab(S,R,listquo,coeffP,coeffR,t01,1,contextptr); 257 if (n01!=n1 && n01!=n2) 258 break; 259 if (n01==n1) 260 T0=t01; 261 else 262 T1=t01; 263 } 264 if (csturm_realroots(S,R,listquo,coeffP,coeffR,a,b,T0,t01,realroots,eps,contextptr)==-2) 265 return -2; 266 if (csturm_realroots(S,R,listquo,coeffP,coeffR,a,b,t01,T1,realroots,eps,contextptr)==-2) 267 return -2; 268 return n; 269 } 270 271 // Find complex sturm sequence for P(a+(b-a)*x) 272 // If P is "pseudo"-real on [a,b] and eps>0 put roots in [a,b] 273 // at precision eps inside realroots 274 // returns a,b,R,S,g,listquo,coeffP,coeffR,typeseq 275 // with typeseq=0 (complex Sturm) or 1 (limit) 276 // If b-a is real and horiz_sturm is not empty, it tries to replace 277 // the variable by im(a)*i in horiz_sturm and if no quotient in horiz_sturm 278 // has a leading 0 coefficient, 279 // it returns im(a)*i,im(a)*i+1,R,S,g,listquo,coeffP,coeffR,typeseq 280 // If b-a is pure imaginary and vert_sturm is not empty, it tries to replace 281 // the variable by re(a) and returns re(a),re(a)+i,R,S,g,listquo,coeffP,coeffR,typeseq csturm_segment_seq(const modpoly & P,const gen & a,const gen & b,vecteur & realroots,double eps,vecteur & horiz_sturm,vecteur & vert_sturm,GIAC_CONTEXT)282 static vecteur csturm_segment_seq(const modpoly & P,const gen & a,const gen & b,vecteur & realroots,double eps,vecteur & horiz_sturm,vecteur & vert_sturm,GIAC_CONTEXT){ 283 // try with horiz_sturm and vert_sturm 284 gen ab(b-a); 285 /* // Optimization fails for sturmab(x^3-1,-1-i,1+i) 286 if (is_zero(re(ab,contextptr))){ // b-a is pure imaginary 287 if (vert_sturm.empty()){ 288 gen A=gen(makevecteur(1,0),_POLY1__VECT); 289 vert_sturm.push_back(undef); 290 vecteur tmp; 291 vert_sturm=csturm_segment_seq(P,A,A+cst_i,tmp,eps,horiz_sturm,vert_sturm,contextptr); 292 if (is_undef(vert_sturm)) 293 return vert_sturm; 294 } 295 if (vert_sturm.size()==9){ 296 vecteur res(vert_sturm); 297 gen A=re(a,contextptr); 298 res[0]=A; // re(a) 299 res[1]=A+cst_i; // re(a)+i 300 res[2]=apply1st(res[2],A,horner); // R 301 res[3]=apply1st(res[3],A,horner); // S 302 res[4]=horner(res[4],A); // g 303 vecteur tmp(*res[5]._VECTptr); 304 int tmps=tmp.size(); 305 for (int j=0;j<tmps;++j) 306 tmp[j]=apply1st(tmp[j],A,horner); 307 res[5]=tmp; // listquo 308 res[6]=apply1st(res[6],A,horner); // coeffP 309 res[7]=apply1st(res[7],A,horner); // coeffR 310 if (res[6].type==_VECT && !equalposcomp(*res[6]._VECTptr,0)) 311 return res; 312 else 313 CERR << "list of quotients is not regular" << '\n'; 314 } 315 } 316 */ 317 modpoly Q(taylor(P,a)); 318 change_scale(Q,b-a); 319 // now x is in 0..1 320 gen gtmp=apply(Q,re,contextptr); 321 if (gtmp.type!=_VECT) 322 return vecteur(1,gensizeerr(contextptr)); 323 modpoly R=trim(*gtmp._VECTptr,0); 324 gtmp=apply(Q,im,contextptr); 325 if (gtmp.type!=_VECT) 326 return vecteur(1,gensizeerr(contextptr)); 327 modpoly S=trim(*gtmp._VECTptr,0); 328 modpoly listquo,coeffP,coeffR; 329 gen g=csturm_seq(S,R,listquo,coeffP,coeffR,contextptr); 330 int typeseq=-1; 331 if (debug_infolevel) 332 *logptr(contextptr) << "segment " << a << ".." << b << ", im/re:" << S << "|" << R << ", gcd:" << g << '\n'; 333 if (g.type==_VECT && g._VECTptr->size()==P.size()){ 334 // if g==P (up to a constant), use real Sturm sequences 335 if (debug_infolevel) 336 *logptr(contextptr) << "Real-kind roots: " << g << '\n'; 337 R=*g._VECTptr; 338 S=derivative(R); 339 g=csturm_seq(S,R,listquo,coeffP,coeffR,contextptr); 340 typeseq=csturm_realroots(S,R,listquo,coeffP,coeffR,a,b-a,0,1,realroots,eps,contextptr); 341 if (typeseq==-2) 342 return realroots; 343 } 344 if (g.type==_VECT) 345 g=inv_linear_changevar(*g._VECTptr,b-a,a); 346 vecteur res= makevecteur(a,b,R,S,g,listquo,coeffP,coeffR,typeseq); 347 return res; 348 } 349 350 // index for segment a,b (2* number of roots when summed over a closed 351 // polygon). Note that if S=ImP along the segment is 0 we remove 352 // the roots on [a,b] using real Sturm sequences 353 // If S=0 at a or b, this is simply ignored 354 // Indeed the computed index is then the same as if S was of the 355 // sign of R, and since R!=0 if S is 0 this is a property of the vertex 356 // not of the segment (note that contrary to counting real roots 357 // on an interval, S can vanish as many times as long as R keeps 358 // the same sign, without modifying the algebraic number of Im=0 359 // cuts if S has the same sign on both end) csturm_segment(const vecteur & seq,const gen & a,const gen & b,GIAC_CONTEXT)360 static int csturm_segment(const vecteur & seq,const gen & a,const gen & b,GIAC_CONTEXT){ 361 gen t0,t1; 362 if (seq.size()!=9) 363 return -(RAND_MAX/2); 364 gen aseq=seq[0]; 365 gen bseq=seq[1]; 366 gen directeur=(b-a)/(bseq-aseq); 367 t0=(a-aseq)/(bseq-aseq); 368 if ( !is_zero(im(directeur,contextptr)) || !is_zero(im(t0,contextptr)) ) 369 return -(RAND_MAX/2); 370 t0=re(t0,contextptr); // t0=normal(t0); 371 t1=re(t0+directeur,contextptr); // t1=normal(t0+directeur); 372 int signe=1; 373 if (is_strictly_greater(t0,t1,contextptr)){ 374 signe=-1; 375 swapgen(t0,t1); 376 } 377 const modpoly & R=*seq[2]._VECTptr; 378 const modpoly & S=*seq[3]._VECTptr; 379 gen g=seq[4]; 380 const modpoly & listquo=*seq[5]._VECTptr; 381 const modpoly & coeffP=*seq[6]._VECTptr; 382 const modpoly & coeffR=*seq[7]._VECTptr; 383 int debut=(seq[8].val==-1)?0:1; 384 int tmp = csturm_vertex_ab(S,R,listquo,coeffP,coeffR,t0,debut,contextptr); 385 int res = tmp; 386 tmp = csturm_vertex_ab(S,R,listquo,coeffP,coeffR,t1,debut,contextptr); 387 res -= tmp; 388 // tmp = (-csturm_vertex_a(S,R,t0,1,contextptr)+csturm_vertex_a(S,R,t1,-1,contextptr)); 389 // res += tmp; 390 res=(debut?1:signe)*res; 391 if (debug_infolevel) 392 *logptr(contextptr) << "segment " << a << ".." << b << " index contribution " << res << '\n'; 393 return res; 394 } 395 csturm_square_seq(const modpoly & P,const gen & a0,const gen & b0,const gen & a1,const gen & b1,gen & pgcd,vecteur & realroots,double eps,vecteur & seq1,vecteur & seq2,vecteur & seq3,vecteur & seq4,vecteur & horiz_sturm,vecteur & vert_sturm,GIAC_CONTEXT)396 static bool csturm_square_seq(const modpoly & P,const gen & a0,const gen & b0,const gen & a1,const gen & b1,gen & pgcd,vecteur & realroots,double eps,vecteur & seq1,vecteur & seq2,vecteur & seq3,vecteur & seq4,vecteur & horiz_sturm,vecteur & vert_sturm,GIAC_CONTEXT){ 397 gen A=a0+cst_i*b0,B=a1+cst_i*b0; 398 vecteur rroots; 399 seq1=csturm_segment_seq(P,A,B,rroots,eps,horiz_sturm,vert_sturm,contextptr); 400 if (is_undef(seq1)) 401 return false; 402 pgcd=seq1[4]; 403 if (!is_one(pgcd)){ 404 return false; 405 } 406 A=a1+cst_i*b0; B=a1+cst_i*b1; 407 seq2=csturm_segment_seq(P,A,B,rroots,eps,horiz_sturm,vert_sturm,contextptr); 408 if (is_undef(seq2)) 409 return false; 410 pgcd=seq2[4]; 411 if (!is_one(pgcd)){ 412 return false; 413 } 414 A=a1+cst_i*b1; B=a0+cst_i*b1; 415 seq3=csturm_segment_seq(P,A,B,rroots,eps,horiz_sturm,vert_sturm,contextptr); 416 if (is_undef(seq3)) 417 return false; 418 pgcd=seq3[4]; 419 if (!is_one(pgcd)){ 420 return false; 421 } 422 A=a0+cst_i*b1; B=a0+cst_i*b0; 423 seq4=csturm_segment_seq(P,A,B,rroots,eps,horiz_sturm,vert_sturm,contextptr); 424 if (is_undef(seq4)) 425 return false; 426 pgcd=seq4[4]; 427 if (!is_one(pgcd)){ 428 return false; 429 } 430 realroots=mergevecteur(realroots,rroots); 431 return true; 432 } 433 434 // find 2* number of roots of P inside the square of vertex of affixes a,b 435 // roots on the square are not counted. P must not vanish at the vertices. 436 // The complex Sturm sequences must be known 437 // returns -1 on error csturm_square(const modpoly & P,const gen & a0,const gen & b0,const gen & a1,const gen & b1,const vecteur & seq1,const vecteur & seq2,const vecteur & seq3,const vecteur & seq4,GIAC_CONTEXT)438 static int csturm_square(const modpoly & P,const gen & a0,const gen & b0,const gen & a1,const gen & b1,const vecteur & seq1,const vecteur & seq2,const vecteur & seq3,const vecteur & seq4,GIAC_CONTEXT){ 439 int ind,tmp; 440 ind = 0; 441 gen A=a0+cst_i*b0,B=a1+cst_i*b0; 442 tmp = csturm_segment(seq1,A,B,contextptr); 443 if (tmp==-(RAND_MAX/2)) 444 return -1; 445 ind += tmp; 446 A=a1+cst_i*b0; B=a1+cst_i*b1; 447 tmp = csturm_segment(seq2,A,B,contextptr); 448 if (tmp==-(RAND_MAX/2)) 449 return -1; 450 ind += tmp; 451 A=a1+cst_i*b1; B=a0+cst_i*b1; 452 tmp = csturm_segment(seq3,A,B,contextptr); 453 if (tmp==-(RAND_MAX/2)) 454 return -1; 455 ind += tmp; 456 A=a0+cst_i*b1; B=a0+cst_i*b0; 457 tmp = csturm_segment(seq4,A,B,contextptr); 458 if (tmp==-(RAND_MAX/2)) 459 return -1; 460 ind += tmp; 461 return ind; 462 } 463 csturm_normalize(modpoly & p,const gen & a0,const gen & b0,const gen & a1,const gen & b1,vecteur & roots)464 static void csturm_normalize(modpoly & p,const gen & a0,const gen & b0,const gen & a1,const gen & b1,vecteur & roots){ 465 int n=int(p.size())-1; 466 // Make sure that x->a+i*x does not return a multiple 467 // of a real polynomial with the multiple non real 468 // If degree of p is even the multiple will be a real (because of lcoeff) 469 if (n%2){ 470 // If degree is odd then look at q=p(x-a_{n-1}/n*an) 471 // it has the same property 472 // if its cst coeff is zero remove 473 gen an=p.front(); 474 gen b=p[1]; 475 gen shift=-b/n/an; 476 modpoly q(taylor(p,shift)); 477 gen q0; 478 // remove valuation 479 int qs=int(q.size()); 480 int n1=0; 481 for (;qs>0;--qs,++n1){ 482 if (!is_zero(q0=q[qs-1])) 483 break; 484 } 485 if (is_zero(re(q0,context0))){ 486 q=cst_i*q; 487 p=cst_i*p; 488 } 489 if (n1){ 490 q=modpoly(q.begin(),q.begin()+qs); 491 gen a=re(shift,context0),b=im(shift,context0); 492 if (is_greater(a,a0,context0) && is_greater(b,b0,context0) && is_greater(a1,a,context0) && is_greater(b1,b,context0)) 493 roots.push_back(makevecteur(shift,n1)); 494 p=taylor(q,-shift); 495 } 496 } 497 } 498 ab2a0b0a1b1(const gen & a,const gen & b,gen & a0,gen & b0,gen & a1,gen & b1,GIAC_CONTEXT)499 void ab2a0b0a1b1(const gen & a,const gen & b,gen & a0,gen & b0,gen & a1,gen & b1,GIAC_CONTEXT){ 500 a0=re(a,contextptr); b0=im(a,contextptr); 501 a1=re(b,contextptr); b1=im(b,contextptr); 502 if (ck_is_greater(a0,a1,contextptr)) swapgen(a0,a1); 503 if (ck_is_greater(b0,b1,contextptr)) swapgen(b0,b1); 504 } 505 506 // find 2* number of roots of P inside the square of vertex of affixes a,b 507 // excluding those on the square 508 // returns -1 on error csturm_square(const gen & p,const gen & a,const gen & b,gen & pgcd,GIAC_CONTEXT)509 int csturm_square(const gen & p,const gen & a,const gen & b,gen& pgcd,GIAC_CONTEXT){ 510 if (p.type==_POLY){ 511 int res=0; 512 factorization f(sqff(*p._POLYptr)); 513 factorization::const_iterator it=f.begin(),itend=f.end(); 514 for (;it!=itend;++it){ 515 int tmp=csturm_square(polynome2poly1(it->fact),a,b,pgcd,contextptr); 516 if (tmp==-1) 517 return -1; 518 res += it->mult*tmp; 519 } 520 return res; 521 } 522 if (p.type!=_VECT) 523 return 0; 524 modpoly P=*p._VECTptr; 525 vecteur realroots; 526 gen a0,b0,a1,b1; 527 ab2a0b0a1b1(a,b,a0,b0,a1,b1,contextptr); 528 csturm_normalize(P,a0,b0,a1,b1,realroots); 529 int evident=0; 530 if (!realroots.empty()){ 531 gen r=realroots.front(); 532 if (r.type==_VECT && r._VECTptr->size()==2) 533 r=r._VECTptr->front(); 534 gen rx=re(r,contextptr),ry=im(r,contextptr); 535 if ( ( is_zero(ry) && (rx==a0 || rx==a1) ) || 536 ( is_zero(rx) && (ry==b0 || ry==b1) ) ) 537 ; 538 else 539 evident=1; 540 } 541 if (P.size()<2) 542 return evident; 543 vecteur seq1,seq2,seq3,seq4,horiz_seq,vert_seq; 544 if (!csturm_square_seq(P,a0,b0,a1,b1,pgcd,realroots,0.0,seq1,seq2,seq3,seq4,horiz_seq,vert_seq,contextptr)){ 545 if (pgcd.type!=_VECT) 546 return -1; 547 modpoly g=(*pgcd._VECTptr)/pgcd[0]; 548 // true factorization found, restart with each factor 549 modpoly p1=P/g; 550 int n1=csturm_square(p1,a,b,pgcd,contextptr); 551 if (n1==-1) 552 return -1; 553 int n2=csturm_square(g,a,b,pgcd,contextptr); 554 if (n2==-1) 555 return -1; 556 return evident+n1+n2; 557 } 558 int tmp=csturm_square(P,a0,b0,a1,b1,seq1,seq2,seq3,seq4,contextptr); 559 if (tmp==-1) 560 return tmp; 561 return evident+tmp; 562 } 563 564 static void complex_roots(const modpoly & P,const gen & a0,const gen & b0,const gen & a1,const gen & b1,vecteur & realroots,vecteur & complexroots,double eps); 565 complex_roots_split(const modpoly & P,const gen & pgcd,const gen & a0,const gen & b0,const gen & a1,const gen & b1,vecteur & realroots,vecteur & complexroots,double eps)566 static bool complex_roots_split(const modpoly & P,const gen & pgcd,const gen & a0,const gen & b0,const gen & a1,const gen & b1,vecteur & realroots,vecteur & complexroots,double eps){ 567 if (pgcd.type!=_VECT) 568 return false; 569 modpoly g=(*pgcd._VECTptr)/pgcd[0]; 570 // true factorization found, restart with each factor 571 modpoly p1=P/g; 572 csturm_normalize(p1,a0,b0,a1,b1,realroots); 573 csturm_normalize(g,a0,b0,a1,b1,realroots); 574 complex_roots(p1,a0,b0,a1,b1,realroots,complexroots,eps); 575 complex_roots(g,a0,b0,a1,b1,realroots,complexroots,eps); 576 return true; 577 } 578 579 #if 0 580 // check that arg is >=pi/8 (assumes im(g)>=0) 581 static bool arg_geq_pi_8(const gen & g){ 582 gen gr=re(g,context0),gi=im(g,context0); 583 if (is_positive(-gr,context0)) 584 return true; 585 // ? gi/gr>=sqrt(2)-1 586 gen r=gi/gr+1; 587 if (is_positive(r*r-2,context0)) 588 return true; 589 return false; 590 } 591 592 // is im(b/a)>=0, tested without quotient 593 static bool arg_in_0_pi(const gen & a,const gen & b){ 594 gen A(a),B(b); 595 if (A.type==_FRAC && is_integer(A._FRACptr->den) && is_positive(A._FRACptr->den,context0)) 596 A=A._FRACptr->num; 597 if (B.type==_FRAC && is_integer(B._FRACptr->den) && is_positive(B._FRACptr->den,context0)) 598 B=B._FRACptr->num; 599 gen c=re(A,context0)*im(B,context0)+re(B,context0)*im(A,context0); 600 return is_positive(c,context0); 601 } 602 603 static gen hornerarg(const modpoly & p,const gen & x){ 604 if (p.empty()) 605 return 0; 606 if (x.type!=_FRAC || !is_integer(x._FRACptr->den)) 607 return horner(p,x); 608 fraction & f =*x._FRACptr; 609 gen num=f.num,den=f.den,d=den; 610 if (is_positive(-f.den,context0)){ 611 num=-num; den=-den; d=den; 612 } 613 modpoly::const_iterator it=p.begin(),itend=p.end(); 614 gen res(*it); 615 ++it; 616 if (it==itend) 617 return res; 618 for (;;){ 619 res=res*num+(*it)*d; 620 ++it; 621 if (it==itend) 622 break; 623 d=d*den; 624 } 625 return res; 626 } 627 628 // Find one complex root inside a0,b0->a1,b1, return false if not found 629 static bool complex_1root(const modpoly & P,const gen & a0,const gen & b0,const gen & a1,const gen & b1,vecteur & complexroots,double eps){ 630 return false; // disabled since it is not faster!! 631 int step,nstep =int(evalf_double(log(max(a1-a0,b1-b0,context0)/eps,context0)/log(2.0,context0),1,context0)._DOUBLE_val+0.5); 632 if (nstep<4) 633 return false; 634 // First compute P at the 4 vertex and check whether P[vertex_n+1]/P[vertex_n] is in C^+ 635 gen P0=hornerarg(P,a0+cst_i*b0),P2=hornerarg(P,a1+cst_i*b0), 636 P4=hornerarg(P,a1+cst_i*b1),P6=hornerarg(P,a0+cst_i*b1); 637 if (!(arg_in_0_pi(P0,P2) && arg_in_0_pi(P2,P4) && arg_in_0_pi(P4,P6) && arg_in_0_pi(P6,P0))) 638 return false; 639 gen A0(a0),A2(a1),B0(b0),B2(b1),A1,B1; 640 for (step=0;step<2*nstep;step++){ 641 A1=(A0+A2)/2; 642 B1=(B0+B2)/2; 643 gen P1=hornerarg(P,A1+cst_i*B0),P7=hornerarg(P,A0+cst_i*B1),P8=hornerarg(P,A1+cst_i*B1),P3,P5; 644 bool found=false; 645 /* 646 P6(A0,B2) - P5(A1,B2) - P4(A2,B2) 647 | | | 648 P7(A0,B1) - P8(A1,B1) - P3(A2,B1) 649 | | | 650 P0(A0,B0) - P1(A1,B0) - P2(A2,B0) 651 */ 652 // ? P0, P1, P8, P7 653 if (arg_in_0_pi(P0,P1) && arg_in_0_pi(P1,P8) && arg_in_0_pi(P8,P7) && arg_in_0_pi(P7,P0)){ 654 A2=A1; 655 B2=B1; 656 P2=P1; 657 P4=P8; 658 P6=P7; 659 if (step<nstep) 660 continue; 661 found=true; 662 } 663 if (!found){ 664 P3=hornerarg(P,A2+cst_i*B1); 665 // ? P1, P2, P3, P8 666 if (arg_in_0_pi(P1,P2) && arg_in_0_pi(P2,P3) && arg_in_0_pi(P3,P8) && arg_in_0_pi(P8,P1)){ 667 A0=A1; 668 B2=B1; 669 P0=P1; 670 P4=P3; 671 P6=P8; 672 if (step<nstep) 673 continue; 674 found=true; 675 } 676 } 677 if (!found){ 678 // P8, P3, P4, P5 679 P5=hornerarg(P,A1+cst_i*B2); 680 if (arg_in_0_pi(P8,P3) && arg_in_0_pi(P3,P4) && arg_in_0_pi(P4,P5) && arg_in_0_pi(P5,P8)){ 681 A0=A1; 682 B0=B1; 683 P0=P8; 684 P2=P3; 685 P6=P5; 686 if (step<nstep) 687 continue; 688 found=true; 689 } 690 } 691 if (!found){ 692 // P7 P8 P5 P6 693 if (arg_in_0_pi(P7,P8) && arg_in_0_pi(P8,P5) && arg_in_0_pi(P5,P6) && arg_in_0_pi(P6,P7)){ 694 A2=A1; 695 B0=B1; 696 P0=P7; 697 P2=P8; 698 P4=P5; 699 if (step<nstep) 700 continue; 701 found=true; 702 } 703 } 704 if (!found) 705 return false; 706 // Final check that there is indeed a root inside rectangle 707 // args must be >= pi/8 and degree of (P)*max square length/distance to original square <= pi/8 708 gen dist=min(min(A0-a0,a1-A2,context0),min(B0-b0,b1-B2,context0),context0); 709 if (is_zero(dist)) 710 continue; 711 gen max_sq=max(A2-A0,B2-B0,context0); 712 if (is_greater((int(P.size())-2)*max_sq/dist,cst_pi/8,context0)) 713 continue; 714 gen r1=P2/P0, r2=P4/P2, r3=P6/P4, r4=P0/P6; 715 if (arg_geq_pi_8(r1) && arg_geq_pi_8(r2) && arg_geq_pi_8(r3) && arg_geq_pi_8(r4)){ 716 complexroots.push_back(makevecteur(makevecteur(A0+cst_i*B0,A2+cst_i*B2),1)); 717 return true; 718 } 719 } 720 return false; 721 } 722 #endif 723 round2util(const gen & num,const gen & den,int n)724 static gen round2util(const gen & num,const gen & den,int n){ 725 if (num.type==_CPLX){ 726 gen r=round2util(*num._CPLXptr,den,n); 727 gen i=round2util(*(num._CPLXptr+1),den,n); 728 return r+cst_i*i; 729 } 730 // num must be a _ZINT 731 mpz_t tmp1,tmp2; 732 mpz_init_set(tmp1,*num._ZINTptr); 733 mpz_mul_2exp(tmp1,tmp1,n+1); // tmp1=2^(n+1)*num 734 mpz_add(tmp1,tmp1,*den._ZINTptr); // + den 735 mpz_init_set(tmp2,*den._ZINTptr); 736 mpz_mul_ui(tmp2,tmp2,2); // tmp2=2*den 737 mpz_fdiv_q(tmp1,tmp1,tmp2); 738 gen res=tmp1; 739 mpz_clear(tmp1); mpz_clear(tmp2); 740 return res; 741 } 742 in_round2(gen & x,const gen & deuxn,int n)743 void in_round2(gen & x,const gen & deuxn, int n){ 744 if (x.type==_INT_ || x.type==_ZINT) 745 return ; 746 if (x.type==_FRAC && x._FRACptr->den.type==_CPLX) 747 x=fraction(x._FRACptr->num*conj(x._FRACptr->den,context0),x._FRACptr->den.squarenorm(context0)); 748 if (x.type==_FRAC && x._FRACptr->den.type==_ZINT && 749 (x._FRACptr->num.type==_ZINT || 750 (x._FRACptr->num.type==_CPLX && x._FRACptr->num._CPLXptr->type==_ZINT && (x._FRACptr->num._CPLXptr+1)->type==_ZINT)) 751 ){ 752 gen num=x._FRACptr->num,d=x._FRACptr->den; 753 x=round2util(num,d,n); 754 x=x/deuxn; 755 return; 756 } 757 x=_floor(x*deuxn+plus_one_half,context0)/deuxn; 758 } 759 round2(gen & x,int n)760 void round2(gen & x,int n){ 761 if (x.type==_INT_ || x.type==_ZINT) 762 return ; 763 gen deuxn; 764 if (n<30) 765 deuxn = (1<<n); 766 else { 767 mpz_t tmp; 768 mpz_init_set_si(tmp,1); 769 mpz_mul_2exp(tmp,tmp,n); 770 deuxn=tmp; 771 mpz_clear(tmp); 772 } 773 in_round2(x,deuxn,n); 774 } 775 round2(gen & x,const gen & deuxn,GIAC_CONTEXT)776 void round2(gen & x,const gen & deuxn,GIAC_CONTEXT){ 777 if (x.type==_INT_ || x.type==_ZINT) 778 return; 779 if (x.type!=_FRAC) 780 x=_floor(x*deuxn+plus_one_half,context0)/deuxn; 781 else { 782 gen n=x._FRACptr->num,d=x._FRACptr->den; 783 if (d.type==_INT_){ 784 int di=d.val,ni=1; 785 while (di>1){ di=di>>1; ni=ni<<1;} 786 if (ni==d.val) 787 return; 788 } 789 n=2*n*deuxn+d; 790 x=iquo(n,2*d)/deuxn; 791 } 792 } 793 794 // Find one complex root inside a0,b0->a1,b1, return false if not found 795 // algo: Newton method in exact mode starting from center newton_complex_1root(const modpoly & P,const gen & a0,const gen & b0,const gen & a1,const gen & b1,vecteur & complexroots,double eps)796 bool newton_complex_1root(const modpoly & P,const gen & a0,const gen & b0,const gen & a1,const gen & b1,vecteur & complexroots,double eps){ 797 if (is_positive(a1-a0-0.01,context0) || 798 is_positive(b1-b0-0.01,context0)) 799 return false; 800 gen x0=(a0+a1)/2+cst_i*(b0+b1)/2; 801 modpoly Pprime=derivative(P); 802 int n=int(-std::log(eps)/std::log(2.0)+.5); // for rounding 803 gen eps2=pow(2,-(n+1),context0); 804 for (int ii=0;ii<n;ii++){ 805 gen Pprimex0=horner(Pprime,x0,0,false); 806 if (is_zero(Pprimex0,context0)) 807 return false; 808 gen dx=horner(P,x0,0,false)/Pprimex0; 809 gen absdx=dx*conj(dx,context0); 810 x0=x0-dx; 811 gen r=re(x0,context0),i=im(x0,context0); 812 if (is_positive(a0-r,context0) || is_positive(r-a1,context0) || 813 is_positive(b0-i,context0) || is_positive(i-b1,context0)) 814 return false; 815 round2(r,n+4); 816 round2(i,n+4); 817 x0=r+cst_i*i; 818 if (is_positive(absdx-eps2*eps2,context0)) 819 continue; 820 // make a small square around x0 821 // and check that there is indeed a root inside 822 gen A0=r-eps2; 823 gen A1=r+eps2; 824 gen B0=i-eps2; 825 gen B1=i+eps2; 826 gen tmp; 827 if (csturm_square(P,A0+cst_i*B0,A1+cst_i*B1,tmp,context0)==2){ 828 complexroots.push_back(makevecteur(makevecteur(A0+cst_i*B0,A1+cst_i*B1),1)); 829 return true; 830 } 831 } 832 return false; 833 } 834 835 // Find complex roots of P in a0,b0 -> a1,b1 complex_roots(const modpoly & P,const gen & a0,const gen & b0,const gen & a1,const gen & b1,const vecteur & seq1,const vecteur & seq2,const vecteur & seq3,const vecteur & seq4,vecteur & realroots,vecteur & complexroots,double eps,vecteur & horiz_sturm,vecteur & vert_sturm)836 static int complex_roots(const modpoly & P,const gen & a0,const gen & b0,const gen & a1,const gen & b1,const vecteur & seq1,const vecteur & seq2,const vecteur & seq3,const vecteur & seq4,vecteur & realroots,vecteur & complexroots,double eps,vecteur & horiz_sturm,vecteur & vert_sturm){ 837 int n=csturm_square(P,a0,b0,a1,b1,seq1,seq2,seq3,seq4,context0); 838 if (debug_infolevel && n) 839 CERR << a0 << "," << b0 << ".." << a1 << "," << b1 << ":" << n/2 << '\n'; 840 if (n<=0) 841 return 2*n; 842 if (eps<=0){ 843 *logptr(context0) << gettext("Bad precision, using 1e-12 instead of ")+print_DOUBLE_(eps,14) << '\n'; 844 eps=1e-12; 845 } 846 if (is_strictly_greater(eps,a1-a0,context0) && is_strictly_greater(eps,b1-b0,context0)){ 847 gen r(makevecteur(a0+cst_i*b0,a1+cst_i*b1)); 848 complexroots.push_back(makevecteur(r,gen(n)/2)); 849 return n; 850 } 851 if (n==2 && newton_complex_1root(P,a0,b0,a1,b1,complexroots,eps)) 852 return n; 853 gen a01=(a0+a1)/2,b01=(b0+b1)/2,pgcd; 854 vecteur seqvert,seqhoriz; 855 gen A=a0+cst_i*b01,B=a1+cst_i*b01; 856 seqhoriz=csturm_segment_seq(P,A,B,realroots,eps,horiz_sturm,vert_sturm,context0); 857 if (is_undef(seqhoriz)){ 858 realroots=seqhoriz; 859 return -2; 860 } 861 pgcd=seqhoriz[4]; 862 if (is_one(pgcd)){ 863 A=a01+cst_i*b0; B=a01+cst_i*b1; 864 seqvert=csturm_segment_seq(P,A,B,realroots,eps,horiz_sturm,vert_sturm,context0); 865 if (is_undef(seqvert)){ 866 realroots=seqvert; 867 return -2; 868 } 869 pgcd=seqvert[4]; 870 } 871 if (!is_one(pgcd)){ 872 complex_roots_split(P,pgcd,a0,b0,a1,b1,realroots,complexroots,eps); 873 return n; 874 } 875 /* 876 (a0,b1) - (a01,b1) - (a1,b1) seq3 seq3 877 | n4 | n3 | seq4 n4 seqvert n3 seq2 878 (a0,b01) - (a01,b01) - (a1,b01) seqhoriz seqhoriz 879 | n1 | n2 | seq4 n1 seqvert n2 seq2 880 (a0,b0) - (a01,b0) - (a1,b0) seq1 seq1 881 */ 882 int n1=complex_roots(P,a0,b0,a01,b01,seq1,seqvert,seqhoriz,seq4,realroots,complexroots,eps,horiz_sturm,vert_sturm),nadd; 883 if (n1==-2) 884 return -2; 885 if (n1==n) 886 return n; 887 n1 += (nadd=complex_roots(P,a01,b0,a1,b01,seq1,seq2,seqhoriz,seqvert,realroots,complexroots,eps,horiz_sturm,vert_sturm)); 888 if (nadd==-2) 889 return -2; 890 if (n1==n) 891 return n; 892 n1 += (nadd=complex_roots(P,a01,b01,a1,b1,seqhoriz,seq2,seq3,seqvert,realroots,complexroots,eps,horiz_sturm,vert_sturm)); 893 if (nadd==-2) 894 return -2; 895 if (n1==n) 896 return n; 897 n1 += (nadd=complex_roots(P,a0,b01,a01,b1,seqhoriz,seqvert,seq3,seq4,realroots,complexroots,eps,horiz_sturm,vert_sturm)); 898 if (nadd==-2) 899 return -2; 900 return n; 901 } 902 903 // Find complex roots of P in a0,b0 -> a1,b1 complex_roots(const modpoly & P,const gen & a0,const gen & b0,const gen & a1,const gen & b1,vecteur & realroots,vecteur & complexroots,double eps)904 static void complex_roots(const modpoly & P,const gen & a0,const gen & b0,const gen & a1,const gen & b1,vecteur & realroots,vecteur & complexroots,double eps){ 905 if (P.size()<2) 906 return; 907 vecteur Seq1,Seq2,Seq3,Seq4,horiz_sturm,vert_sturm; 908 gen pgcd; 909 if (!csturm_square_seq(P,a0,b0,a1,b1,pgcd,realroots,eps,Seq1,Seq2,Seq3,Seq4,horiz_sturm,vert_sturm,context0)) 910 complex_roots_split(P,pgcd,a0,b0,a1,b1,realroots,complexroots,eps); 911 else 912 complex_roots(P,a0,b0,a1,b1,Seq1,Seq2,Seq3,Seq4,realroots,complexroots,eps,horiz_sturm,vert_sturm); 913 } 914 915 // Find complex roots of P in a0,b0 -> a1,b1 complex_roots(const modpoly & P,const gen & a0,const gen & b0,const gen & a1,const gen & b1,gen & pgcd,vecteur & roots,double eps)916 bool complex_roots(const modpoly & P,const gen & a0,const gen & b0,const gen & a1,const gen & b1,gen & pgcd,vecteur & roots,double eps){ 917 vecteur realroots,complexroots; 918 complex_roots(P,a0,b0,a1,b1,realroots,complexroots,eps); 919 if (is_undef(realroots)) 920 return false; 921 roots=mergevecteur(roots,mergevecteur(realroots,complexroots)); 922 return true; 923 } 924 crationalroot(polynome & p,bool complexe)925 vecteur crationalroot(polynome & p,bool complexe){ 926 vectpoly v; 927 int i=1; 928 polynome qrem; 929 environment * env= new environment; 930 env->complexe=complexe || !is_zero(im(p,context0)); 931 vecteur w; 932 if (!do_linearfind(p,env,qrem,v,w,i)) 933 w.clear(); 934 delete env; 935 p=qrem; 936 return w; 937 } 938 keep_in_rectangle(const vecteur & croots,const gen A0,const gen & B0,const gen & A1,const gen & B1,bool embed,GIAC_CONTEXT)939 vecteur keep_in_rectangle(const vecteur & croots,const gen A0,const gen & B0,const gen & A1,const gen & B1,bool embed,GIAC_CONTEXT){ 940 vecteur roots; 941 const_iterateur it=croots.begin(),itend=croots.end(); 942 for (;it!=itend;++it){ 943 gen a=re(*it,contextptr),b=im(*it,contextptr); 944 if (is_greater(a,A0,contextptr)&&is_greater(A1,a,contextptr)&&is_greater(b,B0,contextptr)&&is_greater(B1,b,contextptr)) 945 roots.push_back(embed?makevecteur(*it,1):*it); 946 } 947 return roots; 948 } 949 square_modulus(const gen & g,GIAC_CONTEXT)950 gen square_modulus(const gen & g,GIAC_CONTEXT){ 951 return g.squarenorm(contextptr); 952 } 953 954 // P is the polynomial, P1 derivative, v list of approx roots 955 // (initially should have at least n bits precision), 956 // epsn is the target number of bits precision int(std::log(eps)/std::log(2.)-.5); 957 // epsg2surdeg2 is eps^2/degree(P)^2 as a gen, epsg is the target precision 958 // v[i] is set by newton_improve to be at distance at most vradius[i] of a root 959 // kmax is the maximal number of Newton iterations newton_improve(const vecteur & P,const vecteur & P1,bool Preal,vecteur & v,vecteur & vradius,int i,int kmax,int n,int epsn,const gen & epsg2surdeg2,const gen & epsg)960 bool newton_improve(const vecteur & P,const vecteur & P1,bool Preal,vecteur & v,vecteur & vradius,int i,int kmax,int n,int epsn,const gen & epsg2surdeg2,const gen & epsg){ 961 gen r=v[i]; 962 bool nextconj=false; 963 if (Preal && i+1<int(v.size())) 964 nextconj=is_exactly_zero(r-conj(v[i+1],context0)); 965 if (r.type==_FRAC || is_cinteger(r)) 966 return true; 967 // find nearest root from v 968 gen deltar=plus_inf,delta; 969 for (unsigned j=0;j<v.size();++j){ 970 if (int(j)==i) 971 continue; 972 gen tmp=abs(r-v[j],context0); 973 if (is_strictly_greater(deltar,tmp,context0)) 974 deltar=tmp; 975 } 976 if (is_zero(deltar)) 977 return false; 978 deltar=deltar/3; 979 gen sumdr2=0; 980 int N=n; // effective value of number of bits for computation 981 #ifdef HAVE_LIBMPFR 982 if (r.type==_REAL && mpfr_get_prec(r._REALptr->inf)>N) 983 N=mpfr_get_prec(r._REALptr->inf); 984 if (r.type==_CPLX && r._CPLXptr->type==_REAL && mpfr_get_prec(r._CPLXptr->_REALptr->inf)>N) 985 N=mpfr_get_prec(r._CPLXptr->_REALptr->inf); 986 #endif 987 #if 0 // def HAVE_LIBMPFI 988 gen deuxN=pow(2,N,context0); 989 gen rr,ri,dr,di; 990 reim(r,rr,ri,context0); 991 if (Preal && !nextconj) 992 r=eval(gen(makevecteur(rr-plus_one/deuxN,rr+plus_one/deuxN),_INTERVAL__VECT),1,context0); 993 else 994 r=eval(gen(makevecteur(rr-plus_one/deuxN,rr+plus_one/deuxN),_INTERVAL__VECT),1,context0)+cst_i*eval(gen(makevecteur(ri-plus_one/deuxN,ri+plus_one/deuxN),_INTERVAL__VECT),1,context0); 995 for (int k=0;k<kmax;++k){ 996 // check if root precision is ok 997 // otherwise try to improve root precision with Newton method 998 gen P1r=horner(P1,r,0,false); 999 if (is_exactly_zero(P1r)){ 1000 delta=plus_inf; 1001 break; 1002 } 1003 gen Pr=horner(P,r,0,false); 1004 delta=Pr/P1r; 1005 bool test; 1006 if (Preal && ! nextconj){ 1007 test=contains(delta,dr); 1008 delta=abs(delta,context0); 1009 } 1010 else { 1011 reim(delta,dr,di,context0); 1012 test= contains(rr,dr) && contains(ri,di); 1013 delta=square_modulus(delta,context0); 1014 } 1015 if (test){ 1016 v[i]=r; // we can certify there is a root in r by Brouwer fixed thm 1017 vradius[i]=-1; 1018 if (nextconj){ 1019 v[i+1]=conj(r,context0); 1020 vradius[i+1]=vradius[i]; 1021 ++i; 1022 } 1023 break; 1024 } 1025 if (delta.type==_REAL){ 1026 if (real_interval * ptr=dynamic_cast<real_interval *>(delta._REALptr)){ 1027 mpfr_t tmp; mpfr_init(tmp); 1028 mpfi_get_right(tmp,ptr->infsup); 1029 delta=real_object(tmp); 1030 mpfr_clear(tmp); 1031 } 1032 } 1033 sumdr2 += delta; 1034 if (!is_greater(deltar*deltar,sumdr2,context0)){ 1035 CERR << "Unable to certify " << v[i] << '\n' ; 1036 return false; 1037 } 1038 if (N<P.size()-epsn){ 1039 // add 10 bits of precision or double it 1040 if (N<-epsn/2){ 1041 deuxN=deuxN*deuxN; 1042 N*=2; 1043 } 1044 else { 1045 deuxN=1024*deuxN; 1046 N+=10; 1047 } 1048 } 1049 r -= Pr/P1r; 1050 // change precision to N 1051 reim(r,rr,ri,context0); 1052 if (rr.type==_REAL){ 1053 if (real_interval * ptr=dynamic_cast<real_interval *>(rr._REALptr)) 1054 mpfi_set_prec(ptr->infsup,N); 1055 } 1056 if (ri.type==_REAL){ 1057 if (real_interval * ptr=dynamic_cast<real_interval *>(ri._REALptr)) 1058 mpfi_set_prec(ptr->infsup,N); 1059 } 1060 r=rr+cst_i*ri; 1061 } // end for k 1062 #else 1063 if (N>int(P.size())/4-epsn/2) 1064 N=int(P.size())/4-epsn/2; 1065 gen deuxN=pow(2,N,context0); 1066 for (int k=0;k<kmax;++k){ 1067 in_round2(r,deuxN,N); // round2(r,deuxN,context0); 1068 // check if root precision is ok 1069 // otherwise try to improve root precision with Newton method 1070 gen P1r=horner(P1,r,0,false); 1071 if (is_exactly_zero(P1r)){ 1072 delta=plus_inf; 1073 break; 1074 } 1075 gen Pr=horner(P,r,0,false); 1076 delta=square_modulus(Pr,context0)/square_modulus(P1r,context0); 1077 if (is_greater(epsg2surdeg2,delta,context0)){ 1078 v[i]=r; // we can certify there is a root at distance <= eps from r 1079 if (is_exactly_zero(Pr)) 1080 vradius[i]=0; 1081 else 1082 vradius[i]=n*sqrt(accurate_evalf(delta,100),context0); 1083 // problem with double underflow 1084 if (!is_exactly_zero(vradius[i])) 1085 vradius[i]=min(epsg,pow(plus_two,int(evalf_double(ln(vradius[i],context0),1,context0)._DOUBLE_val/std::log(2.))+1),context0); 1086 if (debug_infolevel) 1087 CERR << CLOCK() << " isolated " << r << " radius " << vradius[i] << '\n'; 1088 if (nextconj){ 1089 v[i+1]=conj(r,context0); 1090 vradius[i+1]=vradius[i]; 1091 ++i; 1092 } 1093 break; 1094 } 1095 sumdr2 += delta; 1096 if (!is_greater(deltar*deltar,sumdr2,context0)){ 1097 CERR << "Unable to certify " << v[i] << '\n' ; 1098 return false; 1099 } 1100 if (N<int(P.size())-epsn){ 1101 // add 10 bits of precision or double it 1102 if (N<-epsn){ 1103 deuxN=deuxN*deuxN; 1104 N*=2; 1105 } 1106 else { 1107 deuxN=1024*deuxN; 1108 N+=10; 1109 } 1110 } 1111 in_round2(Pr,deuxN,N); in_round2(P1r,deuxN,N); // round2(Pr,deuxN,context0); round2(P1r,deuxN,context0); 1112 r -= Pr/P1r; 1113 } // end for k 1114 #endif 1115 if (!is_greater(epsg*epsg,delta,context0)) 1116 return false; 1117 return true; 1118 } 1119 1120 // find roots of polynomial P at precision eps using proot or 1121 // complex Sturm sequences 1122 // P must have numeric coefficients, in Q[i] complex_roots(const modpoly & P,const gen & a0,const gen & b0,const gen & a1,const gen & b1,bool complexe,double eps,bool use_proot)1123 vecteur complex_roots(const modpoly & P,const gen & a0,const gen & b0,const gen & a1,const gen & b1,bool complexe,double eps,bool use_proot){ 1124 if (P.empty()) 1125 return P; 1126 eps=absdouble(eps); 1127 if (eps>1e-6) 1128 eps=1e-6; 1129 if (use_proot){ 1130 int epsn=int(std::log(eps)/std::log(2.)-.5); 1131 gen epsg=pow(2,epsn,context0); 1132 gen epsg2surdeg2=(epsg*epsg)/int((P.size()+1)*(P.size()+1)); 1133 // first try proot with double precision, if it's not enough increase 1134 int n=45; 1135 bool Preal=is_zero(im(P,context0)); 1136 modpoly P1=derivative(P); 1137 for (;n<400;n*=2){ 1138 double cureps=std::pow(2.0,-n); 1139 if (debug_infolevel) 1140 CERR << CLOCK() << " proot at precision " << cureps << '\n'; 1141 vecteur v=proot(P,cureps,n); 1142 if (debug_infolevel) 1143 CERR << CLOCK() << " proot end at precision " << cureps << '\n'; 1144 vecteur vradius(v.size()); 1145 unsigned i=0; 1146 int kmax=int(std::log(eps)/std::log(cureps))+4; 1147 for (;i<v.size();++i){ 1148 newton_improve(P,P1,Preal,v,vradius,i,kmax,n,epsn,epsg2surdeg2,epsg); 1149 } // end for i 1150 if (i==v.size()){ 1151 vecteur res; 1152 for (unsigned j=0;j<v.size();++j){ 1153 if (Preal && is_zero(im(v[j],context0))){ 1154 if (is_exactly_zero(vradius[j]) || vradius[j]==-1){ 1155 if (is_greater(v[j],a0,context0) && is_greater(a1,v[j],context0) && is_greater(0,b0,context0) && is_greater(b1,0,context0)) 1156 res.push_back(makevecteur(v[j],1)); 1157 continue; 1158 } 1159 gen P1=horner(P,v[j]-vradius[j],0,false),P2=horner(P,v[j]+vradius[j],0,false); 1160 if (P1.type==_FRAC) P1=P1._FRACptr->num; 1161 if (P2.type==_FRAC) P2=P2._FRACptr->num; 1162 if (is_strictly_positive(-P1*P2,context0)){ 1163 if (is_greater(v[j],a0,context0) && is_greater(a1,v[j],context0) && is_greater(0,b0,context0) && is_greater(b1,0,context0)) 1164 res.push_back(makevecteur(eval(change_subtype(makevecteur(v[j]-vradius[j],v[j]+vradius[j]),_INTERVAL__VECT),1,context0),1)); 1165 continue; 1166 } 1167 } 1168 gen R,I; 1169 reim(v[j],R,I,context0); 1170 if (is_greater(R,a0,context0) && is_greater(a1,R,context0) && is_greater(I,b0,context0) && is_greater(b1,I,context0)){ 1171 if (is_exactly_zero(vradius[j])) 1172 res.push_back(makevecteur(v[j],1)); 1173 else { 1174 #ifdef HAVE_LIBMPFI 1175 gen a,b; 1176 reim(v[j],a,b,context0); 1177 res.push_back(makevecteur(eval(change_subtype(makevecteur(a-vradius[j],a+vradius[j]),_INTERVAL__VECT),1,context0)+cst_i*eval(change_subtype(makevecteur(b-vradius[j],b+vradius[j]),_INTERVAL__VECT),1,context0),1)); 1178 #else 1179 res.push_back(makevecteur(makevecteur(ratnormal(v[j]-vradius[j]*(1+cst_i)),ratnormal(v[j]+vradius[j]*(1+cst_i))),1)); 1180 #endif 1181 } 1182 } 1183 } 1184 return res; 1185 } // end if i==v.size() 1186 } // end for n 1187 CERR << "proot isolation did not work, trying complex Sturm sequences" << '\n'; 1188 } 1189 bool aplati=(a0==a1) && (b0==b1); 1190 if (!aplati && complexe && (a0==a1 || b0==b1) ) 1191 return vecteur(1,gensizeerr(gettext("Square is flat!"))); 1192 gen A0(a0),B0(b0),A1(a1),B1(b1); 1193 { 1194 // initial rectangle: |roots|< 1+ max(|a_i|)/|a_n| 1195 gen maxai=_max(*apply(P,abs,context0)._VECTptr,context0); 1196 gen tmp=1+maxai/abs(P.front(),context0); 1197 if (aplati){ 1198 A0=-tmp; 1199 B0=-tmp; 1200 A1=tmp; 1201 B1=tmp; 1202 } 1203 if (is_inf(A0)) A0=-tmp; 1204 if (is_inf(B0)) B0=-tmp; 1205 if (is_inf(A1)) A1=tmp; 1206 if (is_inf(B1)) B1=tmp; 1207 } 1208 gen tmp; 1209 modpoly p(*apply(P,exact,context0)._VECTptr); 1210 lcmdeno(p,tmp,context0); 1211 polynome pp(poly12polynome(p)); 1212 if (!complexe){ 1213 gen tmp=gcd(re(pp,context0),im(pp,context0)); 1214 if (tmp.type!=_POLY) 1215 return vecteur(0); 1216 pp=*tmp._POLYptr; 1217 } 1218 vecteur croots=crationalroot(pp,complexe); 1219 vecteur roots=keep_in_rectangle(croots,A0,B0,A1,B1,true,context0); 1220 p=polynome2poly1(pp); 1221 gen an=p.front(); 1222 if (!is_zero(im(an,context0))) 1223 p=conj(p.front(),context0)*p; 1224 if (!complexe){ // real root isolation 1225 modpoly R=p; 1226 modpoly S=derivative(R); 1227 vecteur listquo,coeffP,coeffR; 1228 csturm_seq(S,R,listquo,coeffP,coeffR,context0); 1229 // sparse polynomial patch 1230 if (pp.coord.size()<p.size()/10.){ 1231 R=vecteur(1,poly12polynome(R,1,1)); 1232 S=vecteur(1,poly12polynome(S,1,1)); 1233 } 1234 csturm_realroots(S,R,listquo,coeffP,coeffR,0,1,A0,A1,roots,eps,context0); 1235 return roots; 1236 } 1237 csturm_normalize(p,A0,B0,A1,B1,roots); 1238 gen pgcd; 1239 if (!complex_roots(p,A0,B0,A1,B1,pgcd,roots,eps)) 1240 return vecteur(1,gensizeerr(context0)); 1241 return roots; 1242 } 1243 1244 complexroot(const gen & g,bool complexe,GIAC_CONTEXT)1245 gen complexroot(const gen & g,bool complexe,GIAC_CONTEXT){ 1246 vecteur v=gen2vecteur(g); 1247 bool use_vas=!complexe ,use_proot=true; 1248 #ifndef HAVE_LIBMPFR 1249 use_proot=false; 1250 #endif 1251 bool isolation=false; 1252 if (!v.empty() && v[0]==at_sturm){ 1253 use_vas=false; 1254 use_proot=false; 1255 v.erase(v.begin()); 1256 } 1257 if (v.empty()) 1258 return gensizeerr(contextptr); 1259 if (v.size()<2){ 1260 isolation=true; 1261 v.push_back(epsilon(contextptr)); 1262 } 1263 if (v.size()==3) 1264 v.insert(v.begin()+1,epsilon(contextptr)); 1265 gen p=v.front(),prec=evalf_double(v[1],1,contextptr); 1266 if (prec.type!=_DOUBLE_) 1267 return gentypeerr(contextptr); 1268 double eps=prec._DOUBLE_val; 1269 if (eps>=1.0) 1270 eps=std::pow(10.,-eps); 1271 if (eps<=0) 1272 eps=epsilon(contextptr); 1273 if (eps<=0) 1274 eps=1e-12; 1275 if (v[0].type==_VECT && has_num_coeff(v[0])){ 1276 v=proot(*v[0]._VECTptr,eps); 1277 vecteur w; 1278 for (unsigned i=0;i<v.size();++i){ 1279 if (is_zero(im(v[i],contextptr))) 1280 w.push_back(makevecteur(v[i],1)); 1281 } 1282 return w; 1283 } 1284 unsigned vs=unsigned(v.size()); 1285 gen A(0),B(0),a0(minus_inf),b0(minus_inf),a1(plus_inf),b1(plus_inf); 1286 if (vs>3){ 1287 A=v[2]; 1288 B=v[3]; 1289 a0=re(A,contextptr); b0=im(A,contextptr); 1290 a1=re(B,contextptr);b1=im(B,contextptr); 1291 } 1292 if (is_greater(a0,a1,contextptr)) 1293 swapgen(a0,a1); 1294 if (is_greater(b0,b1,contextptr)) 1295 swapgen(b0,b1); 1296 vecteur vas_res; 1297 if (p.type==_VECT){ 1298 if (use_vas && vas(*p._VECTptr,a0,a1,isolation?1e300:eps,vas_res,true,contextptr)) 1299 return vas_res; 1300 return complex_roots(*p._VECTptr,a0,b0,a1,b1,complexe,eps,use_proot); 1301 } 1302 if (use_vas && vas(symb2poly_num(v[0],contextptr),a0,a1,isolation?1e300:eps,vas_res,true,contextptr)) 1303 return vas_res; 1304 vecteur l,l0; 1305 lidnt(p,l0,false); 1306 if (l0.size()!=1) 1307 return gentypeerr(contextptr); 1308 l=alg_lvar(p); 1309 gen px=_e2r(makesequence(p,l),contextptr); 1310 if (px.type==_FRAC) 1311 px=px._FRACptr->num; 1312 if (px.type!=_POLY) 1313 return vecteur(0); 1314 factorization f(sqff(*px._POLYptr)); 1315 factorization::const_iterator it=f.begin(),itend=f.end(); 1316 vecteur res; 1317 for (;it!=itend;++it){ 1318 gen P=_poly2symb(makesequence(it->fact,l),contextptr); 1319 P=_e2r(makesequence(P,l0.front()),contextptr); 1320 if (P.type!=_VECT) 1321 continue; 1322 vecteur tmp=complex_roots(*P._VECTptr,a0,b0,a1,b1,complexe,eps,use_proot); 1323 if (is_undef(tmp)) 1324 return tmp; 1325 iterateur jt=tmp.begin(),jtend=tmp.end(); 1326 for (;jt!=jtend;++jt){ 1327 if (jt->type==_VECT && jt->_VECTptr->size()==2) 1328 jt->_VECTptr->back()=it->mult*jt->_VECTptr->back(); 1329 } 1330 res=mergevecteur(res,tmp); 1331 } 1332 return res; 1333 } 1334 _complexroot(const gen & g,GIAC_CONTEXT)1335 gen _complexroot(const gen & g,GIAC_CONTEXT){ 1336 if ( g.type==_STRNG && g.subtype==-1) return g; 1337 gen res=complexroot(g,true,contextptr); 1338 if (res.type==_VECT) 1339 gen_sort_f_context(res._VECTptr->begin(),res._VECTptr->end(),complex_sort,contextptr); 1340 return res; 1341 // return _sorta(complexroot(g,true,contextptr),contextptr); 1342 } 1343 static const char _complexroot_s []="complexroot"; 1344 static define_unary_function_eval (__complexroot,&_complexroot,_complexroot_s); 1345 define_unary_function_ptr5( at_complexroot ,alias_at_complexroot,&__complexroot,0,true); 1346 _realroot(const gen & g,GIAC_CONTEXT)1347 gen _realroot(const gen & g,GIAC_CONTEXT){ 1348 if ( g.type==_STRNG && g.subtype==-1) return g; 1349 gen res; bool evalf_after=false; 1350 if (g.type==_VECT && !g._VECTptr->empty() && g._VECTptr->back()==at_evalf){ 1351 res=complexroot(gen(vecteur(g._VECTptr->begin(),g._VECTptr->end()-1),g.subtype),false,contextptr); 1352 evalf_after=true; 1353 } 1354 else 1355 res=complexroot(g,false,contextptr); 1356 if (res.type!=_VECT) 1357 return res; 1358 vecteur v=*res._VECTptr; 1359 for (unsigned i=0;i<v.size();++i){ 1360 if (v[i].type==_VECT && v[i]._VECTptr->size()==2){ 1361 gen a=v[i]._VECTptr->front(),b=v[i]._VECTptr->back(); 1362 if (a.type==_VECT && a.subtype==_INTERVAL__VECT){ 1363 if (evalf_after) 1364 v[i]=evalf((a._VECTptr->front()+a._VECTptr->back())/2,1,contextptr); 1365 else { 1366 a=eval(a,1,contextptr); 1367 v[i]=makevecteur(a,b); 1368 } 1369 } 1370 else { 1371 if (evalf_after) 1372 v[i]=evalf(a,1,contextptr); 1373 } 1374 } 1375 } 1376 return v; 1377 } 1378 static const char _realroot_s []="realroot"; 1379 static define_unary_function_eval (__realroot,&_realroot,_realroot_s); 1380 define_unary_function_ptr5( at_realroot ,alias_at_realroot,&__realroot,0,true); 1381 crationalroot(const gen & g0,bool complexe)1382 static vecteur crationalroot(const gen & g0,bool complexe){ 1383 gen g(g0),a,b; 1384 if (g.type==_VECT){ 1385 if (g.subtype==_SEQ__VECT){ 1386 vecteur & tmp=*g._VECTptr; 1387 if (tmp.size()!=3) 1388 return vecteur(1,gendimerr(context0)); 1389 g=tmp[0]; 1390 a=tmp[1]; 1391 b=tmp[2]; 1392 } 1393 else { 1394 g=poly12polynome(*g._VECTptr); 1395 } 1396 } 1397 gen a0,b0,a1,b1; 1398 ab2a0b0a1b1(a,b,a0,b0,a1,b1,context0); 1399 vecteur l; 1400 lvar(g,l); 1401 if (l.empty()) 1402 return vecteur(0); 1403 if (l.size()!=1) 1404 return vecteur(1,gentypeerr(context0)); 1405 gen px=_e2r(makevecteur(g,l),context0); 1406 if (px.type==_FRAC) 1407 px=px._FRACptr->num; 1408 if (px.type!=_POLY) 1409 return vecteur(0); 1410 factorization f(sqff(*px._POLYptr)); 1411 factorization::const_iterator it=f.begin(),itend=f.end(); 1412 vecteur res; 1413 for (;it!=itend;++it){ 1414 polynome p=it->fact; 1415 vecteur tmp=crationalroot(p,complexe); 1416 res=mergevecteur(res,tmp); 1417 } 1418 if (a0!=a1 || b0!=b1) 1419 res=keep_in_rectangle(res,a0,b0,a1,b1,false,context0); 1420 return res; 1421 } _crationalroot(const gen & g,GIAC_CONTEXT)1422 gen _crationalroot(const gen & g,GIAC_CONTEXT){ 1423 if ( g.type==_STRNG && g.subtype==-1) return g; 1424 return crationalroot(g,true); 1425 } 1426 static const char _crationalroot_s []="crationalroot"; 1427 static define_unary_function_eval (__crationalroot,&_crationalroot,_crationalroot_s); 1428 define_unary_function_ptr5( at_crationalroot ,alias_at_crationalroot,&__crationalroot,0,true); 1429 _rationalroot(const gen & g,GIAC_CONTEXT)1430 gen _rationalroot(const gen & g,GIAC_CONTEXT){ 1431 if ( g.type==_STRNG && g.subtype==-1) return g; 1432 return crationalroot(g,false); 1433 } 1434 static const char _rationalroot_s []="rationalroot"; 1435 static define_unary_function_eval (__rationalroot,&_rationalroot,_rationalroot_s); 1436 define_unary_function_ptr5( at_rationalroot ,alias_at_rationalroot,&__rationalroot,0,true); 1437 1438 // convert numerator of g to a list symb2poly_num(const gen & g_,GIAC_CONTEXT)1439 vecteur symb2poly_num(const gen & g_,GIAC_CONTEXT){ 1440 gen g(g_); 1441 if (g.type!=_VECT) 1442 g=makesequence(g,ggb_var(g)); 1443 gen tmp=_symb2poly(g,contextptr); 1444 if (tmp.type==_FRAC) 1445 tmp=tmp._FRACptr->num; 1446 if (tmp.type!=_VECT) 1447 return vecteur(1,gensizeerr(contextptr)); 1448 return *tmp._VECTptr; 1449 } 1450 // VAS implementation. Based on Xcas implementation by Alkiviadis G. Akritas, 1451 // A first C++ implementation was written by Spyros Kehagias and others 1452 // but it was too close to the Xcas code 1453 // This implementation is much faster, using basic data structures of giac 1454 // number of sign changes of the coefficients of P, returns -1 on error variations(const modpoly & P,GIAC_CONTEXT)1455 int variations(const modpoly & P,GIAC_CONTEXT){ 1456 int res=0,n=int(P.size()); 1457 if (!n) 1458 return -1; 1459 int previous=fastsign(P.front(),contextptr); 1460 if (previous==0) 1461 return -1; 1462 for (int i=1;i<n;i++){ 1463 if (is_exactly_zero(P[i])) 1464 continue; 1465 int current=fastsign(P[i],contextptr); 1466 if (!current) 1467 return -1; 1468 if (current!=previous){ 1469 ++res; 1470 previous=current; 1471 } 1472 } 1473 return res; 1474 } 1475 1476 #ifndef M_LN2 1477 #define M_LN2 0.6931471805599454 1478 #endif 1479 1480 // like (ln(n/d)+shift*ln2)/expo, but faster for large integers LMQ_evalf(const gen & n,const gen & d,double shift,int expo,GIAC_CONTEXT)1481 gen LMQ_evalf(const gen & n,const gen & d,double shift,int expo,GIAC_CONTEXT){ 1482 #ifndef USE_GMP_REPLACEMENTS 1483 if (is_integer(n) && is_integer(d)){ 1484 long int nexp=0,dexp=0; 1485 double nmant,dmant; 1486 if (n.type==_INT_) 1487 nmant=n.val; 1488 else 1489 nmant=mpz_get_d_2exp (&nexp,*n._ZINTptr); 1490 if (d.type==_INT_) 1491 dmant=d.val; 1492 else 1493 dmant=mpz_get_d_2exp (&dexp,*d._ZINTptr); 1494 return ( std::log(-nmant/dmant) + (nexp-dexp+shift)*M_LN2 )/expo; 1495 } 1496 #endif 1497 return ( ln(evalf(-n/d,1,contextptr),contextptr) + shift*M_LN2 )/gen(expo); 1498 } 1499 compute_lnabsmantexpo(const vecteur & cl,vector<double> & cllnabsmant,vector<long int> & clexpo,vector<short int> & clsign,GIAC_CONTEXT)1500 static bool compute_lnabsmantexpo(const vecteur & cl,vector<double> & cllnabsmant,vector<long int> & clexpo,vector<short int> & clsign,GIAC_CONTEXT){ 1501 int k=int(cl.size()); 1502 cllnabsmant.resize(k); 1503 clexpo.resize(k); 1504 clsign.resize(k); 1505 for (int i=0;i<k;++i){ 1506 gen tmp=sign(cl[i],contextptr); 1507 if (tmp.type!=_INT_) 1508 return false; 1509 clsign[i]=tmp.val; 1510 double mant; 1511 long int expo=0; 1512 if (is_integer(cl[i])){ 1513 if (cl[i].type==_ZINT){ 1514 #ifdef USE_GMP_REPLACEMENTS 1515 mant=evalf_double(cl[i],1,contextptr)._DOUBLE_val; 1516 #else 1517 mant=mpz_get_d_2exp (&expo,*cl[i]._ZINTptr); 1518 #endif 1519 } 1520 else mant=cl[i].val; 1521 } 1522 else { 1523 if (cl[i].type==_DOUBLE_){ 1524 mant=cl[i]._DOUBLE_val; 1525 } 1526 else { 1527 mant=evalf_double(cl[i],1,contextptr)._DOUBLE_val; 1528 } 1529 } 1530 mant=std::log(absdouble(mant)); 1531 cllnabsmant[i]=mant; 1532 clexpo[i]=expo; 1533 } 1534 return true; 1535 } 1536 1537 posubLMQ(const modpoly & P,GIAC_CONTEXT)1538 gen posubLMQ(const modpoly & P,GIAC_CONTEXT){ 1539 //---implements the Local_Max_Quadratic method (LMQ) to compute an 1540 //---upper bound on the values of the POSITIVE roots of p(x). 1541 1542 //---Reference:"Linear and Quadratic Complexity Bounds on the Values of the 1543 //---Positive Roots of Polynomials" by Alkiviadis G. Akritas. 1544 //---Journal of Universal Computer Science, Vol. 15, No. 3, 523-537, 2009. 1545 int k=int(P.size()); 1546 if (k<=1) 1547 return 0; 1548 vecteur cl; 1549 if (is_strictly_positive(P.front(),contextptr)) 1550 cl=P; 1551 else 1552 cl=-P; 1553 reverse(cl.begin(),cl.end()); 1554 vector<double> cllnabsmant; 1555 vector<long int> clexpo; 1556 vector<short int> clsign; 1557 if (!compute_lnabsmantexpo(cl,cllnabsmant,clexpo,clsign,contextptr)) 1558 return gensizeerr(contextptr); 1559 gen tempmax=minus_inf; 1560 vector<int> timesused(k+1,1); 1561 for (int m=k-1;m>=1;--m){ 1562 if (clsign[m-1]==-1){ // is_strictly_positive(-cl[m-1],contextptr) 1563 gen tempmin=plus_inf; 1564 for (int n=k;n>m;--n){ 1565 if (clsign[n-1]==1){ // is_strictly_positive(cl[n-1],contextptr) 1566 gen temp= (cllnabsmant[m-1]-cllnabsmant[n-1] + (clexpo[m-1]-clexpo[n-1]+timesused[n-1])*M_LN2)/(n-m);// LMQ_evalf(cl[m-1],cl[n-1],timesused[n-1],n-m,contextptr); 1567 // gen temp=pow(-cl[m-1]/cl[n-1]*pow(plus_two,timesused[n-1]),inv(n-m,contextptr),contextptr); 1568 // temp=evalf(temp,1,contextptr); 1569 ++timesused[n-1]; 1570 if (is_strictly_greater(tempmin,temp,contextptr)) 1571 tempmin=temp; 1572 } 1573 } 1574 if (is_strictly_greater(tempmin,tempmax,contextptr)) 1575 tempmax=tempmin; 1576 } 1577 } 1578 return _ceil(65*exp(tempmax,contextptr)/64,contextptr); // small margin 1579 } 1580 _posubLMQ(const gen & g,GIAC_CONTEXT)1581 gen _posubLMQ(const gen & g,GIAC_CONTEXT){ 1582 if ( g.type==_STRNG && g.subtype==-1) return g; 1583 vecteur v; 1584 if (g.type==_VECT && g.subtype!=_SEQ__VECT) 1585 v=*g._VECTptr; 1586 else 1587 v=symb2poly_num(g,contextptr); 1588 return posubLMQ(v,contextptr); 1589 } 1590 static const char _posubLMQ_s []="posubLMQ"; 1591 static define_unary_function_eval (__posubLMQ,&_posubLMQ,_posubLMQ_s); 1592 define_unary_function_ptr5( at_posubLMQ ,alias_at_posubLMQ,&__posubLMQ,0,true); 1593 poslbdLMQ(const modpoly & P,GIAC_CONTEXT)1594 gen poslbdLMQ(const modpoly & P,GIAC_CONTEXT){ 1595 //---implements the Local_Max_Quadratic method (LMQ) to compute a 1596 //---lower bound on the values of the POSITIVE roots of p(x). 1597 1598 //---Reference:"Linear and Quadratic Complexity Bounds on the Values of the 1599 //---Positive Roots of Polynomials" by Alkiviadis G. Akritas. 1600 //---Journal of Universal Computer Science, Vol. 15, No. 3, 523-537, 2009. 1601 int k=int(P.size()); 1602 if (k<=1) 1603 return 0; 1604 vecteur cl(P); 1605 reverse(cl.begin(),cl.end()); 1606 if (is_strictly_positive(-cl.front(),contextptr)) 1607 cl=-cl; 1608 vector<double> cllnabsmant; 1609 vector<long int> clexpo; 1610 vector<short int> clsign; 1611 if (!compute_lnabsmantexpo(cl,cllnabsmant,clexpo,clsign,contextptr)) 1612 return gensizeerr(contextptr); 1613 gen tempmax=minus_inf; 1614 vector<int> timesused(k,1); 1615 for (int m=1;m<k;++m){ 1616 if (clsign[m]==-1){ // is_strictly_positive(-cl[m],contextptr) 1617 gen tempmin=plus_inf; 1618 for (int n=0;n<m;++n){ 1619 if (clsign[n]==1){ // is_strictly_positive(cl[n],contextptr) 1620 // gen temp=pow(-cl[m]/cl[n]*pow(plus_two,timesused[n]),inv(m-n,contextptr),contextptr); 1621 // temp=evalf(temp,1,contextptr); 1622 gen temp= (cllnabsmant[m]-cllnabsmant[n] + (clexpo[m]-clexpo[n]+timesused[n])*M_LN2)/(m-n);// LMQ_evalf(cl[m],cl[n],timesused[n],m-n,contextptr); 1623 ++timesused[n]; 1624 if (is_strictly_greater(tempmin,temp,contextptr)) 1625 tempmin=temp; 1626 } 1627 } 1628 if (is_strictly_greater(tempmin,tempmax,contextptr)) 1629 tempmax=tempmin; 1630 } 1631 } 1632 tempmax=tempmax/M_LN2; 1633 tempmax=-_ceil(tempmax,contextptr); 1634 tempmax=pow(plus_two,tempmax,contextptr); 1635 return tempmax; 1636 } 1637 _poslbdLMQ(const gen & g,GIAC_CONTEXT)1638 gen _poslbdLMQ(const gen & g,GIAC_CONTEXT){ 1639 if ( g.type==_STRNG && g.subtype==-1) return g; 1640 vecteur v; 1641 if (g.type==_VECT && g.subtype!=_SEQ__VECT) 1642 v=*g._VECTptr; 1643 else 1644 v=symb2poly_num(g,contextptr); 1645 return poslbdLMQ(v,contextptr); 1646 } 1647 static const char _poslbdLMQ_s []="poslbdLMQ"; 1648 static define_unary_function_eval (__poslbdLMQ,&_poslbdLMQ,_poslbdLMQ_s); 1649 define_unary_function_ptr5( at_poslbdLMQ ,alias_at_poslbdLMQ,&__poslbdLMQ,0,true); 1650 makeinterval(const gen & a,const gen & b)1651 vecteur makeinterval(const gen & a,const gen & b){ 1652 if (is_strictly_greater(a,b,context0)) 1653 return makevecteur(b,a); 1654 else 1655 return makevecteur(a,b); 1656 } 1657 vas_sort(const gen & a,const gen & b)1658 bool vas_sort(const gen & a,const gen &b){ 1659 gen a1(a),b1(b); 1660 if (a.type==_VECT && a._VECTptr->size()==2) 1661 a1=a._VECTptr->front(); 1662 if (b.type==_VECT && b._VECTptr->size()==2) 1663 b1=b._VECTptr->front(); 1664 return is_strictly_greater(b1,a1,context0); 1665 } 1666 1667 // P is assumed to be squarefree and without rational roots 1668 // find roots of P((ax+b)/(cx+d)) VAS_positive_roots(const modpoly & P,const gen & ap,const gen & bp,const gen & cp,const gen & dp,GIAC_CONTEXT)1669 vecteur VAS_positive_roots(const modpoly & P,const gen & ap,const gen & bp,const gen & cp,const gen & dp,GIAC_CONTEXT){ 1670 //---The steps below correspond to the steps described in the reference below. 1671 1672 //---Reference: "A Comparative Study of Two Real Root Isolation Methods" 1673 //---by Alkiviadis G. Akritas and Adam W. Strzebonski. 1674 //---Nonlinear Analysis: Modelling and Control, Vol. 10, No. 4, 297-304, 2005. 1675 vecteur res; // root isolation intervals 1676 vecteur intervals_to_process; 1677 // STEP 1 1678 int v0=variations(P,contextptr); 1679 if (!v0) 1680 return res; 1681 gen ub=posubLMQ(P,contextptr); 1682 if (v0==1) 1683 return vecteur(1,makeinterval(0,ap*ub)); 1684 intervals_to_process.push_back(makevecteur(ap, bp, cp, dp, P,v0)); 1685 1686 // STEP 2 1687 while (!intervals_to_process.empty()){ 1688 gen tmp=intervals_to_process.back(); 1689 intervals_to_process.pop_back(); 1690 if (tmp.type!=_VECT || tmp._VECTptr->size()!=6) 1691 return vecteur(1,gensizeerr("VAS interval"+tmp.print())); 1692 vecteur & tmpv=*tmp._VECTptr; 1693 gen a=tmpv[0],b=tmpv[1],c=tmpv[2],d=tmpv[3], genf=tmpv[4],genv=tmpv[5]; 1694 if (genf.type!=_VECT || genv.type!=_INT_) 1695 return vecteur(1,gensizeerr("VAS interval"+tmp.print())); 1696 int v=genv.val; 1697 modpoly f = *genf._VECTptr; 1698 1699 // STEP 3 1700 gen lb=poslbdLMQ(f,contextptr); 1701 1702 // STEP 4 1703 if (is_strictly_greater(lb,16,contextptr)){ 1704 change_scale(f,lb); 1705 a=lb*a; c=lb*c; lb=1; 1706 } 1707 1708 // STEP 5 1709 if (is_greater(lb,1,contextptr)){ 1710 // f=taylor(f,lb); 1711 change_scale(f,lb); 1712 f=taylor(f,1); 1713 back_change_scale(f,lb); 1714 b = lb*a + b; d = lb*c + d; 1715 if (is_zero(f.back())){ 1716 res.push_back(b/d); 1717 f.pop_back(); 1718 } 1719 v=variations(f,contextptr); 1720 if (!v) 1721 continue; 1722 if (v==1){ 1723 if (!is_zero(c)) 1724 res.push_back(makeinterval(a/c,b/d)); 1725 else 1726 res.push_back(makeinterval(b,b+a*posubLMQ(f,contextptr))); 1727 continue; 1728 } 1729 } 1730 1731 // STEP 6 1732 modpoly f1=taylor(f,1),f2; 1733 gen a1=a, b1=a+b, c1=c, d1=c+d; 1734 int r=0; 1735 if (is_zero(f1.back())){ 1736 f1.pop_back(); 1737 res.push_back(b1/d1); 1738 r=1; 1739 } 1740 int v1=variations(f1,contextptr); 1741 int v2=v-v1-r; 1742 gen a2=b, b2=a+b, c2=d, d2=c+d; 1743 1744 // STEP 7 1745 if (v2>1){ 1746 f2=f; 1747 reverse(f2.begin(),f2.end()); 1748 f2=taylor(f2,1); 1749 if (is_zero(f2.back())) 1750 f2.pop_back(); 1751 v2=variations(f2,contextptr); 1752 } 1753 1754 // STEP 8 1755 if (v1<v2){ 1756 swapgen(a1,a2); 1757 swapgen(b1,b2); 1758 swapgen(c1,c2); 1759 swapgen(d1,d2); 1760 swap(f1,f2); 1761 int i=v1; v1=v2; v2=i; 1762 } 1763 1764 // STEP 9 1765 if (v1==0) continue; 1766 if (v1==1){ 1767 if (!is_zero(c1)) 1768 res.push_back(makeinterval(a1/c1,b1/d1)); 1769 else 1770 res.push_back(makeinterval(b1,b1 + a1*posubLMQ(f1,contextptr))); 1771 } 1772 else 1773 intervals_to_process.push_back(makevecteur(a1,b1,c1,d1,f1,v1)); 1774 1775 // STEP 10 1776 if (v2==0) continue; 1777 if (v2==1){ 1778 if (!is_zero(c2)) 1779 res.push_back(makeinterval(a2/c2,b2/d2)); 1780 else 1781 res.push_back(makeinterval(b2,b2 + a2*posubLMQ(f2,contextptr))); 1782 } 1783 else 1784 intervals_to_process.push_back(makevecteur(a2,b2,c2,d2,f2,v2)); 1785 } 1786 gen_sort_f(res.begin(),res.end(),vas_sort); 1787 return res; 1788 } 1789 _VAS_positive(const gen & g,GIAC_CONTEXT)1790 gen _VAS_positive(const gen & g,GIAC_CONTEXT){ 1791 if ( g.type==_STRNG && g.subtype==-1) return g; 1792 vecteur v; 1793 if (g.type==_VECT && g.subtype!=_SEQ__VECT) 1794 v=*g._VECTptr; 1795 else 1796 v=symb2poly_num(g,contextptr); 1797 return VAS_positive_roots(v,1,0,0,1,contextptr); 1798 } 1799 static const char _VAS_positive_s []="VAS_positive"; 1800 static define_unary_function_eval (__VAS_positive,&_VAS_positive,_VAS_positive_s); 1801 define_unary_function_ptr5( at_VAS_positive ,alias_at_VAS_positive,&__VAS_positive,0,true); 1802 1803 // square-free factorization of p, then remove all exponents 1804 // optionally remove factors with even multiplicities remove_multiplicities(const modpoly & p,factorization & f,bool odd_only,GIAC_CONTEXT)1805 modpoly remove_multiplicities(const modpoly & p,factorization & f,bool odd_only,GIAC_CONTEXT){ 1806 vecteur res(1,1),tmp; 1807 polynome P; 1808 poly12polynome(p,1,P,1); 1809 P=P/lgcd(P); 1810 f=sqff(P); 1811 factorization::const_iterator it=f.begin(),itend=f.end(); 1812 for (;it!=itend;++it){ 1813 if (odd_only && it->mult%2==0) 1814 continue; 1815 polynome2poly1(it->fact,1,tmp); 1816 res=operator_times(res,tmp,0); 1817 } 1818 return res; 1819 } 1820 vas(const modpoly & p,GIAC_CONTEXT)1821 gen vas(const modpoly & p,GIAC_CONTEXT){ 1822 vecteur v(p); 1823 vecteur res; 1824 bool has_zero=false; 1825 if (is_zero(v.back())){ 1826 has_zero=true; 1827 v.pop_back(); 1828 } 1829 res=VAS_positive_roots(v,1,0,0,1,contextptr); 1830 vecteur w(v); 1831 change_scale(w,-1); 1832 if (w.size()%2==0) 1833 w=-w; 1834 if (w==v){ 1835 v=-res; 1836 reverse(v.begin(),v.end()); 1837 iterateur it=v.begin(),itend=v.end(); 1838 for (;it!=itend;++it){ 1839 if (it->type==_VECT) 1840 reverse(it->_VECTptr->begin(),it->_VECTptr->end()); 1841 } 1842 if (has_zero) 1843 v.push_back(0); 1844 res=mergevecteur(v,res); 1845 } 1846 else { 1847 v=VAS_positive_roots(w,-1,0,0,1,contextptr); 1848 if (has_zero) 1849 v.push_back(0); 1850 res=mergevecteur(v,res); 1851 } 1852 return res; 1853 } _VAS(const gen & g,GIAC_CONTEXT)1854 gen _VAS(const gen & g,GIAC_CONTEXT){ 1855 if ( g.type==_STRNG && g.subtype==-1) return g; 1856 vecteur v; 1857 if (g.type==_VECT && g.subtype!=_SEQ__VECT) 1858 v=*g._VECTptr; 1859 else 1860 v=symb2poly_num(g,contextptr); 1861 factorization f; 1862 v=remove_multiplicities(v,f,false,contextptr); 1863 return vas(v,contextptr); 1864 } 1865 static const char _VAS_s []="VAS"; 1866 static define_unary_function_eval (__VAS,&_VAS,_VAS_s); 1867 define_unary_function_ptr5( at_VAS ,alias_at_VAS,&__VAS,0,true); 1868 1869 static const double bisection_newton_eps=1e-3; 1870 1871 // returns true if a root of p is found by Newton method, such that res-eps>a0 1872 // res+eps<b0 and sign of p changes between res-eps and res+eps bisection_newton(const modpoly & P,const modpoly & Pprime,const gen & a0,const gen & a1,gen & x0,double eps,gen & eps2,GIAC_CONTEXT)1873 static bool bisection_newton(const modpoly & P,const modpoly & Pprime,const gen & a0,const gen & a1,gen & x0,double eps,gen & eps2,GIAC_CONTEXT){ 1874 if (is_greater(bisection_newton_eps,x0-a0,contextptr) || is_greater(bisection_newton_eps,a1-x0,contextptr)) 1875 return false; // bisection is faster if the root is too close to the isolation interval boundaries 1876 int n=int(-std::log(eps)/std::log(2.0)+.5); // for rounding 1877 eps2=pow(2,-(n+1),contextptr); 1878 gen deuxn4(pow(2,n+4,contextptr)); 1879 in_round2(x0,deuxn4,n+4); // round2(x0,deuxn4,contextptr); 1880 for (int ii=0;ii<n;ii++){ 1881 gen Pprimex0=horner(Pprime,x0,0,false); 1882 if (is_zero(Pprimex0,contextptr)) 1883 return false; 1884 gen Px0=horner(P,x0,0,false); 1885 in_round2(Px0,deuxn4,n+4); // round2(Px0,deuxn4,contextptr); 1886 in_round2(Pprimex0,deuxn4,n+4); // round2(Pprimex0,deuxn4,contextptr); 1887 gen dx=Px0/Pprimex0; 1888 in_round2(dx,deuxn4,n+4); // round2(dx,deuxn4,contextptr); 1889 x0=x0-dx; 1890 if (is_positive(a0-x0,contextptr) || is_positive(x0-a1,contextptr)) 1891 return false; 1892 if (is_greater(abs(dx,contextptr),eps2,contextptr)) 1893 continue; 1894 if (is_positive(-horner(P,x0-eps2,0,false)*horner(P,x0+eps2,0,false),contextptr)) 1895 return true; 1896 } 1897 return false; 1898 1899 } 1900 bisection(const modpoly & p,const gen & a0,const gen & b0,double eps,GIAC_CONTEXT)1901 gen bisection(const modpoly & p,const gen & a0,const gen &b0,double eps,GIAC_CONTEXT){ 1902 int nsteps=int(std::ceil(std::log(evalf_double(b0-a0,1,contextptr)._DOUBLE_val/eps)/M_LN2)); 1903 int trynewtonstep=int(nsteps-std::log(bisection_newton_eps/eps)/M_LN2); 1904 gen a(a0),b(b0),m,eps2; 1905 modpoly dp=derivative(p); 1906 int s1=fastsign(horner(p,a,0,false),contextptr); 1907 if (s1==0) 1908 s1=fastsign(horner(dp,a,0,false),contextptr); 1909 int s2=fastsign(horner(p,b,0,false),contextptr); 1910 if (s2==0) 1911 s2=-fastsign(horner(dp,b,0,false),contextptr); 1912 if (s1==s2) 1913 return undef; 1914 for (int i=0;i<nsteps;++i){ 1915 m=(a+b)/2; 1916 if (i==trynewtonstep && bisection_newton(p,dp,a0,b0,m,eps,eps2,contextptr)) 1917 return makevecteur(m-eps2,m+eps2); 1918 int s=fastsign(horner(p,m,0,false),contextptr); 1919 if (s==0) 1920 return m; 1921 if (s==s1) 1922 a=m; 1923 else 1924 b=m; 1925 } 1926 return makevecteur(a,b); 1927 } 1928 multiplicity(const factorization & f,const gen & interval,GIAC_CONTEXT)1929 int multiplicity(const factorization & f,const gen & interval,GIAC_CONTEXT){ 1930 factorization::const_iterator it=f.begin(),itend=f.end(); 1931 if (interval.type==_VECT && interval._VECTptr->size()==2){ 1932 for (;it!=itend;++it){ 1933 if (is_strictly_positive(-it->fact(interval._VECTptr->front())*it->fact(interval._VECTptr->back()),contextptr)) 1934 return it->mult; 1935 } 1936 for (it=f.begin();it!=itend;++it){ 1937 if (is_positive(-it->fact(interval._VECTptr->front())*it->fact(interval._VECTptr->back()),contextptr)) 1938 return it->mult; 1939 } 1940 } 1941 else { 1942 for (;it!=itend;++it){ 1943 if (is_zero(it->fact(interval))) 1944 return it->mult; 1945 } 1946 } 1947 return 0; 1948 } 1949 add_vasres(vecteur & vasres,const gen & a,const gen & a0,const gen & b0,int mult,bool with_mult,GIAC_CONTEXT)1950 static void add_vasres(vecteur & vasres,const gen & a,const gen & a0,const gen & b0,int mult,bool with_mult,GIAC_CONTEXT){ 1951 if (a0==b0 || (is_greater(a,a0,contextptr) && is_greater(b0,a,contextptr)) ) 1952 vasres.push_back(with_mult?gen(makevecteur(a,mult)):a); 1953 } 1954 1955 // isolate and find real roots of P at precision eps between a and b 1956 // returns a list of intervals or of rationals vas(const modpoly & P,const gen & a0,const gen & b0,double eps,vecteur & vasres,bool with_mult,GIAC_CONTEXT)1957 bool vas(const modpoly & P,const gen & a0,const gen &b0,double eps,vecteur & vasres,bool with_mult,GIAC_CONTEXT){ 1958 if (P.size()<=3){ 1959 if (P.size()<2) 1960 return true; 1961 gen a(P[0]),b(P[1]); 1962 if (P.size()==2){ 1963 a=-b/a; 1964 add_vasres(vasres,a,a0,b0,1,with_mult,contextptr); 1965 return true; 1966 } 1967 gen c(P[2]); 1968 gen delta=b*b-4*a*c; 1969 if (is_zero(delta)){ 1970 a=-b/a/2; 1971 add_vasres(vasres,a,a0,b0,2,with_mult,contextptr); 1972 return true; 1973 } 1974 if (is_positive(delta,contextptr)){ 1975 delta=sqrt(delta,contextptr)/a/2; 1976 c=-b/a/2; 1977 add_vasres(vasres,c-delta,a0,b0,1,with_mult,contextptr); 1978 add_vasres(vasres,c+delta,a0,b0,1,with_mult,contextptr); 1979 } 1980 return true; 1981 } 1982 gen a(a0),b(b0); 1983 if (a==b){ 1984 a=minus_inf; 1985 b=plus_inf; 1986 } 1987 // check and convert coeffs of P 1988 modpoly p(P); 1989 iterateur it=p.begin(),itend=p.end(); 1990 for (;it!=itend;++it){ 1991 *it=exact(*it,contextptr); 1992 } 1993 gen tmp; 1994 lcmdeno(p,tmp,contextptr); 1995 for (it=p.begin();it!=itend;++it){ 1996 if (!is_integer(*it)) 1997 return false; 1998 } 1999 p=divvecteur(p,lgcd(p)); 2000 factorization f; 2001 p=remove_multiplicities(p,f,false,contextptr); 2002 tmp=vas(p,contextptr); 2003 if (tmp.type!=_VECT) 2004 return false; 2005 vecteur v=*tmp._VECTptr; 2006 // now improve precision by bisection 2007 it=v.begin(),itend=v.end(); 2008 for (;it!=itend;++it){ 2009 if (it->type!=_VECT){ 2010 if (is_greater(*it,a,contextptr) && is_greater(b,*it,contextptr)){ 2011 if (with_mult){ 2012 int n=multiplicity(f,*it,contextptr); 2013 vasres.push_back(makevecteur(*it,n)); 2014 } 2015 else 2016 vasres.push_back(*it); 2017 } 2018 continue; 2019 } 2020 if (it->_VECTptr->size()!=2) 2021 return false; 2022 gen A=it->_VECTptr->front(),B=it->_VECTptr->back(); 2023 if (is_strictly_greater(a,B,contextptr) || is_strictly_greater(A,b,contextptr)) 2024 continue; 2025 gen interval=bisection(p,max(A,a,contextptr),min(B,b,contextptr),eps,contextptr); 2026 if (is_undef(interval)) 2027 continue; 2028 if (interval.type==_VECT) 2029 interval.subtype=_INTERVAL__VECT; 2030 if (with_mult) 2031 vasres.push_back(makevecteur(interval,multiplicity(f,interval,contextptr))); 2032 else 2033 vasres.push_back( (interval.type==_VECT && interval._VECTptr->size()==2)?evalf((interval._VECTptr->front()+interval._VECTptr->back())/2,1,contextptr):interval); 2034 } 2035 return true; 2036 } 2037 2038 #ifndef NO_NAMESPACE_GIAC 2039 } // namespace giac 2040 #endif // ndef NO_NAMESPACE_GIAC 2041