tests/rtest3.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 /* rtest3 */
  10 kill(all);
  11 done;
  12 for a from -3 step 7 thru 26 do ldisplay(a);
  13 done$
  14 s:0;
  15 0$
  16 for i while i <= 10 do s:s+i;
  17 done$
  18 s;
  19 55$
  20 series:1;
  21 1$
  22 term:exp(sin(x));
  23 %e^sin(x)$
  24 for p unless p > 7 do
  25     (term:diff(term,x)/p,series:series+subst(x = 0,term)*x^p);
  26 done$
  27 series;
  28 x^7/90-x^6/240-x^5/15-x^4/8+x^2/2+x+1$
  29 poly:0;
  30 0$
  31 for i thru 5 do (for j from i step -1 thru 1 do poly:poly+i*x^j);
  32 done$
  33 poly;
  34 5*x^5+9*x^4+12*x^3+14*x^2+15*x$
  35 guess:-3.0;
  36 -3.0$
  37 for i thru 10 do
  38     (guess:subst(guess,x,0.5*(x+10/x)),
  39      if abs(guess^2-10) < 5.0e-5 then return(guess));
  40 -3.162280701754386;
  41 /* -3.1622806$ */
  42 for count from 2 next 3*count thru 20 do ldisplay(count);
  43 done$
  44 x:1000;
  45 1000$
  46 thru 10 while x # 0 do x:0.5*(x+5/x);
  47 done$
  48 x;
  49 2.282429035887867$
  50 remvalue(x);
  51 [x]$
  52 newton(f,guess):=block([numer,y],local(f,df,x,guess),numer:true,
  53        define(df(x),diff(f(x),x)),
  54        do (y:df(guess),if y = 0 then error("derivative at",guess,"is zero"),
  55            guess:guess-f(guess)/y,
  56            if abs(f(guess)) < 5.0e-6 then return(guess)));
  57 newton(f,guess):=block([numer,y],local(f,df,x,guess),numer:true,
  58        define(df(x),diff(f(x),x)),
  59        do (y:df(guess),if y = 0 then error("derivative at",guess,"is zero"),
  60            guess:guess-f(guess)/y,
  61            if abs(f(guess)) < 5.0e-6 then return(guess)))$
  62 sqr(x):=x^2-5.0;
  63 sqr(x):=x^2-5.0$
  64 newton(sqr,1000);
  65 2.236068027062195;
  66 for f in [log,rho,atan] do ldisp(f(1.0));
  67 done$
  68 ev(concat(e,linenum-1),numer);
  69 e10$
  70 kill(functions,values,arrays);
  71 done$
  72 done;
  73 done$
  74 exp:diff(x*f(x),x);
  75 x*'diff(f(x),x,1)+f(x)$
  76 f(x):=sin(x);
  77 f(x):=sin(x)$
  78 ev(exp,diff);
  79 sin(x)+x*cos(x)$
  80 x;
  81 x$
  82 x:3;
  83 3$
  84 x;
  85 3$
  86 'x;
  87 x$
  88 f(x):=x^2;
  89 f(x):=x^2$
  90 'f(2);
  91 'f(2)$
  92 ev(%,f);
  93 4$
  94 '(f(2));
  95 f(2)$
  96 f(2);
  97 4$
  98 sum(i!,i,1,4);
  99 33$
 100 'sum(i!,i,1,4);
 101 'sum(i!,i,1,4)$
 102 remvalue(x);
 103 [x]$
 104 'integrate(f(x),x,a,b);
 105 'integrate(x^2,x,a,b)$
 106 for i thru 5 do s:s+i^2;
 107 done$
 108 exp:s;
 109 s+55$
 110 ev(%,s:0);
 111 55$
 112 ev(exp);
 113 s+110$
 114 exp:'sum(g(i),i,0,n);
 115 'sum(g(i),i,0,n)$
 116 z*%e^z;
 117 z*%e^z$
 118 ev(%,z:x^2);
 119 x^2*%e^x^2$
 120 subst(x^2,z,exp);
 121 'sum(g(i),i,0,n)$
 122 a:%;
 123 'sum(g(i),i,0,n)$
 124 a+1;
 125 'sum(g(i),i,0,n)+1$
 126 kill(a,y);
 127 done$
 128 a;
 129 a$
 130 declare(integrate,noun);
 131 done$
 132 integrate(y^2,y);
 133 integrate(y^2,y)$
 134 ''integrate(y^2,y);
 135 y^3/3$
 136 f(y):=diff(y*log(y),y,2);
 137 f(y):=diff(y*log(y),y,2)$
 138 f(y):=1/y;
 139 f(y):=1/y$
 140 c10;
 141 c10$
 142 (x+y)^3;
 143 (y+x)^3$
 144 diff(%,x);
 145 3*(y+x)^2$
 146 y:x^2+1;
 147 x^2+1$
 148
 149 /* begin fix */
 150 kill(all);
 151 done;
 152  ev(%e^x*sin(x)^2,exponentialize);
 153  -%e^x*(%e^(%i*x)-%e^-(%i*x))^2/4;
 154   integrate(%,x);
 155 -((%e^((2*%i+1)*x)/(2*%i+1)+%e^((1-2*%i)*x)/(1-2*%i)-2*%e^x)/4);
 156  ev(%,demoivre);
 157  -((%e^x*(%i*sin(2*x)+cos(2*x))/(2*%i+1)
 158       +%e^x*(cos(2*x)-%i*sin(2*x))/(1-2*%i)-2*%e^x)
 159       /4);
 160  ans:ev(%,ratexpand);
 161  -%e^x*sin(2*x)/5-%e^x*cos(2*x)/10+%e^x/2;
 162  ev(ans,x:1,numer)-ev(ans,x:0,numer);
 163  0.5779160182042402;
 164  (fpprec : 35, 0);
 165  0;
 166  ev(ans,x:1,bfloat)-ev(ans,x:0,bfloat);
 167  5.7791601820424019599988308251707781b-1;
 168  integrate(%e^x*sin(x)^2,x);
 169  -(((2*%e^x*sin(2*x)+%e^x*cos(2*x)-5*%e^x)/10));
 170  trigreduce(%);
 171  -((2*%e^x*sin(2*x)+%e^x*cos(2*x)-5*%e^x)/10);
 172  % - ans,ratsimp;
 173  0 ;
 174  reset (fpprec);
 175  [fpprec];
 176
 177 /* end fix*/
 178
 179 ev(sin(x),%emode);
 180 sin(x)$
 181 sin(%pi/12)+tan(%pi/6);
 182 sin(%pi/12)+1/sqrt(3)$
 183 ev(%,numer);
 184 0.8361693142921465;
 185 /* tops 20 : 0.83616931$ */
 186 sin(1);
 187 sin(1)$
 188 ev(sin(1),numer);
 189 0.8414709848078965$
 190 beta(1/2,2/5);
 191 beta(1/2,2/5)$
 192 ev(%,numer);
 193 3.679093980405881;
 194 /* tops 20: 3.67909265$ */
 195 diff(atanh(sqrt(x)),x);
 196 1/(2*(1-x)*sqrt(x))$
 197 fpprec:25;
 198 25$
 199 sin(5.0b-1);
 200 4.794255386042030002732879b-1$
 201 (reset (fpprec), 0);
 202 0;
 203 /*begin fix */
 204  exp:cos(x)^2-sin(x)^2;
 205  cos(x)^2-sin(x)^2$
 206  ev(%,x:%pi/3);
 207  -1/2$
 208  diff(exp,x);
 209  -4*cos(x)*sin(x)$
 210  integrate(exp,x);
 211  (sin(2*x)/2+x)/2-(x-sin(2*x)/2)/2$
 212  expand(%);
 213  sin(2*x)/2$
 214  trigexpand(%);
 215  cos(x)*sin(x)$
 216  trigreduce(%);
 217  sin(2*x)/2$
 218  diff(%,x);
 219  cos(2*x)$
 220  %-exp,trigreduce,ratsimp;
 221   0;
 222 /*end fix*/
 223 sech(x)^2*sinh(x)*tanh(x)/coth(x)^2+cosh(x)^2*sech(x)^2*tanh(x)/coth(x)^2
 224                                    +sech(x)^2*tanh(x)/coth(x)^2;
 225 sech(x)^2*sinh(x)*tanh(x)/coth(x)^2+cosh(x)^2*sech(x)^2*tanh(x)/coth(x)^2
 226                                    +sech(x)^2*tanh(x)/coth(x)^2$
 227 trigsimp(%);
 228 (sinh(x)^5+sinh(x)^4+2*sinh(x)^3)/cosh(x)^5$
 229 /* These are from the trgsmp.dem file.
 230  * I (rtoy) hand-verified these results (using maxima, of course)
 231  */
 232 (1-sin(x)^2)*cos(x)/cos(x)^2+tan(x)*sec(x)^2;
 233 (1-sin(x)^2)*cos(x)/cos(x)^2+tan(x)*sec(x)^2$
 234 trigsimp(%);
 235 (sin(x)+cos(x)^4)/cos(x)^3$
 236
 237 tan(x)^2+sec(x)^2/(1-tan(x)*sec(x));
 238 tan(x)^2+sec(x)^2/(1-tan(x)*sec(x))$
 239 trigsimp(%);
 240 (sin(x)^4+sin(x)^3-1)/(cos(x)^2*sin(x)-cos(x)^4)$
 241
 242 (sin(x)^4-6*cos(x)^2*sin(x)^2+4*(cos(x)^2-sin(x)^2)+8*sin(x)+cos(x)^4+3)/(8*cos(x)^3);
 243 (sin(x)^4-6*cos(x)^2*sin(x)^2+4*(cos(x)^2-sin(x)^2)+8*sin(x)+cos(x)^4+3)/(8*cos(x)^3)$
 244 trigsimp(%);
 245 (sin(x)+cos(x)^4)/cos(x)^3$
 246
 247
 248 sech(x)^2*sinh(x)*tanh(x)/coth(x)^2+cosh(x)^2*sech(x)^2*tanh(x)/coth(x)^2+sech(x)^2*tanh(x)/coth(x)^2;
 249 sech(x)^2*sinh(x)*tanh(x)/coth(x)^2+cosh(x)^2*sech(x)^2*tanh(x)/coth(x)^2+sech(x)^2*tanh(x)/coth(x)^2$
 250 trigsimp(%);
 251 (sinh(x)^5+sinh(x)^4+2*sinh(x)^3)/cosh(x)^5$
 252
 253 -sech(x)^5*(sinh(x)^5+2*(sinh(x)^4+6*cosh(x)^2*sinh(x)^2+cosh(x)^4)-13*(sinh(x)^3+3*cosh(x)^2*sinh(x))+10*cosh(x)^2*sinh(x)^3-8*(sinh(x)^2+cosh(x)^2)+5*cosh(x)^4*sinh(x)+34*sinh(x)+6)/16;
 254 -sech(x)^5*(sinh(x)^5+2*(sinh(x)^4+6*cosh(x)^2*sinh(x)^2+cosh(x)^4)-13*(sinh(x)^3+3*cosh(x)^2*sinh(x))+10*cosh(x)^2*sinh(x)^3-8*(sinh(x)^2+cosh(x)^2)+5*cosh(x)^4*sinh(x)+34*sinh(x)+6)/16$
 255 trigsimp(%);
 256 -((sinh(x)^5+sinh(x)^4-2*sinh(x)^3)/cosh(x)^5)$
 257
 258 cos(x)*(sec(x)^2*tan(x)+1)-sec(x)^2*sin(x)-cos(x);
 259 cos(x)*(sec(x)^2*tan(x)+1)-sec(x)^2*sin(x)-cos(x)$
 260 trigsimp(%);
 261 0$
 262
 263 v*cos(x)*sec(x)^2*tan(x)+(-v*sec(x)^2-2*'diff(v,x))*sin(x)+'diff(v,x)*cos(x)*sec(x)+'diff(v,x,2)*cos(x);
 264 v*cos(x)*sec(x)^2*tan(x)+(-v*sec(x)^2-2*'diff(v,x))*sin(x)+'diff(v,x)*cos(x)*sec(x)+'diff(v,x,2)*cos(x)$
 265 trigsimp(%);
 266 -2*'diff(v,x,1)*sin(x)+'diff(v,x,2)*cos(x)+'diff(v,x,1)$
 267
 268 triginverses : all;
 269 all;
 270
 271 sinh(acosh(x));
 272 sqrt(x-1)*sqrt(x+1);
 273
 274 sinh(atanh(x));
 275 x/(sqrt(1-x)*sqrt(x+1));
 276
 277 cosh(asinh(x));
 278 sqrt(x^2+1);
 279
 280 cosh(atanh(x));
 281 1/(sqrt(1-x)*sqrt(x+1));
 282
 283 tanh(asinh(x));
 284 x/sqrt(x^2+1);
 285
 286 tanh(acosh(x));
 287 sqrt(x-1)*sqrt(x+1)/x;
 288
 289 /* A few checks to see that triginverses false disables the above transformations */
 290 triginverses: false;
 291 false;
 292
 293 cos(acosh(x));
 294 cos(acosh(x));
 295
 296 triginverses : all;
 297 all;
 298
 299 /* SF bug # 1981518, Calling desolve inside a "for...do" makes it loop endlessly
 300  * (protect against endless loop by throw--catch in case bug is triggered)
 301  */
 302 catch (block ([foo:1],
 303  for i thru 3 do (ilt (1/s, s, t),
 304  if foo > 3 then throw ('i = i) else foo : foo + 1)));
 305 done;
 306
 307 /* bug reported to mailing list 2009-05-09
 308  * unexpected behavior in for loop with variable step
 309  */
 310
 311 block ([L : []], for r:0 thru 7 step +2 do L : cons (r, L), L);
 312 [6, 4, 2, 0];
 313
 314 block ([L : []], for r:7 thru 0 step -2 do L : cons (r, L), L);
 315 [1, 3, 5, 7];
 316
 317 block ([L : [], r0 : 0, r1 : 7, s : +2], for r:r0 thru r1 step s do L : cons (r, L), L);
 318 [6, 4, 2, 0];
 319
 320 block ([L : [], r0 : 7, r1 : 0, s : -2], for r:r0 thru r1 step s do L : cons (r, L), L);
 321 [1, 3, 5, 7];
 322
 323 /* step is evaluated once at start of loop, so these loops are defined */
 324
 325 block ([L : [], s : +2], for i:1 thru 10 step s do L : cons (s : -s, L), L);
 326 [-2, 2, -2, 2, -2];
 327
 328 block ([L : [], s : -2], for i:10 thru 1 step s do L : cons (s : -s, L), L);
 329 [2, -2, 2, -2, 2];
 330
 331 /* bug reported to mailing list 2009-05-13 "reset ( radexpand,  domain )"
 332  *
 333  * display2d is a resetable option variable. We save the value of display2d
 334  * and restore it after the reset. This allows to run the testsuite in both
 335  * display modes.
 336  */
 337 (save:display2d, done);
 338 done$
 339 (reset (), [radexpand, domain]);
 340 [true, real];
 341 (display2d:save, done);
 342 done$
 343
 344 [radexpand, domain] : [all, complex];
 345 [all, complex];
 346
 347 reset (radexpand, domain);
 348 [radexpand, domain];
 349
 350 [radexpand, domain];
 351 [true, real];
 352
 353 ([foo, bar, baz] : [1, 2, 3],
 354  /* should ignore these non-defmvar's */
 355  reset (foo, bar, baz));
 356 [];
 357
 358 /* verify that ORDFNA can handle CRE.
 359  */
 360 (kill (a, b), [doallmxops, doscmxops] : [false, false], 0);
 361 0;
 362
 363 b*matrix([rat(a)]);
 364 b*matrix([''(rat(a))]);
 365
 366 (reset (doallmxops, doscmxops), 0);
 367 0;
 368
 369 /* SF bug #2936: stack overflow in integrate */
 370
 371 kill (x, A, B, MU, SIGMA);
 372 done;
 373
 374 trigsimp (gamma_incomplete (1, log (x)));
 375 gamma_incomplete (1, log (x));
 376
 377 trigsimp ((%i*gamma_incomplete(1,(1-2*log(x))^2/4)*(1-2*log(x))^2)
 378            /(2*log(x)-1)^2);
 379 %i*gamma_incomplete(1,(4*log(x)^2-4*log(x)+1)/4)$
 380
 381 trigsimp (integrate (%e^((-log(x)^2)-1)*log(x),x));
 382 -(%e^-(3/4)*(2*gamma_incomplete(1,(4*log(x)^2-4*log(x)+1)/4)*abs(2*log(x)-1)
 383             +2*gamma_incomplete(1/2,(4*log(x)^2-4*log(x)+1)/4)*log(x)
 384             -gamma_incomplete(1/2,(4*log(x)^2-4*log(x)+1)/4)))
 385  /(4*abs(2*log(x)-1))$
 386
 387 /* throw away results of integrate, just make sure it runs without crashing */
 388 block ([foo, bar, ctxt],
 389   foo : exp(-(log(x) - MU)*(log(x) - MU)/(2*SIGMA*SIGMA))/(x*SIGMA*sqrt(2*%pi)),
 390   bar : (log(B) - log(x*SIGMA) + ((x-A)*(x-A)/(2*B*B) - (log(x) -MU)*(log(x) -MU)/(2*SIGMA*SIGMA))),
 391   [foo, bar] : subst ([A=2, MU=3], [foo, bar]),
 392   ctxt : newcontext (),
 393   assume (SIGMA > 0, B > 0),
 394   integrate (expand (foo*bar), x),
 395   integrate (expand (foo*bar), x, 2, inf),
 396   killcontext (ctxt),
 397   0);
 398 0;
 399
 400 /* mailing list 2015-09-07: How can I catch this error? "errcatch" doesn't do the trick. */
 401
 402 /* result for this test will change if ever "quotient by zero" bug is fixed */
 403 errcatch (integrate(ev(ratsimp(1/(x^(5/2)+3*x^(1/3))),algebraic),x));
 404 [];
 405
 406 (kill (foo, bar, baz), foo () := bar (), bar () := (errcatch (integrate(ev(ratsimp(1/(x^(5/2)+3*x^(1/3))),algebraic),x)), throw ('baz)), catch (foo ()));
 407 baz;
 408
 409 /* result for this test will change if ever "quotient by zero" bug is fixed */
 410 errcatch (integrate (exp(2^(3/2)*x^2 - sqrt(2)*x^2), x));
 411 [];
 412
 413 (kill (quux), bar () := (errcatch (integrate (exp(2^(3/2)*x^2 - sqrt(2)*x^2), x)), throw ('quux)), catch (foo ()));
 414 quux;
 415
 416 /* result for this test will change if ever "quotient by zero" bug is fixed */
 417 errcatch (taylor(coth(x), x, %i*%pi, 0));
 418 [];
 419
 420 (kill (mumble), bar () := (errcatch (taylor(coth(x), x, %i*%pi, 0)), throw ('mumble)), catch (foo ()));
 421 mumble;