tests/rtest15.mac

   1 /* This file contains tests added since April 2002 */
   2
   3 kill(all);
   4 done$
   5
   6 /* This file assumes fpprec has its default value of 16 */
   7 fpprec:16;
   8 16;
   9
  10 /* apropos function added 7 april 2002
  11 */
  12 apropos("tr_optimize_max_loop");
  13 [tr_optimize_max_loop]$
  14
  15 apropos(tr_optimize_max_loop);
  16 [tr_optimize_max_loop]$
  17
  18 /* test introduced 7 april 2002. bug fix incorporated 9 june 2002 */
  19 integrate(3^log(x),x);
  20 3^((1/log(3)+1)*log(x))/((1/log(3)+1)*log(3))$
  21
  22 /* wester[1995] problem 84
  23    Bug 541030 fixed 2006-02-27, rev 1.14 of src/sin.lisp */
  24 integrate(sqrt(x+1/x-2),x,0,1);
  25 4/3$
  26
  27 /* bug reported by kevin ellwood.  fixed mar 11 2004 */
  28 integrate(exp(-k*t)/sqrt(k*t),t);
  29 sqrt(%pi)*erf(sqrt(k*t))/k$
  30
  31 /* a bug in chebyf.  fixed 20 apr 2004. */
  32
  33 kill(t,k);
  34 done$
  35 integrate(sqrt(k*t)*t,t);
  36 2*t^2*sqrt(k*t)/5$
  37 integrate(sqrt(k*t)*t^(1/3),t);
  38 6*t^(4/3)*sqrt(k*t)/11$
  39 integrate(sqrt(k*t)/t^(3/2),t);
  40 sqrt(k*t)*log(t)/sqrt(t)$
  41 integrate(sqrt(k*t)/sqrt(t),t);
  42 sqrt(t)*sqrt(k*t)$
  43 kill(t,k);
  44 done$
  45
  46 /* lisp error observed by stavros macrakis (#956730).  fixed 2004-05-20. */
  47
  48 block([context:'foobar],
  49   assume(n+1<0),
  50   declare(n,integer),
  51   errcatch(integrate(x^n,x,0,inf)));
  52 []$
  53
  54 error;
  55 ["defint: integral is divergent."];
  56
  57 killcontext('foobar);
  58 done$
  59
  60 /* second thru tenth code added 9 june 2002 */
  61 l : [1, 2, 3, 4, 5, 6, 7, 8, 9, 10];
  62 [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]$
  63 first(l);
  64 1$
  65 second(l);
  66 2$
  67 third(l);
  68 3$
  69 fourth(l);
  70 4$
  71 fifth(l);
  72 5$
  73 sixth(l);
  74 6$
  75 seventh(l);
  76 7$
  77 eighth(l);
  78 8$
  79 ninth(l);
  80 9$
  81 tenth(l);
  82 10$
  83
  84 /* assoc tests */
  85 assoc(1,[x=2]);
  86 false$
  87 assoc(x,[x=2]);
  88 2$
  89 assoc(1,[x=2],foobar);
  90 foobar$
  91
  92 kill(foo, f, x, g, y, z, h, u, xx, yy, none);
  93 done;
  94 assoc (f(x), foo(g(x) = y, f(x) = z + 1, h(x) = 2*u));
  95 z + 1;
  96 assoc (yy, [xx = 111, yy = 222, yy = 333, yy = 444]);
  97 222;
  98 assoc (abc, [[x, 111], [y, 222], [z, 333]], none);
  99 none;
 100 assoc (abc, [[x, 111], [y, 222], [z, 333]]);
 101 false;
 102
 103 kill(bar, baz, ww, zz);
 104 done;
 105 assoc (f(x, y), foo(baz(g(x), y), baz(f(x, y), z + 1), baz(h(x), 2*u)));
 106 z + 1;
 107 assoc (yy, [xx^11, yy^22, yy^33, yy^44]);
 108 22;
 109 assoc (yy, [xx . yy, yy . zz, yy . ww, yy . xx]);
 110 zz;
 111 assoc (foo(baz(xx), yy), bar(), none);
 112 none;
 113 assoc (baz(foo(ww)), []);
 114 false;
 115
 116 /* some older versions of gcl had a bug that resulted in a failure of the */
 117 /* following. */
 118 primep(80630964769);
 119 true$
 120
 121
 122 /* maxima has a bug causing an incorrect answer for the second integral. */
 123 integrate(diff(f(x,y),x,1,y,1),x);
 124 'diff(f(x,y),y,1)$
 125 integrate(diff(f(x,y),x,1,y,1),y);
 126 'diff(f(x,y),x,1)$
 127 /* the same bug causes a bug with the second integral in this set. */
 128 h(x,y):=x*y*diff(f(x,y),x,1,y,1);
 129 h(x,y):=x*y*diff(f(x,y),x,1,y,1)$
 130 integrate(h(x,y),x);
 131 'integrate(x*'diff(f(x,y),x,1,y,1)*y,x)$
 132 integrate(h(x,y),y);
 133 'integrate(x*'diff(f(x,y),x,1,y,1)*y,y)$
 134
 135 /* trigrat example from manual - fixed june 2002 */
 136 trigrat(sin(3*a)/sin(a+%pi/3));
 137 sqrt(3)*sin(2*a)+cos(2*a)-1;
 138
 139 /* another trigrat from manual
 140  * call ratsimp because result is not simplified same as expected result ... strange
 141  */
 142
 143 (sin(a)^2*sin(3*b+3*a)^2/sin(b+a)^2
 144     -2*sin(a)*sin(3*a)*cos(b)*sin(b+a-%pi/3)*sin(3*b+3*a)
 145     /(sin(a-%pi/3)*sin(b+a))
 146     +sin(3*a)^2*sin(b+a-%pi/3)^2/sin(a-%pi/3)^2,
 147  result : trigrat (%%),
 148  expected :
 149      (- (sqrt(3)*sin(4*b + 4*a) - cos(4*b + 4*a)
 150       - 2*sqrt(3)*sin(4*b + 2*a) + 2*cos(4*b + 2*a)
 151       - 2*sqrt(3)*sin(2*b + 4*a) + 2*cos(2*b + 4*a)
 152       + 4*sqrt(3)*sin(2*b + 2*a) - 8*cos(2*b + 2*a) - 4*cos(2*b - 2*a)
 153       + sqrt(3)*sin(4*b) - cos(4*b) - 2*sqrt(3)*sin(2*b) + 10*cos(2*b)
 154       + sqrt(3)*sin(4*a) - cos(4*a) - 2*sqrt(3)*sin(2*a) + 10*cos(2*a)
 155       - 9)/4),
 156  ratsimp (result - expected));
 157 0;
 158
 159 /* verify that trigrat distributes over composite objects
 160  * resolves SF bug reports # 1732315 and 1562340.
 161  */
 162
 163 trigrat (matrix ([1, 2], [3, 4]));
 164 ''(matrix ([trigrat(1), trigrat(2)], [trigrat(3), trigrat(4)]));
 165
 166 trigrat ({1, 2, 3});
 167 ''({trigrat(1), trigrat(2), trigrat(3)});
 168
 169 trigrat ([1, 2, 3]);
 170 ''([trigrat(1), trigrat(2), trigrat(3)]);
 171
 172 (kill (a, b), trigrat ([a < b, a <= b, a = b, a # b, a >= b, a > b]));
 173 ''([trigrat(a) < trigrat(b),
 174     trigrat(a) <= trigrat(b),
 175     trigrat(a) = trigrat(b),
 176     trigrat(a) # trigrat(b),
 177     trigrat(a) >= trigrat(b),
 178     trigrat(a) > trigrat(b)]);
 179
 180 /* Bug ID: 742909 - trigrat(sin(x/2)) makes a mess
 181  * Bug ID: 2999635 - trigrat(sin(1)) makes mess
 182  */
 183 trigrat(sin(x/2));
 184 sin(x/2);
 185 trigrat(sin(1));
 186 sin(1);
 187
 188 /* -----------------------------------------------------------------------------
 189  * Bug ID: 3398047 - trigrat() causes an error
 190  * -------------------------------------------------------------------------- */
 191
 192 trigrat((cos(x)+(sqrt(3)+2)*sin(x)-4*sin(5*%pi/12)*cos(x-(7*%pi)/12))/cos(x));
 193 2;
 194
 195 /* compile() will fail with gcl if gcc not installed */
 196 f(x):=x+2;
 197 f(x):=x+2;
 198 f(2);
 199 4;
 200 compile(f);
 201 [f];
 202 f(2);
 203 4;
 204
 205 /* gcl 2.6.8 sets small floats to 0.0 in compiled code */
 206 (f():=1e-6,compile(f),is(f()>0));
 207 true;
 208
 209 /* some tests for lambda expressions.
 210   we test for both translate and compile because there are some compiler
 211   macros for translated code */
 212 define_variable(qwerty,1,fixnum);
 213 1;
 214 /* m-tlambda */
 215 f():=apply(lambda([u],u+qwerty),[1]);
 216 f():=apply(lambda([u],u+qwerty),[1]);
 217 f();
 218 2;
 219 translate(f);
 220 [f];
 221 f();
 222 2;
 223 compile(f);
 224 [f];
 225 f();
 226 2;
 227 /* m-tlambda&env */
 228 f(x):=apply(lambda([u],u+x),[1]);
 229 f(x):=apply(lambda([u],u+x),[1]);
 230 f(1);
 231 2;
 232 translate(f);
 233 [f];
 234 f(1);
 235 2;
 236 compile(f);
 237 [f];
 238 f(1);
 239 2;
 240 /* m-tlambda&env */
 241 f(x,qwerty):=apply(lambda([u],u+x+qwerty),[1]);
 242 f(x,qwerty):=apply(lambda([u],u+x+qwerty),[1]);
 243 f(1,0);
 244 2;
 245 translate(f);
 246 [f];
 247 f(1,0);
 248 2;
 249 compile(f);
 250 [f];
 251 f(1,0);
 252 2;
 253 /* m-tlambda&env (outer) and m-tlambda (inner) */
 254 f(x):=apply(lambda([u],x+apply(lambda([v],v+qwerty),[u])),[-1]);
 255 f(x):=apply(lambda([u],x+apply(lambda([v],v+qwerty),[u])),[-1]);
 256 f(2);
 257 2;
 258 translate(f);
 259 [f];
 260 f(2);
 261 2;
 262 compile(f);
 263 [f];
 264 f(2);
 265 2;
 266 /* m-tlambda& */
 267 f():=apply(lambda([u,[v]],[u+qwerty,v]),[0,2,3]);
 268 f():=apply(lambda([u,[v]],[u+qwerty,v]),[0,2,3]);
 269 f();
 270 [1,[2,3]];
 271 translate(f);
 272 [f];
 273 f();
 274 [1,[2,3]];
 275 compile(f);
 276 [f];
 277 f();
 278 [1,[2,3]];
 279 /* m-tlambda&env& */
 280 f(x):=apply(lambda([u,[v]],[u+x,v]),[0,2,3]);
 281 f(x):=apply(lambda([u,[v]],[u+x,v]),[0,2,3]);
 282 f(1);
 283 [1,[2,3]];
 284 translate(f);
 285 [f];
 286 f(1);
 287 [1,[2,3]];
 288 compile(f);
 289 [f];
 290 f(1);
 291 [1,[2,3]];
 292 /* m-tlambda&env (from the sum), the inner lambda currently remains
 293   untranslated.  this is really a test for fungen&env-for-mevalsumarg. */
 294 f(n):=sum(apply(lambda([x],i+x),[i]),i,1,n);
 295 f(n):=sum(apply(lambda([x],i+x),[i]),i,1,n);
 296 f(3);
 297 12;
 298 translate(f);
 299 [f];
 300 f(3);
 301 12;
 302 compile(f);
 303 [f];
 304 f(3);
 305 12;
 306 kill(qwerty,f);
 307 done;
 308 /* this should kill f, but doesn't which upsets subsequent tests
 309    redefining it then killing it does the right thing */
 310 f(x):=x+3;
 311 f(x):=x+3;
 312 kill(f);
 313 done;
 314
 315 /* trignometric and hyperbolic functions of complex arguments */
 316 sin(%i*z);
 317 %i*sinh(z);
 318 cos(%i*z);
 319 cosh(z);
 320 tan(%i*z);
 321 %i*tanh(z);
 322 csc(%i*z);
 323 -%i*csch(z);
 324 sec(%i*z);
 325 sech(z);
 326 cot(%i*z);
 327 -%i*coth(z);
 328 sinh(%i*z);
 329 %i*sin(z);
 330 cosh(%i*z);
 331 cos(z);
 332 tanh(%i*z);
 333 %i*tan(z);
 334 csch(%i*z);
 335 -%i*csc(z);
 336 sech(%i*z);
 337 sec(z);
 338 coth(%i*z);
 339 -%i*cot(z);
 340
 341 /* trigreduce(sinh(x)^n)) wrong for some cases.  fixed 4 oct 2003 */
 342 trigreduce(sin(x)^2);
 343 (1-cos(2*x))/2;
 344 trigreduce(sin(x)^3);
 345 (3*sin(x)-sin(3*x))/4;
 346 trigreduce(sin(x)^4);
 347 (cos(4*x)-4*cos(2*x)+3)/8;
 348 trigreduce(sin(x)^5);
 349 (sin(5*x)-5*sin(3*x)+10*sin(x))/16;
 350 trigreduce(cos(x)^2);
 351 (cos(2*x)+1)/2;
 352 trigreduce(cos(x)^3);
 353 (cos(3*x)+3*cos(x))/4;
 354 trigreduce(cos(x)^4);
 355 (cos(4*x)+4*cos(2*x)+3)/8;
 356 trigreduce(cos(x)^5);
 357 (cos(5*x)+5*cos(3*x)+10*cos(x))/16;
 358 trigreduce(sinh(x)^2);
 359 (cosh(2*x)-1)/2;
 360 trigreduce(sinh(x)^3);
 361 (sinh(3*x)-3*sinh(x))/4;
 362 trigreduce(sinh(x)^4);
 363 (cosh(4*x)-4*cosh(2*x)+3)/8;
 364 trigreduce(sinh(x)^5);
 365 (sinh(5*x)-5*sinh(3*x)+10*sinh(x))/16;
 366 trigreduce(cosh(x)^2);
 367 (cosh(2*x)+1)/2;
 368 trigreduce(cosh(x)^3);
 369 (cosh(3*x)+3*cosh(x))/4;
 370 trigreduce(cosh(x)^4);
 371 (cosh(4*x)+4*cosh(2*x)+3)/8;
 372 trigreduce(cosh(x)^5);
 373 (cosh(5*x)+5*cosh(3*x)+10*cosh(x))/16;
 374
 375 /* de moivre's theorem - abramowitz & stegun 4.3.48, 4.5.53 */
 376 expand(trigreduce(expand((cos(x)+%i*sin(x))^2)));
 377 %i*sin(2*x)+cos(2*x);
 378 expand(trigreduce(expand((cos(x)+%i*sin(x))^3)));
 379 %i*sin(3*x)+cos(3*x);
 380 expand(trigreduce(expand((cos(x)+%i*sin(x))^4)));
 381 %i*sin(4*x)+cos(4*x);
 382 expand(trigreduce(expand((cos(x)+%i*sin(x))^5)));
 383 %i*sin(5*x)+cos(5*x);
 384 expand(trigreduce(expand((cos(x)+%i*sin(x))^6)));
 385 %i*sin(6*x)+cos(6*x);
 386 expand(trigreduce(expand((cos(x)+%i*sin(x))^7)));
 387 %i*sin(7*x)+cos(7*x);
 388 expand(trigreduce(expand((cos(x)+%i*sin(x))^8)));
 389 %i*sin(8*x)+cos(8*x);
 390 expand(trigreduce(expand((cos(x)+%i*sin(x))^9)));
 391 %i*sin(9*x)+cos(9*x);
 392 expand(trigreduce(expand((cosh(x)+sinh(x))^2)));
 393 sinh(2*x)+cosh(2*x);
 394 expand(trigreduce(expand((cosh(x)+sinh(x))^3)));
 395 sinh(3*x)+cosh(3*x);
 396 expand(trigreduce(expand((cosh(x)+sinh(x))^4)));
 397 sinh(4*x)+cosh(4*x);
 398 expand(trigreduce(expand((cosh(x)+sinh(x))^5)));
 399 sinh(5*x)+cosh(5*x);
 400 expand(trigreduce(expand((cosh(x)+sinh(x))^6)));
 401 sinh(6*x)+cosh(6*x);
 402 expand(trigreduce(expand((cosh(x)+sinh(x))^7)));
 403 sinh(7*x)+cosh(7*x);
 404 expand(trigreduce(expand((cosh(x)+sinh(x))^8)));
 405 sinh(8*x)+cosh(8*x);
 406 expand(trigreduce(expand((cosh(x)+sinh(x))^9)));
 407 sinh(9*x)+cosh(9*x);
 408
 409 /* bug 835287 */
 410 solve('diff(y,x),'diff(y,x));
 411 ['diff(y,x,1)=0];
 412
 413 /* multiplicity bug in solve.lisp 1.3 - 13 may 2004 */
 414 eigenvalues(matrix([3,1,0,0],[-4,-1,0,0],[7,1,2,1],[-17,-6,-1,0]));
 415 [[1],[4]];
 416
 417 /* Bug #2183 */
 418 eigenvalues(matrix([0,0,0,1,0,1,1,1,0,1],[1,1,0,0,1,0,1,0,1,0],[1,1,1,1,0,0,0,1,0,1],[0,1,1,0,0,1,1,1,1,0],[0,0,1,0,1,0,1,0,1,0],[0,1,0,0,0,0,0,1,0,0],[1,1,0,0,0,0,0,0,1,1],[1,0,1,1,1,0,1,0,0,0],[0,1,1,0,0,1,0,1,0,1],[1,0,0,1,1,0,0,1,0,0]));
 419 [[0],[1]];
 420
 421 eigenvalues(matrix([0,0,0,0,0,0,0,0,0,0,0],[0,-5276*A/5,18*A,18*A,18*A,0,0,0,0,0,0],[0,18*A,-4648*A/5,18*A,18*A,18*A,0,0,0,0,0],[0,18*A,18*A,-804*A,18*A,18*A,18*A,0,0,0,0],[0,18*A,18*A,18*A,-3392*A/5,18*A,18*A,18*A,0,0,0],[0,0,18*A,18*A,18*A,-3392*A/5,18*A,18*A,18*A,0,0],[0,0,0,18*A,18*A,18*A,-3392*A/5,18*A,18*A,18*A,0],[0,0,0,0,18*A,18*A,18*A,-3392*A/5,18*A,18*A,18*A],[0,0,0,0,0,18*A,18*A,18*A,-804*A,18*A,18*A],[0,0,0,0,0,0,18*A,18*A,18*A,-4648*A/5,18*A],[0,0,0,0,0,0,0,18*A,18*A,18*A,-5276*A/5]));
 422 [[0],[1]];
 423
 424 is(equal(eigenvalues(matrix([2,0,-1,0,0,0,-1,0,0,0,0,0,0],[0,2,0,-1,-1,0,0,0,0,0,0,0,0],
 425     [-1,0,3,-1,0,0,0,0,0,0,0,0,-1],[0,-1,-1,2,0,0,0,0,0,0,0,0,0],
 426     [0,-1,0,0,1,0,0,0,0,0,0,0,0],[0,0,0,0,0,1,0,-1,0,0,0,0,0],
 427     [-1,0,0,0,0,0,2,0,0,0,0,-1,0],[0,0,0,0,0,-1,0,1,0,0,0,0,0],
 428     [0,0,0,0,0,0,0,0,0,0,0,0,0],[0,0,0,0,0,0,0,0,0,2,-1,-1,0],
 429     [0,0,0,0,0,0,0,0,0,-1,1,0,0],[0,0,0,0,0,0,-1,0,0,-1,0,2,0],
 430     [0,0,-1,0,0,0,0,0,0,0,0,0,1])),
 431 [[2,-(sqrt(5)-3)/2,(sqrt(5)+3)/2,-(sqrt(5)-5)/2,(sqrt(5)+5)/2,0],
 432 [1,1,1,1,1,3]]));
 433 true;
 434
 435 is(equal(eigenvectors(matrix([2,0,-1,0,0,0,-1,0,0,0,0,0,0],[0,2,0,-1,-1,0,0,0,0,0,0,0,0],
 436     [-1,0,3,-1,0,0,0,0,0,0,0,0,-1],[0,-1,-1,2,0,0,0,0,0,0,0,0,0],
 437     [0,-1,0,0,1,0,0,0,0,0,0,0,0],[0,0,0,0,0,1,0,-1,0,0,0,0,0],
 438     [-1,0,0,0,0,0,2,0,0,0,0,-1,0],[0,0,0,0,0,-1,0,1,0,0,0,0,0],
 439     [0,0,0,0,0,0,0,0,0,0,0,0,0],[0,0,0,0,0,0,0,0,0,2,-1,-1,0],
 440     [0,0,0,0,0,0,0,0,0,-1,1,0,0],[0,0,0,0,0,0,-1,0,0,-1,0,2,0],
 441     [0,0,-1,0,0,0,0,0,0,0,0,0,1])),
 442 [[[2,-(sqrt(5)-3)/2,(sqrt(5)+3)/2,-(sqrt(5)-5)/2,(sqrt(5)+5)/2,0],
 443         [1,1,1,1,1,3]],
 444        [[[0,0,0,0,0,1,0,-1,0,0,0,0,0]],
 445         [[1,-1,1,0,-(sqrt(5)+1)/2,0,(sqrt(5)-1)/2,0,0,-(sqrt(5)-1)/2,-1,0,
 446           (sqrt(5)+1)/2]],
 447         [[1,-1,1,0,(sqrt(5)-1)/2,0,-(sqrt(5)+1)/2,0,0,(sqrt(5)+1)/2,-1,0,
 448           -(sqrt(5)-1)/2]],
 449         [[1,1,1,sqrt(5)+1,-(sqrt(5)+3)/2,0,(sqrt(5)-3)/2,0,0,(sqrt(5)-3)/2,1,
 450           1-sqrt(5),-(sqrt(5)+3)/2]],
 451         [[1,1,1,1-sqrt(5),(sqrt(5)-3)/2,0,-(sqrt(5)+3)/2,0,0,-(sqrt(5)+3)/2,1,
 452           sqrt(5)+1,(sqrt(5)-3)/2]],
 453         [[1,1,1,1,1,0,1,0,0,1,1,1,1],[0,0,0,0,0,1,0,1,0,0,0,0,0],
 454         [0,0,0,0,0,0,0,0,1,0,0,0,0]]]]));
 455     true;
 456
 457 /* Bug 1369669 */
 458 trigexpand(csc(3*x));
 459 csc(x)^3*sec(x)^3/(3*csc(x)^2*sec(x)-sec(x)^3);
 460
 461 /* Bug 1369451 */
 462 trigexpand(sec(x+y));
 463 csc(x)*sec(x)*csc(y)*sec(y)/(csc(x)*csc(y)-sec(x)*sec(y));
 464
 465 /* Bug 1309377 */
 466 (fpred(x) := (prederror:'false, is(x <= -1)),0);
 467 0$
 468
 469 fpred(xyzzy);
 470 unknown$
 471 translate(fpred);
 472 [fpred]$
 473 fpred(xyzzy);
 474 unknown$
 475
 476 /* Bug 1281737 */
 477 limit(atan(x)/(1/exp(1)-exp(-(1+x)^2)),x,inf);
 478 %e*%pi/2;
 479
 480 /* Bug 1363411 */
 481 'sum(1+f(k),k,1,2),simpsum;
 482 f(2) + f(1) + 2;
 483
 484 /* Bug 1403415, fixed in clmacs.lisp rev 1.20 */
 485 float(0.0b0);
 486 0.0;
 487
 488 /* Bug 1404754, fixed in float.lisp rev 1.23 */
 489 1.0b0*x;
 490 1.0b0*x;
 491
 492 /* Bug 1374704.  See also 1418010 because result isn't simplified. */
 493 (assume(cos(x)>0),0);
 494 0;
 495 integrate(sin(x)/cos(x)^2,x,0,%pi/3);
 496 1$
 497 kill(all);
 498 done;
 499
 500 /*
 501  * Test compiled array access
 502  */
 503 use_fast_arrays:false;
 504 false;
 505
 506 (ar:make_array('fixnum,3,2),0);
 507 0;
 508 (aref(array,row,col):=array[row,col],0);
 509 0;
 510 ar[2,1]:1;
 511 1;
 512 aref(ar,2,1);
 513 1;
 514 compile(aref);
 515 [aref];
 516 aref(ar,2,1);
 517 1;
 518 /* A old-style maxima hashed array */
 519 har[2,1]:1;
 520 1;
 521 aref(har,2,1);
 522 1;
 523
 524 use_fast_arrays:true;
 525 true;
 526 har2[3,4]:12;
 527 12;
 528 aref(har2,3,4);
 529 12;
 530
 531 /*
 532  * This test from Fabrizio Caruso, maxima list, 2006/02/14
 533  *
 534  * We're basically trying to test that the compiled version of foo
 535  * returns a fixnum array, just like the interpreted one.
 536  */
 537 (foo(x) := block([r], r : make_array('fixnum,2), r[0]:x,r[1]:x+1, return(r)),0);
 538 0$
 539
 540 (expr : foo(7), arrayinfo(expr));
 541 [declared,1,[1]]$
 542
 543 compile(foo);
 544 [foo];
 545
 546 (expr : foo(7), arrayinfo(expr));
 547 [declared,1,[1]]$
 548
 549 use_fast_arrays:false;
 550 false;
 551
 552 /* Bug 1405931: integrate(log(x)/(x^2+1)^2,x,0,inf)
 553  */
 554 integrate(log(x)/(x^2+1)^2,x,0,inf);
 555 -%pi/4;
 556
 557 /*
 558  * Bug 941457: integrate(1/x^5,x,1,2^(1/78))
 559  * Fixed 2006/03/13.  Solution needs work, though.
 560  */
 561 integrate(1/x^5,x,1,2^(1/78));
 562 1/4-(1/8)*2^(37/39);
 563
 564 integrate(1/x^3, x, 1, inf);
 565 1/2;
 566
 567 /*
 568  * Bug 938235: integrate((1/2)*u^2-1/u^5,u,1,sqrt(2))
 569  * Fixed 2006/03/13.  Solution needs work, though.
 570  */
 571 integrate((1/2)*u^2-1/u^5,u,1,sqrt(2));
 572 sqrt(2)/3-17/48;
 573
 574 /*
 575  * These next two cases are regression tests.  Rev 1.11 of sin.lisp
 576  * caused different (much messier) results to be returned.
 577  */
 578 integrate(sin(2*x)/sec(x),x);
 579 -2*cos(x)^3/3;
 580
 581 integrate(-a*sin(2*x)/(a*sin(x)^2+b),x);
 582 -log(a*sin(x)^2+b);
 583
 584 /*
 585  * Bug 1471861: limit(abs(sin(x)/x),x,0);
 586  * Fixed 2006/05/05, limit.lisp, rev 1.20.
 587  */
 588 limit(abs(sin(x)/x),x,0);
 589 1;
 590
 591 /*
 592  * Bug 1482843: subscripted variable causes trouble in integrate
 593  *
 594  * Fixed in defint.lisp, rev 1.25, 2006/05/06.
 595  */
 596 integrate(exp(-(x-mu[1])^2),x,0,inf);
 597 sqrt(%pi)*erf(mu[1])/2+sqrt(%pi)/2;
 598
 599 /*
 600  * Bug 1487703: integrate((sqrt(x^4-6*x^2+1)-x^2+1)/(2*x),x) fails
 601  *
 602  * Fixed in sin.lisp, 2006/05/14.  We don't loop forever anymore.
 603  */
 604 integrate((sqrt(x^4-6*x^2+1)-x^2+1)/(2*x),x);
 605 ('integrate((sqrt(x^2+2*x-1)*sqrt(x^2-2*x-1))/x,x)+log(x)-x^2/2)/2;
 606
 607 /*
 608  * Bug mentioned in mailing-list, 2006-05-31, "problem with bigfloats"
 609  */
 610 exp(1.0b-1);
 611 1.105170918075648b0;
 612
 613 /* sqrt(0b0) loops endlessly (bug report # 1515703)
 614  * asin(1b0) triggers same bug (it eventually calls sqrt(0b0))
 615  * Problem was in FPROOT.
 616  */
 617
 618 sqrt(0b0);
 619 0b0;
 620
 621 /* 0b0^(1/n) is not routed through FPROOT. Test it anyway. */
 622 [0b0^(1/2), 0b0^(1/3), 0b0^(1/17)];
 623 [0b0, 0b0, 0b0];
 624
 625 asin(1b0);
 626 ''(bfloat(%pi/2));
 627
 628 asin(-1b0);
 629 ''(-bfloat(%pi/2));
 630
 631 /* li[2](1.0) stack overflow (bug report # 1514861)
 632  */
 633
 634 li[2](1.0);
 635 ''(ev (%pi^2/6, numer));
 636
 637 /* bug report [ 782046 ] limit(abs(x),x,0) fails
 638  * appears to be fixed in 5.9.3cvs -- verify
 639  */
 640
 641 limit (abs(x), x, 0);
 642 0;
 643
 644 /* Related to Bug [1044318] defint(1/(sin(x)^2+1),x,0,3*%pi)
 645  *
 646  * integrate(1/(sin(x)^2+1),x,0,n*2*%pi) should be
 647  * n*integrate(1/(sin(x)^2+1),x,0,2*%pi), for positive integer n.  We
 648  * were just returning the same value.
 649  */
 650 factor(integrate(1/(sin(x)^2+1),x,0,2*%pi));
 651 sqrt(2)*%pi;
 652
 653 integrate(1/(sin(x)^2+1),x,0,4*%pi);
 654 2*sqrt(2)*%pi;
 655
 656 integrate(1/(sin(x)^2+1),x,0,20*%pi);
 657 10*sqrt(2)*%pi;
 658
 659 /*
 660  * Bug  1547769:  integrate(sqrt(x^3/(2*a-x)),x,0,2*a);
 661  *
 662  * We don't get an internal error, and we should be able to evaluate
 663  * this.  defint.lisp, rev 1.27
 664  */
 665 (assume(a > 0), 0);
 666 0;
 667
 668 integrate(sqrt(x^3/(2*a-x)),x,0,2*a);
 669 3*%pi*a^2/2;
 670
 671 /*
 672  * Bug [1044318] defint(1/(sin(x)^2+1),x,0,3*%pi)
 673  *
 674  * But there are other bugs related to this one.
 675  */
 676 factor(integrate(1/(sin(x)^2+1),x,0,3*%pi));
 677 (3 / 2) * sqrt(2)*%pi;
 678
 679 /*
 680  * Bug [1504505] integrate( 1/(x^8-1),x,0,1/2) => internal error
 681  *
 682  * Fixed in residu.lisp, 1.5.
 683  *
 684  * The call to factor is to get rid of some numerical factors.
 685  */
 686 ratsimp(ev(logcontract((integrate(1/(x^8-1),x,0,1/2))),algebraic) +
 687 (sqrt(2)*log((2^(3/2)+5)/(2^(3/2)-5)/-1)
 688       +2*sqrt(2)*atan((sqrt(2)+2)/2)
 689       +2*sqrt(2)*atan((sqrt(2)-2)/2)+log(9)
 690       +4*atan(1/2))
 691  /16);
 692 0;
 693
 694 /*
 695  * Bug [ 1582625 ] integrate(t^2*log(t)/((t^2-1)*(t^4+1)), t, 0, 1) wrong?
 696  *
 697  * Fixed in defint.lisp
 698  *
 699  */
 700 integrate(t^2*log(t)/((t^2-1)*(t^4+1)), t, 0, 1);
 701 /* The following result has changed to an equivalent result,
 702  * because of a change in simp.lisp revision 14.11.2011.
 703  * -((sqrt(2)-2)*%pi^2/32);
 704  */
 705 (sqrt(2)-1)*%pi^2/2^(9/2)$
 706
 707 /*
 708  * Bug [ 1073338 ] integrate yields incorrect result on rational function
 709  *
 710  * Here's one of the tests listed in that bug report.
 711  */
 712 logcontract(ratsimp(factor(integrate (1/((x-3)^4+1), x,0,1)))),algebraic;
 713 (sqrt(2)*log((5*sqrt(2)-38)/(5*sqrt(2)+38)/-1)
 714  +2^(3/2)*atan((7*sqrt(2)+5)/73)+2^(3/2)*atan((7*sqrt(2)-5)/73))
 715
 716  /8$
 717 /*
 718  * Work around in residu.lisp, rev 1.9
 719  */
 720 factor(expand(sqrtdenest(integrate (1/((x-3)^4+1/2), x,0,1))))
 721 /* we factor the result and subtract it */
 722 -factor(-(2*atan((2^(13/4)+2^(5/2)+2^(3/4))/(-2^(13/4)+2^(3/4)-98))
 723  +2*atan((-2^(13/4)+2^(5/2)-2^(3/4))/(-2^(13/4)+2^(3/4)+98))
 724  +log((3*2^(9/4)-2^(3/4)+sqrt(2)+73)/33)
 725  -log((-3*2^(9/4)+2^(3/4)+sqrt(2)+73)/33))
 726  /2^(7/4));
 727 0;
 728
 729 /*
 730 (2*atan((sqrt(2)-4*2^(1/4)+8)/(49*2^(3/4)+sqrt(2)-8))
 731  -log((2^(3/4)+12*sqrt(2)+73*2^(1/4)-2)/(33*2^(1/4)))
 732  +log((2^(3/4)-12*sqrt(2)+73*2^(1/4)+2)/(33*2^(1/4)))
 733  -2*atan((sqrt(2)+4*2^(1/4)+8)/(-49*2^(3/4)+sqrt(2)-8)))
 734  /(2*2^(3/4))$
 735  */
 736
 737 /*
 738  * Bug [ 1607567 ] trigreduce([atan(sin(a)/cos(a))]) => [ atan(tan(a)) ] (FIX)
 739  */
 740 (assume(a>-%pi/2, a<%pi/2),0);
 741 0;
 742
 743 trigreduce(atan(sin(a)/cos(a)));
 744 a;
 745 trigreduce([atan(sin(a)/cos(a))]);
 746 [a];
 747
 748 (forget(a>-%pi/2, a<%pi/2),0);
 749 0;
 750
 751 trigreduce([sin(x)^2]);
 752 [(1-cos(2*x))/2];
 753
 754 /*
 755  * Bug [ 1620977 ] limit(5^n/(2^n*3^n),n,inf) is wrong
 756  */
 757 limit(5^n/(2^n*3^n),n,inf);
 758 0;
 759
 760 /*
 761  * Bug [ 1370433 ] trigsimp(sqrt(%i2)) != sqrt(trigsimp(%i2))
 762  */
 763 trigsimp(sqrt(2*(cos(x)^2-sin(x)^2)+2));
 764 ''(sqrt(trigsimp(2*(cos(x)^2-sin(x)^2)+2)));
 765
 766 /* Tests for COERCE-FLOAT-FUN via QUADPACK functions.
 767  * COERCE-FLOAT-FUN is also an important element in plotting.
 768  *
 769  * Following cases are tested here:
 770  *  EXPR is a symbol
 771  *    name of a Lisp function
 772  *    name of a Maxima function
 773  *    name of a Maxima macro
 774  *    a string which is the name of a Maxima operator (e.g., "!")
 775  *    name of a simplifying function
 776  *  EXPR is a Maxima lambda expression
 777  *  EXPR is a general Maxima expression with atomic variable
 778  *  EXPR is a general Maxima expression with subscripted variable
 779  *
 780  * The one case not tested: EXPR is the name of a DEFMSPEC function
 781  */
 782
 783 (kill(F1, F2, F3),
 784  F1(foo) := sin(foo),
 785  F2(bar) ::= sin(bar),
 786  F3 : lambda([baz], sin(baz)),
 787  postfix ("F4"),
 788  "F4"(zorg) := sin(zorg),
 789  0);
 790 0;
 791
 792 (tol : 1e-8,
 793  expected : quad_qags (sin(quux), quux, 0, %pi, 'epsrel=tol),
 794  [first (expected),
 795   is (second (expected) < tol * first (expected)),
 796   is (abs (first (expected) - 2.0) < tol * 2.0),
 797   last (expected)]);
 798 [2.0, true, true, 0];
 799
 800 quad_qags (F1, blurf, 0, %pi, 'epsrel=tol);
 801 ''(expected);
 802
 803 quad_qags (F2, mumble, 0, %pi, 'epsrel=tol);
 804 ''(expected);
 805
 806 quad_qags (F3, gronk, 0, %pi, 'epsrel=tol);
 807 ''(expected);
 808
 809 quad_qags ("F4", flopt, 0, %pi, 'epsrel=tol);
 810 ''(expected);
 811
 812 quad_qags (sin, foobar, 0, %pi, 'epsrel=tol);
 813 ''(expected);
 814
 815 (translate (F1), ?not(?null(?fboundp (F1))));
 816 true;
 817
 818 quad_qags (F1, glump, 0, %pi, 'epsrel=tol);
 819 ''(expected);
 820
 821 quad_qags (sin(y[1]), y[1], 0, %pi, 'epsrel=tol);
 822 ''(expected);
 823
 824 /*
 825  * [ 1654183 ] integrate(x^2 / (1+x^6)^(3/2),x);
 826  */
 827 integrate(x^2/(1+x^6)^(3/2),x);
 828 x^3/(3*sqrt(x^6+1));
 829
 830 integrate(x^2/(1+x^6)^(3/2),x,-1,1);
 831 sqrt(2)/3;
 832
 833 /*
 834  * A few more tests of CHEBYF that handles integrals of the form
 835  *
 836  * x^r1*(c1+c2*x^q)^r2
 837  */
 838 /* (r1-q+1)/q a positive integer */
 839 integrate(x^7/(1+x^4)^(3/2),x);
 840 sqrt(x^4+1)/2+1/(2*sqrt(x^4+1));
 841
 842 /* r2 an integer */
 843 integrate(sqrt(x)/(1+x^(3/2))^2,x);
 844 -2/(3*(x^(3/2)+1));
 845
 846 /* (r1-q+1)/q a negative integer */
 847 integrate(x^(-1)/(1+x^4)^(3/2),x);
 848 -log(sqrt(x^4+1)+1)/4+log(sqrt(x^4+1)-1)/4+1/(2*sqrt(x^4+1));
 849
 850 /* (r1-q+1)/q + r2 an integer */
 851 integrate(sqrt(x)*(1+x^2)^(1/4),x);
 852 log((x^2+1)^(1/4)/sqrt(x)+1)/8-log((x^2+1)^(1/4)/sqrt(x)-1)/8
 853                               +atan((x^2+1)^(1/4)/sqrt(x))/4
 854                               +(x^2+1)^(1/4)/(sqrt(x)*(2*(x^2+1)/x^2-2))$
 855
 856 /*
 857  * Bug [ 1552789 ] integrate(1/(sin(x)^2+1),x,1,1+%pi)
 858  *
 859  * This bug said it was slow.  That's no longer true, but the result
 860  * was wrong.
 861  *
 862  * This has been fixed for this particular case.
 863  */
 864 integrate(1/(sin(x)^2+1),x,1,1+%pi);
 865 %pi/sqrt(2);
 866
 867 /* A few more related tests.  I think these are right */
 868 integrate(1/(sin(x)^2+1),x,1,1+4*%pi);
 869 2*sqrt(2)*%pi;
 870
 871 /*
 872  * Bug [ 1690374 ] asin(1 / sqrt(2))
 873  *
 874  * We return %pi/4 now.
 875  */
 876 asin(1/sqrt(2)),%piargs;
 877 %pi/4;
 878 asin(-1/sqrt(2)),%piargs;
 879 -%pi/4;
 880 acos(1/sqrt(2)),%piargs;
 881 %pi/4;
 882 acos(-1/sqrt(2)),%piargs;
 883 3*%pi/4;
 884
 885 /* Bug [ 1045287 ] */
 886 floatnump(float(exp(exp(2))));
 887 true$
 888
 889 /*
 890  * Bug [ 1778796 ] integrate( (x^3+1)/(x^4 + 4*x + 1), 0, 1)
 891  *
 892  * Not really a bug, but maxima takes way to long to compute the
 893  * integral, and the result is extremely long.  The fix is it try
 894  * the antiderivative first before trying the keyhole contour.
 895  */
 896 integrate( (x^3+1)/(x^4 + 4*x + 1), x, 0, 1);
 897 log(6)/4;
 898
 899 /* [ 1884711 ] bug when adding fractions involving square roots */
 900 factor(sqrt(2)/6-2*sqrt(2)/6);
 901 -sqrt(2)/6;
 902
 903 sqrt(3)/12 - 5*sqrt(3)/12;
 904 -1/sqrt(3);
 905
 906 /* Bug reported to mailing list 2008-03-23
 907  * Bug causes a Lisp error in FREEL (eventually called by $DEFINT).
 908  */
 909
 910 block ([foo],
 911     apply (forget, facts ()),
 912     assume (equal (a, 0)),
 913     foo : integrate (cos(a*x)/(1 + x^2), x, 0, inf),
 914     forget (equal (a, 0)),
 915     foo);
 916 %e^-a * (%pi*%e^(2*a) + %pi)/4;
 917
 918 /* verify that verbified math functions are not evaluated to numbers
 919  * bug reported to mailing list 2007-11-19
 920  */
 921
 922 (kill (x, t),
 923  'integrate (sqrt (9*x^2 + 37), x, 0, 2),
 924  changevar (%%, 3*x = sqrt(37)*sinh(t), t, x),
 925  ev (%%, nouns));
 926 37*%e^-(2*asinh(6/sqrt(37)))
 927   *(%e^(4*asinh(6/sqrt(37)))+4*asinh(6/sqrt(37))
 928                               *%e^(2*asinh(6/sqrt(37)))-1) /24;
 929
 930 /* considering the verbification stuff in more detail:
 931  * (1) verify that foo(non-float), nouns => non-numeric
 932  * (2) verify that foo(non-float), numer => float
 933  * (3) verify that foo(float) => float
 934  */
 935
 936 /* ignore last few digits; not important in this context */
 937
 938 (kill (all), float_approx_equal_tolerance : 1e-12, 0);
 939 0;
 940
 941 /* LIST OF MATH FUNCTIONS FOR WHICH THERE ARE FLOAT FUNCTIONS
 942  * (1) FUNCTIONS OF 1 ARGUMENT
 943  */
 944 (F1 :
 945  [erf, airy_ai, airy_bi, airy_dai, airy_dbi, elliptic_ec,
 946   elliptic_kc, exp, log, factorial, gamma, acosh, acoth, acsch,
 947   asech, asinh, atanh, cosh, coth, csch, sech, sinh, tanh, sqrt,
 948   acos, acot, acsc, asec, asin, atan, cos, cot, csc, sec, sin, tan,
 949   li[1], li[2], li[3], psi[0], psi[1], psi[2]],
 950  Y1 : makelist (foo(1/7), foo, F1),
 951  ev (Y1, nouns));
 952 [erf(1/7), 'airy_ai(1/7), 'airy_bi(1/7), 'airy_dai(1/7),
 953  'airy_dbi(1/7), 'elliptic_ec(1/7), 'elliptic_kc(1/7), %e^(1/7),
 954  -log(7), (1/7)!, gamma(1/7), acosh(1/7), acoth(1/7), acsch(1/7),
 955  asech(1/7), asinh(1/7), atanh(1/7), cosh(1/7), coth(1/7), csch(1/7),
 956  sech(1/7), sinh(1/7), tanh(1/7), 1/sqrt(7), acos(1/7), acot(1/7),
 957  acsc(1/7), asec(1/7), asin(1/7), atan(1/7), cos(1/7), cot(1/7),
 958  csc(1/7), sec(1/7), sin(1/7), tan(1/7), -log(6/7), li[2](1/7),
 959  li[3](1/7), psi[0](1/7), psi[1](1/7), psi[2](1/7)];
 960
 961 ev (Y1, numer);
 962 [0.16010712672873, 0.31821739764818, 0.67928220872456,
 963 - 0.25544752010866, 0.45500004582164, 1.513096652187913,
 964 1.631906796078438, 1.153564994895108, - 1.945910149055313,
 965 0.93543756289255, 6.548062940247834, 1.427448757889531*%i,
 966 0.14384103622589 - 1.570796326794897*%i, 2.644120761058629,
 967 2.633915793849633, 0.1423756431678, 0.14384103622589,
 968 1.010221447322645, 7.047554385466551, 6.976247043798608,
 969 0.98988197355171, 0.14334354757246, 0.14189319376693,
 970 0.37796447300923, 1.427448757889531, 1.428899272190733,
 971 1.570796326794897 - 2.633915793849634*%i, 2.633915793849634*%i,
 972 0.14334756890537, 0.14189705460416, 0.98981326044662,
 973 6.952316038379697, 7.023866335396166, 1.010291577169605,
 974 0.14237172979226, 0.14383695943619, 0.15415067982726,
 975 0.1483117974987926, 0.1455231699304894, - 7.363980242224349,
 976 50.3574714369117, - 687.6815220686585];
 977
 978 makelist (foo(float(1/7)), foo, F1);
 979 [0.16010712672873, 0.31821739764818, 0.67928220872456,
 980 - 0.25544752010866, 0.45500004582164, 1.513096652187913,
 981 1.631906796078438, 1.153564994895108, - 1.945910149055313,
 982 0.93543756289255, 6.548062940247834, 1.427448757889531*%i,
 983 0.14384103622589 - 1.570796326794897*%i, 2.644120761058629,
 984 2.633915793849633, 0.1423756431678, 0.14384103622589,
 985 1.010221447322645, 7.047554385466551, 6.976247043798608,
 986 0.98988197355171, 0.14334354757246, 0.14189319376693,
 987 0.37796447300923, 1.427448757889531, 1.428899272190733,
 988 1.570796326794897 - 2.633915793849634*%i, 2.633915793849634*%i,
 989 0.14334756890537, 0.14189705460416, 0.98981326044662,
 990 6.952316038379697, 7.023866335396166, 1.010291577169605,
 991 0.14237172979226, 0.14383695943619, 0.15415067982726,
 992 0.1483117974987926, 0.1455231699304894, - 7.363980242224349,
 993 50.3574714369117, - 687.6815220686585];
 994
 995 /* (2) FUNCTIONS OF 2 ARGUMENTS
 996  */
 997 (F2 :
 998  [atan2, bessel_i, bessel_j, bessel_k, bessel_y,
 999 jacobi_cd, jacobi_cn,
1000   jacobi_cs, jacobi_dc, jacobi_dn, jacobi_ds, jacobi_nc, jacobi_nd,
1001   jacobi_ns, jacobi_sc, jacobi_sd, jacobi_sn, elliptic_e,
1002   elliptic_eu, elliptic_f, beta],
1003  Y2 : makelist (foo(1/7, 2/7), foo, F2),
1004  ev (Y2, nouns));
1005 [atan(1/2), 'bessel_i(1/7, 2/7), 'bessel_j(1/7, 2/7),
1006  'bessel_k(1/7, 2/7), 'bessel_y(1/7, 2/7),
1007  'jacobi_cd(1/7, 2/7), 'jacobi_cn(1/7, 2/7), 'jacobi_cs(1/7, 2/7),
1008  'jacobi_dc(1/7, 2/7), 'jacobi_dn(1/7, 2/7), 'jacobi_ds(1/7, 2/7),
1009  'jacobi_nc(1/7, 2/7), 'jacobi_nd(1/7, 2/7), 'jacobi_ns(1/7, 2/7),
1010  'jacobi_sc(1/7, 2/7), 'jacobi_sd(1/7, 2/7), 'jacobi_sn(1/7, 2/7),
1011  elliptic_e(1/7, 2/7), elliptic_eu(1/7, 2/7), elliptic_f(1/7, 2/7),
1012  beta(1/7, 2/7)];
1013
1014 ev (Y2, numer);
1015 [0.46364760900081, 0.82410033570153, 0.79518665715578,
1016 1.44325571710677, - 1.035878262667884, 0.99270611823031,
1017 0.98983293044762, 6.959141929885592, 1.007347473371774,
1018 0.99710570154659, 7.010274039905824, 1.010271500613521,
1019 1.002902699732754, 7.030622760487995, 0.14369587660018,
1020 0.14264777586547, 0.1422349106284, 0.14271875687,
1021 0.14258093144727, 0.1429957703483, 9.973633393356877];
1022
1023 makelist (foo (float (1/7), float (2/7)), foo, F2);
1024 [0.46364760900081, 0.82410033570153, 0.79518665715578,
1025 1.44325571710677, - 1.035878262667884, 0.99270611823031,
1026 0.98983293044762, 6.959141929885592, 1.007347473371774,
1027 0.99710570154659, 7.010274039905824, 1.010271500613521,
1028 1.002902699732754, 7.030622760487995, 0.14369587660018,
1029 0.14264777586547, 0.1422349106284, 0.14271875687,
1030 0.14258093144727, 0.1429957703483, 9.973633393356877];
1031
1032 (F2 :
1033  [inverse_jacobi_cd, inverse_jacobi_cn, inverse_jacobi_cs,
1034   inverse_jacobi_sc, inverse_jacobi_sd, inverse_jacobi_sn],
1035  Y2 : makelist (foo(1/7, 2/7), foo, F2),
1036  ev (Y2, nouns));
1037 ['inverse_jacobi_cd(1/7,2/7),'inverse_jacobi_cn(1/7,2/7),
1038  'inverse_jacobi_cs(1/7,2/7),'inverse_jacobi_sc(1/7,2/7),
1039  'inverse_jacobi_sd(1/7,2/7),'inverse_jacobi_sn(1/7,2/7)];
1040
1041 ev (Y2, numer);
1042 [1.562139916774403, 1.536246964132837, 1.537956342530595,
1043   .1420329085375161, .1430674391291362, 0.143487627597992];
1044
1045 makelist (foo (float (1/7), float (2/7)), foo, F2);
1046 [1.562139916774403, 1.536246964132837, 1.537956342530595,
1047   .1420329085375161, .1430674391291362, 0.143487627597992];
1048
1049 (kill (F2),
1050  Y2 : [inverse_jacobi_dc (8/7, 1/7),
1051   inverse_jacobi_dn (1/7, 8/7),
1052   inverse_jacobi_ds (1/7, 8/7),
1053   inverse_jacobi_nc (8/7, 9/7),
1054   inverse_jacobi_nd (8/7, 9/7),
1055   inverse_jacobi_ns (8/7, 9/7)],
1056  ev (%%, nouns));
1057 ['inverse_jacobi_dc(8/7,1/7),'inverse_jacobi_dn(1/7,8/7),
1058  'inverse_jacobi_ds(1/7,8/7),'inverse_jacobi_nc(8/7,9/7),
1059  'inverse_jacobi_nd(8/7,9/7),'inverse_jacobi_ns(8/7,9/7)];
1060
1061 ev (Y2, numer);
1062 [.5422366033071744, 1.9430926302695156045465694,
1063 1.969205044088928940255018036, 0.53608536161539688257016016,
1064 0.461019342719434980049828553, 1.7162052237576627569674777881];
1065
1066 [inverse_jacobi_dc (float (8/7), float (1/7)),
1067  inverse_jacobi_dn (float (1/7), float (8/7)),
1068  inverse_jacobi_ds (float (1/7), float (8/7)),
1069  inverse_jacobi_nc (float (8/7), float (9/7)),
1070  inverse_jacobi_nd (float (8/7), float (9/7)),
1071  inverse_jacobi_ns (float (8/7), float (9/7))];
1072 [.5422366033071744, 1.9430926302695156045465694,
1073 1.969205044088928940255018036, 0.53608536161539688257016016,
1074 0.461019342719434980049828553, 1.7162052237576627569674777881];
1075
1076 /* (3) FUNCTIONS OF 3 ARGUMENTS */
1077 (Y3 : elliptic_pi (1/2, 1/3, 1/4),
1078  ev (Y3, nouns));
1079 elliptic_pi (1/2, 1/3, 1/4);
1080
1081 ev (Y3, numer);
1082 .3411527670552568;
1083
1084 elliptic_pi (float (1/2), float (1/3), float (1/4));
1085 .3411527670552568;
1086
1087 /* (4) FUNCTIONS WHICH ARE DON'T SEEM TO
1088  * GENERALLY YIELD NOUNS WHICH CAN BE EVALUATED TO NUMBERS
1089  * AND WHICH ARE THEREFORE EXCLUDED FROM THESE TESTS
1090  *  hermite chebyshev_t chebyshev_u scaled_bessel_i scaled_bessel_i0
1091  *  stirling jacobi_p laguerre legendre_p legendre_q assoc_legendre_p
1092  *  assoc_legendre_q spherical_bessel_j spherical_bessel_y
1093  *  spherical_hankel1 spherical_hankel2 spherical_harmonic
1094  *  ultraspherical pochhammer zeta genfact gen_laguerre
1095  */
1096
1097 (reset (float_approx_equal_tolerance), 0);
1098 0;
1099
1100 /* disallow bigfloat conversion for floating point infinity and not-a-number.
1101  * SF bug [ 2013654 ] bfloat(NaN) => finite number
1102  */
1103
1104 (foo : ?most\-positive\-double\-float,
1105  bar : ?most\-negative\-double\-float,
1106  block ([ratepsilon : 0], [bfloat (foo), bfloat (bar)]));
1107 /* might need to replace the expected result w/ a rat-ified value to ensure equality ... */
1108 [1.797693134862316b308, - 1.797693134862316b308];
1109
1110 errcatch (bfloat (2*foo));
1111 [];
1112
1113 errcatch (bfloat (2*bar));
1114 [];
1115
1116 errcatch (bfloat (2*foo + 3*bar));
1117 [];
1118
1119 /* Test that sinh(-x) for large x doesn't crash */
1120 sinh(-100b0)+sinh(100b0);
1121 0b0;
1122
1123
1124 /* tests for push and pop */
1125 errcatch(push());
1126 []$
1127
1128 errcatch(push(x,5));
1129 []$
1130
1131 errcatch(pop(2014));
1132 []$
1133
1134 errcatch(push(2014));
1135 []$
1136
1137 errcatch(push(a,b,c));
1138 []$
1139
1140 errcatch(push(a,[1,2]));
1141 []$
1142
1143 errcatch(pop());
1144 []$
1145
1146 errcatch(pop(2014));
1147 []$
1148
1149 errcatch(pop([1,2]));
1150 []$
1151
1152 (l : [1,2],0);
1153 0$
1154
1155 errcatch(pop(l,12));
1156 []$
1157
1158 errcatch(push(x,l, 2014));
1159 []$
1160
1161 (remvalue(k),a[k] : [a], a[k+1] : [b],  push(x, a[k : k+1]), [k,a[k], a[k+1]]);
1162 [k+1,[b],[x,b]]$
1163
1164 (remvalue(k),0);
1165 0$
1166
1167 (l : [], push(1,l), push(2,l), push(3,l),l);
1168 [3,2,1]$
1169
1170 (l : [l : [1]], push(x,l),l);
1171 [x,[1]]$
1172
1173 (l : [5,6], push(l : [],l),l);
1174 [[]]$
1175
1176 (l : [1,2], push(x,l),l);
1177 [x,1,2]$
1178
1179 [pop(l), l];
1180 [x, [1,2]]$
1181
1182 (l : [false], [pop(l), l]);
1183 [false,[]]$
1184
1185 (l : [1,2], a : push(0,l), pop(l), [a,l]);
1186 [[0,1,2],[1,2]]$
1187
1188 (l : [1,2], push(l,l));
1189 [[1,2],1,2]$
1190
1191 (remvalue(a), a[%pi] : [0,1,[8]], 0);
1192 0$
1193
1194 (push(2014, a[%pi]), a[%pi]);
1195 [2014,0,1,[8]]$
1196
1197 (pop(a[%pi]), a[%pi]);
1198 [0,1,[8]]$
1199
1200 (l1 : [0], l2 : [0,1], l3 : [0,1,2],0);
1201 0$
1202
1203 map('pop, '[l1,l2,l3]);
1204 [0,0,0]$
1205
1206 [l1,l2,l3];
1207 [[],[1],[1,2]]$
1208
1209 (a : b, b : c, c : d, d : e, e : f,0);
1210 0$
1211
1212 (l : [a,b,c,d,e], push(x,l), pop(l),l);
1213 [b,c,d,e,f]$
1214
1215 [pop(l),pop(l),pop(l),pop(l),pop(l)];
1216 [b,c,d,e,f]$
1217
1218 (declare("[", symmetric), l : [3,2,1], push(4,l));
1219 [1,2,3,4]$
1220
1221 (remove("[", symmetric),0);
1222 0$
1223
1224 /* ensure that remove("[", symmetric) doesn't interfere with other stuff --
1225  * bug previously caused rtest_levin to fail if executed after this one.
1226  */
1227
1228 ?mopp (?mlist);
1229 true;
1230
1231 constantp ([1.0, 2.0]);
1232 true;
1233
1234 member ("[", props);
1235 false;
1236
1237 /* resume other tests */
1238
1239 (aa[1] : [1], push(aa[1],aa[1]), [aa[1],pop(aa[1]), aa[1]]);
1240 [[[1],1],[1],[1]]$
1241
1242
1243 (remvalue(a,l,l1,l2,l3,b,c,d,e,f),remarray(aa), remarray(a), 0);
1244 0$
1245
1246 /* end push and pop tests */
1247
1248 /* tests for cons and endcons--especially tests for on nonlists */
1249
1250 errcatch(cons(a,a));
1251 []$
1252
1253 errcatch(cons(false,false));
1254 []$
1255
1256 errcatch(cons(x,a^b));
1257 []$
1258
1259 cons(false,[]);
1260 [false]$
1261
1262 cons(a,[a]);
1263 [a,a]$
1264
1265 cons(a,[true,false]);
1266 [a,true,false]$
1267
1268 cons(b, cons(a,[b]));
1269 [b,a,b]$
1270
1271 (l : [1,2,3], ll : cons(x,l), l[3] : 42, ll);
1272 [x,1,2,42]$
1273
1274 (remvalue(l,ll),0);
1275 0$
1276
1277 cons(a,[[]]);
1278 [a,[]]$
1279
1280 cons(2,a[2014]);
1281 a[2,2014]$
1282
1283 cons(x,a+b);
1284 a+b+x$
1285
1286 cons(x,a*b);
1287 a*b*x$
1288
1289 block([inflag : true],cons(x,a/b));
1290 a*x/b$
1291
1292 block([inflag : false],errcatch(cons(x,a/b)));
1293 []$
1294
1295 block([inflag : true], cons(x,-a));
1296 -a*x$
1297
1298 block([inflag : false], cons(x,-a));
1299 x-a$
1300
1301 errcatch(endcons(a,a));
1302 []$
1303
1304 errcatch(endcons(false,false));
1305 []$
1306
1307 errcatch(endcons(x,a^b));
1308 []$
1309
1310 endcons(false,[]);
1311 [false]$
1312
1313 endcons(a,[a]);
1314 [a,a]$
1315
1316 endcons(a,[true,false]);
1317 [true,false,a]$
1318
1319 endcons(k,endcons(n,[u]));
1320 [u,n,k]$
1321
1322 (remvalue(l,ll),0);
1323 0$
1324
1325 endcons(a,[[]]);
1326 [[],a]$
1327
1328 endcons(2,a[2014]);
1329 a[2014,2]$
1330
1331 endcons(x,a+b);
1332 a+b+x$
1333
1334 endcons(x,a*b);
1335 a*b*x$
1336
1337 block([inflag : true],endcons(x,a/b));
1338 x*a/b$
1339
1340 block([inflag : false],errcatch(endcons(x,a/b)));
1341 []$
1342
1343 block([inflag : true], endcons(x,-a));
1344 -a*x$
1345
1346 block([inflag : false], endcons(x,-a));
1347 a-x$
1348
1349 block([partswitch : true], rest([]));
1350 end$
1351
1352 /* end of tests for cons and endcons--especially tests for on nonlists */
1353
1354 /* SF bug #2812: lambda doesn't work,but %lambda does work. */
1355
1356 (kill (x, lambda), freeof (x, 1 + x * lambda[2]));
1357 false;
1358
1359 /* original bug trigger for #2812 -- result probably dependent on simplification details so don't bother
1360  *
1361 (kill (X, X11, X12, X21, X22, alpha, c, delta, f, eta, b, d),
1362  X11 : alpha[1]*c[1]^2/3 + c[3]^3 - log(delta[1]) + f^5,
1363  X12 : exp( eta[1] + b*eta[2] ) - log(d),
1364  X21 : log( eta[1] + eta[2]/b ) + exp(d),
1365  X22 : log(alpha[2]*c[2]^5/5 - c[1]/5) + exp(delta[2]) - lambda[2],
1366  X : matrix ([X11, X12], [X21, X22]),
1367  eigenvalues (X));
1368  *
1369  */
1370
1371 /* SF bug #2913: trigrat crashes with variable name "e" */
1372
1373 trigrat(c^x * sin(x));
1374 %i*abs(c)^x*sin(x)*sin(atan2(0, c)*x) + abs(c)^x*sin(x)*cos(atan2(0, c)*x);
1375
1376 trigrat(d^x * sin(x));
1377 %i*abs(d)^x*sin(x)*sin(atan2(0, d)*x) + abs(d)^x*sin(x)*cos(atan2(0, d)*x);
1378
1379 trigrat(e^x * sin(x));
1380 %i*abs(e)^x*sin(x)*sin(atan2(0, e)*x) + abs(e)^x*sin(x)*cos(atan2(0, e)*x);
1381
1382 trigrat(f^x * sin(x));
1383 %i*abs(f)^x*sin(x)*sin(atan2(0, f)*x) + abs(f)^x*sin(x)*cos(atan2(0, f)*x);
1384
1385 (check_eigenvectors (M) := block ([vals, vecs, mults, eqns],
1386   [[vals, mults], vecs] : eigenvectors (M),
1387   eqns : [lsum (m, m, mults) = length (M),
1388           expand (apply ("*", map (lambda ([v, m], v^m), vals, mults)) = determinant (M)),
1389           length (vecs) = length (vals),
1390           every (lambda ([mults1, vecs1], length (vecs1) <= mults1), mults, vecs),
1391           every (lambda ([val, vecs1], every (lambda ([vec], expand (transpose (M . vec) = val * matrix (vec))), vecs1)), vals, vecs)],
1392   if apply ("and", eqns)
1393   then true
1394   else [M, eqns]),
1395  0);
1396 0;
1397
1398 /* SF bug #3008: "Eigenvectors are missing" */
1399
1400 eigenvectors (matrix([1, 1, 0], [0, 1, 0], [0, 0, 2]));
1401 [[[1, 2], [2, 1]], [[[1, 0, 0]], [[0, 0, 1]]]];
1402
1403 check_eigenvectors (matrix([1, 1, 0], [0, 1, 0], [0, 0, 2]));
1404 true;
1405
1406 /* SF bug #3085: "missed eigenvectors" */
1407
1408 eigenvectors (matrix([3,1,0,0,0], [1,3,0,0,0], [0,0,2,1,1], [0,0,1,2,1], [0,0,1,1,2]));
1409 [[[2, 1, 4], [1, 2, 2]],
1410  [[[1, - 1, 0, 0, 0]],
1411   [[0, 0, 1, 0, - 1], [0, 0, 0, 1, - 1]],
1412   [[1, 1, 0, 0, 0], [0, 0, 1, 1, 1]]]];
1413
1414 check_eigenvectors (matrix([3,1,0,0,0], [1,3,0,0,0], [0,0,2,1,1], [0,0,1,2,1], [0,0,1,1,2]));
1415 true;
1416
1417 /* SF bug #3149: "eigenvectors does not show eigenvectors with eigenvalue zero" */
1418
1419 eigenvectors (matrix([2,3,1],[4,4,2],[-4,-8,-2]));
1420 [[[2, 0], [2, 1]], [[[1, 1, - 3]], [[1, 0, - 2]]]];
1421
1422 check_eigenvectors (matrix([2,3,1],[4,4,2],[-4,-8,-2]));
1423 true;
1424
1425 /* mailing list 2016-05-13: "Eigenvectors problem from a teacher" */
1426
1427 eigenvectors (matrix ([-1,1,0,1], [1,-1,1,0], [0,1,-1,1], [1,0,1,-1]));
1428 [[[1, - 1, - 3], [1, 2, 1]],
1429  [[[1, 1, 1, 1]], [[1, 0, - 1, 0], [0, 1, 0, - 1]], [[1, - 1, 1, - 1]]]];
1430
1431 check_eigenvectors (matrix ([-1,1,0,1], [1,-1,1,0], [0,1,-1,1], [1,0,1,-1]));
1432 true;
1433
1434 /* examples from reference manual */
1435
1436 eigenvectors (matrix ([11, -1], [1, 7]));
1437 [[[9 - sqrt(3), sqrt(3) + 9], [1, 1]],
1438  [[[1, sqrt(3) + 2]], [[1, 2 - sqrt(3)]]]];
1439
1440 check_eigenvectors (matrix ([11, -1], [1, 7]));
1441 true;
1442
1443 eigenvectors (matrix ([0, 1, 0, 0], [0, 0, 0, 0], [0, 0, 2, 0], [0, 0, 0, 2]));
1444 [[[0, 2], [2, 2]],
1445  [[[1, 0, 0, 0]], [[0, 0, 1, 0], [0, 0, 0, 1]]]];
1446
1447 check_eigenvectors (matrix ([0, 1, 0, 0], [0, 0, 0, 0], [0, 0, 2, 0], [0, 0, 0, 2]));
1448 true;
1449
1450 /* Test changevar for summation */
1451 changevar(sum(a[N-1-r]*x^r, r, 0, N-1), s = N-1-r, s, r);
1452 sum(a[s]*x^(N-1-s), s, 0, N-1);
1453
1454 /* SF bug report #3170: "Error in simtran" */
1455
1456 (A:matrix([-2,1,1],[1,-2,1],[1,1,-2]),
1457  simtran(A))$
1458 [[[-3,0],[2,1]],
1459  [[[1/sqrt(2),0,-1/sqrt(2)],[0,1/sqrt(2),-1/sqrt(2)]],
1460   [[1/sqrt(3),1/sqrt(3),1/sqrt(3)]]]]$
1461
1462 leftmatrix . A . rightmatrix $
1463 matrix([-3,0,0],[0,-3,0],[0,0,0])$
1464
1465 /* following trigsimp tests comprise the examples in demo(trgsmp),
1466  * with expected output = output produced by Maxima 5.39.0.
1467  */
1468
1469 (check_points(a, b) :=
1470  /* keep check points away from multiples of pi/2 */
1471  (makelist (0.05 + random (float(%pi/2) - 0.1), 100),
1472   map (lambda ([x1], (random(4) - 2)*float (%pi/2) + x1), %%),
1473   lmax (map (lambda ([x1], ev (abs (a - b), x = x1)), %%)),
1474   if %% < 1e-8 then true else %%),
1475  0);
1476 0;
1477
1478 (foo:((1-sin(x)^2)*cos(x))/cos(x)^2+tan(x)*sec(x)^2,
1479  bar:trigsimp(%%));
1480 (sin(x)+cos(x)^4)/cos(x)^3;
1481
1482 check_points (foo, bar);
1483 true;
1484
1485 (foo:tan(x)^2+sec(x)^2/(1-tan(x)*sec(x)),
1486  bar:trigsimp(%%));
1487 (sin(x)^4+sin(x)^3-1)/(cos(x)^2*sin(x)-cos(x)^4);
1488
1489 check_points (foo, bar);
1490 true;
1491
1492 (foo:(sin(x)^4-6*cos(x)^2*sin(x)^2+4*(cos(x)^2-sin(x)^2)+8*sin(x)+cos(x)^4+3)
1493  /(8*cos(x)^3),
1494  bar:trigsimp(%%));
1495 (sin(x)+cos(x)^4)/cos(x)^3;
1496
1497 check_points (foo, bar);
1498 true;
1499
1500 (foo:(sech(x)^2*sinh(x)*tanh(x))/coth(x)^2
1501  +(cosh(x)^2*sech(x)^2*tanh(x))/coth(x)^2+(sech(x)^2*tanh(x))/coth(x)^2,
1502  bar:trigsimp(%%));
1503 (sinh(x)^5+sinh(x)^4+2*sinh(x)^3)/cosh(x)^5;
1504
1505 check_points (foo, bar);
1506 true;
1507
1508 (foo:((-sech(x)^5)*(sinh(x)^5+2*(sinh(x)^4+6*cosh(x)^2*sinh(x)^2+cosh(x)^4)
1509                              +(-13)*(sinh(x)^3+3*cosh(x)^2*sinh(x))
1510                              +10*cosh(x)^2*sinh(x)^3
1511                              +(-8)*(sinh(x)^2+cosh(x)^2)+5*cosh(x)^4*sinh(x)
1512                              +34*sinh(x)+6)) /16,
1513  bar:trigsimp(%%));
1514 -((sinh(x)^5+sinh(x)^4-2*sinh(x)^3)/cosh(x)^5);
1515
1516 check_points (foo, bar);
1517 true;
1518
1519 (foo:cos(x)*(sec(x)^2*tan(x)+1)-sec(x)^2*sin(x)-cos(x),
1520  bar:trigsimp(%%));
1521 0;
1522
1523 check_points (foo, bar);
1524 true;
1525
1526 (foo:v*cos(x)*sec(x)^2*tan(x)+((-v)*sec(x)^2-2*'diff(v,x))*sin(x)
1527                              +'diff(v,x)*cos(x)*sec(x)+'diff(v,x,2)*cos(x),
1528  bar:trigsimp(%%));
1529 (-2*'diff(v,x,1)*sin(x))+'diff(v,x,2)*cos(x)+'diff(v,x,1);
1530
1531 /* bug reported to mailing list 2017-06-03: "trigsimp() bug?" */
1532
1533 trigsimp (1 - cos(x)^2);
1534 sin(x)^2;
1535
1536 trigsimp (cos(x)^2 - 1);
1537 -sin(x)^2;
1538
1539 trigsimp (1 - sin(x)^2);
1540 cos(x)^2;
1541
1542 trigsimp (sin(x)^2 - 1);
1543 -cos(x)^2;
1544
1545 /* Bug 3375: missing eigenvectors
1546    Results are ugly.  Just check that no eigenvectors are []
1547  */
1548 ([[vals,mult],vects]:eigenvectors(matrix([0,20,50],[1/20,0,0],[0,1/10,0])),
1549  some(emptyp,vects));
1550 false;
1551
1552 ([[vals,mult],vects]:eigenvectors(matrix([0,1,0],[1,0,1],[2,0,1])),
1553  some(emptyp,vects));
1554 false;
1555
1556 /* Bug 1820: missing eigenvectors */
1557 ([[vals,mult],vects]:eigenvectors(matrix([-1.8890437,-0.3],[1,0])),
1558  some(emptyp,vects));
1559 false;
1560
1561 kill(vals,mult,vects);
1562 done;
1563
1564 (reset(use_fast_arrays),1);
1565 1;
1566
1567 /* SF bug #3399: "trigreduce gives malformed property list error" */
1568
1569 trigreduce (atan (tan (x + y + z)));
1570 x + y + z;
1571
1572 /* this is the original example from the bug report;
1573  * verify that it doesn't trigger an error, but do not inspect the result itself
1574  */
1575 block ([expr, A, B, a, b, x, y],
1576   expr:-(2*(B*atan(((A*b+A*a-B)*(cos(B*x)*sin(A*y)+sin(B*x)*cos(A*y)))/(sqrt(A^2*b^2-2*A*B*b-A^2*a^2+B^2)*(sin(B*x)*sin(A*y)-cos(B*x)*cos(A*y)-1)))
1577     +sqrt(A^2*b^2-2*A*B*b-A^2*a^2+B^2)*atan((cos(B*x)*sin(A*y)+sin(B*x)*cos(A*y))/(sin(B*x)*sin(A*y)-cos(B*x)*cos(A*y)-1))))
1578      /(A*B*sqrt(A^2*b^2-2*A*B*b-A^2*a^2+B^2)),
1579   errcatch (trigreduce(expr)),
1580   if %% = [] then 'failed('trigreduce(expr)) else true);
1581 true;
1582
1583 /* SF bug #3413: "false in definite integral of rational"
1584  * SF bug #2644: "integrate(1/(1+s^7),s,0,%pi) includes a 'false' term"
1585  * obviously we can update these next two results if someday Maxima can compute these integrals
1586  */
1587
1588 integrate(1/(x^8+x+1),x,0,1);
1589 'integrate(1/(x^8+x+1),x,0,1);
1590
1591 integrate(1/(1+s^7),s,0,%pi);
1592 'integrate(1/(1+s^7),s,0,%pi);