1 // -*- mode:C++ ; compile-command: "g++ -I.. -I../include -DHAVE_CONFIG_H -DIN_GIAC -DGIAC_GENERIC_CONSTANTS -fno-strict-aliasing -g -c alg_ext.cc -Wall" -*- 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 21 using namespace std; 22 #include <cmath> 23 #include <cstdlib> 24 #include <stdexcept> 25 #include <errno.h> 26 #include <map> 27 #include "gen.h" 28 #include "gausspol.h" 29 #include "identificateur.h" 30 #include "poly.h" 31 #include "usual.h" 32 #include "sym2poly.h" 33 #include "vecteur.h" 34 #include "modpoly.h" 35 #include "alg_ext.h" 36 #include "vecteur.h" 37 #include "solve.h" 38 #include "subst.h" 39 #include "plot.h" 40 #include "derive.h" 41 #include "ezgcd.h" 42 #include "prog.h" 43 #include "intg.h" 44 #include "csturm.h" 45 #include "lin.h" 46 #include "ti89.h" 47 #include "giacintl.h" 48 49 #ifndef NO_NAMESPACE_GIAC 50 namespace giac { 51 #endif // ndef NO_NAMESPACE_GIAC 52 islesscomplex(const gen & a,const gen & b)53 bool islesscomplex(const gen & a,const gen & b){ 54 if (a==b) 55 return false; 56 return a.islesscomplexthan(b); 57 } 58 59 // symbolic_rootof_list() protected with a mutex in multi-thread environment operator ()(const gen & a,const gen & b) const60 bool comparegen::operator ()(const gen & a,const gen & b) const { 61 if (a.type==_INT_ && b.type==_INT_) 62 return a.val<b.val; 63 gen A1,A2,B1,B2; 64 if (a.type==_VECT && a._VECTptr->size()==2 && (A1=a._VECTptr->front()).type==_INT_ && (A2=a._VECTptr->back()).type==_INT_ && b.type==_VECT && b._VECTptr->size()==2 && (B1=b._VECTptr->front()).type==_INT_ && (B2=b._VECTptr->back()).type==_INT_){ 65 return (A1.val!=B1.val)?A1.val<B1.val:A2.val<B2.val; 66 } 67 return a.islesscomplexthan(b); 68 } 69 symbolic_rootof_list()70 rootmap & symbolic_rootof_list(){ 71 static rootmap * ans= 0; 72 if (!ans) ans=new rootmap; 73 return *ans; 74 } proot_list()75 rootmap & proot_list(){ 76 static rootmap * ans= 0; 77 if (!ans) ans=new rootmap; 78 return *ans; 79 } 80 galoisconj_list()81 rootmap & galoisconj_list(){ 82 static rootmap * ans= 0; 83 if (!ans) ans=new rootmap; 84 return *ans; 85 } 86 87 #ifdef HAVE_LIBPTHREAD 88 static pthread_mutex_t rootof_mutex = PTHREAD_MUTEX_INITIALIZER; rootof_trylock()89 static int rootof_trylock(){ 90 return pthread_mutex_trylock(&rootof_mutex); 91 } rootof_unlock()92 static void rootof_unlock(){ 93 pthread_mutex_unlock(&rootof_mutex); 94 } 95 #else rootof_trylock()96 static int rootof_trylock(){ return 0; } rootof_unlock()97 static void rootof_unlock(){ } 98 99 #endif 100 // get Galois conjugates in the same number field from cache galoisconj_cached(const vecteur & v,vecteur & res)101 bool galoisconj_cached(const vecteur & v,vecteur & res){ 102 if (rootof_trylock()) 103 return false; 104 res.clear(); 105 rootmap::iterator ritend=galoisconj_list().end(),rit=galoisconj_list().find(v); 106 if (rit!=ritend && rit->second.type==_VECT) 107 res=*rit->second._VECTptr; 108 rootof_unlock(); 109 return !res.empty(); 110 } 111 112 // cache list of Galois conjugates galoisconj_cache(const vecteur & v,const vecteur & res)113 bool galoisconj_cache(const vecteur & v,const vecteur & res){ 114 if (rootof_trylock()) 115 return false; 116 rootmap::iterator ritend=galoisconj_list().end(),rit=galoisconj_list().find(v); 117 if (rit==ritend) 118 galoisconj_list()[v]=res; 119 rootof_unlock(); 120 return true; 121 } 122 galoisconj(const vecteur & v,GIAC_CONTEXT)123 vecteur galoisconj(const vecteur & v,GIAC_CONTEXT){ 124 vecteur res; 125 if (galoisconj_cached(v,res)) 126 return res; 127 gen g=symb_horner(v,vx_var); 128 if (pari_galoisconj(g,res,contextptr)) 129 return res; 130 if (int(v.size())>MAX_COMMON_ALG_EXT_ORDER_SIZE) return res; 131 // factor v over rootof(v) if degree is small 132 g=_factors(makesequence(g,rootof(g,contextptr)),contextptr); 133 if (g.type!=_VECT) return res; 134 vecteur w=*g._VECTptr; 135 for (int i=0;i<int(w.size())-1;i+=2){ 136 gen a,b; 137 if (is_linear_wrt(w[i],vx_var,a,b,contextptr) && !is_zero(a)){ 138 res.push_back(-b/a); 139 } 140 } 141 galoisconj_cache(v,res); 142 return res; 143 } 144 _galoisconj(const gen & args,GIAC_CONTEXT)145 gen _galoisconj(const gen & args,GIAC_CONTEXT){ 146 gen g=args; 147 if (g.type==_SYMB) 148 g=_symb2poly(args,contextptr); 149 if (g.type!=_VECT) return gensizeerr(contextptr); 150 return galoisconj(*g._VECTptr,contextptr); 151 } 152 static const char _galoisconj_s []="galoisconj"; 153 static define_unary_function_eval (__galoisconj,&_galoisconj,_galoisconj_s); 154 define_unary_function_ptr5( at_galoisconj ,alias_at_galoisconj,&__galoisconj,0,true); 155 156 // if true, g is a rootof such that conj(rootof(w))=g conj_in_nf(const vecteur & w,gen & g,GIAC_CONTEXT)157 bool conj_in_nf(const vecteur & w,gen & g,GIAC_CONTEXT){ 158 gen r1=rootof(w,contextptr); 159 vecteur c=galoisconj(w,contextptr); 160 gen pow10=pow(10,14,contextptr); 161 int maxdigits=1000; 162 if (c.size()<w.size()-1) 163 maxdigits=50; 164 #ifndef HAVE_LIBMPFR 165 maxdigits=14; 166 #endif 167 gen borne=100; 168 for (int ndigits=14;ndigits<=maxdigits;ndigits*=2){ 169 gen R1=conj(_evalf(makesequence(r1,ndigits),contextptr),contextptr); 170 for (int i=0;i<int(c.size());++i){ 171 gen r2=c[i]; 172 gen R2=_evalf(makesequence(r2,ndigits),contextptr); 173 if (is_greater(borne*abs(R1,contextptr),abs(R1-R2,contextptr)*pow10,contextptr)){ 174 g=r2; 175 return true; 176 } 177 } 178 pow10=pow10*pow10; 179 borne=borne*borne; 180 } 181 return false; 182 } 183 proot_cached(const vecteur & v,double eps,vecteur & res)184 bool proot_cached(const vecteur & v,double eps,vecteur & res){ 185 if (rootof_trylock()) 186 return false; 187 res.clear(); 188 double oldeps=1e300; 189 rootmap::iterator ritend=proot_list().end(),rit=proot_list().find(v); 190 if (rit!=ritend && rit->second.type==_VECT){ 191 res=*rit->second._VECTptr; 192 if (res.size()==2 && res.front().type==_VECT && res.back().type==_DOUBLE_){ 193 oldeps=res.back()._DOUBLE_val; 194 res=*res.front()._VECTptr; 195 } 196 else 197 res.clear(); 198 } 199 rootof_unlock(); 200 return !res.empty() && oldeps<=eps; 201 } 202 proot_cache(const vecteur & v,double eps,const vecteur & res)203 bool proot_cache(const vecteur & v,double eps,const vecteur & res){ 204 if (rootof_trylock()) 205 return false; 206 rootmap::iterator ritend=proot_list().end(),rit=proot_list().find(v); 207 if (rit!=ritend){ 208 if (rit->second.type!=_VECT || rit->second._VECTptr->size()!=2 || rit->second._VECTptr->front().type!=_VECT || rit->second._VECTptr->back().type!=_DOUBLE_ || rit->second._VECTptr->back()._DOUBLE_val>eps) 209 rit->second=makevecteur(res,eps); 210 } 211 else 212 proot_list()[v]=makevecteur(res,eps); 213 rootof_unlock(); 214 return true; 215 } 216 algebraic_EXTension(const gen & a_,const gen & v)217 gen algebraic_EXTension(const gen & a_,const gen & v){ 218 gen a(a_); 219 if (a.type==_VECT && !a._VECTptr->empty() && is_zero(a._VECTptr->front())){ 220 a=trim(*a._VECTptr,0); 221 } 222 if (is_zero(a) ) 223 return 0; 224 if (a.type==_VECT){ 225 if (a._VECTptr->empty()) 226 return zero; 227 if (a._VECTptr->size()==1) 228 return a._VECTptr->front(); 229 // a.subtype=_POLY1__VECT; 230 } 231 gen res; 232 #ifdef SMARTPTR64 233 * ((ulonglong * ) &res) = ulonglong(new ref_algext) << 16; 234 #else 235 res.__EXTptr=new ref_algext; 236 #endif 237 res.type=_EXT; 238 *(res._EXTptr+1) = v; 239 // if (v.type==_VECT) (res._EXTptr+1)->subtype=_POLY1__VECT; 240 if (a.type==_FRAC){ 241 *res._EXTptr = a._FRACptr->num; 242 return fraction(res,a._FRACptr->den); 243 } 244 *res._EXTptr = a; 245 return res; 246 } 247 in_select_root(const vecteur & a,bool reel,GIAC_CONTEXT,double eps)248 gen in_select_root(const vecteur & a,bool reel,GIAC_CONTEXT,double eps){ 249 if (a.empty() || is_undef(a)) 250 return undef; 251 gen current(a.front()); 252 double max_re(re(current,contextptr).evalf_double(1,contextptr)._DOUBLE_val),max_im(im(current,contextptr).evalf_double(1,contextptr)._DOUBLE_val); 253 const_iterateur it=a.begin(),itend=a.end(); 254 for (;it!=itend;++it){ 255 double cur_re(re(*it,contextptr).evalf_double(1,contextptr)._DOUBLE_val),cur_im(im(*it,contextptr).evalf_double(1,contextptr)._DOUBLE_val); 256 if (cur_re > (1+eps)*max_re ){ 257 current=*it; 258 max_re=cur_re; 259 max_im=cur_im; 260 } 261 else { // same argument 262 if ( absdouble(cur_re-max_re)<eps*max_re && (cur_im>max_im) ){ 263 current=*it; 264 max_im=cur_im; 265 } 266 } 267 } 268 if (reel && is_strictly_positive(-im(current,contextptr),contextptr)) 269 current=conj(current,contextptr); 270 return current; 271 } 272 select_root(const vecteur & v,GIAC_CONTEXT)273 gen select_root(const vecteur & v,GIAC_CONTEXT){ 274 int n=decimal_digits(contextptr); 275 if (n<12) n=12; 276 if (n>307) n=307; 277 double eps=std::pow(0.1,n); 278 int rprec=int(n*3.3); 279 vecteur a=proot(v,eps,rprec); 280 gen r=in_select_root(a,is_real(v,contextptr),contextptr); 281 return r; 282 } 283 alg_evalf(const gen & a,const gen & b,GIAC_CONTEXT)284 gen alg_evalf(const gen & a,const gen &b,GIAC_CONTEXT){ 285 if (a.type==_FRAC) 286 return rdiv(alg_evalf(a._FRACptr->num,b,contextptr),alg_evalf(a._FRACptr->den,b,contextptr),contextptr); 287 gen a1=a.evalf(1,contextptr),b1=b.evalf(1,contextptr); 288 if (a1.type!=_VECT) 289 return a1; 290 if (b1.type!=_VECT) 291 return algebraic_EXTension(a1,b1); 292 gen r(select_root(*b1._VECTptr,contextptr)); 293 if (is_undef(r)) 294 return algebraic_EXTension(a1,b1); 295 return horner(*a1._VECTptr,r); 296 } 297 ext_reduce(const gen & a,const gen & v)298 gen ext_reduce(const gen & a, const gen & v){ 299 if (a.type==_FRAC) 300 return fraction(ext_reduce(a._FRACptr->num,v),ext_reduce(a._FRACptr->den,v)); 301 if (a.type!=_VECT) 302 return a;// algebraic_EXTension(a,v); 303 if (a._VECTptr->empty()) 304 return zero; 305 if (a._VECTptr->size()==1) 306 return a._VECTptr->front(); 307 if (v.type==_VECT){ 308 if (a._VECTptr->size()<v._VECTptr->size()) 309 return algebraic_EXTension(a,v); 310 #if 1 311 // special code for quadratic extension, if v=[1,0,-a] 312 if (v._VECTptr->size()==3 && v[0]==1 && v[1]==0){ 313 gen x=-v[2],r1,r0; 314 if (a._VECTptr->size()==3){ 315 r0=a._VECTptr->front()*x+a._VECTptr->back(); 316 r1=(*a._VECTptr)[1]; 317 } 318 else { 319 const_iterateur it=a._VECTptr->begin(),itend=a._VECTptr->end()-1; 320 for (;it<itend;){ 321 r0 = r0*x+(*it); 322 ++it; 323 r1 = r1*x+(*it); 324 ++it; 325 } 326 if (it==itend) 327 r0 = r0*x+(*it); 328 else swapgen(r0,r1); 329 } 330 if (r1==0) 331 return r0; 332 gen c=new ref_vecteur(2); 333 c._VECTptr->front()=r1; 334 c._VECTptr->back()=r0; 335 return algebraic_EXTension(c,v); 336 } 337 gen c=new ref_vecteur; 338 vecteur & rem=*c._VECTptr; 339 modpoly quo; 340 environment env; 341 DivRem(*a._VECTptr,*v._VECTptr,0,quo,rem); 342 if (rem.empty()) return 0; 343 if (rem.size()==1) return rem.front(); 344 return algebraic_EXTension(c,v); 345 #endif 346 return algebraic_EXTension((*a._VECTptr) % (*v._VECTptr),v); 347 } 348 if (v.type==_FRAC) 349 return horner(*a._VECTptr,*v._FRACptr,true); 350 if (v.type!=_EXT) 351 return gentypeerr(gettext("ext_reduce")); 352 gen va=*v._EXTptr,vb=*(v._EXTptr+1); 353 if (va.type==_FRAC) 354 return ext_reduce(horner(*a._VECTptr,*va._FRACptr,true),vb); 355 if (va.type!=_VECT){ 356 if (vb.type!=_VECT) 357 return gensizeerr(gettext("alg_ext.cc/ext_reduce")); 358 return algebraic_EXTension( (*a._VECTptr) % (*vb._VECTptr),v); 359 } 360 return ext_reduce(horner(*a._VECTptr,gen(*va._VECTptr,_POLY1__VECT)),vb); 361 } 362 ext_reduce(const gen & e)363 gen ext_reduce(const gen & e){ 364 #ifdef DEBUG_SUPPORT 365 if (e.type!=_EXT){ 366 gensizeerr(gettext("alg_ext.cc/ext_reduce")); 367 CERR << gettext("alg_ext.cc/ext_reduce"); 368 return e; 369 } 370 #endif 371 if ( (e._EXTptr->type==_VECT) && ((e._EXTptr+1)->type==_VECT) && 372 (e._EXTptr->_VECTptr->size()<(e._EXTptr+1)->_VECTptr->size()) ) 373 return e; 374 return ext_reduce(*(e._EXTptr),*(e._EXTptr+1)); 375 } 376 polynome2vecteur(const polynome & p,int na,int nb,vecteur & v)377 static bool polynome2vecteur(const polynome & p,int na,int nb,vecteur & v){ 378 v=vecteur(na*nb,zero); 379 int i,j; 380 if (p.dim!=2){ 381 #ifdef NO_STDEXCEPT 382 return false; 383 #else 384 setsizeerr(gettext("alg_ext.cc/polynome2vecteur")); 385 return false; 386 #endif 387 } 388 vector< monomial<gen> >::const_iterator it=p.coord.begin(),itend=p.coord.end(); 389 for (;it!=itend;++it){ 390 i=it->index.front(); 391 j=it->index.back(); 392 // cerr << nb*(na-i-1)+nb-j-1 << " " << na*nb << '\n'; 393 v[nb*(na-i-1)+nb-j-1]=it->value; 394 } 395 return true; 396 } 397 is_known_rootof(const vecteur & v,gen & symroot,GIAC_CONTEXT)398 bool is_known_rootof(const vecteur & v,gen & symroot,GIAC_CONTEXT){ 399 const_iterateur it=v.begin(),itend=v.end(); 400 for (;it!=itend;++it){ 401 if (it->type!=_INT_) 402 return false; 403 } 404 if (rootof_trylock()) 405 return false; 406 rootmap::iterator ritend=symbolic_rootof_list().end(),rit=symbolic_rootof_list().find(v); 407 if (rit!=ritend) 408 symroot=rit->second; 409 rootof_unlock(); 410 if (rit!=ritend) 411 return true; 412 if (v.size()==3){ 413 vecteur w; 414 identificateur x(" x"); 415 in_solve(symb_horner(v,x),x,w,0,contextptr); 416 if (w.empty()) 417 return false; 418 symroot=w.front(); 419 return true; 420 } 421 return false; 422 } 423 424 // replace _EXT == to ext by g in v replace_ext(const vecteur & v,const vecteur & ext,const gen & g,GIAC_CONTEXT)425 static vecteur replace_ext(const vecteur & v,const vecteur &ext,const gen & g,GIAC_CONTEXT){ 426 vecteur res; 427 const_iterateur it=v.begin(),itend=v.end(); 428 res.reserve(int(itend-it)); 429 for (;it!=itend;++it){ 430 gen numtmp=*it,dentmp=1; 431 if (it->type==_FRAC){ 432 numtmp=it->_FRACptr->num; 433 dentmp=it->_FRACptr->den; 434 } 435 // if numtmp is an ext, it must be the same ext as a 436 if (numtmp.type==_EXT){ 437 if (*(numtmp._EXTptr+1)!=ext) 438 return vecteur(1,gensizeerr(gettext("Invalid _EXT in replace_ext"))); 439 res.push_back(horner(*numtmp._EXTptr,g)/dentmp); 440 } 441 else 442 res.push_back(evalf_double(*it,1,contextptr)); 443 } 444 return res; 445 } 446 447 // given theta1 and theta2 with minimal poly va and vb (inside gen ga and gb) 448 // find k / Q[theta1+ k*theta2 ] contains theta1 and theta2 449 // return the minimal poly of theta=theta1+k*theta2 450 // and return in a and b theta1 and theta2 as ext (in terms of theta) common_minimal_POLY(const gen & ga,const gen & gb,gen & a,gen & b,int & k,GIAC_CONTEXT)451 gen common_minimal_POLY(const gen & ga,const gen & gb, gen & a,gen & b,int & k,GIAC_CONTEXT){ 452 const vecteur & va=*ga._VECTptr; 453 const vecteur & vb=*gb._VECTptr; 454 int na=int(va.size()-1),nb=int(vb.size()-1); 455 if (nb==1){ 456 k=0; 457 vecteur un(2,zero); 458 un[0]=plus_one; 459 gen vag(va); 460 a=algebraic_EXTension(un,vag); 461 gen tmp=-vb[1]; 462 if (tmp.type!=_POLY) 463 b=tmp; 464 else { 465 if (tmp._POLYptr->coord.empty()) 466 b=zero; 467 else 468 b=tmp._POLYptr->coord.front().value; 469 } 470 return vag; 471 } 472 // create minimal polynomial of theta1/theta2 as 2-d polynomials 473 // with main variable respectively a and b 474 // (since pb is used for reduction after var reordering of p) 475 polynome pa(2),pb(2); 476 const_iterateur it=va.begin(),itend=va.end(); 477 for (int d=na;it!=itend;++it,--d){ 478 if (!is_zero(*it)) 479 pa.coord.push_back(monomial<gen>(*it,d,1,2)); // deg=d, var=1, dim=2 480 } 481 it=vb.begin(),itend=vb.end(); 482 int k_init=0; 483 for (int d=nb;it!=itend;++it,--d){ 484 if (!is_zero(*it)){ 485 gen numtmp=*it,dentmp=1; 486 if (it->type==_FRAC){ 487 numtmp=it->_FRACptr->num; 488 dentmp=it->_FRACptr->den; 489 } 490 polynome pbadd(pb.dim); 491 // if numtmp is an ext, it must be the same ext as a 492 if (numtmp.type==_EXT){ 493 pbadd=poly12polynome(*(numtmp._EXTptr->_VECTptr),1,1).untrunc1(d); 494 k_init=1; 495 } 496 else 497 pbadd.coord.push_back(monomial<gen>(numtmp,d,1,2)); 498 pb = pb + pbadd/dentmp; 499 } 500 } 501 if (k_init){ 502 vecteur v1=*evalf_double(va,1,contextptr)._VECTptr; 503 if (is_fully_numeric(v1)){ 504 // when theta2 depends on theta1, theta1+k*theta2 is not necessarily 505 // the largest root, because the numeric value of v2 depends 506 // on the selected root of v1 507 // 508 // we should compute k*theta1+theta2 for a sufficiently large 509 // value of k to insure largest root, e.g. 510 // this implies computing approx value of theta1 and theta2 511 // 512 vecteur rac=real_proot(v1,1e-12,contextptr); 513 if (rac.empty()){ 514 vecteur rac1=proot(v1,1e-12); 515 gen theta1=in_select_root(rac1,is_real(v1,contextptr),contextptr); 516 // replace _EXT in vb by r1 and evaluate numerically 517 vecteur v2=replace_ext(vb,va,theta1,contextptr); 518 if (!v2.empty() && is_undef(v2)) 519 return v2.front(); 520 // find theta2 521 if (is_fully_numeric(v2)){ 522 vecteur rac2=proot(v2,1e-12); 523 if (!rac2.empty() && !is_undef(rac2)){ 524 gen theta2=in_select_root(rac2,is_real(v2,contextptr),contextptr); 525 int racs=int(rac1.size()); 526 for (int i=0;i<racs;++i){ 527 gen r1=rac1[i],K; 528 if (r1==theta1) 529 continue; 530 v2=replace_ext(vb,va,r1,contextptr); 531 if (!v2.empty() && is_undef(v2)) 532 return v2.front(); 533 #ifndef NO_STDEXCEPT 534 try { 535 #endif 536 gen r2=select_root(v2,contextptr); 537 K=(r2-theta2)/(theta1-r1); // must be <= k 538 if (is_undef(K)) 539 K=0; 540 #ifndef NO_STDEXCEPT 541 } 542 catch (std::runtime_error & ){ 543 last_evaled_argptr(contextptr)=NULL; 544 K=0; 545 } 546 #endif 547 // so that r2-theta2 <= k*theta1-k*r1 548 // or k*r1+r2 <= k*theta1+theta2 549 K=_floor(re(K,contextptr),contextptr)+1; 550 if (is_positive(K,contextptr) && K.type!=_INT_) 551 return gensizeerr(gettext("Unable to find common minimal polynomial")); 552 k_init=std::max(k_init,K.val); 553 } 554 } // !rac2.empty 555 } // is_fully_numeric(v2) 556 } 557 if (!rac.empty() && !is_undef(rac)){ 558 gen theta1=_max(rac,contextptr); 559 // replace _EXT in vb by r1 and evaluate numerically 560 vecteur v2=replace_ext(vb,va,theta1,contextptr); 561 if (!v2.empty() && is_undef(v2)) 562 return v2.front(); 563 // find largest root (i.e. theta2) 564 if (is_fully_numeric(v2)){ 565 vecteur rac2=real_proot(v2,1e-12,contextptr); 566 if (!rac2.empty() && !is_undef(rac2)){ 567 gen theta2=_max(rac2,contextptr); 568 int racs=int(rac.size()); 569 for (int i=0;i<racs;++i){ 570 gen r1=rac[i],K; 571 if (r1==theta1) 572 continue; 573 v2=replace_ext(vb,va,r1,contextptr); 574 if (!v2.empty() && is_undef(v2)) 575 return v2.front(); 576 #ifndef NO_STDEXCEPT 577 try { 578 #endif 579 gen r2=_max(real_proot(v2,1e-12,contextptr),contextptr); 580 K=(r2-theta2)/(theta1-r1); // must be <= k 581 if (is_undef(K)) 582 K=0; 583 #ifndef NO_STDEXCEPT 584 } 585 catch (std::runtime_error & ){ 586 last_evaled_argptr(contextptr)=NULL; 587 K=0; 588 } 589 #endif 590 // so that r2-theta2 <= k*theta1-k*r1 591 // or k*r1+r2 <= k*theta1+theta2 592 K=_floor(K,0)+1; 593 if (is_positive(K,contextptr) && K.type!=_INT_) 594 return gensizeerr(gettext("Unable to find common minimal polynomial")); 595 k_init=std::max(k_init,K.val); 596 } 597 } // !rac2.empty 598 } 599 } // if !rac.empty() 600 } // fully_numeric 601 } 602 else 603 ++k_init; // start with k=1 if theta2 does not depend on theta1 604 matrice m; 605 m.reserve(na*nb); 606 for (k=k_init;;++k){ 607 polynome p(2); 608 p.coord.push_back(monomial<gen>(1,2)); 609 polynome q(2),tmpq(2),tmpr(2); 610 q.coord.push_back(monomial<gen>(k,1,1,2)); // k*a: deg=1, var=1, dim=2 611 q.coord.push_back(monomial<gen>(1,1,2,2)); // b: deg=1, var=2 612 // create the matrix 613 // lines are 1, k*a+b, ..., (k*a+b)^(na*nb) 614 // in terms of (columns) 615 // a^(na-1)*b^(nb-1) ... a^(na-1) ... ab^(nb-1) ... ab a b^(nb-1) ... b 1 616 m.clear(); 617 vecteur ligne; 618 for (int j=0;j<=na*nb;++j){ 619 if (!polynome2vecteur(p,na,nb,ligne)) 620 return gensizeerr(gettext("alg_ext.cc/polynome2vecteur")); 621 // ligne.push_back(pow(theta,j)); 622 m.push_back(ligne); 623 p=p*q; 624 // permutation of indices order before making division by pb 625 p.reorder(transposition(0,1,2)); 626 p.TDivRem(pb,tmpq,tmpr,true); p.coord.swap(tmpr.coord); // p=p%pb; // 627 p.reorder(transposition(0,1,2)); 628 // division by a after because b might depend on a 629 p.TDivRem(pa,tmpq,tmpr,true); p.coord.swap(tmpr.coord); // p=p%pa; // 630 } 631 // Add the lines corresponding to b and a (i.e. theta2, theta1) 632 ligne=vecteur(na*nb); 633 ligne[na*nb-2]=plus_one; 634 m.push_back(ligne); 635 ligne=vecteur(na*nb); 636 ligne[na*nb-nb-1]=plus_one; 637 m.push_back(ligne); 638 // Transpose matrix 639 // then we have the na*nb+3 columns 1, theta, ..., theta^(na*nb), b, a 640 // in terms of a basis (with na*nb coordinates) 641 m=mtran(m); 642 // reduce the matrix m to echelon form and test rank=na*nb 643 // if ok break, else try another value of k 644 matrice m_red; 645 vecteur pivots; 646 gen det; 647 int st=step_infolevel(contextptr); 648 step_infolevel(contextptr)=0; 649 if (!mrref(m,m_red,pivots,det,0,na*nb,0,na*nb+3, 650 /* fullreduction */1,0,true,1,0, 651 contextptr)){ 652 step_infolevel(contextptr)=st; 653 return gensizeerr(contextptr); 654 } 655 step_infolevel(contextptr)=st; 656 m=m_red; 657 // the reduced matrix m should have the form 658 // * 0 ... 0 * * * 659 // 0 * ... 0 * * * 660 // ... 661 // 0 0 ... 0 * * * 662 // 0 0 ... 0 * * * 663 // 0 0 ... ? * * * 664 // with ? != 0, we check ?, if it is zero we try another value k 665 vecteur v(m[na*nb-1]._VECTptr->begin(),m[na*nb-1]._VECTptr->end()-1); 666 if (!is_zero__VECT(v,contextptr)) 667 break; 668 } 669 mdividebypivot(m); 670 // add a -1 at the end of column na*nb (C convention, index starting at 0) 671 // to get the min poly 672 vecteur v(na*nb+1); 673 for (int i=0;i<na*nb;++i) 674 v[i]=-m[i][na*nb]; 675 v[na*nb]=plus_one; 676 reverse(v.begin(),v.end()); 677 // remove denominators 678 gen e; 679 lcmdeno_converted(v,e,contextptr); 680 // use column na*nb+1 to find b=theta2 in terms of theta 681 vecteur w(na*nb); 682 for (int i=0;i<na*nb;++i) 683 w[i]=m[i][na*nb+1]; 684 reverse(w.begin(),w.end()); 685 w=trim(w,0); 686 lcmdeno_converted(w,e,contextptr); 687 b=fraction(w,e); 688 // to get a=theta1 we use column na*nb+2 689 w=vecteur(na*nb); 690 for (int i=0;i<na*nb;++i) 691 w[i]=m[i][na*nb+2]; 692 reverse(w.begin(),w.end()); 693 w=trim(w,0); 694 lcmdeno_converted(w,e,contextptr); 695 a=fraction(w,e); 696 // convert to algebraic extensions 697 gen vg(v); 698 b=algebraic_EXTension(b,vg); 699 a=algebraic_EXTension(a,vg); 700 // add v to the rootof_list 701 if (!rootof_trylock()){ 702 rootmap::iterator ritend=symbolic_rootof_list().end(),rit=symbolic_rootof_list().find(v); 703 rootof_unlock(); 704 if (rit==ritend){ 705 // should first check that va/vb are solvable poly 706 gen gaa,gbb; 707 if (is_known_rootof(va,gaa,contextptr) && is_known_rootof(vb,gbb,contextptr)){ 708 if (!rootof_trylock()){ 709 symbolic_rootof_list()[v]=gaa +k*gbb; 710 rootof_unlock(); 711 } 712 } 713 } 714 } 715 return vg; 716 } 717 718 // assuming a is the extptr+1 of an ext, return the min pol of 719 // theta generating the algebraic extension min_pol(gen & a)720 vecteur min_pol(gen & a){ 721 if (a.type==_VECT) 722 return *a._VECTptr; 723 else { 724 if ( (a.type!=_EXT) || ((a._EXTptr+1)->type!=_VECT) ) 725 return vecteur(1,gensizeerr(gettext("alg_ext.cc/min_pol"))); 726 return *((a._EXTptr+1)->_VECTptr); 727 } 728 } 729 730 // Find an evaluation point for p at b where pb=p[b] is squarefree find_good_eval(const polynome & F,polynome & Fb,vecteur & b)731 bool find_good_eval(const polynome & F,polynome & Fb,vecteur & b){ 732 int Fdeg=F.lexsorted_degree(),nvars=int(b.size()); 733 gen Fg; 734 int essai=0; 735 for (;;++essai){ 736 Fb=peval_1(F,b,0); 737 if (Fb.lexsorted_degree()==Fdeg && gcd(Fb,Fb.derivative()).lexsorted_degree()==0 ){ 738 return true; 739 } 740 b=vranm(nvars,0,0); // find another random point 741 } 742 } 743 744 static void clean(gen & g); clean(polynome & p)745 static void clean(polynome & p){ 746 vector< monomial<gen> >::iterator it=p.coord.begin(),itend=p.coord.end(); 747 for (;it!=itend;++it) 748 clean(it->value); 749 } 750 clean(gen & g)751 static void clean(gen & g){ 752 if (g.is_symb_of_sommet(at_neg) && is_integer(g._SYMBptr->feuille)) 753 g=-g._SYMBptr->feuille; 754 if (g.type==_POLY){ 755 clean(*g._POLYptr); 756 return; 757 } 758 if (g.type==_VECT){ 759 iterateur it=g._VECTptr->begin(),itend=g._VECTptr->end(); 760 for (;it!=itend;++it) 761 clean(*it); 762 return; 763 } 764 if (g.type==_EXT){ 765 clean(*g._EXTptr); 766 clean(*(g._EXTptr+1)); 767 } 768 } 769 770 // in-place reduction of algebraic extensions clean_ext_reduce(vecteur & v)771 void clean_ext_reduce(vecteur & v){ 772 iterateur it=v.begin(),itend=v.end(); 773 for (;it!=itend;++it) 774 clean_ext_reduce(*it); 775 } clean_ext_reduce(gen & g)776 void clean_ext_reduce(gen & g){ 777 if (g.type==_EXT){ 778 g=ext_reduce(g); 779 return; 780 } 781 if (g.type==_VECT){ 782 clean_ext_reduce(*g._VECTptr); 783 return; 784 } 785 if (g.type==_POLY){ 786 vector< monomial<gen> >::iterator it=g._POLYptr->coord.begin(),itend=g._POLYptr->coord.end(); 787 for (;it!=itend;++it) 788 clean_ext_reduce(it->value); 789 return; 790 } 791 if (g.type==_FRAC) 792 clean_ext_reduce(g._FRACptr->num); 793 } 794 795 // a and b are supposed to be *(_EXTptr+1) of some algebraic extension 796 // common_EXT will return a new algebraic extension 797 // (suitable to be an extptr+1) 798 // and will modify a and b to be ext of the returned common_EXT common_EXT(gen & a,gen & b,const vecteur * l,GIAC_CONTEXT)799 gen common_EXT(gen & a,gen & b,const vecteur * l,GIAC_CONTEXT){ 800 if (a==b) 801 return a; 802 if (a.type==_FRAC) 803 return common_EXT(a._FRACptr->num,b,l,contextptr); 804 if (b.type==_FRAC) 805 return common_EXT(a,b._FRACptr->num,l,contextptr); 806 // extract minimal polynomials 807 gen a_orig(a),b_orig(b); 808 gen a__VECT,b__VECT; 809 if (a.type==_VECT) 810 a__VECT=a; 811 else { 812 if ( (a.type!=_EXT) || ((a._EXTptr+1)->type!=_VECT) ) 813 return gensizeerr(gettext("alg_ext.cc/common_EXT")); 814 a__VECT=*(a._EXTptr+1); 815 } 816 if (b.type==_VECT) 817 b__VECT=b; 818 else { 819 if ( (b.type!=_EXT) || ((b._EXTptr+1)->type!=_VECT) ) 820 return gensizeerr(gettext("alg_ext.cc/common_EXT")); 821 b__VECT=*(b._EXTptr+1); 822 } 823 int as=int(a__VECT._VECTptr->size()),bs=int(b__VECT._VECTptr->size()); 824 if (bs>as) 825 return common_EXT(b,a,l,contextptr); 826 if (as==3 && bs==3 && is_one(a[0]) && is_one(b[0]) && is_zero(a[1]) && is_zero(b[1]) && a[2]==-b[2]){ // sqrt(X) and sqrt(-X) 827 b=algebraic_EXTension(makevecteur(cst_i,0),a); 828 gen tmp=a; 829 a=algebraic_EXTension(makevecteur(1,0),a); 830 return tmp; 831 } 832 // special handling if both extensions are cyclotomic 833 int ac=is_cyclotomic(*a__VECT._VECTptr,epsilon(contextptr)),bc; 834 if (ac && (bc=is_cyclotomic(*b__VECT._VECTptr,epsilon(contextptr))) ){ 835 int cc=ac*bc/gcd(ac,bc); 836 gen res=gen(cyclotomic(cc),_POLY1__VECT); 837 a__VECT=gen(vecteur(cc/ac+1),_POLY1__VECT); 838 a__VECT._VECTptr->front()=1; 839 a=algebraic_EXTension(a__VECT,res); 840 b__VECT=gen(vecteur(cc/bc+1),_POLY1__VECT); 841 b__VECT._VECTptr->front()=1; 842 b=algebraic_EXTension(b__VECT,res); 843 return res; 844 } 845 // reduce extension degree by factorizing b__VECT over Q[a] 846 polynome p(poly12polynome(*b__VECT._VECTptr)); 847 polynome p_content(p.dim); 848 factorization f; 849 gen an,extra_div; 850 ext_factor(p,algebraic_EXTension(a__VECT,a__VECT),an,p_content,f,false,extra_div); 851 // now choose in the factorization which factor is relevant for b 852 // this is done by approximation if possible 853 // or by choosing the factor of lowest degree 854 // this way we update b__VECT 855 int min_deg=int(b__VECT._VECTptr->size()); 856 factorization::const_iterator f_it=f.begin(),f_itend=f.end(); 857 bool trouve=false; 858 if (f_itend-f_it==1) 859 trouve=true; 860 vecteur racines; 861 vector<double> real_racines; 862 int innerdim=0; 863 const_iterateur b_it=b__VECT._VECTptr->begin(),b_itend=b__VECT._VECTptr->end(); 864 for (;b_it!=b_itend;++b_it){ 865 if (b_it->type==_POLY) 866 innerdim=b_it->_POLYptr->dim; 867 } 868 vecteur vb(innerdim); 869 gen racine_max=undef; 870 bool deep_emb=false; // marker for deep embedding 871 if (!trouve){ 872 // Change for multivariate polynomials p, added evaluation 873 if (innerdim){ 874 gen params; 875 polynome pb(1),px(unsplitmultivarpoly(p,innerdim)); 876 *logptr(contextptr) << gettext("Warning, need to choose a branch for the root of a polynomial with parameters. This might be wrong.") << '\n'; 877 if (l && l->size()>=2){ 878 for (int i=1;i<l->size();++i){ 879 params=(*l)[i]; 880 if (params.type==_VECT && !params._VECTptr->empty()) 881 break; 882 } 883 // IMPROVE: using context and *l look for assumptions 884 if (params.type==_VECT){ 885 vecteur paramv=*params._VECTptr; 886 int nessais=3*paramv.size(); 887 for (int essai=0;essai<nessais;++essai){ 888 for (unsigned j=0;j<paramv.size() && j<vb.size();++j){ 889 gen p=paramv[j]; 890 if (p.type!=_IDNT) 891 continue; 892 if (p==cst_pi){ 893 vb[j]=p; 894 continue; 895 } 896 gen g,g2=p._IDNTptr->eval(1,g,contextptr); 897 if ((g2.type==_VECT) && (g2.subtype==_ASSUME__VECT)){ 898 vecteur V=*g2._VECTptr; 899 if (V.size()==2) 900 vb[j]=V[1]; 901 if ( V.size()==3 && V[1].type==_VECT && V[2].type==_VECT){ 902 for (unsigned i=0;i<V[1]._VECTptr->size();++i){ 903 gen tmp=(*V[1]._VECTptr)[i]; 904 if (tmp.type==_VECT && tmp._VECTptr->size()==2){ 905 gen a=tmp._VECTptr->front(),b=tmp._VECTptr->back(); 906 int decal=1; 907 if (1 || essai) 908 decal += int((giac_rand(contextptr)*100.0)/rand_max2); 909 if (a==minus_inf) 910 vb[j]=b-decal; 911 else { 912 if (b==plus_inf) 913 vb[j]=a+decal; 914 else { 915 if (a+b==0) 916 vb[j]=(decal%2?a:b)/(decal+1); 917 else 918 vb[j]=(decal*a+b)/(decal+1); 919 } 920 } 921 } 922 } 923 } // end if V.size()==3 924 } // end g2 assume_vect 925 } // end for j 926 vecteur vb0=vb; 927 find_good_eval(px,pb,vb); 928 if (vb0==vb) 929 break; 930 } // end trying to find a good eval point satisfying assumptions 931 } // end params.type==_VECT 932 } 933 vecteur vb0=vb; 934 find_good_eval(px,pb,vb); // find_good_eval does not take care of assumptions, but vb should be ok (loop above) 935 if (vb==vb0) 936 *logptr(contextptr) << gettext("The choice was done assuming ") << params << "=" << vb << '\n'; 937 else 938 *logptr(contextptr) << gettext("Non regular value ") << vb0 << gettext(" was discarded and replaced randomly by ") << params << "=" << vb << '\n'; 939 // checking for embedded polynomial coefficients 940 vector< monomial<gen> >::const_iterator it=pb.coord.begin(),itend=pb.coord.end(); 941 for (;0 && it!=itend;++it){ // disabled, computations would be too complex 942 if (it->value.type==_POLY){ 943 deep_emb=true; 944 break; 945 } 946 } 947 if (!deep_emb) 948 racines=proot(gen2vecteur(evalf(polynome2poly1(pb),1,contextptr))); 949 } 950 else 951 racines=proot(*evalf(b__VECT,1,contextptr)._VECTptr); // evalf to avoid recursion if computing exact roots of b__VECT 952 if (is_undef(racines)) return gensizeerr(contextptr); 953 // racines= list of approx roots if b__VECT is numeric 954 // empty if not numeric 955 racine_max=in_select_root(racines,is_real(b__VECT,contextptr),contextptr); 956 } // if (!trouve) 957 if (!deep_emb && !trouve && !is_undef(racine_max)){ // select root for b 958 // now eval each factor over racine_max and choose the one with 959 // minimal absolute value 960 double min_abs=0,racine_max_d=evalf_double(abs(racine_max,contextptr),1,contextptr)._DOUBLE_val; 961 int ndig=14; 962 while (!trouve){ 963 for (;f_it!=f_itend;++f_it){ 964 vecteur vtmp(polynome2poly1(f_it->fact)); 965 gen tmp; 966 lcmdeno_converted(vtmp,tmp,contextptr); 967 int maxsave=max_sum_sqrt(contextptr); 968 max_sum_sqrt(0,contextptr); 969 if (innerdim) 970 tmp=r2sym(vtmp,vecteur(1,vb),contextptr); 971 else 972 tmp=r2sym(vtmp,vecteur(1,vecteur(0)),contextptr); 973 max_sum_sqrt(maxsave,contextptr); 974 #if defined HAVE_LIBMPFR && defined HAVE_LIBPARI // change for Martin Deraux big extensions 975 if (ndig<15) 976 tmp=evalf(tmp,1,contextptr); 977 else 978 tmp=_evalf(makesequence(tmp,ndig),contextptr); 979 #else 980 tmp=evalf(tmp,1,contextptr); 981 #endif 982 if (tmp.type==_VECT && !tmp._VECTptr->empty()) 983 tmp=tmp/tmp._VECTptr->front(); 984 gen f_racine_max(evalf_double(abs(horner(tmp,racine_max),contextptr),1,contextptr)); 985 if (f_racine_max.type!=_DOUBLE_) 986 continue; 987 double current_evaluation=fabs(f_racine_max._DOUBLE_val); 988 if (!trouve){ 989 trouve=true; 990 min_abs=current_evaluation; 991 p=f_it->fact; 992 } 993 else { 994 if (min_abs>current_evaluation){ 995 min_abs=current_evaluation; 996 p=f_it->fact; 997 } 998 } 999 } // end for on f_it 1000 if (min_abs>1e-4*racine_max_d){ 1001 *logptr(contextptr) << "Precision problem choosing root in common_EXT, current precision " << ndig << '\n'; 1002 trouve=false; 1003 ndig=2*ndig; 1004 f_it=f.begin(); 1005 #if defined HAVE_LIBMPFR && defined HAVE_LIBPARI 1006 if (ndig>1000) 1007 #endif 1008 break; 1009 } 1010 } // end while !trouve 1011 } // end racine_max defined 1012 if (!trouve) { 1013 for (;f_it!=f_itend;++f_it){ 1014 if ( (b.type==_EXT) && is_zero(horner(polynome2poly1(f_it->fact,1),*b._EXTptr)) ){ 1015 p=f_it->fact; 1016 break; 1017 } 1018 int d=f_it->fact.lexsorted_degree(); 1019 if (d && (d<=min_deg)){ 1020 p=f_it->fact; 1021 min_deg=d; 1022 } 1023 } 1024 } // end choose by degree 1025 clean(p); 1026 b__VECT=polynome2poly1(p/p.coord.front().value); // p must be monic (?) 1027 // compute new minimal polynomial 1028 int k; 1029 gen res1=common_minimal_POLY(a__VECT,b__VECT,a,b,k,contextptr); 1030 if ((a_orig.type==_EXT) && (b_orig.type==_EXT) && !is_undef(res1)) 1031 return algebraic_EXTension(a_orig+gen(k)*b_orig,res1); 1032 else 1033 return res1; 1034 } 1035 ext_add(const gen & aa,const gen & bb,GIAC_CONTEXT)1036 gen ext_add(const gen & aa,const gen & bb,GIAC_CONTEXT){ 1037 gen a(ext_reduce(aa)),b(ext_reduce(bb)); 1038 if ( (a.type!=_EXT) || (b.type!=_EXT) ) 1039 return a+b; 1040 if (*(a._EXTptr+1)==*(b._EXTptr+1)){ 1041 if ( (a._EXTptr->type==_VECT) && (b._EXTptr->type==_VECT)){ 1042 gen c=new ref_vecteur; 1043 addmodpoly(*a._EXTptr->_VECTptr,*b._EXTptr->_VECTptr,*c._VECTptr); 1044 return ext_reduce(c,*(a._EXTptr+1)); 1045 return ext_reduce(*(a._EXTptr->_VECTptr)+ *(b._EXTptr->_VECTptr),*(a._EXTptr+1)); 1046 } 1047 else 1048 return ext_reduce(*a._EXTptr+*b._EXTptr,*(a._EXTptr+1)); 1049 } 1050 gen c=common_EXT(*(a._EXTptr+1),*(b._EXTptr+1),0,contextptr); 1051 if (is_undef(c)) return c; 1052 // if c.type==_INT_/_ZINT, call ichinrem on a.extptr,b.extptr,... 1053 return ext_reduce(a)+ext_reduce(b); 1054 } 1055 ext_sub(const gen & a,const gen & b,GIAC_CONTEXT)1056 gen ext_sub(const gen & a,const gen & b,GIAC_CONTEXT){ 1057 if (*(a._EXTptr+1)==*(b._EXTptr+1)){ 1058 if ( (a._EXTptr->type==_VECT) && (b._EXTptr->type==_VECT)){ 1059 #if 1 1060 gen c=new ref_vecteur; 1061 submodpoly(*a._EXTptr->_VECTptr,*b._EXTptr->_VECTptr,*c._VECTptr); 1062 return ext_reduce(c,*(a._EXTptr+1)); 1063 #endif 1064 return ext_reduce(*(a._EXTptr->_VECTptr)- *(b._EXTptr->_VECTptr),*(a._EXTptr+1)); 1065 } 1066 else 1067 return ext_reduce(*a._EXTptr-*b._EXTptr,*(a._EXTptr+1)); 1068 } 1069 return ext_add(a,-b,contextptr); 1070 } 1071 ext_mul(const gen & aa,const gen & bb,GIAC_CONTEXT)1072 gen ext_mul(const gen & aa,const gen & bb,GIAC_CONTEXT){ 1073 gen a(ext_reduce(aa)),b(ext_reduce(bb)); 1074 if ( (a.type!=_EXT) || (b.type!=_EXT) ) 1075 return a*b; 1076 if (*(a._EXTptr+1)==*(b._EXTptr+1)){ 1077 if ((a._EXTptr->type==_VECT) && (b._EXTptr->type==_VECT)){ 1078 #if 1 1079 gen c=new ref_vecteur; 1080 operator_times(*a._EXTptr->_VECTptr,*b._EXTptr->_VECTptr,0,*c._VECTptr); 1081 return ext_reduce(c,*(a._EXTptr+1)); 1082 #endif 1083 return ext_reduce( *(a._EXTptr->_VECTptr) * *(b._EXTptr->_VECTptr),*(a._EXTptr+1)); 1084 } 1085 else 1086 return ext_reduce((*a._EXTptr)*(*b._EXTptr),*(a._EXTptr+1)); 1087 } 1088 gen c=common_EXT(*(a._EXTptr+1),*(b._EXTptr+1),0,contextptr); 1089 if (is_undef(c)) return c; 1090 // if c.type==_INT_/_ZINT, call ichinrem on a._EXTptr,b._EXTptr,... 1091 return ext_reduce(a)*ext_reduce(b); 1092 } 1093 inv_EXT(const gen & aa)1094 gen inv_EXT(const gen & aa){ 1095 if (aa.type!=_EXT) 1096 return inv(aa,context0); 1097 gen a(ext_reduce(aa)); 1098 if (a.type==_FRAC){ 1099 return a._FRACptr->den*inv_EXT(a._FRACptr->num); 1100 } 1101 if (a.type!=_EXT) 1102 return inv(a,context0); 1103 if (a._EXTptr->type==_VECT){ 1104 vecteur u,v,d; 1105 egcd(*(a._EXTptr->_VECTptr),*((a._EXTptr+1)->_VECTptr),0,u,v,d); 1106 if (d.size()!=1) 1107 return gensizeerr(gettext("inv_EXT")); 1108 gen de=d.front(),du=u; 1109 simplify(du,de); 1110 return fraction(algebraic_EXTension(du,*(a._EXTptr+1)),de); 1111 } 1112 return gentypeerr(gettext("inv_EXT")); 1113 } 1114 horner_rootof(const vecteur & p,const gen & g,GIAC_CONTEXT)1115 gen horner_rootof(const vecteur & p,const gen & g,GIAC_CONTEXT){ 1116 if (g.type==_SYMB && g._SYMBptr->feuille.type==_VECT && 1117 // false 1118 int(g._SYMBptr->feuille._VECTptr->size())>max_sum_sqrt(contextptr) 1119 ) 1120 return symb_horner(p,g); 1121 const_iterateur it=p.begin(),itend=p.end(); 1122 gen res; 1123 for (;it!=itend;++it){ 1124 res=ratnormal(res*g+*it,contextptr); 1125 } 1126 return ratnormal(res,contextptr); 1127 } 1128 has_rootof_value(const gen & Pmin,gen & value,GIAC_CONTEXT)1129 bool has_rootof_value(const gen & Pmin,gen & value,GIAC_CONTEXT){ 1130 value=undef; 1131 if (contextptr && contextptr->globalcontextptr->rootofs){ 1132 const vecteur & r=*contextptr->globalcontextptr->rootofs; 1133 for (unsigned i=0;i<r.size();++i){ 1134 gen ri=r[i]; 1135 if (ri.type==_VECT && ri._VECTptr->size()==2 && Pmin.type==_VECT && ri._VECTptr->front().type==_VECT && *Pmin._VECTptr==*ri._VECTptr->front()._VECTptr){ 1136 value=ri._VECTptr->back(); 1137 return true; 1138 } 1139 } 1140 } 1141 return !is_undef(value); 1142 } 1143 printasrootof(const gen & g,const char * s,GIAC_CONTEXT)1144 static string printasrootof(const gen & g,const char * s,GIAC_CONTEXT){ 1145 if (contextptr && g.type==_VECT && g._VECTptr->size()==2){ 1146 gen value; 1147 if (g._VECTptr->front().type==_VECT && has_rootof_value(g._VECTptr->back(),value,contextptr)){ 1148 value=horner_rootof(*g._VECTptr->front()._VECTptr,value,contextptr); 1149 string res=value.print(contextptr); 1150 if (need_parenthesis(value)) 1151 res=("("+res)+')'; 1152 return res; 1153 } 1154 } 1155 string res(s); 1156 res+='('; 1157 res+=g.print(contextptr); 1158 res+=')'; 1159 return res; 1160 } 1161 1162 // rootof has 2 args: P(theta) and Pmin(theta) symb_rootof(const gen & p,const gen & pmin,GIAC_CONTEXT)1163 gen symb_rootof(const gen & p,const gen &pmin,GIAC_CONTEXT){ 1164 if (p.type!=_VECT) 1165 return p; 1166 // first check that pmin is in the list of known rootof 1167 gen value(undef); 1168 if (!rootof_trylock()){ 1169 rootmap::iterator it=symbolic_rootof_list().find(pmin),itend=symbolic_rootof_list().end(); 1170 if (it!=itend) 1171 value=it->second; 1172 rootof_unlock(); 1173 } 1174 if (is_undef(value)) 1175 return symbolic(at_rootof,makevecteur(p,pmin)); 1176 return horner_rootof(*p._VECTptr,value,contextptr); 1177 // return ratnormal(ratnormal(symb_horner(*p._VECTptr,it->second))); 1178 } rootof(const gen & e,GIAC_CONTEXT)1179 gen rootof(const gen & e,GIAC_CONTEXT){ 1180 if (e.type!=_VECT){ 1181 vecteur v=lidnt(e); 1182 if (v.size()==1) 1183 return rootof(_symb2poly(makesequence(e,v.front()),contextptr),contextptr); 1184 return gentypeerr(gettext("rootof")); 1185 } 1186 if (e.type==_VECT && *e._VECTptr==makevecteur(1,0,1)){ 1187 *logptr(contextptr) << "rootof([1,0,1]) was converted to i" << '\n'; 1188 return cst_i; 1189 } 1190 if (e._VECTptr->size()==2 && e._VECTptr->front().type!=_VECT){ 1191 vecteur v=lidnt(e); 1192 if (v.size()!=1) 1193 return gentypeerr(gettext("rootof")); 1194 return rootof(makesequence(_symb2poly(makesequence(e._VECTptr->front(),v.front()),contextptr),_symb2poly(makesequence(e._VECTptr->back(),v.front()),contextptr)),contextptr); 1195 } 1196 if (e._VECTptr->size()!=2 || e._VECTptr->back().type!=_VECT) 1197 return rootof(makesequence(makevecteur(1,0),e),contextptr); 1198 if (has_num_coeff(e)) 1199 return approx_rootof(e,contextptr); 1200 if (!lop(lvar(e),at_pow).empty()){ 1201 *logptr(contextptr) << gettext("Algebraic extensions not allowed in a rootof")<<'\n'; 1202 return approx_rootof(e,contextptr); 1203 } 1204 // should call factor before returning unevaluated rootof 1205 if (e.type==_VECT && e._VECTptr->size()==2 && e._VECTptr->back().type==_VECT){ 1206 vecteur v2=*e._VECTptr->back()._VECTptr; 1207 gen g(1); 1208 lcmdeno(v2,g,contextptr); 1209 if (is_minus_one(v2[0])) 1210 v2=-v2; 1211 if (!is_one(v2[0])) 1212 return gensizeerr("rootof minimal polynomial must be unitary"); 1213 return symbolic(at_rootof,gen(makevecteur(e._VECTptr->front(),gen(v2,e._VECTptr->back().subtype)),e.subtype)); 1214 } 1215 return symbolic(at_rootof,e); 1216 } approx_rootof(const gen & e,GIAC_CONTEXT)1217 gen approx_rootof(const gen & e,GIAC_CONTEXT){ 1218 if ( (e.type!=_VECT) || (e._VECTptr->size()!=2) ) 1219 return gensizeerr(contextptr); 1220 if (!lidnt(e).empty()) 1221 return symbolic(at_rootof,e); 1222 gen a=e._VECTptr->front(),b=e._VECTptr->back(); 1223 return alg_evalf(a,b,contextptr); 1224 } 1225 /* statically in derive.cc 1226 static gen d1_rootof(const gen & args,GIAC_CONTEXT){ 1227 return gentypeerr(contextptr); 1228 return zero; 1229 } 1230 static gen d2_rootof(const gen & args,GIAC_CONTEXT){ 1231 return gentypeerr(contextptr); 1232 return zero; 1233 } 1234 define_unary_function_ptr( D1_rootof,alias_D1_rootof,new unary_function_eval(0,&d1_rootof,"")); 1235 define_unary_function_ptr( D2_rootof,alias_D2_rootof,new unary_function_eval(0,&d2_rootof,"")); 1236 static unary_function_ptr d_rootof(int i){ 1237 if (i==1) 1238 return D1_rootof; 1239 if (i==2) 1240 return D2_rootof; 1241 return gensizeerr(contextptr); 1242 return 0; 1243 } 1244 partial_derivative_multiargs D_rootof(&d_rootof); 1245 */ 1246 static const char _rootof_s []="rootof"; 1247 static define_unary_function_eval2 (__rootof,&rootof,_rootof_s,&printasrootof); 1248 define_unary_function_ptr5( at_rootof ,alias_at_rootof,&__rootof,0,true); 1249 max_algext(const gen & args,GIAC_CONTEXT)1250 gen max_algext(const gen & args,GIAC_CONTEXT){ 1251 gen g=args; 1252 if (!is_integral(g) || g.type!=_INT_ || g.val<3) 1253 return gensizeerr(contextptr); 1254 return MAX_ALG_EXT_ORDER_SIZE=MAX_COMMON_ALG_EXT_ORDER_SIZE=g.val; 1255 } 1256 static const char _max_algext_s []="max_algext"; 1257 static define_unary_function_eval (__max_algext,&max_algext,_max_algext_s); 1258 define_unary_function_ptr5( at_max_algext ,alias_at_max_algext,&__max_algext,0,true); 1259 sturm(const gen & g)1260 static vecteur sturm(const gen & g){ 1261 if (g.type!=_POLY) 1262 return vecteur(1,g); 1263 polynome p(*g._POLYptr); 1264 polynome pl(lgcd(p)); 1265 polynome pp=p/pl; 1266 polynome cont(p.dim); 1267 factorization f(sqff(pp)); 1268 factorization::const_iterator it=f.begin(),itend=f.end(); 1269 gen a=p.coord.front().value; 1270 for (;it!=itend;++it){ 1271 if (it->mult %2) 1272 a=a/it->fact.coord.front().value; 1273 } 1274 vecteur v(1,pl.coord.empty()?a:a/pl.coord.front().value*pl); 1275 for (it=f.begin();it!=itend;++it){ 1276 if (it->mult %2) 1277 v.push_back(sturm_seq(it->fact,cont)); 1278 } 1279 return v; 1280 } sturm(const gen & g,const gen & x,GIAC_CONTEXT)1281 vecteur sturm(const gen &g,const gen & x,GIAC_CONTEXT){ 1282 if (g.type==_VECT) 1283 return vecteur(1,gensizeerr(contextptr)); 1284 vecteur l; 1285 if (!is_zero(x)) 1286 l.push_back(x); 1287 lvar(g,l); 1288 fraction fa(e2r(exact(g,contextptr),l,contextptr)); 1289 gen n,d; 1290 fxnd(fa,n,d); 1291 vecteur v=mergevecteur(sturm(n),sturm(d)); 1292 vecteur res,tmp,ll=cdr_VECT(l); 1293 const_iterateur it=v.begin(),itend=v.end(); 1294 for (;it!=itend;++it){ 1295 if (it->type==_VECT){ 1296 const_iterateur jt=it->_VECTptr->begin(),jtend=it->_VECTptr->end(); 1297 vecteur tmpres; 1298 tmpres.reserve(int(jtend-jt)); 1299 for (;jt!=jtend;++jt){ 1300 if (jt->type==_POLY){ 1301 tmp=polynome2poly1(*(jt->_POLYptr),1); 1302 tmpres.push_back(r2e(tmp,ll,contextptr)); 1303 } 1304 else 1305 tmpres.push_back(*jt); 1306 } 1307 res.push_back(tmpres); 1308 } 1309 else { // it->type != _VECT but we must convert anyway the cst coeff! 1310 if (it->type==_POLY){ 1311 gen tmpg=polynome2poly1(*(it->_POLYptr),1).front(); 1312 res.push_back(r2e(tmpg,ll,contextptr)); 1313 } 1314 else 1315 res.push_back(*it); 1316 } 1317 } 1318 return res; 1319 } 1320 // v is a sequence of dense polynomials 1321 // each poly is evaluated at a, then we count # of sign changes 1322 // ignoring zeros 1323 // The function modifies a sign variable according to the sign first 1324 // non-zero element of v number_of_sign_changes(const vecteur & v,const gen & a0,int & global_sign,GIAC_CONTEXT)1325 static int number_of_sign_changes(const vecteur & v,const gen & a0,int & global_sign,GIAC_CONTEXT){ 1326 gen a=exact(a0,contextptr); 1327 gen w=normal(apply1st(v,a,horner),contextptr); 1328 int previous_sign=0,current_sign,res=0; 1329 const_iterateur it=w._VECTptr->begin(),itend=w._VECTptr->end(); 1330 for (;it!=itend;++it){ 1331 if (is_exactly_zero(*it)) 1332 continue; 1333 if (ck_is_strictly_positive(*it,contextptr)) 1334 current_sign=1; 1335 else 1336 current_sign=-1; 1337 if (!previous_sign) {// assign first non-zero sign 1338 previous_sign=current_sign; 1339 global_sign = global_sign *current_sign; 1340 } 1341 if (previous_sign==current_sign) 1342 continue; 1343 ++res; 1344 previous_sign=current_sign; 1345 } 1346 return res; 1347 } sturmab(const gen & g,const gen & x,const gen & a,const gen & b,bool remove_b_root,GIAC_CONTEXT)1348 static int sturmab(const gen & g,const gen & x,const gen & a,const gen & b,bool remove_b_root,GIAC_CONTEXT){ 1349 if (g.type==_VECT){ 1350 #ifdef NO_STDEXCEPT 1351 return -2; 1352 #else 1353 setsizeerr(contextptr); 1354 #endif 1355 } 1356 if (ck_is_strictly_greater(a,b,contextptr)) 1357 return sturmab(g,x,b,a,contextptr); 1358 if (a==b){ 1359 gen tmp; 1360 if (is_inf(a) && x.type==_IDNT) 1361 tmp=limit(g,*x._IDNTptr,a,0,contextptr); 1362 else 1363 tmp=subst(g,x,a,false,contextptr); 1364 int s=fastsign(tmp,contextptr); 1365 if (s==1 || s==-1) 1366 return (s-1)/2; 1367 } 1368 #ifdef NO_STDEXCEPT 1369 vecteur lvarg(lvar(g)); 1370 if (!lvarg.empty() && lvarg!=vecteur(1,x)) 1371 return -2; 1372 #endif 1373 int res=0,dontcare,global_sign=1; 1374 vecteur v=sturm(g,x,contextptr); 1375 if (is_undef(v)) 1376 return -2; 1377 const_iterateur it=v.begin(),itend=v.end(); 1378 for (;it!=itend;++it){ 1379 if (it->type==_VECT){ 1380 res += number_of_sign_changes(*it->_VECTptr,a,global_sign,contextptr)-number_of_sign_changes(*it->_VECTptr,b,dontcare,contextptr); 1381 if (remove_b_root && is_zero(horner(it->_VECTptr->front(),b))) 1382 --res; 1383 } 1384 else { 1385 if (!ck_is_positive(*it,contextptr)) 1386 global_sign = -global_sign; 1387 } 1388 } 1389 if (res) 1390 return res; 1391 return (global_sign-1)/2; 1392 } sturmab(const gen & g,const gen & x,const gen & a,const gen & b,GIAC_CONTEXT)1393 int sturmab(const gen & g,const gen & x,const gen & a,const gen & b,GIAC_CONTEXT){ 1394 return sturmab(g,x,a,b,false,contextptr); 1395 } _sturmab(const gen & g_orig,GIAC_CONTEXT)1396 gen _sturmab(const gen & g_orig,GIAC_CONTEXT){ 1397 if ( g_orig.type==_STRNG && g_orig.subtype==-1) return g_orig; 1398 if ( g_orig.type!=_VECT || g_orig._VECTptr->size()<3 ) 1399 return gensizeerr(contextptr); 1400 vecteur v(*g_orig._VECTptr); 1401 int s=int(v.size()); 1402 gen P(v[0]),x(vx_var),a,b; 1403 if (s==3){ a=v[1]; b=v[2]; } 1404 else { 1405 x=v[1]; a=v[2]; b=v[3]; 1406 if (P.type==_VECT) 1407 *logptr(contextptr) << gettext("Warning: variable name ignored: ") << x << '\n'; 1408 } 1409 gen ai=im(a,contextptr); 1410 gen bi=im(b,contextptr); 1411 if (!is_zero(ai) || !is_zero(bi)){ 1412 gen p=_e2r(gen(makevecteur(P,vecteur(1,x)),_SEQ__VECT),contextptr),n,d,g1,g2; 1413 if (is_undef(p)) return p; 1414 fxnd(p,n,d); 1415 vecteur nr; 1416 int n1; 1417 #if 0 // replace by 1 if you want to count complex rational root on the edges 1418 if (n.type==_POLY && n._POLYptr->dim==1){ 1419 polynome nrp=*n._POLYptr; 1420 nr=crationalroot(nrp,true); 1421 n1=csturm_square(nrp,a,b,g1,contextptr); 1422 } 1423 else 1424 n1=csturm_square(n,a,b,g1,contextptr); 1425 #else 1426 n1=csturm_square(n,a,b,g1,contextptr); 1427 #endif 1428 int d1=csturm_square(d,a,b,g2,contextptr); 1429 if (n1==-1 || d1==-1) 1430 return gensizeerr(contextptr); 1431 return int(nr.size())+gen(n1)/2+cst_i*gen(d1)/2; 1432 } 1433 if (s==5 && v[4].type==_INT_) 1434 return sturmab(P,x,a,b,v[4].val!=0,contextptr); 1435 return sturmab(P,x,a,b,contextptr); 1436 } 1437 static const char _sturmab_s []="sturmab"; 1438 static define_unary_function_eval (__sturmab,&_sturmab,_sturmab_s); 1439 define_unary_function_ptr5( at_sturmab ,alias_at_sturmab,&__sturmab,0,true); 1440 _sturm(const gen & g,GIAC_CONTEXT)1441 gen _sturm(const gen & g,GIAC_CONTEXT){ 1442 if ( g.type==_STRNG && g.subtype==-1) return g; 1443 if (g.type!=_VECT || (g.type==_VECT && g.subtype!=_SEQ__VECT) ) 1444 return sturm(g,zero,contextptr); 1445 vecteur & v = *g._VECTptr; 1446 int s=int(v.size()); 1447 if (s==2) 1448 return sturm(v.front(),v.back(),contextptr); 1449 if (s==4) 1450 return _sturmab(g,contextptr); 1451 if (s==3){ 1452 if (v[2].type!=_IDNT) 1453 return gensizeerr(contextptr); 1454 gen S=_e2r(gen(makevecteur(v[0],v[2]),_SEQ__VECT),contextptr); 1455 if (is_undef(S)) return S; 1456 gen R=_e2r(gen(makevecteur(v[1],v[2]),_SEQ__VECT),contextptr); 1457 if (is_undef(R)) return R; 1458 if (S.type==_FRAC) 1459 S=S._FRACptr->num; 1460 if (R.type==_FRAC) 1461 R=R._FRACptr->num; 1462 modpoly r0(gen2vecteur(S)),r1(gen2vecteur(R)); 1463 vecteur listquo,coeffP,coeffR; 1464 gen pgcd=csturm_seq(r0,r1,listquo,coeffP,coeffR,contextptr); 1465 return makevecteur(r0,r1,pgcd,listquo,coeffP,coeffR); 1466 } 1467 return gendimerr(contextptr); 1468 } 1469 static const char _sturm_s []="sturm"; 1470 static define_unary_function_eval (__sturm,&_sturm,_sturm_s); 1471 define_unary_function_ptr5( at_sturm ,alias_at_sturm,&__sturm,0,true); 1472 _sturmseq(const gen & g,GIAC_CONTEXT)1473 gen _sturmseq(const gen & g,GIAC_CONTEXT){ 1474 if ( g.type==_STRNG && g.subtype==-1) return g; 1475 // should return in the same format as maple 1476 return _sturm(g,contextptr); 1477 } 1478 static const char _sturmseq_s []="sturmseq"; 1479 static define_unary_function_eval (__sturmseq,&_sturmseq,_sturmseq_s); 1480 define_unary_function_ptr5( at_sturmseq ,alias_at_sturmseq,&__sturmseq,0,true); 1481 recompute_minmax(const vecteur & w,const vecteur & range,const gen & expr,const gen & var,gen & resmin,gen & resmax,vecteur & xmin,vecteur & xmax,int direction,GIAC_CONTEXT)1482 void recompute_minmax(const vecteur & w,const vecteur & range,const gen & expr,const gen & var,gen & resmin,gen & resmax,vecteur & xmin,vecteur & xmax,int direction,GIAC_CONTEXT){ 1483 const_iterateur it=w.begin(),itend=w.end(); 1484 for (;it!=itend;++it){ 1485 if (ck_is_strictly_greater(*it,range[1],contextptr) || ck_is_strictly_greater(range[0],*it,contextptr)) 1486 continue; 1487 #ifdef NO_STDEXCEPT 1488 gen tmp=limit(expr,*var._IDNTptr,*it,direction,contextptr); 1489 #else 1490 gen tmp; 1491 try { 1492 tmp=limit(expr,*var._IDNTptr,*it,direction,contextptr); 1493 } catch (std::runtime_error & err){ 1494 last_evaled_argptr(contextptr)=NULL; 1495 tmp=undef; 1496 } 1497 #endif 1498 if (is_undef(tmp) || tmp==unsigned_inf) 1499 continue; 1500 if (tmp==resmax && !equalposcomp(xmax,*it)) 1501 xmax.push_back(*it); 1502 else { 1503 if (ck_is_strictly_greater(tmp,resmax,contextptr)){ 1504 resmax=tmp; 1505 xmax=vecteur(1,*it); 1506 } 1507 } 1508 if (tmp==resmin && !equalposcomp(xmin,*it)) 1509 xmin.push_back(*it); 1510 else { 1511 if (ck_is_strictly_greater(resmin,tmp,contextptr)){ 1512 resmin=tmp; 1513 xmin=vecteur(1,*it); 1514 } 1515 } 1516 } 1517 } 1518 1519 // minmax=0 both 1 min, 2 max, /3 =1 if return x instead of f(x) fminmax(const gen & g,int minmax,GIAC_CONTEXT)1520 gen fminmax(const gen & g,int minmax,GIAC_CONTEXT){ 1521 gen expr,var; 1522 vecteur v(gen2vecteur(g)); 1523 if (v.size()==1) 1524 v.push_back(vx_var); 1525 if (v.size()!=2) 1526 return gensizeerr(contextptr); 1527 expr=v[0]; 1528 var=v[1]; 1529 // avoid inf recursion like g0(x):=ln(abs(ln(x))); 1530 // g1(x,xp):=x/(ln(x))^(xp);g0(g1(x,.3)); 1531 gen varev=eval(var,1,contextptr); 1532 if (varev!=var && contains(varev,var)) 1533 return undef; 1534 if (expr.type==_SYMB){ 1535 unary_function_ptr & u=expr._SYMBptr->sommet; 1536 if (u==at_exp || u==at_ln || u==at_atan || u==at_abs){ 1537 gen tmp=fminmax(makevecteur(expr._SYMBptr->feuille,var),minmax,contextptr); 1538 if (is_undef(tmp)) 1539 return tmp; 1540 if (u==at_abs && tmp.type==_VECT && tmp._VECTptr->size()==2 ){ 1541 gen t1=tmp._VECTptr->front(); 1542 gen t2=tmp._VECTptr->back(); 1543 if (is_positive(t1,contextptr)) 1544 return tmp; 1545 if (is_positive(-t2,contextptr)){ 1546 return gen(makevecteur(-t2,-t1),_LINE__VECT); 1547 } 1548 // t1<=0 t2>=0 1549 if (is_greater(-t1,t2,contextptr)) 1550 return gen(makevecteur(0,-t1),_LINE__VECT); 1551 else 1552 return gen(makevecteur(0,t2),_LINE__VECT); 1553 } 1554 if (u==at_ln && tmp.type==_VECT && tmp._VECTptr->size()==2 && is_positive(-tmp._VECTptr->front(),contextptr) ) 1555 tmp._VECTptr->front()=zero; 1556 if (minmax/3) 1557 return tmp; 1558 else 1559 return u(tmp,contextptr); 1560 } 1561 } 1562 bool do_find_range=true; 1563 vecteur range; 1564 if (is_equal(var)){ 1565 gen tmp=var._SYMBptr->feuille; 1566 if (tmp.type==_VECT && tmp._VECTptr->size()==2){ 1567 gen varminmax=tmp._VECTptr->back(); 1568 var=tmp._VECTptr->front(); 1569 if (varminmax.is_symb_of_sommet(at_interval) && varminmax._SYMBptr->feuille.type==_VECT){ 1570 range=*varminmax._SYMBptr->feuille._VECTptr; 1571 do_find_range=false; 1572 } 1573 } 1574 } 1575 // gensizeerr replaced by undef because otherwise abs(sin(exp(x))) fails on emcc 1576 if (var.type!=_IDNT) 1577 return undef; // gensizeerr(contextptr); 1578 if (do_find_range){ 1579 find_range(var,range,contextptr); 1580 if (range.size()!=1 || range.front().type!=_VECT) 1581 return gensizeerr(gettext("Or condition not implemented")); 1582 range=*range.front()._VECTptr; 1583 } 1584 if (range.size()!=2) 1585 return gensizeerr(gettext("fminmax, range ")+gen(range).print(contextptr)); 1586 if (range[0]==minus_inf || range[1]==plus_inf){ 1587 // periodic function? 1588 vecteur w=lvarx(trig2exp(expr,contextptr),var); 1589 gen period=0; 1590 for (unsigned i=0;i<w.size();++i){ 1591 if (!w[i].is_symb_of_sommet(at_exp)){ 1592 period=0; 1593 break; 1594 } 1595 gen tmp=w[i]._SYMBptr->feuille,a,b; 1596 if (!is_linear_wrt(tmp,var,a,b,contextptr) || !is_zero(re(a,contextptr))){ 1597 period=0; 1598 break; 1599 } 1600 if (is_zero(a)) 1601 continue; 1602 a=ratnormal(cst_two_pi/im(a,contextptr),contextptr); // current period 1603 if (is_zero(period)) 1604 period=a; 1605 else { // find common period (if it exists) 1606 b=ratnormal(period/a,contextptr); 1607 if (b.type!=_INT_ && b.type!=_FRAC){ 1608 period=0; 1609 break; 1610 } 1611 if (b.type==_FRAC) 1612 period=period*b._FRACptr->den; 1613 } 1614 } 1615 if (!is_zero(period)){ 1616 if (w.size()>1) 1617 expr=simplify(expr,contextptr); 1618 if (range[0]==minus_inf){ 1619 if (range[1]==plus_inf){ 1620 range[1]=period/2; 1621 range[0]=-range[1]; 1622 } 1623 else 1624 range[0]=range[1]-period; 1625 } 1626 else 1627 range[1]=range[0]+period; 1628 } 1629 } 1630 gen df(derive(expr,var,contextptr)); 1631 if (is_undef(df)) 1632 return df; 1633 vecteur w; 1634 if (range==makevecteur(minus_inf,plus_inf)) 1635 w=solve(df,var,2,contextptr); 1636 else { 1637 // FIXME: check if var is quoted, otherwise it will be erased 1638 gen savevar=var; 1639 var._IDNTptr->in_eval(1,var,savevar,contextptr); 1640 // if (var._IDNTptr->in_eval(1,var,savevar,contextptr)) 1; 1641 giac_assume(symbolic(at_and,makevecteur(symb_superieur_egal(var,range[0]),symb_inferieur_egal(var,range[1]))),contextptr); 1642 w=solve(df,var,2,contextptr); 1643 if (savevar==var) 1644 purgenoassume(var,contextptr); 1645 else 1646 sto(savevar,var,contextptr); 1647 } 1648 if (w.empty() && debug_infolevel) 1649 *logptr(contextptr) << gettext("Warning: ") << df << gettext("=0: no solution found") << '\n'; 1650 vecteur wvar=makevecteur(cst_pi); 1651 lidnt(w,wvar,false); 1652 if (wvar.size()>1) 1653 return undef; 1654 gen resmin=plus_inf; 1655 gen resmax=minus_inf; 1656 vecteur xmin,xmax; 1657 // Extrema 1658 recompute_minmax(w,range,expr,var,resmin,resmax,xmin,xmax,0,contextptr); 1659 // Limits at begin and end of range 1660 recompute_minmax(vecteur(1,range[0]),range,expr,var,resmin,resmax,xmin,xmax,1,contextptr); 1661 recompute_minmax(vecteur(1,range[1]),range,expr,var,resmin,resmax,xmin,xmax,-1,contextptr); 1662 // Singularities 1663 vecteur ws=find_singularities(expr,*var._IDNTptr,0,contextptr); 1664 int wss=int(ws.size()); 1665 w.clear(); 1666 for (int i=0;i<wss;++i){ 1667 if (ws[i]!=range[0] && ws[i]!=range[1]) 1668 w.push_back(ws[i]); 1669 } 1670 recompute_minmax(w,range,expr,var,resmin,resmax,xmin,xmax,1,contextptr); 1671 recompute_minmax(w,range,expr,var,resmin,resmax,xmin,xmax,-1,contextptr); 1672 if (minmax/3){ 1673 if (minmax %3 ==1) 1674 return xmin; 1675 if (minmax %3 ==2) 1676 return xmax; 1677 return gen(makevecteur(xmin,xmax),_LINE__VECT); 1678 } 1679 else { 1680 if (minmax %3 ==1) 1681 return resmin; 1682 if (minmax %3 ==2) 1683 return resmax; 1684 return gen(makevecteur(resmin,resmax),_LINE__VECT); 1685 } 1686 } 1687 bool is_constant_idnt(const gen & g); // FIXME -> prog.h 1688 // find extremals values of g 1689 // should be improved (currently return -1..1 for sin and cos find_range(const gen & g,vecteur & a,GIAC_CONTEXT)1690 int find_range(const gen & g,vecteur & a,GIAC_CONTEXT){ 1691 if (g.type==_IDNT){ 1692 gen g2=g._IDNTptr->eval(1,g,contextptr); 1693 if ((g2.type==_VECT) && (g2.subtype==_ASSUME__VECT)){ 1694 vecteur v=*g2._VECTptr; 1695 if ( (v.size()==3) && (v.front()==vecteur(0) || v.front()==_DOUBLE_ || v.front()==_ZINT || v.front()==_SYMB || v.front()==0) && (v[1].type==_VECT)){ 1696 a=*v[1]._VECTptr; 1697 if (v.front()==_ZINT) return 3; 1698 return 1; 1699 } 1700 if (v.size()==1 && v.front()==_ZINT) 1701 return 2; 1702 } 1703 } 1704 if (g.type==_SYMB){ 1705 #ifndef NO_STDEXCEPT 1706 try { 1707 #endif 1708 if (g._SYMBptr->feuille.type==_SPOL1) 1709 return 0; 1710 vecteur lv0(lvar(g._SYMBptr->feuille)),lv; // remove cst idnt 1711 for (unsigned i=0;i<lv0.size();++i){ 1712 if (evalf_double(lv0[i],1,contextptr).type!=_DOUBLE_) // if (!is_constant_idnt(lv0[i])) 1713 lv.push_back(lv0[i]); 1714 } 1715 if (!lv.empty()){ 1716 gen res=fminmax(makevecteur(g,lv[0]),0,contextptr); 1717 if (is_undef(res)) 1718 return 0; 1719 a=vecteur(1,res); 1720 return 1; 1721 } 1722 #ifndef NO_STDEXCEPT 1723 } 1724 catch (std::runtime_error & ){ 1725 last_evaled_argptr(contextptr)=NULL; 1726 } 1727 #endif 1728 unary_function_ptr s(g._SYMBptr->sommet); 1729 if ( (s==at_sin) || (s==at_cos) ){ 1730 a=vecteur(1,gen(makevecteur(minus_one,plus_one),_LINE__VECT)); 1731 return 1; 1732 } 1733 } 1734 a=vecteur(1,gen(makevecteur(minus_inf,plus_inf),_LINE__VECT)); 1735 return 1; 1736 } 1737 is_sqrt(const gen & a,gen & arg)1738 bool is_sqrt(const gen & a,gen & arg){ 1739 if (a.is_symb_of_sommet(at_sqrt)){ 1740 arg=a._SYMBptr->feuille; 1741 return true; 1742 } 1743 if (!a.is_symb_of_sommet(at_pow)) 1744 return false; 1745 gen & f = a._SYMBptr->feuille; 1746 if (f.type!=_VECT || f._VECTptr->size()!=2) 1747 return false; 1748 arg = f._VECTptr->front(); 1749 gen & expo = f._VECTptr->back(); 1750 if (expo.type!=_FRAC || !is_one(expo._FRACptr->num)) 1751 return false; 1752 gen & d =expo._FRACptr->den; 1753 if (d.type!=_INT_ || d.val!=2) 1754 return false; 1755 return true; 1756 } 1757 insturmsign1(const gen & g0,bool strict,GIAC_CONTEXT)1758 static int insturmsign1(const gen & g0,bool strict,GIAC_CONTEXT){ 1759 gen g=recursive_normal(exact(g0,contextptr),contextptr); 1760 if (has_i(g)) 1761 return 0; 1762 vecteur v(lvar(g)); 1763 // search for a sqrt inside v: sign(a+b*sqrt(c))= 1764 // = sign(a) if a^2-c*b^2 > 0, 1765 // = sign(b) if a^2-c*b^2 < 0 1766 int s=int(v.size()); 1767 if (!s 1768 #ifdef EMCC 1769 || s>1 1770 #endif 1771 ) 1772 return fastsign(g,contextptr); 1773 gen v0(v[0]); 1774 for (int i=0;i<s;++i){ // replace by first idnt with an assumption 1775 if (v[i].type==_IDNT && v[i]._IDNTptr->eval(1,v[i],contextptr).type!=_IDNT){ 1776 v0=v[i]; 1777 } 1778 gen a,b,c; 1779 if (is_sqrt(v[i],c)){ 1780 identificateur x(" x"); 1781 gen g1=subst(g,v[i],x,false,contextptr); 1782 if (is_linear_wrt(g1,x,b,a,contextptr)){ 1783 gen s=sign(a*b,contextptr); 1784 if (is_one(s) && (s=sign(a,contextptr)).type==_INT_) 1785 return s.val; 1786 s=sign(a*a-c*b*b,contextptr); 1787 if (s.type!=_INT_ || is_zero(s.val)) 1788 return 0; 1789 s=(is_one(s))?sign(a,contextptr):sign(b,contextptr); 1790 if (is_one(s)) 1791 return 1; 1792 if (is_minus_one(s)) 1793 return -1; 1794 return 0; 1795 } 1796 } 1797 } 1798 vecteur a; 1799 if (!find_range(v0,a,contextptr)) 1800 return -2; 1801 int previous_sign=2,current_sign=0; 1802 #ifndef NO_STDEXCEPT 1803 try { 1804 #endif 1805 const_iterateur ita=a.begin(),itaend=a.end(); 1806 for (;ita!=itaend;++ita){ 1807 if ( (ita->type!=_VECT) || (ita->subtype!=_LINE__VECT) || (ita->_VECTptr->size()!=2) ) 1808 return 0; 1809 gen last(ita->_VECTptr->back()); 1810 gen gg(g); 1811 identificateur idnttmp("t"); 1812 gen testg(subst(g,v0,idnttmp,false,contextptr)); 1813 if (is_zero(limit(testg,idnttmp,last,-1,contextptr))){ 1814 if (strict && (v0.is_symb_of_sommet(at_sin) || v0.is_symb_of_sommet(at_cos))) 1815 return 0; 1816 gen tmp=_fxnd(gg,contextptr); 1817 if (tmp.type!=_VECT || tmp._VECTptr->size()!=2){ 1818 #ifdef NO_STDEXCEPT 1819 return -2; 1820 #else 1821 setsizeerr(contextptr); 1822 #endif 1823 } 1824 gen num=tmp._VECTptr->front(),den=tmp._VECTptr->back(),tmpden; 1825 tmp=_e2r(makevecteur(num,v0),contextptr); 1826 tmpden=_e2r(makevecteur(den,v0),contextptr); 1827 if (is_undef(tmp) || is_undef(tmpden)) 1828 return -2; 1829 if (is_inf(last) && tmpden.type==_VECT) 1830 den=den/pow(v0,2*int(tmpden._VECTptr->size()/2)); 1831 tmp=gen2vecteur(tmp); 1832 modpoly p(*tmp._VECTptr),q; 1833 if (!is_inf(last)) { 1834 while (is_zero(horner(p,last,0,q))) 1835 p=-q; 1836 } 1837 gg=_r2e(gen(makevecteur(p,v0),_SEQ__VECT),contextptr)/den; 1838 } 1839 current_sign=sturmab(gg,v0,ita->_VECTptr->front(),last,true,contextptr); 1840 if (current_sign>0 || current_sign==-2) 1841 return 0; 1842 if (previous_sign==2) 1843 previous_sign=current_sign; 1844 if (previous_sign!=current_sign) 1845 return 0; 1846 } 1847 #ifndef NO_STDEXCEPT 1848 } 1849 catch (std::runtime_error & ){ 1850 last_evaled_argptr(contextptr)=NULL; 1851 return 0; 1852 } 1853 #endif 1854 return 2*current_sign+1; 1855 } 1856 insturmsign(const gen & g0,bool strict,GIAC_CONTEXT)1857 static int insturmsign(const gen & g0,bool strict,GIAC_CONTEXT){ 1858 //bool absb=eval_abs(contextptr); 1859 //eval_abs(false,contextptr); 1860 int res=insturmsign1(g0,strict,contextptr); 1861 return res; 1862 //eval_abs(absb,contextptr); 1863 } 1864 sturmsign(const gen & g0,bool strict,GIAC_CONTEXT)1865 int sturmsign(const gen & g0,bool strict,GIAC_CONTEXT){ 1866 int fs=fastsign(g0,contextptr); 1867 if (fs) return fs; 1868 gen g=simplifier(g0,contextptr); 1869 // first check some operators inv, *, exp, sqrt 1870 int tmp; 1871 if (g.is_symb_of_sommet(at_neg)){ 1872 tmp=sturmsign(g._SYMBptr->feuille,strict,contextptr); 1873 return tmp==-2?tmp:-tmp; 1874 } 1875 if (g.is_symb_of_sommet(at_inv)){ 1876 tmp=sturmsign(g._SYMBptr->feuille,strict,contextptr); 1877 return tmp; 1878 } 1879 if (g.is_symb_of_sommet(at_exp)) 1880 return 1; 1881 /* if (g.is_symb_of_sommet(at_pow) && g._SYMBptr->feuille[1]==plus_one_half) 1882 return 1; */ 1883 if (g.is_symb_of_sommet(at_prod)){ 1884 gen &f=g._SYMBptr->feuille; 1885 vecteur v(gen2vecteur(f)); 1886 int s=int(v.size()); 1887 vecteur w; 1888 int res=1,currentsign; 1889 // remove cst coeffs and exp/ 1890 for (int i=0;i<s;++i){ 1891 if (v[i].is_symb_of_sommet(at_sqrt) && sturmsign(v[i]._SYMBptr->feuille,strict,contextptr)==1) 1892 continue; 1893 if ( (currentsign=fastsign(v[i],contextptr)) ) 1894 res *= currentsign; 1895 else 1896 w.push_back(v[i]); 1897 } 1898 switch (w.size()){ 1899 case 0: 1900 return res; 1901 case 1: 1902 tmp=insturmsign(w.front(),strict,contextptr); return tmp==-2?-2:res*tmp; 1903 default: 1904 tmp=insturmsign(symbolic(at_prod,w),strict,contextptr); return tmp==-2?-2:res*tmp; 1905 } 1906 } 1907 return insturmsign(g,strict,contextptr); 1908 } 1909 1910 #ifndef NO_NAMESPACE_GIAC 1911 } // namespace giac 1912 #endif // ndef NO_NAMESPACE_GIAC 1913