tests/rtest8.mac

   1 /*************** -*- Mode: MACSYMA; Package: MAXIMA -*-  ******************/
   2 /***************************************************************************
   3 ***                                                                    *****
   4 ***     Copyright (c) 1984 by William Schelter,University of Texas     *****
   5 ***     All rights reserved                                            *****
   6 ***************************************************************************/
   7
   8
   9 block([],kill(all),%rnum:0);
  10 0$
  11 x^10-2*x^4+1/2;
  12 x^10-2*x^4+1/2$
  13 nroots(%,-6,9.1);
  14 4$
  15 realroots(x^5-x-1,5.0e-6);
  16 [x = 612003/524288]$
  17 ev(%[1],float);
  18 x = 1.1673030853271484$
  19 ev(x^5-x-1,%);
  20 -7.396496210176906e-6;
  21
  22 /* This failed when a gcl bug set rootsepsilon to 0.0, 2013-10-01 */
  23 realroots(x^3+2*x^2-2*x-1);
  24 [x = -87846643/33554432,x = -12816653/33554432,x = 1];
  25
  26 (2*x+1)^3 = 13.5*(x^5+1);
  27 (2*x+1)^3 = 13.5*(x^5+1)$
  28 sort(allroots(%));
  29 [x = -1.0157555438281209,x = 0.82967499021293611,x = 1.0,
  30        x = -0.96596251521963683*%i-0.40695972319240747,
  31        x = 0.96596251521963683*%i-0.40695972319240747];
  32
  33 /* SF bug [ 1951128 ] curious warning from allroots(x=0) */
  34
  35 allroots (x = 0);
  36 [x = 0.0];
  37
  38 allroots (17*x = 0);
  39 [x = 0.0];
  40
  41 allroots (19*x^4 = 0);
  42 [x = 0.0, x = 0.0, x = 0.0, x = 0.0];
  43
  44 allroots (x^3 * (x^4 - 1) = 0);
  45 [x = 0.0, x = 0.0, x = 0.0,
  46  x = 1.0, x = - 1.0, x = 1.0*%i, x = - 1.0*%i];
  47
  48 allroots (%i*x^5 = 0); /* this one goes through CPOLY-SL */
  49 [x = 0.0, x = 0.0, x = 0.0, x = 0.0, x = 0.0];
  50
  51 /* additional tests for allroots */
  52
  53 allroots (x = 1);
  54 [x = 1.0];
  55
  56 allroots (8*u = 1);
  57 [u = 0.125];
  58
  59 allroots (u^2 - 2*u = 35);
  60 [u = -5.0, u = 7.0];
  61
  62 (complex_float_approx_equal (a, b) :=
  63     if listp (a) and listp (b)
  64         then apply ("and", map (complex_float_approx_equal, a, b))
  65         elseif equationp (a) and equationp (b)
  66             then is (lhs (a) = lhs (b))
  67                 and complex_float_approx_equal (rhs (a), rhs (b))
  68             else
  69                 my_float_approx_equal (realpart (a), realpart (b))
  70                     and my_float_approx_equal (imagpart (a), imagpart (b)),
  71
  72  equationp (e) := not atom (e) and op (e) = "=",
  73
  74  my_float_approx_equal (x, y) :=
  75     if equal (y, 0)
  76         then is (abs (x) <= float_approx_equal_tolerance)
  77         else float_approx_equal (x, y),
  78
  79  float_approx_equal_tolerance : 1e-12,
  80
  81  0);
  82 0;
  83
  84 /* (u - 5/4)*(u - 7/4)*(u + 1/4)*(u^2 - 2*u + 5/4) which has roots
  85  * 5/4, 7/4, -1/4, and 1 + %i/2, 1 - %i/2.
  86  */
  87 complex_float_approx_equal
  88    (allroots (u^5 + 131*u^3/16 + 45*u/64 + 175/256 = 19*u^4/4 + 369*u^2/64),
  89     [u = -0.25, u = 0.5*%i + 1, u = 1 - 0.5*%i, u = 1.25, u = 1.75]);
  90 true;
  91
  92 /* (v - 5/4)*(v - 7/4*%i)*(v + 1/4*%i)*(v^2 - 2*v + 5/4) which has roots
  93  * 5/4, 7/4*%i, -1/4*%i, and 1 + %i/2, 1 - %i/2.
  94  */
  95 complex_float_approx_equal
  96    (allroots (expand ((v - 5/4)*(v - 7/4*%i)*(v + 1/4*%i)*(v^2 - 2*v + 5/4))),
  97     [v = - 0.25*%i,
  98      v = 1.25,
  99      v = 1.0 - 0.5*%i,
 100      v = 0.5*%i + 1.0,
 101      v = 1.75*%i]);
 102 true;
 103
 104 reset (float_approx_equal_tolerance);
 105 [float_approx_equal_tolerance];
 106
 107 exp1:x+z = y;
 108 z+x = y$
 109 exp:2*a*x-y = 2*a^2;
 110 2*a*x-y = 2*a^2$
 111 y-2*z = 2;
 112 y-2*z = 2$
 113 ev(linsolve([exp,exp1,%],[x,y,z]),globalsolve);
 114 [x = a+1,y = 2*a,z = a-1]$
 115
 116 /* see http://trac.sagemath.org/sage_trac/ticket/8731 and
 117 https://sourceforge.net/tracker/index.php?func=detail&aid=2990307&group_id=4933&atid=104933 */
 118 algsys([8*x=1],[x])$
 119 [[x = 1/8]]$
 120 block([realonly:true],algsys([8*x=1],[x]))$
 121 [[x = 1/8]]$
 122 algsys([x^2+1],[x])$
 123 [[x = %i], [x = - %i]]$
 124 block([realonly:true],algsys([x^2+1],[x]))$
 125 []$
 126 /* see https://sourceforge.net/tracker/index.php?func=detail&aid=2786017&group_id=4933&atid=104933 */
 127 algsys([x^4+1],[x])$
 128 [[x = (-1)^(1/4)],[x = -(-1)^(1/4)*%i],[x = -(-1)^(1/4)],[x = sqrt(-%i)]]$
 129 block([realonly:true],algsys([x^4+1],[x]))$
 130 []$
 131 algsys([x^4-1],[x])$
 132 [[x = 1],[x = -1],[x = %i],[x = -%i]]$
 133 algsys([x^4-1],[x]),realonly=true $
 134 [[x = 1],[x = -1]]$
 135
 136 block(
 137   [f1:2*x*(1-l1)-2*(x-1)*l2,
 138   f2:l2-l1,
 139   f3:l1*(1-x^2-y),
 140   f4:l2*(y-(x-1)^2)],
 141   algsys([f1,f2,f3,f4],[x,y,l1,l2]))$
 142 [[x = 0,y = %r1,l1 = 0,l2 = 0],[x = 1,y = 0,l1 = 1,l2 = 1]]$
 143 block(
 144   [f1:x^2-y^2,
 145   f2:x^2-x+2*y^2-y-1],
 146   algsys([f1,f2],[x,y]))$
 147 [[x = -1/sqrt(3),y = 1/sqrt(3)],[x = 1/sqrt(3),y = -1/sqrt(3)],
 148  [x = -1/3,y = -1/3],[x = 1,y = 1]]$
 149 solve(asin(cos(3*x))*(f(x)-1),x);
 150 [x = %pi/6,f(x) = 1]$
 151 ev(solve(5^f(x) = 125,f(x)),solveradcan:true);
 152 [f(x) = 3]$
 153
 154
 155 (float_approx_equal_tolerance : 1e-12, 0);
 156 0$
 157
 158 solve([4*x^2-y^2 = 12,x*y-x = 2], [x,y]);
 159 [[x = 2,y = 2],
 160  [x = 0.5202594388652008*%i-0.1331240357358706,
 161   y = 0.07678378523787777-3.608003221870287*%i],
 162  [x = -0.5202594388652008*%i-0.1331240357358706,
 163   y = 3.608003221870287*%i+0.07678378523787777],
 164  [x = -1.733751846381093,y = -0.1535675710019696]]$
 165
 166 (reset (float_approx_equal_tolerance), 0);
 167 0$
 168
 169 (eq :1+a*x+x^3, sol : solve(eq,x), makelist(ratsimp(subst(s,eq)),s, sol));
 170 [0,0,0]$
 171
 172 (eq : x^3-1, sol : solve(eq,x), makelist(ratsimp(subst(s,eq)),s,sol));
 173 [0,0,0]$
 174
 175 sol : solve(x^6-1);
 176 [x = (sqrt(3)*%i+1)/2,x = (sqrt(3)*%i-1)/2,x = -1,x = -((sqrt(3)*%i+1)/2), x = -((sqrt(3)*%i-1)/2),x = 1]$
 177
 178 ratsimp(makelist(subst(s, x^6-1), s, sol));
 179 [0,0,0,0,0,0]$
 180
 181 (remvalue(eq, sol),0);
 182 0$
 183
 184 exp:x^2-1;
 185 x^2-1$
 186 solve(%,x);
 187 [x = -1,x = 1]$
 188 ev(exp,%[1]);
 189 0$
 190 h[i,j]:=1/(i+j-1);
 191 h[i,j]:=1/(i+j-1)$
 192 genmatrix(h,3,3);
 193 matrix([1,1/2,1/3],[1/2,1/3,1/4],[1/3,1/4,1/5])$
 194 [2*x-(a-1)*y = 5*b,a*x+b*y+c = 0];
 195 [2*x-(a-1)*y = 5*b,b*y+a*x+c = 0]$
 196 augcoefmatrix(%,[x,y]);
 197 matrix([2,1-a,-5*b],[a,b,c])$
 198 matrix([2,1-a,-5*b],[a,b,c]);
 199 matrix([2,1-a,-5*b],[a,b,c])$
 200
 201 echelon(%);
 202 matrix([1,-((a-1)/2),-(5*b/2)],[0,1,(2*c+5*a*b)/(2*b+a^2-a)])$
 203
 204 matrix([2,1-a,-5*b],[a,b,c]);
 205 matrix([2,1-a,-5*b],[a,b,c])$
 206 triangularize(%);
 207 matrix([2,1-a,-5*b],[0,2*b+a^2-a,2*c+5*a*b])$
 208 matrix([2,1-a,-5*b],[a,b,c]);
 209 matrix([2,1-a,-5*b],[a,b,c])$
 210 rank(%);
 211 2$
 212 a:matrix([3,1],[2,4]);
 213 matrix([3,1],[2,4])$
 214 expand(charpoly(a,lambda));
 215 lambda^2-7*lambda+10$
 216 exp:(programmode:true,solve(%));
 217 [lambda = 5,lambda = 2]$
 218 matrix([x1],[x2]);
 219 matrix([x1],[x2])$
 220 ev(a . %-lambda*%,exp[1]);
 221 matrix([x2-2*x1],[2*x1-x2])$
 222 exp:%[1,1] = 0;
 223 x2-2*x1 = 0$
 224 x1^2+x2^2 = 1;
 225 x2^2+x1^2 = 1$
 226 solve([exp,%],[x1,x2]);
 227 [[x1 = -1/sqrt(5),x2 = -2/sqrt(5)],[x1 = 1/sqrt(5),x2 = 2/sqrt(5)]]$
 228
 229 /* verify that find_root is happy with %e
 230  * (problem reported on mailing list 2007/01/25)
 231  */
 232
 233 find_root (2*x = -(log((4 + %e)/(2*%pi)))*(((4 + %e)/(2*%pi))^x), x, -1, 0);
 234 -0.03340289826874122;
 235
 236 find_root (2*x = cos((%e + %pi)*x), x, 0, 1);
 237 0.1984210505656873;
 238
 239 /* verify that find_root evaluates its first argument
 240  * (problem reported to mailing list 2007/06/07)
 241  */
 242 (expr : x^2 - 5, find_root (expr, x, 0, 10));
 243 sqrt (5.0);
 244
 245 /* other tests for find_root:
 246  * verify that find_root expression is returned for non-numeric expression or bounds
 247  */
 248 (kill (a, b), find_root (x^a - 5, x, 0, 10));
 249 find_root (x^a - 5, x, 0.0, 10.0);
 250
 251 find_root (x^3 - 5, x, a, b);
 252 find_root (x^3 - 5, x, a, b);
 253
 254 quad_closeto(qest,qtru,qtol) :=
 255   map(lambda([est, tru, tol],
 256              block([numer:true, abse],
 257                    abse:abs(est-tru),
 258                    if (abse <= tol) then true else abse)),
 259       qest, qtru, qtol);
 260 quad_closeto(qest,qtru,qtol) :=
 261   map(lambda([est, tru, tol],
 262              block([numer:true, abse],
 263                    abse:abs(est-tru),
 264                    if (abse <= tol) then true else abse)),
 265       qest, qtru, qtol);
 266
 267 /* verify that find_root nested inside another function call is OK */
 268 quad_closeto(quad_qags (find_root (x^a - 5, x, 0, 10), a, 1, 3),
 269              quad_qags (5^(1/a), a, 1, 3),
 270              [1d-15, 5d-15, 0, 0]);
 271 [true,true,true,true];
 272
 273 find_root (find_root (a^2 = x, a, 0, x) = 7, x, 0, 100); /* inner find_root returns sqrt(x) */
 274 49.0;
 275
 276 /* verify that symbolic function name is OK */
 277 (foo (a) := 3^a - 5, bar : foo, find_root (bar, 0, 10));
 278 log(5.0) / log(3.0);
 279
 280 /* example from mailing list 2006/12/01 */
 281 (expr : t = (297 * exp ((1000000 * t) / 33) - 330) / 10000000, find_root (expr, t, 1e-9, 0.003));
 282 1.7549783076857198E-5;
 283
 284 /* example from mailing list 2007/06/07 */
 285 (expr : 6096 * tan((2 * atan(c/(2 * fl))) / r) / (tan((1/60) * (%pi/180))),
 286  ev (find_root (expr=6096, fl, 1, 10), c=7.176, r=3264));
 287 6.98149293282488;
 288
 289 /* adapted from mailing list 2007/01/13 */
 290
 291 (g (a) := find_root (f (x, a), x, 0, 200),
 292  f (x, a) := x^a - 5,
 293  0);
 294 0;
 295
 296 g (0.5);
 297 25.0;
 298
 299 expr : g (z + z);
 300 find_root (x^(2 * z) - 5, x, 0.0, 200.0);
 301
 302 ''(at (expr, z=0.25));
 303 25.0;
 304
 305 quad_closeto(quad_qags (g (z), z, 1, 3),
 306              quad_qags (5^(1/z), z, 1, 3),
 307              [1d-15, 5d-15, 0, 0]);
 308 [true,true,true,true];
 309
 310 /* adapted from the reference manual */
 311
 312 (f(x) := sin(x) - x/2, 0);
 313 0;
 314
 315 [find_root (sin(x) - x/2, x, 0.1, %pi),
 316  find_root (sin(x) = x/2, x, 0.1, %pi),
 317  find_root (f(x), x, 0.1, %pi),
 318  find_root (f, 0.1, %pi)];
 319 [1.895494267033981, 1.895494267033981, 1.895494267033981, 1.895494267033981];
 320
 321 [find_root (f, 1/(%pi*%e), 2*%pi*sin(%e)),
 322  find_root (f, log(%pi), %e^%pi),
 323  find_root (f, exp(1/5), exp(cos(%e + %pi))),
 324  find_root (f, cos(exp(2))/10, 10*cos(exp(2)))];
 325 [1.895494267033981, 1.895494267033981, 1.895494267033981, 1.895494267033981];
 326
 327 /* adapted from the mailing list 2007/06/10
 328  * charfun2 copied from the interpol share package
 329  */
 330
 331 block ([expr],
 332  charfun2 (z, l1, l2) := charfun (l1 <= z and z < l2),
 333  expr : (-.329*x^3+.494*x^2 +.559*x+.117) *charfun2(x,minf,1.0)
 334     +(.215*x^3-1.94*x^2 +4.85*x-2.77) *charfun2(x,2.5,inf) +(.0933*x^3-1.02*x^2
 335     +2.56*x-.866) *charfun2(x,2.0,2.5) +(.0195*x^3-.581*x^2
 336     +1.67*x-.275) *charfun2(x,1.5,2.0) +(.00117*x^3-.498*x^2 +1.55*x -.213)
 337     *charfun2(x,1.0,1.5),
 338  block ([float_approx_equal_tolerance : 1e-12],
 339     float_approx_equal (find_root (expr, x, 0, 4), 3.127271310605426)));
 340 true;
 341
 342 /* SF bug report [ 607079 ] solve with repeated variable
 343  */
 344 solve ('[x = 1], '[x, x]);
 345 [x = 1];
 346
 347 solve ('[u = v], '[u, u, u, u]);
 348 [u = v];
 349
 350 /* verify that quadpack functions return partially-evaluated expressions
 351  * instead of barfing on non-numeric values in limits or integrand.
 352  */
 353
 354 (kill (foo, u, au, bu, cu, omega, trig, alfa, vita, wfn), 0);
 355 0;
 356
 357 e1 : quad_qag (foo (u), u, au, bu, 3);
 358 quad_qag (foo (u), u, au, bu, 3, epsrel = 1.0E-8, epsabs = 0.0, limit = 200);
 359
 360 e1 : ev (e1, foo(u)=u^3, au=1);
 361 quad_qag (u^3, u, 1, bu, 3, epsrel=1e-8, epsabs=0.0, limit=200);
 362
 363 ev (e1, bu=4);
 364 [63.75, 7.077671781985375E-13, 31, 0];
 365
 366 e2 : quad_qags (foo (u), u, au, bu, epsrel=1e-4, epsabs=1e-4);
 367 quad_qags (foo (u), u, au, bu, epsrel=1e-4, epsabs=1e-4, limit=200);
 368
 369 e2 : ev (e2, au= -1, bu=1);
 370 quad_qags (foo (u), u, -1, 1, epsrel=1e-4, epsabs=1e-4, limit=200);
 371
 372 ev (e2, foo(u)=u^4);
 373 [0.4, 4.440892098500628E-15, 21, 0];
 374
 375 e3 : quad_qagi (foo (u), u, minf, au, epsabs=2e-3);
 376 quad_qagi (foo (u), u, minf, au, epsrel=1.0E-8, epsabs=2e-3, limit=200);
 377
 378 e3 : ev (e3, foo(u)=1/u^3);
 379 quad_qagi (1/u^3, u, minf, au, epsrel=1.0E-8, epsabs=2e-3, limit=200);
 380
 381 ev (e3, au= -1);
 382 [- 0.5, 5.551115123125784E-15, 15, 0];
 383
 384 e4 : quad_qawc (foo (u), u, cu, au, bu, limit=16);
 385 quad_qawc (foo (u), u, cu, au, bu, epsrel=1.0E-8, epsabs=0.0, limit=16);
 386
 387 e4 : ev (e4, cu=1, au=0, bu=2);
 388 quad_qawc (foo (u), u, 1, 0, 2, epsrel=1.0E-8, epsabs=0.0, limit=16);
 389
 390 ev (e4, foo(u)=u);
 391 [1.999999999999999, 2.220446049250313E-16, 25, 0];
 392
 393 e5 : quad_qawf (foo (u), u, au, omega, sin, limit=32);
 394 quad_qawf (foo (u), u, au, omega, sin, epsabs=1e-10, limit=32, maxp1=100, limlst=10);
 395
 396 e5 : ev (e5, foo(u)=exp(-u));
 397 quad_qawf (exp (- u), u, au, omega, sin, epsabs=1e-10, limit=32, maxp1=100, limlst=10);
 398
 399 ev (e5, au=0, omega=2);
 400 [.4000000000000001, 2.216570948815925E-11, 175, 0];
 401
 402 e6 : quad_qawo (foo (u), u, au, bu, omega, cos, limit=64);
 403 quad_qawo (foo (u), u, au, bu, omega, cos, epsrel=1e-8, epsabs=0.0, limit=64, maxp1=100);
 404
 405 e6 : ev (e6, au=0, bu=%pi/2);
 406 quad_qawo (foo (u), u, 0, %pi/2, omega, cos, epsrel=1e-8, epsabs=0.0, limit=64, maxp1=100);
 407
 408 ev (e6, foo(u)=1, omega=1);
 409 [1.0, 1.110223024625157E-14, 15, 0];
 410
 411 e7 : quad_qaws (foo (u), u, au, bu, alfa, vita, wfn, limit=48);
 412 quad_qaws (foo (u), u, au, bu, alfa, vita, wfn, epsrel=1e-8, epsabs=0.0, limit=48);
 413
 414 e7 : ev (e7, foo(u)=1/u, au=1, bu=2, wfn=1);
 415 quad_qaws (1/u, u, 1, 2, alfa, vita, 1, epsrel=1e-8, epsabs=0.0, limit=48);
 416
 417 /* expect [.05296102778655729, 5.551115123125782E-17, 50, 0] */
 418 (ev (e7, alfa=2, vita=1),
 419  [float_approx_equal (%%[1], .05296102778655729),
 420   /* checking relative error is problematic when expected value is close to zero; check absolute error instead */
 421   is (abs (%%[2] - 5.551115123125782E-17) < float_approx_equal_tolerance),
 422   is (%%[3] = 50),
 423   is (%%[4] = 0)]);
 424 [true, true, true, true];
 425
 426 /* Tests for bfallroots.  Same as the allroots tests above */
 427 bfallroots(x=0);
 428 [x = 0b0];
 429
 430 bfallroots (17*x = 0);
 431 [x = 0b0];
 432
 433 bfallroots (19*x^4 = 0);
 434 [x = 0b0, x = 0b0, x = 0b0, x = 0b0];
 435
 436 bfallroots (x^3 * (x^4 - 1) = 0);
 437 [x = 0b0, x = 0b0, x = 0b0,
 438  x = 1b0, x = - 1b0, x = 1b0*%i, x = - 1b0*%i];
 439
 440 bfallroots (%i*x^5 = 0); /* this one goes through CPOLY-SL */
 441 [x = 0b0, x = 0b0, x = 0b0, x = 0b0, x = 0b0];
 442
 443 /* additional tests for bfallroots */
 444
 445 bfallroots (x = 1);
 446 [x = 1b0];
 447
 448 bfallroots (8*u = 1);
 449 [u = 0.125b0];
 450
 451 bfallroots (u^2 - 2*u = 35);
 452 [u = -5b0, u = 7b0];
 453
 454 (float_approx_equal_tolerance : 1e-12, 0);
 455 0;
 456
 457 /* (u - 5/4)*(u - 7/4)*(u + 1/4)*(u^2 - 2*u + 5/4) which has roots
 458  * 5/4, 7/4, -1/4, and 1 + %i/2, 1 - %i/2.
 459  */
 460 complex_float_approx_equal
 461    (bfallroots (u^5 + 131*u^3/16 + 45*u/64 + 175/256 = 19*u^4/4 + 369*u^2/64),
 462     [u = -0.25, u = 0.5*%i + 1, u = 1 - 0.5*%i, u = 1.25, u = 1.75]);
 463 true;
 464
 465 /* (v - 5/4)*(v - 7/4*%i)*(v + 1/4*%i)*(v^2 - 2*v + 5/4) which has roots
 466  * 5/4, 7/4*%i, -1/4*%i, and 1 + %i/2, 1 - %i/2.
 467  */
 468 complex_float_approx_equal
 469    (bfallroots (expand ((v - 5/4)*(v - 7/4*%i)*(v + 1/4*%i)*(v^2 - 2*v + 5/4))),
 470     [v = - 0.25*%i,
 471      v = 1.25,
 472      v = 1.0 - 0.5*%i,
 473      v = 0.5*%i + 1.0,
 474      v = 1.75*%i]);
 475 true;
 476
 477 /* [ 940835 ] rectform fails with float/numer flags */
 478 rectform(log(-%i)),float;
 479 -0.5 * %i * %pi;
 480
 481 /* verify that exp(foo) evaluates to a number
 482  * probably should try several variations on this
 483  * adapted from sage mailing list
 484  */
 485
 486 first (quad_qags (sin (%pi * exp (x / 2)), x, 0, 2));
 487 - 0.4373454748252497;
 488
 489 /* verify that nested numerical integral is handled correctly
 490  * adapted from sage mailing list
 491  */
 492
 493 quad_qags (w^2 * quad_qags (1/(s - w), s, 1, 5) [1], w, -5, -1) [1];
 494 25.83639378805382;
 495
 496 /* find_root example from sage mailing list */
 497
 498 find_root (.05^(x + 1) = (x + 1)*10^(-10), x, 5, 100);
 499 6.034992572983213;
 500
 501 /* another nested example, collected from mma user forum */
 502
 503 (f : diff (1/(1 + (1 + (a - b)^2)), a),
 504  g : quad_qags (f*b*(1 - b)^2, b, 0, 1) [1],
 505  find_root (g = 0, a, 0, 1));
 506 0.3978613590133817;
 507
 508 /* a variation -- not sure what g:=... means in mma */
 509
 510 (f : diff (1/(1 + (1 + (a - b)^2)), a),
 511  g : 'quad_qags (f*b*(1 - b)^2, b, 0, 1) [1],
 512  find_root (g = 0, a, 0, 1));
 513 0.3978613590133817;
 514
 515 /* from mailing list 2009-02-18
 516  * "Re: [Maxima] I want to tell maxima (-1)^0.33333333333333=-1, what should i do?"
 517  * see also tests/rtest_plot
 518  */
 519
 520 (foo17(x):=(sqrt(-16*x^4-16*x^3+20*x^2+12*x+23)/(6*sqrt(3))+(16*x^3-12*x^2-6*x-25)/54)^(1/3),
 521  float_approx_equal_tolerance : 1e-12,
 522  0);
 523 0;
 524
 525 first (quad_qags (foo17 (u), u, -1, 0));
 526 - 0.359753467469551;
 527
 528 first (quad_qags (foo17, u, -1, 0));
 529 - 0.359753467469551;
 530
 531 find_root (foo17 (u) = -0.2, u, -1, 0);
 532 - 0.246809031968399;
 533
 534 (bar17 (u) := foo17 (u) + 0.2, find_root (bar17, u, -1, 0));
 535 - 0.246809031968399;
 536
 537 (compile (foo17), first (quad_qags (foo17, u, -1, 0)));
 538 - 0.359753467469551;
 539
 540 find_root (bar17, u, -1, 0);
 541 - 0.246809031968399;
 542
 543 /* SF bug # 2937837 "find_root_error documentation incorrect"
 544  */
 545
 546 (find_root_error : true,
 547  errcatch (find_root (1 + x^2, x, 0, 1)));
 548 [];
 549
 550 (find_root_error : "FOO",
 551  errcatch (find_root (1 + x^2, x, 0, 1)));
 552 ["FOO"];
 553
 554 reset (float_approx_equal_tolerance, find_root_error);
 555 [float_approx_equal_tolerance, find_root_error];
 556
 557 /* Here are some tests of bf_find_root, based on the find_root tests above.
 558    Use larger precision to catch any strangeness.
 559 */
 560 fpprec:32;
 561 32;
 562
 563 bf_find_root (2*x = -(log((4 + %e)/(2*%pi)))*(((4 + %e)/(2*%pi))^x), x, -1, 0);
 564 -3.3402898268741287760799570603459b-2;
 565
 566 bf_find_root (2*x = cos((%e + %pi)*x), x, 0, 1);
 567 1.9842105056568722553872075784746b-1;
 568
 569 /* verify that bf_find_root evaluates its first argument
 570  * (problem reported to mailing list 2007/06/07)
 571  */
 572 (expr : x^2 - 5, bf_find_root (expr, x, 0, 10));
 573 sqrt (5b0);
 574
 575 /* other tests for bf_find_root:
 576  * verify that bf_find_root expression is returned for non-numeric expression or bounds
 577  */
 578 (kill (a, b), bf_find_root (x^a - 5, x, 0, 10));
 579 bf_find_root (x^a - 5, x, 0b0, 10b0);
 580
 581 bf_find_root (x^3 - 5, x, a, b);
 582 bf_find_root (x^3 - 5, x, a, b);
 583
 584 /* verify that bf_find_root nested inside another function call is OK */
 585 quad_closeto(quad_qags (bf_find_root (x^a - 5, x, 0, 10), a, 1, 3),
 586              quad_qags (5^(1/a), a, 1, 3),
 587              [1d-15, 5d-15, 0, 0]);
 588 [true,true,true,true];
 589
 590 bf_find_root (bf_find_root (a^2 = x, a, 0, x) = 7, x, 0, 100); /* inner bf_find_root returns sqrt(x) */
 591 49b0;
 592
 593 /* verify that symbolic function name is OK */
 594 (foo (a) := 3^a - 5, bar : foo, bf_find_root (bar, 0, 10));
 595 log(5b0) / log(3b0);
 596
 597 /* example from mailing list 2006/12/01 */
 598 (expr : t = (297 * exp ((1000000 * t) / 33) - 330) / 10000000, bf_find_root (expr, t, 1e-9, 0.003));
 599 1.7549783076857196664805799825747b-5;
 600
 601 /* example from mailing list 2007/06/07 */
 602 (expr : 6096 * tan((2 * atan(c/(2 * fl))) / r) / (tan((1/60) * (%pi/180))),
 603  ev (bf_find_root (expr=6096, fl, 1, 10), c=7.176, r=3264));
 604 6.9814929328248795062474005396418b0;
 605
 606 /* adapted from mailing list 2007/01/13 */
 607
 608 (g (a) := bf_find_root (f (x, a), x, 0, 200),
 609  f (x, a) := x^a - 5,
 610  0);
 611 0;
 612
 613 g (0.5);
 614 25b0;
 615
 616 expr : g (z + z);
 617 bf_find_root (x^(2 * z) - 5, x, 0b0, 200b0);
 618
 619 ''(at (expr, z=0.25));
 620 25b0;
 621
 622 quad_closeto(quad_qags (g (z), z, 1, 3),
 623              quad_qags (5^(1/z), z, 1, 3),
 624              [1d-15, 5d-15, 0, 0]);
 625 [true,true,true,true];
 626
 627
 628 /* adapted from the reference manual */
 629
 630 (kill(f), f(x) := sin(x) - x/2, 0);
 631 0;
 632
 633 [bf_find_root (sin(x) - x/2, x, 0.1, %pi),
 634  bf_find_root (sin(x) = x/2, x, 0.1, %pi),
 635  bf_find_root (f(x), x, 0.1, %pi),
 636  bf_find_root (f, 0.1, %pi)];
 637 [1.8954942670339809471440357380936b0,
 638  1.8954942670339809471440357380936b0,
 639  1.8954942670339809471440357380936b0,
 640  1.8954942670339809471440357380936b0];
 641
 642 [bf_find_root (f, 1/(%pi*%e), 2*%pi*sin(%e)),
 643  bf_find_root (f, log(%pi), %e^%pi),
 644  bf_find_root (f, exp(1/5), exp(cos(%e + %pi))),
 645  bf_find_root (f, cos(exp(2))/10, 10*cos(exp(2)))];
 646 [1.8954942670339809471440357380936b0,
 647  1.8954942670339809471440357380936b0,
 648  1.8954942670339809471440357380936b0,
 649  1.8954942670339809471440357380936b0];
 650 /* adapted from the mailing list 2007/06/10
 651  * charfun2 copied from the interpol share package
 652  */
 653
 654 block ([expr],
 655  charfun2 (z, l1, l2) := charfun (l1 <= z and z < l2),
 656  expr : (-.329*x^3+.494*x^2 +.559*x+.117) *charfun2(x,minf,1.0)
 657     +(.215*x^3-1.94*x^2 +4.85*x-2.77) *charfun2(x,2.5,inf) +(.0933*x^3-1.02*x^2
 658     +2.56*x-.866) *charfun2(x,2.0,2.5) +(.0195*x^3-.581*x^2
 659     +1.67*x-.275) *charfun2(x,1.5,2.0) +(.00117*x^3-.498*x^2 +1.55*x -.213)
 660     *charfun2(x,1.0,1.5),
 661  block ([tru : 3.12727131060542283643481895355b0,
 662          est : bf_find_root (expr, x, 0, 4),
 663          err],
 664    err : abs(est-tru),
 665    if is(err < 2*10b0^(-fpprec)) then true else err));
 666 true;
 667
 668 /* bf_find_root example from sage mailing list */
 669
 670 bf_find_root (.05^(x + 1) = (x + 1)*10^(-10), x, 5, 100);
 671 6.0349925729832129297929340832397b0;
 672
 673 (find_root_error : true,
 674  errcatch (bf_find_root (1 + x^2, x, 0, 1)));
 675 [];
 676
 677 (find_root_error : "FOO",
 678  errcatch (bf_find_root (1 + x^2, x, 0, 1)));
 679 ["FOO"];
 680
 681 /* From bug 3010567.  Just checking that we don't get
 682    overflows when using bf_find_root */
 683 block([],
 684   F(x,y):=(log(x)/log(y))^x-x^(log(x)/log(y)),
 685   bf_find_root(F(400,z),z,2,1000));
 686 3.6541530643502285043078342270912b2;
 687
 688 reset (find_root_error);
 689 [find_root_error];
 690
 691 /* verify that SF bug #2564 "Regression in solve?" remains fixed */
 692
 693 block ([foo, bar],
 694  foo : [[x = 3,y = 2,z = 1],
 695  [x = .2768050783899193-2.987202528885064*%i,y = 1.478017834441328-1.347391287293138*%i,z = -0.526432162877356*%i-.8502171357296144],
 696  [x = -2.885476929518458*%i-.8209889702162483,y = 1.596034454560479*%i-1.205269272758513,z = .9957341762950345*%i+.09226835946330206],
 697  [x = 1.337215067329613-2.685489874065195*%i,y = 1.052864325754712*%i-1.700434271459228,z = .9324722294043555-.3612416661871523*%i],
 698  [x = -2.394051681840712*%i-1.807903909137758,y = .8914767115530776-1.790326582710134*%i,z = .7390089172206591-.6736956436465571*%i],
 699  [x = 2.217026751662001-2.021086930939692*%i,y = 1.864944458808694-.7224833323742995*%i,z = .9618256431728189*%i-.2736629900720828],
 700  [x = -1.579296488632072*%i-2.550651407188846,y = 1.923651286345638*%i-.5473259801441661,z = -.1837495178165701*%i-.9829730996839015],
 701  [x = 2.797416688213066-1.08372499856146*%i,y = .3674990356331407*%i-1.965946199367804,z = -.7980172272802396*%i-.6026346363792563],
 702  [x = -.5512485534497117*%i-2.948919299051704,y = 0.184536718926604-1.991468352590069*%i,z = .8951632913550623*%i+.4457383557765383],
 703  [x = .5512485534497115*%i-2.948919299051704,y = 1.991468352590069*%i+0.184536718926604,z = .4457383557765383-.8951632913550623*%i],
 704  [x = 1.083724998561459*%i+2.797416688213064,y = -.3674990356331408*%i-1.965946199367804,z = .7980172272802396*%i-.6026346363792563],
 705  [x = 1.57929648863207*%i-2.550651407188845,y = -1.923651286345638*%i-.5473259801441662,z = .1837495178165701*%i-.9829730996839015],
 706  [x = 2.021086930939673*%i+2.217026751661979,y = 0.722483332374306*%i+1.864944458808712,z = -.9618256431728189*%i-.2736629900720828],
 707  [x = 2.394051681840719*%i-1.80790390913777,y = 1.790326582710125*%i+.8914767115530766,z = .6736956436465571*%i+.7390089172206591],
 708  [x = 2.685489874065194*%i+1.337215067329613,y = -1.052864325754712*%i-1.700434271459228,z = .3612416661871523*%i+.9324722294043555],
 709  [x = 2.885476929518458*%i-0.820988970216246,y = -1.59603445456048*%i-1.205269272758512,z = .09226835946330206-.9957341762950345*%i],
 710  [x = 2.9872025288851*%i+.2768050783899063,y = 1.347391287293114*%i+1.478017834441318,z = 0.526432162877356*%i-.8502171357296144]],
 711  bar : sort (solve([x^2*y*z = 18,x*y^3*z = 24,x*y*z^4 = 6],[x,y,z])),
 712  /* some gyrations to round small floats to zero */
 713  matchdeclare (xx, lambda ([x], floatnump (x) and abs (x) < 1e-14)),
 714  defrule (rfoo, xx, 0),
 715  apply1 (abs (foo - bar), rfoo),
 716  apply ("and", map (is, flatten (%))));
 717 true;
 718
 719 /* SF bug #3102: "find_root(x,x,-1e300,1e300) => overflow" */
 720
 721 find_root(x,x,-1e300,1e300);
 722 0.0;
 723
 724 /* SF bug #3145: "solve sometimes gives wrong solution when using ratvars" */
 725
 726 kill(x, c, XXX);
 727 done;
 728
 729 is (solve([x^3+x+c],[x]) = ev(solve([x^3+x+c],[x]), ratvars:[XXX]));
 730 true;
 731
 732 is (solve([x^4+x+c],[x]) = ev(solve([x^4+x+c],[x]), ratvars:[XXX]));
 733 true;
 734
 735 (kill(a0, a1, a2, a3, x1, x2),
 736  eq:(-x^3*a3)+x^2*a2-x*a1+a0,
 737  ratvars:[x1,x2],
 738  map (lambda ([e], radcan (subst (e, eq))), solve (eq, x)));
 739 [0, 0, 0];
 740
 741 (eq:(-x^3*a3)+x^2*a2-x*a1+a0,
 742  ratvars:[x1,x2],
 743  solve(eq, x));
 744 [x = ((-1)/2-(sqrt(3)*%i)/2)*(sqrt(27*a0^2*a3^2+(4*a1^3-18*a0*a1*a2)*a3
 745                                                +4*a0*a2^3-a1^2*a2^2)
 746                              /(2*3^(3/2)*a3^2)
 747                              +((3*a0)/a3-(a1*a2)/a3^2)/6+a2^3/(27*a3^3))
 748                              ^(1/3)
 749    -(((sqrt(3)*%i)/2+(-1)/2)*(a1/(3*a3)+((-1)*a2^2)/(9*a3^2)))
 750     /(sqrt(27*a0^2*a3^2+(4*a1^3-18*a0*a1*a2)*a3+4*a0*a2^3-a1^2*a2^2)
 751      /(2*3^(3/2)*a3^2)
 752      +((3*a0)/a3-(a1*a2)/a3^2)/6+a2^3/(27*a3^3))
 753      ^(1/3)+a2/(3*a3),
 754  x = ((sqrt(3)*%i)/2+(-1)/2)*(sqrt(27*a0^2*a3^2+(4*a1^3-18*a0*a1*a2)*a3
 755                                                +4*a0*a2^3-a1^2*a2^2)
 756                              /(2*3^(3/2)*a3^2)
 757                              +((3*a0)/a3-(a1*a2)/a3^2)/6+a2^3/(27*a3^3))
 758                              ^(1/3)
 759    -(((-1)/2-(sqrt(3)*%i)/2)*(a1/(3*a3)+((-1)*a2^2)/(9*a3^2)))
 760     /(sqrt(27*a0^2*a3^2+(4*a1^3-18*a0*a1*a2)*a3+4*a0*a2^3-a1^2*a2^2)
 761      /(2*3^(3/2)*a3^2)
 762      +((3*a0)/a3-(a1*a2)/a3^2)/6+a2^3/(27*a3^3))
 763      ^(1/3)+a2/(3*a3),
 764  x = (sqrt(27*a0^2*a3^2+(4*a1^3-18*a0*a1*a2)*a3+4*a0*a2^3-a1^2*a2^2)
 765    /(2*3^(3/2)*a3^2)
 766    +((3*a0)/a3-(a1*a2)/a3^2)/6+a2^3/(27*a3^3))
 767    ^(1/3)
 768    -(a1/(3*a3)+((-1)*a2^2)/(9*a3^2))/(sqrt(
 769                                      27*a0^2*a3^2+(4*a1^3-18*a0*a1*a2)*a3
 770                                                  +4*a0*a2^3-a1^2*a2^2)
 771                                      /(2*3^(3/2)*a3^2)
 772                                      +((3*a0)/a3-(a1*a2)/a3^2)/6
 773                                      +a2^3/(27*a3^3))
 774                                      ^(1/3)+a2/(3*a3)]$
 775
 776 /* SF bug #3158: "triangularize gives incorrect result on a matrix containing %i" */
 777
 778 triangularize (matrix([%i,-1,0,0,1,0,0,0],
 779                       [1,%i,0,0,0,1,0,0],
 780                       [1,0,%i,-1,0,0,1,0],
 781                       [0,1,1,%i,0,0,0,1]));
 782 matrix([%i,-1,0,0,1,0,0,0],
 783        [0,1,-1,-%i,-1,0,%i,0],
 784        [0,0,2,2*%i,1,0,-%i,1],
 785        [0,0,0,0,2*%i,2,0,0]);
 786
 787 /* in this next example, one would hope that ev(..., algebraic) isn't necessary,
 788  * however, it is necessary due to ALGPCHK at line 468, src/rat3e.lisp, namely:
 789
 790      ((not $algebraic) nil)
 791
 792  * That line can't be removed because there are a couple of calls to ALGPGET
 793  * in src/nalgfa.lisp which appear to have potentially different behavior
 794  * if that line is removed.
 795  */
 796
 797 (kill(foo),
 798  tellrat(foo^2 + 1),
 799  ev (triangularize (subst (%i=foo, matrix([%i,-1,0,0,1,0,0,0],
 800                                           [1,%i,0,0,0,1,0,0],
 801                                           [1,0,%i,-1,0,0,1,0],
 802                                           [0,1,1,%i,0,0,0,1]))), algebraic));
 803 matrix([foo,-1,0,0,1,0,0,0],
 804        [0,1,-1,-foo,-1,0,foo,0],
 805        [0,0,2,2*foo,1,0,-foo,1],
 806        [0,0,0,0,2*foo,2,0,0]);
 807
 808 /* mailing list 2017-02-27: "Small diff bug?" */
 809
 810 (kill(I_Out, t),
 811  depends(I_Out,t),
 812  diff(I_Out,t),
 813  float(%%));
 814 'diff(I_Out, t, 1);
 815
 816 (kill(f, x, y),
 817  float (diff (f(x, y), x, 1, y, 2)));
 818 'diff (f(x, y), x, 1, y, 2);
 819
 820 float (diff (f(x + 3/2, %pi*y), x, 1, y, 2));
 821 'diff (f(x + 1.5, 3.141592653589793*y), x, 1, y, 2);
 822
 823 /* mailing list 2018-03-22: "bug in quad_qag" */
 824
 825 errcatch (quad_qag(r, r, 1/2, 1, 'epsrel=5d-8));
 826 [];
 827
 828 errcatch (quad_qag(r, r, 1/2, 1, 3, 'epsrel=5d-8));
 829 [[0.3750000000000001, 4.163336342344338e-15, 31, 0]];
 830