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/p/maxima/bugs/1960 */
 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/p/maxima/bugs/1657 */
 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 /*
 478  * Make sure bfallroots has the correct precision when given non-rational bfloat coefficients.
 479  * Use (x-sqrt(2))*(x^2-3) as the test since we know apriori what the roots should be.
 480  */
 481 fpprec:32;
 482 32;
 483 map(lambda([actual, expected],
 484     bfloat_approx_equal(rhs(actual), expected)),
 485   bfallroots(bfloat(expand((x-sqrt(2))*(x^2-3)))),
 486   bfloat([sqrt(2),sqrt(3),-sqrt(3)]));
 487 [true, true, true];
 488
 489 /* [ 940835 ] rectform fails with float/numer flags */
 490 rectform(log(-%i)),float;
 491 -0.5 * %i * %pi;
 492
 493 /* verify that exp(foo) evaluates to a number
 494  * probably should try several variations on this
 495  * adapted from sage mailing list
 496  */
 497
 498 first (quad_qags (sin (%pi * exp (x / 2)), x, 0, 2));
 499 - 0.4373454748252497;
 500
 501 /* verify that nested numerical integral is handled correctly
 502  * adapted from sage mailing list
 503  */
 504
 505 quad_qags (w^2 * quad_qags (1/(s - w), s, 1, 5) [1], w, -5, -1) [1];
 506 25.83639378805382;
 507
 508 /* find_root example from sage mailing list */
 509
 510 find_root (.05^(x + 1) = (x + 1)*10^(-10), x, 5, 100);
 511 6.034992572983213;
 512
 513 /* another nested example, collected from mma user forum */
 514
 515 (f : diff (1/(1 + (1 + (a - b)^2)), a),
 516  g : quad_qags (f*b*(1 - b)^2, b, 0, 1) [1],
 517  find_root (g = 0, a, 0, 1));
 518 0.3978613590133817;
 519
 520 /* a variation -- not sure what g:=... means in mma */
 521
 522 (f : diff (1/(1 + (1 + (a - b)^2)), a),
 523  g : 'quad_qags (f*b*(1 - b)^2, b, 0, 1) [1],
 524  find_root (g = 0, a, 0, 1));
 525 0.3978613590133817;
 526
 527 /* from mailing list 2009-02-18
 528  * "Re: [Maxima] I want to tell maxima (-1)^0.33333333333333=-1, what should i do?"
 529  * see also tests/rtest_plot
 530  */
 531
 532 (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),
 533  float_approx_equal_tolerance : 1e-12,
 534  0);
 535 0;
 536
 537 first (quad_qags (foo17 (u), u, -1, 0));
 538 - 0.359753467469551;
 539
 540 first (quad_qags (foo17, u, -1, 0));
 541 - 0.359753467469551;
 542
 543 find_root (foo17 (u) = -0.2, u, -1, 0);
 544 - 0.246809031968399;
 545
 546 (bar17 (u) := foo17 (u) + 0.2, find_root (bar17, u, -1, 0));
 547 - 0.246809031968399;
 548
 549 (compile (foo17), first (quad_qags (foo17, u, -1, 0)));
 550 - 0.359753467469551;
 551
 552 find_root (bar17, u, -1, 0);
 553 - 0.246809031968399;
 554
 555 /* SF bug # 2937837 "find_root_error documentation incorrect"
 556  */
 557
 558 (find_root_error : true,
 559  errcatch (find_root (1 + x^2, x, 0, 1)));
 560 [];
 561
 562 (find_root_error : "FOO",
 563  errcatch (find_root (1 + x^2, x, 0, 1)));
 564 ["FOO"];
 565
 566 reset (float_approx_equal_tolerance, find_root_error);
 567 [float_approx_equal_tolerance, find_root_error];
 568
 569 /* Here are some tests of bf_find_root, based on the find_root tests above.
 570    Use larger precision to catch any strangeness.
 571 */
 572 fpprec:32;
 573 32;
 574
 575 bf_find_root (2*x = -(log((4 + %e)/(2*%pi)))*(((4 + %e)/(2*%pi))^x), x, -1, 0);
 576 -3.3402898268741287760799570603459b-2;
 577
 578 bf_find_root (2*x = cos((%e + %pi)*x), x, 0, 1);
 579 1.9842105056568722553872075784746b-1;
 580
 581 /* verify that bf_find_root evaluates its first argument
 582  * (problem reported to mailing list 2007/06/07)
 583  */
 584 (expr : x^2 - 5, bf_find_root (expr, x, 0, 10));
 585 sqrt (5b0);
 586
 587 /* other tests for bf_find_root:
 588  * verify that bf_find_root expression is returned for non-numeric expression or bounds
 589  */
 590 (kill (a, b), bf_find_root (x^a - 5, x, 0, 10));
 591 bf_find_root (x^a - 5, x, 0b0, 10b0);
 592
 593 bf_find_root (x^3 - 5, x, a, b);
 594 bf_find_root (x^3 - 5, x, a, b);
 595
 596 /* verify that bf_find_root nested inside another function call is OK */
 597 quad_closeto(quad_qags (bf_find_root (x^a - 5, x, 0, 10), a, 1, 3),
 598              quad_qags (5^(1/a), a, 1, 3),
 599              [1d-15, 5d-15, 0, 0]);
 600 [true,true,true,true];
 601
 602 bf_find_root (bf_find_root (a^2 = x, a, 0, x) = 7, x, 0, 100); /* inner bf_find_root returns sqrt(x) */
 603 49b0;
 604
 605 /* verify that symbolic function name is OK */
 606 (foo (a) := 3^a - 5, bar : foo, bf_find_root (bar, 0, 10));
 607 log(5b0) / log(3b0);
 608
 609 /* example from mailing list 2006/12/01 */
 610 (expr : t = (297 * exp ((1000000 * t) / 33) - 330) / 10000000, bf_find_root (expr, t, 1e-9, 0.003));
 611 1.7549783076857196664805799825747b-5;
 612
 613 /* example from mailing list 2007/06/07 */
 614 (expr : 6096 * tan((2 * atan(c/(2 * fl))) / r) / (tan((1/60) * (%pi/180))),
 615  ev (bf_find_root (expr=6096, fl, 1, 10), c=7.176, r=3264));
 616 6.9814929328248795062474005396418b0;
 617
 618 /* adapted from mailing list 2007/01/13 */
 619
 620 (g (a) := bf_find_root (f (x, a), x, 0, 200),
 621  f (x, a) := x^a - 5,
 622  0);
 623 0;
 624
 625 g (0.5);
 626 25b0;
 627
 628 expr : g (z + z);
 629 bf_find_root (x^(2 * z) - 5, x, 0b0, 200b0);
 630
 631 ''(at (expr, z=0.25));
 632 25b0;
 633
 634 quad_closeto(quad_qags (g (z), z, 1, 3),
 635              quad_qags (5^(1/z), z, 1, 3),
 636              [1d-15, 5d-15, 0, 0]);
 637 [true,true,true,true];
 638
 639
 640 /* adapted from the reference manual */
 641
 642 (kill(f), f(x) := sin(x) - x/2, 0);
 643 0;
 644
 645 [bf_find_root (sin(x) - x/2, x, 0.1, %pi),
 646  bf_find_root (sin(x) = x/2, x, 0.1, %pi),
 647  bf_find_root (f(x), x, 0.1, %pi),
 648  bf_find_root (f, 0.1, %pi)];
 649 [1.8954942670339809471440357380936b0,
 650  1.8954942670339809471440357380936b0,
 651  1.8954942670339809471440357380936b0,
 652  1.8954942670339809471440357380936b0];
 653
 654 [bf_find_root (f, 1/(%pi*%e), 2*%pi*sin(%e)),
 655  bf_find_root (f, log(%pi), %e^%pi),
 656  bf_find_root (f, exp(1/5), exp(cos(%e + %pi))),
 657  bf_find_root (f, cos(exp(2))/10, 10*cos(exp(2)))];
 658 [1.8954942670339809471440357380936b0,
 659  1.8954942670339809471440357380936b0,
 660  1.8954942670339809471440357380936b0,
 661  1.8954942670339809471440357380936b0];
 662 /* adapted from the mailing list 2007/06/10
 663  * charfun2 copied from the interpol share package
 664  */
 665
 666 block ([expr],
 667  charfun2 (z, l1, l2) := charfun (l1 <= z and z < l2),
 668  expr : (-.329*x^3+.494*x^2 +.559*x+.117) *charfun2(x,minf,1.0)
 669     +(.215*x^3-1.94*x^2 +4.85*x-2.77) *charfun2(x,2.5,inf) +(.0933*x^3-1.02*x^2
 670     +2.56*x-.866) *charfun2(x,2.0,2.5) +(.0195*x^3-.581*x^2
 671     +1.67*x-.275) *charfun2(x,1.5,2.0) +(.00117*x^3-.498*x^2 +1.55*x -.213)
 672     *charfun2(x,1.0,1.5),
 673  block ([tru : 3.12727131060542283643481895355b0,
 674          est : bf_find_root (expr, x, 0, 4),
 675          err],
 676    err : abs(est-tru),
 677    if is(err < 2*10b0^(-fpprec)) then true else err));
 678 true;
 679
 680 /* bf_find_root example from sage mailing list */
 681
 682 bf_find_root (.05^(x + 1) = (x + 1)*10^(-10), x, 5, 100);
 683 6.0349925729832129297929340832397b0;
 684
 685 (find_root_error : true,
 686  errcatch (bf_find_root (1 + x^2, x, 0, 1)));
 687 [];
 688
 689 (find_root_error : "FOO",
 690  errcatch (bf_find_root (1 + x^2, x, 0, 1)));
 691 ["FOO"];
 692
 693 /* From bug 3010567.  Just checking that we don't get
 694    overflows when using bf_find_root */
 695 block([],
 696   F(x,y):=(log(x)/log(y))^x-x^(log(x)/log(y)),
 697   bf_find_root(F(400,z),z,2,1000));
 698 3.6541530643502285043078342270912b2;
 699
 700 reset (find_root_error);
 701 [find_root_error];
 702
 703 /* verify that SF bug #2564 "Regression in solve?" remains fixed */
 704
 705 block ([foo, bar],
 706  foo : [[x = 3,y = 2,z = 1],
 707  [x = .2768050783899193-2.987202528885064*%i,y = 1.478017834441328-1.347391287293138*%i,z = -0.526432162877356*%i-.8502171357296144],
 708  [x = -2.885476929518458*%i-.8209889702162483,y = 1.596034454560479*%i-1.205269272758513,z = .9957341762950345*%i+.09226835946330206],
 709  [x = 1.337215067329613-2.685489874065195*%i,y = 1.052864325754712*%i-1.700434271459228,z = .9324722294043555-.3612416661871523*%i],
 710  [x = -2.394051681840712*%i-1.807903909137758,y = .8914767115530776-1.790326582710134*%i,z = .7390089172206591-.6736956436465571*%i],
 711  [x = 2.217026751662001-2.021086930939692*%i,y = 1.864944458808694-.7224833323742995*%i,z = .9618256431728189*%i-.2736629900720828],
 712  [x = -1.579296488632072*%i-2.550651407188846,y = 1.923651286345638*%i-.5473259801441661,z = -.1837495178165701*%i-.9829730996839015],
 713  [x = 2.797416688213066-1.08372499856146*%i,y = .3674990356331407*%i-1.965946199367804,z = -.7980172272802396*%i-.6026346363792563],
 714  [x = -.5512485534497117*%i-2.948919299051704,y = 0.184536718926604-1.991468352590069*%i,z = .8951632913550623*%i+.4457383557765383],
 715  [x = .5512485534497115*%i-2.948919299051704,y = 1.991468352590069*%i+0.184536718926604,z = .4457383557765383-.8951632913550623*%i],
 716  [x = 1.083724998561459*%i+2.797416688213064,y = -.3674990356331408*%i-1.965946199367804,z = .7980172272802396*%i-.6026346363792563],
 717  [x = 1.57929648863207*%i-2.550651407188845,y = -1.923651286345638*%i-.5473259801441662,z = .1837495178165701*%i-.9829730996839015],
 718  [x = 2.021086930939673*%i+2.217026751661979,y = 0.722483332374306*%i+1.864944458808712,z = -.9618256431728189*%i-.2736629900720828],
 719  [x = 2.394051681840719*%i-1.80790390913777,y = 1.790326582710125*%i+.8914767115530766,z = .6736956436465571*%i+.7390089172206591],
 720  [x = 2.685489874065194*%i+1.337215067329613,y = -1.052864325754712*%i-1.700434271459228,z = .3612416661871523*%i+.9324722294043555],
 721  [x = 2.885476929518458*%i-0.820988970216246,y = -1.59603445456048*%i-1.205269272758512,z = .09226835946330206-.9957341762950345*%i],
 722  [x = 2.9872025288851*%i+.2768050783899063,y = 1.347391287293114*%i+1.478017834441318,z = 0.526432162877356*%i-.8502171357296144]],
 723  bar : sort (solve([x^2*y*z = 18,x*y^3*z = 24,x*y*z^4 = 6],[x,y,z])),
 724  /* some gyrations to round small floats to zero */
 725  matchdeclare (xx, lambda ([x], floatnump (x) and abs (x) < 1e-14)),
 726  defrule (rfoo, xx, 0),
 727  apply1 (abs (foo - bar), rfoo),
 728  apply ("and", map (is, flatten (%))));
 729 true;
 730
 731 /* SF bug #3102: "find_root(x,x,-1e300,1e300) => overflow" */
 732
 733 find_root(x,x,-1e300,1e300);
 734 0.0;
 735
 736 /* SF bug #3145: "solve sometimes gives wrong solution when using ratvars" */
 737
 738 kill(x, c, XXX);
 739 done;
 740
 741 is (solve([x^3+x+c],[x]) = ev(solve([x^3+x+c],[x]), ratvars:[XXX]));
 742 true;
 743
 744 is (solve([x^4+x+c],[x]) = ev(solve([x^4+x+c],[x]), ratvars:[XXX]));
 745 true;
 746
 747 (kill(a0, a1, a2, a3, x1, x2),
 748  eq:(-x^3*a3)+x^2*a2-x*a1+a0,
 749  ratvars:[x1,x2],
 750  map (lambda ([e], radcan (subst (e, eq))), solve (eq, x)));
 751 [0, 0, 0];
 752
 753 (eq:(-x^3*a3)+x^2*a2-x*a1+a0,
 754  ratvars:[x1,x2],
 755  solve(eq, x));
 756 [x = ((-1)/2-(sqrt(3)*%i)/2)*(sqrt(27*a0^2*a3^2+(4*a1^3-18*a0*a1*a2)*a3
 757                                                +4*a0*a2^3-a1^2*a2^2)
 758                              /(2*3^(3/2)*a3^2)
 759                              +((3*a0)/a3-(a1*a2)/a3^2)/6+a2^3/(27*a3^3))
 760                              ^(1/3)
 761    -(((sqrt(3)*%i)/2+(-1)/2)*(a1/(3*a3)+((-1)*a2^2)/(9*a3^2)))
 762     /(sqrt(27*a0^2*a3^2+(4*a1^3-18*a0*a1*a2)*a3+4*a0*a2^3-a1^2*a2^2)
 763      /(2*3^(3/2)*a3^2)
 764      +((3*a0)/a3-(a1*a2)/a3^2)/6+a2^3/(27*a3^3))
 765      ^(1/3)+a2/(3*a3),
 766  x = ((sqrt(3)*%i)/2+(-1)/2)*(sqrt(27*a0^2*a3^2+(4*a1^3-18*a0*a1*a2)*a3
 767                                                +4*a0*a2^3-a1^2*a2^2)
 768                              /(2*3^(3/2)*a3^2)
 769                              +((3*a0)/a3-(a1*a2)/a3^2)/6+a2^3/(27*a3^3))
 770                              ^(1/3)
 771    -(((-1)/2-(sqrt(3)*%i)/2)*(a1/(3*a3)+((-1)*a2^2)/(9*a3^2)))
 772     /(sqrt(27*a0^2*a3^2+(4*a1^3-18*a0*a1*a2)*a3+4*a0*a2^3-a1^2*a2^2)
 773      /(2*3^(3/2)*a3^2)
 774      +((3*a0)/a3-(a1*a2)/a3^2)/6+a2^3/(27*a3^3))
 775      ^(1/3)+a2/(3*a3),
 776  x = (sqrt(27*a0^2*a3^2+(4*a1^3-18*a0*a1*a2)*a3+4*a0*a2^3-a1^2*a2^2)
 777    /(2*3^(3/2)*a3^2)
 778    +((3*a0)/a3-(a1*a2)/a3^2)/6+a2^3/(27*a3^3))
 779    ^(1/3)
 780    -(a1/(3*a3)+((-1)*a2^2)/(9*a3^2))/(sqrt(
 781                                      27*a0^2*a3^2+(4*a1^3-18*a0*a1*a2)*a3
 782                                                  +4*a0*a2^3-a1^2*a2^2)
 783                                      /(2*3^(3/2)*a3^2)
 784                                      +((3*a0)/a3-(a1*a2)/a3^2)/6
 785                                      +a2^3/(27*a3^3))
 786                                      ^(1/3)+a2/(3*a3)]$
 787
 788 /* SF bug #3434: "kill(ratvars) should give an error" */
 789
 790 kill(all);
 791 done;
 792
 793 ratvars;
 794 [];
 795
 796 ratdisrep (rat(3));
 797 3;
 798
 799 errcatch (ratvars: 23);
 800 [];
 801
 802 ratdisrep (rat(3));
 803 3;
 804
 805 kill(ratvars);
 806 done;
 807
 808 ratdisrep (rat(3));
 809 3;
 810
 811 ratvars;
 812 [];
 813
 814 /* SF bug #3158: "triangularize gives incorrect result on a matrix containing %i" */
 815
 816 triangularize (matrix([%i,-1,0,0,1,0,0,0],
 817                       [1,%i,0,0,0,1,0,0],
 818                       [1,0,%i,-1,0,0,1,0],
 819                       [0,1,1,%i,0,0,0,1]));
 820 matrix([%i,-1,0,0,1,0,0,0],
 821        [0,1,-1,-%i,-1,0,%i,0],
 822        [0,0,2,2*%i,1,0,-%i,1],
 823        [0,0,0,0,2*%i,2,0,0]);
 824
 825 /* in this next example, one would hope that ev(..., algebraic) isn't necessary,
 826  * however, it is necessary due to ALGPCHK at line 468, src/rat3e.lisp, namely:
 827
 828      ((not $algebraic) nil)
 829
 830  * That line can't be removed because there are a couple of calls to ALGPGET
 831  * in src/nalgfa.lisp which appear to have potentially different behavior
 832  * if that line is removed.
 833  */
 834
 835 (kill(foo),
 836  tellrat(foo^2 + 1),
 837  ev (triangularize (subst (%i=foo, matrix([%i,-1,0,0,1,0,0,0],
 838                                           [1,%i,0,0,0,1,0,0],
 839                                           [1,0,%i,-1,0,0,1,0],
 840                                           [0,1,1,%i,0,0,0,1]))), algebraic));
 841 matrix([foo,-1,0,0,1,0,0,0],
 842        [0,1,-1,-foo,-1,0,foo,0],
 843        [0,0,2,2*foo,1,0,-foo,1],
 844        [0,0,0,0,2*foo,2,0,0]);
 845
 846 /* mailing list 2017-02-27: "Small diff bug?" */
 847
 848 (kill(I_Out, t),
 849  depends(I_Out,t),
 850  diff(I_Out,t),
 851  float(%%));
 852 'diff(I_Out, t, 1);
 853
 854 (kill(f, x, y),
 855  float (diff (f(x, y), x, 1, y, 2)));
 856 'diff (f(x, y), x, 1, y, 2);
 857
 858 float (diff (f(x + 3/2, %pi*y), x, 1, y, 2));
 859 'diff (f(x + 1.5, 3.141592653589793*y), x, 1, y, 2);
 860
 861 /* mailing list 2018-03-22: "bug in quad_qag" */
 862
 863 errcatch (quad_qag(r, r, 1/2, 1, 'epsrel=5d-8));
 864 [];
 865
 866 errcatch (quad_qag(r, r, 1/2, 1, 3, 'epsrel=5d-8));
 867 [[0.3750000000000001, 4.163336342344338e-15, 31, 0]];
 868
 869 /* attempt to verify that float_approx_equal works as advertised
 870  * some discussion on maxima-discuss 2019-01-14: "float_approx_equal too stringent by factor of 1/2?"
 871  */
 872
 873 (reset (float_approx_equal_tolerance),
 874  set_random_state (make_random_state (1234)));
 875 done;
 876
 877 /* following tests assume that foo + something*float_approx_equal_tolerance
 878  * is not equal to foo -- test that assumption
 879  */
 880
 881 block ([a: 1.0, b: 1.0 + 0.75*float_approx_equal_tolerance], [is(b = a), is(b - a > 0.0)]);
 882 [false, true];
 883
 884 block ([a: 1.0, b: 1.0 + 1.25*float_approx_equal_tolerance], [is(b = a), is(b - a > 0.0)]);
 885 [false, true];
 886
 887 sublist (makelist (block ([u: random(0.5) + 0.5, u1],
 888                           u1: u + 0.75*float_approx_equal_tolerance,
 889                           if not float_approx_equal (u, u1)
 890                             then failed_float_approx_equal (u, u1)
 891                             else true),
 892                     50),
 893          lambda ([e], e # true));
 894 [];
 895
 896 sublist (makelist (block ([u: random(0.5) + 0.5, u2],
 897                           u2: u + 1.25*float_approx_equal_tolerance,
 898                           if float_approx_equal (u, u2)
 899                             then failed_not_float_approx_equal (u, u2)
 900                             else false),
 901                     50),
 902          lambda ([e], e # false));
 903 [];
 904
 905 /* SF bug #3273: "quotient by polynomial of higher degree" */
 906
 907 (kill (a, b, c, n1, n2, n3, n4),
 908  reset (%rnum),
 909  solve ([(- a*n2) - a*n1 - 2*a, sqrt(a)*sqrt(- c)*(n2 - n1), b*n4 + n3 + b*n2 + b*n1 + b + 1], [n1,n2,n3,n4]));
 910 [[n1 = - 1, n2 = - 1, n3 = (1 - %r1)*b - 1, n4 = %r1]];
 911
 912 /* GENMATRIX is improved in order to save time when creating a large matrix.
 913  * it is related to SF bug #4056 "Cannot create very large data array"
 914  */
 915
 916 (kill(a), array(a,fixnum,2,2), a[1,1]:%e, a[2,2]:%pi, genmatrix(a,2,2));
 917 matrix([%e,0],[0,%pi]);
 918 genmatrix(lambda([i,j],j-i),3,3);
 919 matrix([0,1,2],[-1,0,1],[-2,-1,0]);
 920 (kill(B),genmatrix(B,2));
 921 matrix([B[1,1],B[1,2]],[B[2,1],B[2,2]]);
 922
 923 /* SF bug #4386: "quad_qagi doesn't recognize -inf as equivalent to minf" */
 924
 925 first( quad_qagi (1/u^3, u, -inf, -1, epsrel=1.0E-8, epsabs=2e-3, limit=200) );
 926 - 0.5;
 927
 928 first( quad_qagi (exp(-x^2), x, -inf, inf*inf, epsrel=1.0E-3, epsabs=2e-3, limit=100) );
 929 1.7724538509067376;
 930