Windows installer: update wxMaxima.
[maxima/cygwin.git] / tests / rtest3.mac
blob7041b5af4d7fba6ec394307694bab352162f5711
1 /*************** -*- Mode: MACSYMA; Package: MAXIMA -*-  ******************/
2 /***************************************************************************
3 ***                                                                    *****
4 ***     Copyright (c) 1984 by William Schelter,University of Texas     *****
5 ***     All rights reserved                                            *****
6 ***************************************************************************/
9 /* rtest3 */                                            
10 kill(all);
11 done;
12 for a from -3 step 7 thru 26 do ldisplay(a);
13 done$
14 s:0;
16 for i while i <= 10 do s:s+i;
17 done$
19 55$
20 series: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;
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$
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)$
82 x:3;
86 'x;
88 f(x):=x^2;
89 f(x):=x^2$
90 'f(2);
91 'f(2)$
92 ev(%,f);
94 '(f(2));
95 f(2)$
96 f(2);
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);
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$
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$
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];
177 /* end fix*/
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;
199 sin(5.0b-1);
200 4.794255386042030002732879b-1$
201 (reset (fpprec), 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$
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)$
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$
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$
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)$
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(%);
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)$
268 triginverses : all;
269 all;
271 sinh(acosh(x));
272 sqrt(x-1)*sqrt(x+1);
274 sinh(atanh(x));
275 x/(sqrt(1-x)*sqrt(x+1));
277 cosh(asinh(x));
278 sqrt(x^2+1);
280 cosh(atanh(x));
281 1/(sqrt(1-x)*sqrt(x+1));
283 tanh(asinh(x));
284 x/sqrt(x^2+1);
286 tanh(acosh(x));
287 sqrt(x-1)*sqrt(x+1)/x;
289 /* A few checks to see that triginverses false disables the above transformations */
290 triginverses: false;
291 false;
293 cos(acosh(x));
294 cos(acosh(x));
296 triginverses : all;
297 all;
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;
307 /* bug reported to mailing list 2009-05-09
308  * unexpected behavior in for loop with variable step
309  */
311 block ([L : []], for r:0 thru 7 step +2 do L : cons (r, L), L);
312 [6, 4, 2, 0];
314 block ([L : []], for r:7 thru 0 step -2 do L : cons (r, L), L);
315 [1, 3, 5, 7];
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];
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];
323 /* step is evaluated once at start of loop, so these loops are defined */
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];
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];
331 /* bug reported to mailing list 2009-05-13 "reset ( radexpand,  domain )"
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$
344 [radexpand, domain] : [all, complex];
345 [all, complex];
347 reset (radexpand, domain);
348 [radexpand, domain];
350 [radexpand, domain];
351 [true, real];
353 ([foo, bar, baz] : [1, 2, 3],
354  /* should ignore these non-defmvar's */
355  reset (foo, bar, baz));
358 /* verify that ORDFNA can handle CRE.
359  */
360 (kill (a, b), [doallmxops, doscmxops] : [false, false], 0);
363 b*matrix([rat(a)]);
364 b*matrix([''(rat(a))]);
366 (reset (doallmxops, doscmxops), 0);
369 /* SF bug #2936: stack overflow in integrate */
371 kill (x, A, B, MU, SIGMA);
372 done;
374 trigsimp (gamma_incomplete (1, log (x)));
375 gamma_incomplete (1, log (x));
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)$
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))$
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);
400 /* mailing list 2015-09-07: How can I catch this error? "errcatch" doesn't do the trick. */
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));
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;
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));
413 (kill (quux), bar () := (errcatch (integrate (exp(2^(3/2)*x^2 - sqrt(2)*x^2), x)), throw ('quux)), catch (foo ()));
414 quux;
416 /* result for this test will change if ever "quotient by zero" bug is fixed */
417 errcatch (taylor(coth(x), x, %i*%pi, 0));
420 (kill (mumble), bar () := (errcatch (taylor(coth(x), x, %i*%pi, 0)), throw ('mumble)), catch (foo ()));
421 mumble;
423 /* Verify that some special variables that were recently given DEFMVAR's are now resettable.
424  * If ever the default values of these variables are changed,
425  * some of the following tests will fail.
426  * In that case just update these tests to use the new default value.
427  */
429 (kill(mydefaults),
430  mydefaults['verbose] : false,
431  mydefaults['exptsubst] : false,
432  mydefaults['partswitch] : false,
433  mydefaults['inflag] : false,
434  mydefaults['derivsubst] : false,
435  mydefaults['opsubst] : true,
436  mydefaults['demoivre] : false,
437  mydefaults['nointegrate] : false,
438  mydefaults['tlimswitch] : true,
439  mydefaults['limsubst] : false,
440  mydefaults['abconvtest] : false,
441  mydefaults['packagefile] : false,
442  mydefaults['factlim] : 100000,
443  mydefaults['cflength] : 1,
444  mydefaults['taylordepth] : 3,
445  mydefaults['maxtaydiff] : 4,
446  mydefaults['lhospitallim] : 4,
447  mydefaults['linel] : 79,
448  0);
451 (reset (),
452  every (lambda ([v], is(mydefaults[v] = ev(v))), flatten (rest (arrayinfo (mydefaults), 2))));
453 true;
455 (myflags : '[verbose, exptsubst, partswitch, inflag, derivsubst, opsubst, demoivre,
456              nointegrate, tlimswitch, limsubst, abconvtest, packagefile],
457  /* "not" causes an extra evaluation of its argument,
458   * so this next line won't work as expected if "not" ever becomes
459   * simplifying operator; update this test case as needed if that
460   * ever comes to pass.
461   */
462  myflags :: map ("not", myflags),
463  ev (myflags));
464 [true, true, true, true, true, false, true, true, false, true, true, true];
466 map (reset, myflags);
467 [[verbose], [exptsubst], [partswitch], [inflag], [derivsubst], [opsubst], [demoivre],
468  [nointegrate], [tlimswitch], [limsubst], [abconvtest], [packagefile]];
470 ev (myflags);
471 [false, false, false, false, false, true, false, false, true, false, false, false];
473 (myvals : '[factlim, cflength, taylordepth, maxtaydiff, lhospitallim, linel],
474  myvals :: makelist (99, v, myvals),
475  ev (myvals));
476 [99, 99, 99, 99, 99, 99];
478 map (reset, myvals);
479 [[factlim], [cflength], [taylordepth], [maxtaydiff], [lhospitallim], [linel]];
481 ev (myvals);
482 [100000, 1, 3, 4, 4, 79];
484 declare_index_properties
485  (A, [],
486   B, [postsubscript],
487   C, [postsubscript, postsuperscript],
488   D, [postsubscript, postsuperscript, presuperscript],
489   E, [postsubscript, postsuperscript, presuperscript, presubscript],
490   F, [postsubscript, postsuperscript, presuperscript, presubscript,
491       postsubscript, postsuperscript, presuperscript, presubscript],
492   G, [postsubscript, postsuperscript, presuperscript, presubscript,
493       postsuperscript, presubscript]);
494 done;
496 map (get_index_properties, '[A, B, C, D, E, F, G]);
497 [[],
498  [postsubscript],
499  [postsubscript, postsuperscript],
500  [postsubscript, postsuperscript, presuperscript],
501  [postsubscript, postsuperscript, presuperscript, presubscript],
502  [postsubscript, postsuperscript, presuperscript, presubscript,
503   postsubscript, postsuperscript, presuperscript, presubscript],
504  [postsubscript, postsuperscript, presuperscript, presubscript,
505   postsuperscript, presubscript]];
507 kill (a, b, c, d, w, x, y, z);
508 done;
510 /* This business about capturing the output to a string and comparing that
511  * is a little fragile in the sense that changing invisible stuff (trailing
512  * spaces, newline character) can make these tests fail.
513  * But the pretty printer code changes very infrequently,
514  * and also I believe that these tests should work the same on Windows as
515  * on Unix-like platforms, because TERPRI (according to CLHS) just outputs #\Newline.
516  * Changing this file from NL to CR-NL line endings will presumably
517  * make these tests fail.
518  */
520 block ([display2d: true, S: make_string_output_stream ()],
521   with_stdout (S, ?terpri (), print (A[x])),
522   get_output_stream_string (S));
524 A  
528 block ([display2d: true, S: make_string_output_stream ()],
529   with_stdout (S, ?terpri (), print (B[x])),
530   get_output_stream_string (S));
532 B  
536 block ([display2d: true, S: make_string_output_stream ()],
537   with_stdout (S, ?terpri (), print (C[x, y])),
538   get_output_stream_string (S));
541 C  
545 block ([display2d: true, S: make_string_output_stream ()],
546   with_stdout (S, ?terpri (), print (D[x, y, z])),
547   get_output_stream_string (S));
549 z y
550  D  
551   x
554 block ([display2d: true, S: make_string_output_stream ()],
555   with_stdout (S, ?terpri (), print (E[w, x, y, z])),
556   get_output_stream_string (S));
558 y x
559  E  
560 z w
563 block ([display2d: true, S: make_string_output_stream ()],
564   with_stdout (S, ?terpri (), print (F[a, b, c, d, w, x, y, z])),
565   get_output_stream_string (S));
567 c, y b, x
568     F     
569 d, z a, w
572 block ([display2d: true, S: make_string_output_stream ()],
573   with_stdout (S, ?terpri (), print (G[a, b, c, d, w, x])),
574   get_output_stream_string (S));
576    c b, w
577     G     
578 d, x a
581 block ([display2d: true, S: make_string_output_stream ()],
582   with_stdout (S, ?terpri (), print (G[a, B[a], C[a, b], D[a, b, c], E[a, b, c, d], F[a, b, c, d, w, x, y, z]])),
583   get_output_stream_string (S));
585              b     c b
586             C  B ,  E
587              a  a  d a
588               G        
589 c b  c, y b, x a
590  D ,     F
591   a  d, z a, w
594 block ([display2d: true, S: make_string_output_stream ()],
595   with_stdout (S, ?terpri (), print (sqrt (G[a, b, c, d, w, x]))),
596   get_output_stream_string (S));
598         c b, w
599 sqrt(    G    ) 
600      d, x a
603 block ([display2d: true, S: make_string_output_stream ()],
604   with_stdout (S, ?terpri (), print ((1 - G[a, b, c, d, w, x])/E[1, 1/2, 2/3, 17/29])),
605   get_output_stream_string (S));
607        c b, w
608 1 -     G
609     d, x a
610 ------------- 
611     2/3 1/2
612        E
613   17/29 1
616 remove_index_properties (A, B, C, D, E, F, G);
617 done;
619 map (get_index_properties, '[A, B, C, D, E, F, G]);
620 [[], [], [], [], [], [], []];