1/*************** -*- Mode: MACSYMA; Package: MAXIMA -*-  ******************/
2/***************************************************************************
3***                                                                    *****
4***     Copyright (c) 1984 by William Schelter,University of Texas     *****
5***     All rights reserved                                            *****
6***************************************************************************/
7
8
9kill(functions,arrays,values);
10done$
11use_fast_arrays:false;
12false;
13a[n]:=n*a[n-1];
14a[n]:=n*a[n-1]$
15a[0]:1;
161$
17a[5];
18120$
19a[n]:=n;
20a[n]:=n$
21a[6];
226$
23a[4];
2424$
25(use_fast_arrays:true,kill(a));
26done;
27lambda([x,y,z],x^2+y^2+z^2);
28lambda([x,y,z],x^2+y^2+z^2)$
29%(1,2,a);
30a^2+5$
311+2+a;
32a+3$
33exp:[x^2,y/3,-2];
34[x^2,y/3,-2]$
35%[1]*x;
36x^3$
37[a,exp,%];
38[a,[x^2,y/3,-2],x^3]$
39m:matrix([a,0],[b,1]);
40matrix([a,0],[b,1])$
41m^2;
42matrix([a^2,0],[b^2,1])$
43exp:m . m;
44matrix([a^2,0],[a*b+b,1])$
45m[1,1]*m;
46matrix([a^2,0],[a*b,a])$
47%-exp+1;
48matrix([1,1],[1-b,a])$
49m^^(-1);
50matrix([1/a,0],[-b/a,1])$
51[x,y] . m;
52matrix([b*y+a*x,y])$
53matrix([a,b,c],[d,e,f],[g,h,i]);
54matrix([a,b,c],[d,e,f],[g,h,i])$
55%^^2;
56matrix([c*g+b*d+a^2,c*h+b*e+a*b,c*i+b*f+a*c],
57       [f*g+d*e+a*d,f*h+e^2+b*d,f*i+e*f+c*d],
58       [g*i+d*h+a*g,h*i+e*h+b*g,i^2+f*h+c*g])$
59exp:x+1 = y^2;
60x+1 = y^2$
61x-1 = 2*y+1;
62x-1 = 2*y+1$
63exp+%;
642*x = y^2+2*y+1$
65exp/y;
66(x+1)/y = y$
671/%;
68y/(x+1) = 1/y$
69fib[n]:=if n = 1 or n = 2 then 1 else fib[n-1]+fib[n-2];
70fib[n]:=if n = 1 or n = 2 then 1 else fib[n-1]+fib[n-2]$
71fib[1]+fib[2];
722$
73fib[3];
742$
75fib[5];
765$
77eta(mu,nu):=if mu = nu then mu else (if mu > nu then mu-nu else mu+nu);
78eta(mu,nu):=if mu = nu then mu else (if mu > nu then mu-nu else mu+nu)$
79eta(5,6);
8011$
81eta(eta(7,7),eta(1,2));
824$
83if not 5 >= 2 and 6 <= 5 or 4+1 > 3 then a else b;
84a$
85kill(f);
86done$
87
88kill(x,y,z);
89done$
90determinant(hessian(x^3-3*a*x*y*z+y^3,[x,y,z]));
91-3*a*y*(9*a^2*x*z+18*a*y^2)-27*a^3*x*y*z-54*a^2*x^3$
92
93subst(1,z,quotient(%,-54*a^2));
94y^3+a*x*y+x^3$
95f(x):=block([a,y],local(a),y:4,a[y]:x,display(a[y]));
96f(x):=block([a,y],local(a),y:4,a[y]:x,display(a[y]))$
97y:2;
982$
99a[y+2]:0;
1000$
101f(9);
102done$
103a[y+2];
1040$
105
106(use_fast_arrays : false, kill(a), 0);
1070$
108
109/* ensure that matrix construction works as advertised */
110(L : makelist ([i], i, 1, 100), apply (matrix, L), [op (%%), args (%%)]);
111[matrix, ''(makelist ([i], i, 1, 100))];
112
113(L : makelist ([i], i, 1, 100), apply (matrix, L), transpose (%%));
114''(matrix (tree_reduce (append, L)));   /* call tree_reduce instead of append because GCL barfs ... */
115
116(matrix (), [op (%%), args (%%)]);
117[matrix, []];
118
119/* construct a matrix of modest size */
120(apply (matrix, makelist ([i], i, 1, 1000)), 0);
1210;
122
123/* construct a matrix of modest size */
124(apply (matrix, makelist ([i], i, 1, 10000)), 0);
1250;
126
127/* verify that arguments are evaluated exactly once */
128block ([a : b, b : c, c: d, d : 1], matrix ([a, b], [c, d]), [op (%%), args (%%)]);
129[matrix, '[[b, c], [d, 1]]];
130
131/* verify that arguments are evaluated exactly once */
132block ([a : b, b : c, c: d, d : 1, L1 : '[a, b], L2 : '[c, d]], matrix (L1, L2), [op (%%), args (%%)]);
133[matrix, '[[a, b], [c, d]]];
134
135/* another evaluation puzzle, derived from discussion on mailing list circa 2013-10-28 */
136
137(kill (q, x),
138 q : '[[x]],
139 x : 3,
140 apply (matrix, q));
141matrix ([x]);
142
143/* a more elaborate version of the preceding evaluation puzzle;
144 * result not checked for correctness
145 */
146
147(kill (all),
148 load (diag),
149 A : matrix ([a, 1], [1, 0]),
150 integer_pow(x) := block ([k], declare (k, integer), x^k),
151 mat_function (integer_pow, A));
152
153matrix([(sqrt(a^2+4)-a)^(k+1)*2^(-k-1)*(-1)^k
154         /((sqrt(a^2+4)+a)/2+(sqrt(a^2+4)-a)/2)
155         +(sqrt(a^2+4)+a)^(k+1)*2^(-k-1)
156          /((sqrt(a^2+4)+a)/2+(sqrt(a^2+4)-a)/2),
157        (sqrt(a^2+4)+a)^k/(((sqrt(a^2+4)+a)/2+(sqrt(a^2+4)-a)/2)*2^k)
158         -(sqrt(a^2+4)-a)^k*(-1)^k/(((sqrt(a^2+4)+a)/2+(sqrt(a^2+4)-a)/2)
159                                   *2^k)],
160       [(sqrt(a^2+4)-a)*(sqrt(a^2+4)+a)^(k+1)*2^(-k-2)
161         /((sqrt(a^2+4)+a)/2+(sqrt(a^2+4)-a)/2)
162         -(sqrt(a^2+4)-a)^(k+1)*(sqrt(a^2+4)+a)*2^(-k-2)*(-1)^k
163          /((sqrt(a^2+4)+a)/2+(sqrt(a^2+4)-a)/2),
164        (sqrt(a^2+4)-a)^k*(sqrt(a^2+4)+a)*2^(-k-1)*(-1)^k
165         /((sqrt(a^2+4)+a)/2+(sqrt(a^2+4)-a)/2)
166         +(sqrt(a^2+4)-a)*(sqrt(a^2+4)+a)^k*2^(-k-1)
167          /((sqrt(a^2+4)+a)/2+(sqrt(a^2+4)-a)/2)]);
168
169kill (all);
170done;
171
172/* should trigger an error */
173errcatch (matrix ([1], [1, 2]));
174[];
175
176/* should trigger an error */
177errcatch (matrix ([1], '(a + b)));
178[];
179
180/* SF bug # 3014545 "submatrix does not work as expected"
181 * works for me, throw in these tests to make sure
182 */
183
184(submatrix (10, 20, zeromatrix (20, 20)), [length (%%), length (%%[1])]);
185[18, 20];
186
187(kill (F), F : 1 + zeromatrix (5, 5), submatrix (2, 5, F, 2, 5));
188matrix ([1, 1, 1], [1, 1, 1], [1, 1, 1]);
189
190submatrix (3, 5, F, 3, 5);
191matrix ([1, 1, 1], [1, 1, 1], [1, 1, 1]);
192
193F;
194matrix ([1, 1, 1, 1, 1], [1, 1, 1, 1, 1], [1, 1, 1, 1, 1], [1, 1, 1, 1, 1], [1, 1, 1, 1, 1]);
195
196(F : matrix ([1, 2, 3, 4], [5, 6, 7, 8], [9, 10, 11, 12]),
197 submatrix (F, 2, 4));
198matrix ([1, 3], [5, 7], [9, 11]);
199
200submatrix (1, 3, F);
201matrix ([5, 6, 7, 8]);
202
203/* next one is mostly just to ensure it doesn't trigger an error */
204submatrix (1, 2, 3, F);
205matrix ();
206
207/* next one is mostly just to ensure it doesn't trigger an error */
208submatrix (F, 1, 2, 3, 4);
209matrix ([], [], []);
210
211F;
212matrix ([1, 2, 3, 4], [5, 6, 7, 8], [9, 10, 11, 12]);
213
214submatrix (F);
215matrix ([1, 2, 3, 4], [5, 6, 7, 8], [9, 10, 11, 12]);
216
217/* name collision with special variables in 1-d output
218 * see mailing list circa 2012-01-09, "invert_by_lu does not work as expected"
219 */
220
221invert_by_lu (matrix ([v [0]]));
222matrix ([1 / v [0]]);
223
224/* additional tests for invert */
225
226/* Attempting to verify the effect of the ratmx and detout ev flags
227 * is quite a mess. ratmx produces CRE but the parser produces
228 * expressions which have a different operator (RAT, versus MRAT for CRE).
229 * detout produces an unsimplified "*" expression, which is quite
230 * readily simplified away; I am reminded of 19th century efforts to
231 * isolate halogens and alkali metals. Anyway, we'll do what we can.
232 */
233
234/* symbolic elements */
235
236(kill (M, M1), M : matrix ([a, b], [c, d]), 0);
2370;
238
239M1 : invert (M), ratsimp;
240matrix([d/(a*d-b*c),-b/(a*d-b*c)],[-c/(a*d-b*c),a/(a*d-b*c)]);
241
242ratsimp ([M1 . M, M . M1]);
243[matrix ([1, 0], [0, 1]), matrix ([1, 0], [0, 1])];
244
245is (invert (M) = M^^-1);
246true;
247
248(M1 : ev (invert (M), detout=true, doscmxops=false, doallmxops=false),
249 block ([inflag:true], [op (M1), ratsimp (args (M1))]));
250["*",[1/(a*d-b*c),matrix([d,-b],[-c,a])]];
251
252is (invert (M) = M^^-1), detout=true, doscmxops=false, doallmxops=false;
253true;
254
255block ([foo : matrix([d/(d*a-c*b),-(b/(d*a-c*b))],[-(c/(d*a-c*b)),a/(d*a-c*b)])],
256 ev (invert (M), ratmx=true), if equal (%%, foo) then true else %%);
257true;
258
259is (invert (M) = M^^-1), ratmx=true;
260true;
261
262block ([foo : ev (invert (M), ratmx=true, detout=true, doscmxops=false, doallmxops=false)],
263  [op (foo), first (foo), second (foo)],
264  if equal (%%, ["/", matrix ([d, -b], [-c, a]), a*d - b*c]) then true else %%);
265true;
266
267is (invert (M) = M^^-1), ratmx=true, detout=true, doscmxops=false, doallmxops=false;
268true;
269
270/* bigfloat elements */
271
272(M : ev (M, a = 1b0, b = 2b0, c = 3b0, d = -4b0), 0);
2730;
274
275invert (M);
276matrix([4.0b-1,2.0b-1],[3.0b-1,-1.0b-1]);
277
278is (invert (M) = M^^-1);
279true;
280
281(M1 : ev (invert (M), detout=true, doscmxops=false, doallmxops=false),
282 ev ([op (M1), args (M1)], simp=false, inflag=true));
283["*", [-0.1b0, matrix([-4.0b0, -2.0b0], [-3.0b0, 1.0b0])]];
284
285is (invert (M) = M^^-1), detout=true, doscmxops=false, doallmxops=false;
286true;
287
288(M1 : ev (invert (M), ratmx=true),
289 if every (ratp, M1) and equal (M1, matrix ([2/5, 1/5], [3/10, -(1/10)])) then true else M1);
290true;
291
292is (invert (M) = M^^-1), ratmx=true;
293true;
294
295(M1 : ev (invert (M), ratmx=true, detout=true, doscmxops=false, doallmxops=false),
296 [o, a] : ev ([op (M1), args (M1)], simp=false, inflag=true),
297 if ?caar (a [1]) = ?rat and every (ratp (a [2])) and equal (%%, ["*", [-1/10, matrix ([-4, -2], [-3, 1])]]) then true else %%);
298true;
299
300is (invert (M) = M^^-1), ratmx=true, detout=true, doscmxops=false, doallmxops=false;
301true;
302
303/* float elements */
304
305(M : float (M), 0);
3060;
307
308invert (M);
309matrix([4.0e-1,2.0e-1],[3.0e-1,-1.0e-1]);
310
311is (invert (M) = M^^-1);
312true;
313
314(M1 : ev (invert (M), detout=true, doscmxops=false, doallmxops=false),
315 ev ([op (M1), args (M1)], simp=false, inflag=true));
316["*", [-0.1e0, matrix([-4.0e0, -2.0e0], [-3.0e0, 1.0e0])]];
317
318is (invert (M) = M^^-1), detout=true, doscmxops=false, doallmxops=false;
319true;
320
321(M1 : ev (invert (M), ratmx=true),
322 if every (ratp, M1) and equal (M1, matrix ([2/5, 1/5], [3/10, -(1/10)])) then true else M1);
323true;
324
325is (invert (M) = M^^-1), ratmx=true;
326true;
327
328(M1 : ev (invert (M), ratmx=true, detout=true, doscmxops=false, doallmxops=false),
329 [o, a] : ev ([op (M1), args (M1)], simp=false, inflag=true),
330 if ?caar (a [1]) = ?rat and every (ratp (a [2])) and equal (%%, ["*", [-1/10, matrix ([-4, -2], [-3, 1])]]) then true else %%);
331true;
332
333is (invert (M) = M^^-1), ratmx=true, detout=true, doscmxops=false, doallmxops=false;
334true;
335
336/* handle detout=true correctly when determinant=1
337 * reported to mailing list 2015-01-22, "Matrix inversion with detout = true?"
338 */
339M : ident (4) $
340matrix ([1,0,0,0], [0,1,0,0], [0,0,1,0], [0,0,0,1]) $
341
342M^^-1, detout=true, doscmxops=false, doallmxops=false;
343matrix ([1,0,0,0], [0,1,0,0], [0,0,1,0], [0,0,0,1]) $
344
345/* test the various matrix inversion functions to make sure they all handle detout correctly */
346
347M^^-1, detout=true, doscmxops=false, doallmxops=false, invert_method='adjoint;
348matrix ([1,0,0,0], [0,1,0,0], [0,0,1,0], [0,0,0,1]) $
349
350M^^-1, detout=true, doscmxops=false, doallmxops=false, invert_method='lu;
351matrix ([1,0,0,0], [0,1,0,0], [0,0,1,0], [0,0,0,1]) $
352
353M^^-1, detout=true, doscmxops=false, doallmxops=false, invert_method='gausselim;
354matrix ([1,0,0,0], [0,1,0,0], [0,0,1,0], [0,0,0,1]) $
355
356M : matrix ([17, 29], [1, 42]);
357matrix ([17, 29], [1, 42]);
358
359(M1 : ev (M^^-1, detout=true, doscmxops=false, doallmxops=false, invert_method='adjoint),
360 block ([inflag:true], [op(M1), args(M1)]));
361["*",[1/685,matrix([42,-29],[-1,17])]] $
362
363(M1 : ev (M^^-1, detout=true, doscmxops=false, doallmxops=false, invert_method='lu),
364 block ([inflag:true], [op(M1), args(M1)]));
365["*",[1/685,matrix([42,-29],[-1,17])]] $
366
367(M1 : ev (M^^-1, detout=true, doscmxops=false, doallmxops=false, invert_method='gausselim),
368 block ([inflag:true], [op(M1), args(M1)]));
369["*",[1/685,matrix([42,-29],[-1,17])]] $
370
371
372/* a matrix of modest size, the subject of bug report #2362 */
373
374(M:matrix([1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1],
375 [0,1,0,-1,0,1,-1,-1,1,3,0,-3,0,3,1,-1,-3,-3,-1,1,3,3,-3,-3,3],
376 [0,0,1,0,-1,1,1,-1,-1,0,3,0,-3,1,3,3,1,-1,-3,-3,-1,3,3,-3,-3],
377 [0,1,0,1,0,1,1,1,1,9,0,9,0,9,1,1,9,9,1,1,9,9,9,9,9],
378 [0,0,1,0,1,1,1,1,1,0,9,0,9,1,9,9,1,1,9,9,1,9,9,9,9],
379 [0,0,0,0,0,1,-1,1,-1,0,0,0,0,3,3,-3,-3,3,3,-3,-3,9,-9,9,-9],
380 [0,1,0,-1,0,1,-1,-1,1,27,0,-27,0,27,1,-1,-27,-27,-1,1,27,27,-27,-27,27],
381 [0,0,0,0,0,1,1,-1,-1,0,0,0,0,9,3,3,9,-9,-3,-3,-9,27,27,-27,-27],
382 [0,0,0,0,0,1,-1,-1,1,0,0,0,0,3,9,-9,-3,-3,-9,9,3,27,-27,-27,27],
383 [0,0,1,0,-1,1,1,-1,-1,0,27,0,-27,1,27,27,1,-1,-27,-27,-1,27,27,-27,-27],
384 [0,1,0,1,0,1,1,1,1,81,0,81,0,81,1,1,81,81,1,1,81,81,81,81,81],
385 [0,0,0,0,0,1,-1,1,-1,0,0,0,0,27,3,-3,-27,27,3,-3,-27,81,-81,81,-81],
386 [0,0,0,0,0,1,1,1,1,0,0,0,0,9,9,9,9,9,9,9,9,81,81,81,81],
387 [0,0,0,0,0,1,-1,1,-1,0,0,0,0,3,27,-27,-3,3,27,-27,-3,81,-81,81,-81],
388 [0,0,1,0,1,1,1,1,1,0,81,0,81,1,81,81,1,1,81,81,1,81,81,81,81],
389 [0,0,0,0,0,1,1,-1,-1,0,0,0,0,81,3,3,81,-81,-3,-3,-81,243,243,-243,-243],
390 [0,0,0,0,0,1,-1,-1,1,0,0,0,0,27,9,-9,-27,-27,-9,9,27,243,-243,-243,243],
391 [0,0,0,0,0,1,1,-1,-1,0,0,0,0,9,27,27,9,-9,-27,-27,-9,243,243,-243,-243],
392 [0,0,0,0,0,1,-1,-1,1,0,0,0,0,3,81,-81,-3,-3,-81,81,3,243,-243,-243,243],
393 [0,0,0,0,0,1,1,1,1,0,0,0,0,81,9,9,81,81,9,9,81,729,729,729,729],
394 [0,0,0,0,0,1,-1,1,-1,0,0,0,0,27,27,-27,-27,27,27,-27,-27,729,-729,729,-729],
395 [0,0,0,0,0,1,1,1,1,0,0,0,0,9,81,81,9,9,81,81,9,729,729,729,729],
396 [0,0,0,0,0,1,1,-1,-1,0,0,0,0,81,27,27,81,-81,-27,-27,-81,2187,2187,-2187,-2187],
397 [0,0,0,0,0,1,-1,-1,1,0,0,0,0,27,81,-81,-27,-27,-81,81,27,2187,-2187,-2187,2187],
398 [0,0,0,0,0,1,1,1,1,0,0,0,0,81,81,81,81,81,81,81,81,6561,6561,6561,6561]),
399 invert (M));
400matrix([1,0,0,-10/9,-10/9,0,0,0,0,0,1/9,0,100/81,0,1/9,0,0,0,0,-10/81,0,-10/81,0,0,1/81],
401       [0,9/16,0,9/16,0,0,-1/16,0,-5/8,0,-1/16,0,-5/8,0,0,0,5/72,0,1/16,5/72,0,1/16,0,-1/144,-1/144],
402       [0,0,9/16,0,9/16,0,0,-5/8,0,-1/16,0,0,-5/8,0,-1/16,1/16,0,5/72,0,1/16,0,5/72,-1/144,0,-1/144],
403       [0,-9/16,0,9/16,0,0,1/16,0,5/8,0,-1/16,0,-5/8,0,0,0,-5/72,0,-1/16,5/72,0,1/16,0,1/144,-1/144],
404       [0,0,-9/16,0,9/16,0,0,5/8,0,1/16,0,0,-5/8,0,-1/16,-1/16,0,-5/72,0,1/16,0,5/72,1/144,0,-1/144],
405       [0,0,0,0,0,81/256,0,81/256,81/256,0,0,-9/256,81/256,-9/256,0,-9/256,-9/256,-9/256,-9/256,-9/256,1/256,-9/256,1/256,1/256,1/256],
406       [0,0,0,0,0,-81/256,0,81/256,-81/256,0,0,9/256,81/256,9/256,0,-9/256,9/256,-9/256,9/256,-9/256,-1/256,-9/256,1/256,-1/256,1/256],
407       [0,0,0,0,0,81/256,0,-81/256,-81/256,0,0,-9/256,81/256,-9/256,0,9/256,9/256,9/256,9/256,-9/256,1/256,-9/256,-1/256,-1/256,1/256],
408       [0,0,0,0,0,-81/256,0,-81/256,81/256,0,0,9/256,81/256,9/256,0,9/256,-9/256,9/256,-9/256,-9/256,-1/256,-9/256,-1/256,1/256,1/256],
409       [0,-1/48,0,-1/144,0,0,1/48,0,5/216,0,1/144,0,5/648,0,0,0,-5/216,0,-1/432,-5/648,0,-1/1296,0,1/432,1/1296],
410       [0,0,-1/48,0,-1/144,0,0,5/216,0,1/48,0,0,5/648,0,1/144,-1/432,0,-5/216,0,-1/1296,0,-5/648,1/432,0,1/1296],
411       [0,1/48,0,-1/144,0,0,-1/48,0,-5/216,0,1/144,0,5/648,0,0,0,5/216,0,1/432,-5/648,0,-1/1296,0,-1/432,1/1296],
412       [0,0,1/48,0,-1/144,0,0,-5/216,0,-1/48,0,0,5/648,0,1/144,1/432,0,5/216,0,-1/1296,0,-5/648,-1/432,0,1/1296],
413       [0,0,0,0,0,-3/256,0,-1/256,-3/256,0,0,3/256,-1/256,1/768,0,1/256,3/256,1/2304,1/768,1/256,-1/768,1/2304,-1/2304,-1/768,-1/2304],
414       [0,0,0,0,0,-3/256,0,-3/256,-1/256,0,0,1/768,-1/256,3/256,0,1/768,1/2304,3/256,1/256,1/2304,-1/768,1/256,-1/768,-1/2304,-1/2304],
415       [0,0,0,0,0,3/256,0,-3/256,1/256,0,0,-1/768,-1/256,-3/256,0,1/768,-1/2304,3/256,-1/256,1/2304,1/768,1/256,-1/768,1/2304,-1/2304],
416       [0,0,0,0,0,3/256,0,-1/256,3/256,0,0,-3/256,-1/256,-1/768,0,1/256,-3/256,1/2304,-1/768,1/256,1/768,1/2304,-1/2304,1/768,-1/2304],
417       [0,0,0,0,0,-3/256,0,1/256,3/256,0,0,3/256,-1/256,1/768,0,-1/256,-3/256,-1/2304,-1/768,1/256,-1/768,1/2304,1/2304,1/768,-1/2304],
418       [0,0,0,0,0,-3/256,0,3/256,1/256,0,0,1/768,-1/256,3/256,0,-1/768,-1/2304,-3/256,-1/256,1/2304,-1/768,1/256,1/768,1/2304,-1/2304],
419       [0,0,0,0,0,3/256,0,3/256,-1/256,0,0,-1/768,-1/256,-3/256,0,-1/768,1/2304,-3/256,1/256,1/2304,1/768,1/256,1/768,-1/2304,-1/2304],
420       [0,0,0,0,0,3/256,0,1/256,-3/256,0,0,-3/256,-1/256,-1/768,0,-1/256,3/256,-1/2304,1/768,1/256,1/768,1/2304,1/2304,-1/768,-1/2304],
421       [0,0,0,0,0,1/2304,0,1/6912,1/6912,0,0,-1/2304,1/20736,-1/2304,0,-1/6912,-1/6912,-1/6912,-1/6912,-1/20736,1/2304,-1/20736,1/6912,1/6912,1/20736],
422       [0,0,0,0,0,-1/2304,0,1/6912,-1/6912,0,0,1/2304,1/20736,1/2304,0,-1/6912,1/6912,-1/6912,1/6912,-1/20736,-1/2304,-1/20736,1/6912,-1/6912,1/20736],
423       [0,0,0,0,0,1/2304,0,-1/6912,-1/6912,0,0,-1/2304,1/20736,-1/2304,0,1/6912,1/6912,1/6912,1/6912,-1/20736,1/2304,-1/20736,-1/6912,-1/6912,1/20736],
424       [0,0,0,0,0,-1/2304,0,-1/6912,1/6912,0,0,1/2304,1/20736,1/2304,0,1/6912,-1/6912,1/6912,-1/6912,-1/20736,-1/2304,-1/20736,-1/6912,1/6912,1/20736])$
425
426/* 16 by 16 example from mailing list 2013-06-27 */
427
428(kill (K, invK),
429 K:matrix([1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],[0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0],[0,0,494.5054945054945,-178.5714285714286,0,0,54.94505494505494,-13.73626373626375,0,-4.578754578754582,109.8901098901099,4.578754578754582,0,59.52380952380952,-73.26007326007327,-59.52380952380952],[0,0,-178.5714285714286,494.5054945054945,0,0,13.73626373626375,-302.1978021978022,0,-73.26007326007327,4.578754578754582,109.8901098901099,0,-73.26007326007327,-4.578754578754582,109.8901098901099],[0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0],[0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0],[0,0,54.94505494505494,13.73626373626375,0,0,494.5054945054945,178.5714285714286,0,59.52380952380952,73.26007326007327,-59.52380952380952,0,-4.578754578754582,-109.8901098901099,4.578754578754582],[0,0,-13.73626373626375,-302.1978021978022,0,0,178.5714285714286,494.5054945054945,0,73.26007326007327,-4.578754578754582,-109.8901098901099,0,73.26007326007327,4.578754578754582,-109.8901098901099],[0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0],[0,0,-4.578754578754582,-73.26007326007327,0,0,59.52380952380952,73.26007326007327,0,110.2389673818245,19.84126984126984,-48.49119134833421,0,7.674864817721962,-19.84126984126984,-22.85016570730857],[0,0,109.8901098901099,4.578754578754582,0,0,73.26007326007327,-4.578754578754582,0,19.84126984126984,110.2389673818245,-19.84126984126984,0,19.84126984126984,-48.49119134833421,-19.84126984126984],[0,0,4.578754578754582,109.8901098901099,0,0,-59.52380952380952,-109.8901098901099,0,-48.49119134833421,-19.84126984126984,110.2389673818245,0,-22.85016570730857,19.84126984126984,7.674864817721962],[0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0],[0,0,59.52380952380952,-73.26007326007327,0,0,-4.578754578754582,73.26007326007327,0,7.674864817721962,19.84126984126984,-22.85016570730857,0,110.2389673818245,-19.84126984126984,-48.49119134833421],[0,0,-73.26007326007327,-4.578754578754582,0,0,-109.8901098901099,4.578754578754582,0,-19.84126984126984,-48.49119134833421,19.84126984126984,0,-19.84126984126984,110.2389673818245,19.84126984126984],[0,0,-59.52380952380952,109.8901098901099,0,0,4.578754578754582,-109.8901098901099,0,-22.85016570730857,-19.84126984126984,7.674864817721962,0,-48.49119134833421,19.84126984126984,110.2389673818245]),
430 invK : invert (K),
431 is (mat_norm (K . invK - ident (16), inf) < 1e-14));
432true;
433
434/* 4 by 4 example from mailing list 2013-04-16 */
435
436(kill (invert_R_from_RealRefraction, R_from_RealRefraction, G, u, a, b, c, d),
437 R_from_RealRefraction: matrix(
438 [   G[a]^2-1,     -G[a]*u[a,1],   -G[a]*u[a,2] ,     -G[a]*u[a,3]   ],
439 [   -G[b]*u[b,1],  1+u[b,1]^2,       u[b,1]*u[b,2],    u[b,1]*u[b,3]  ],
440 [   -G[c]*u[c,2],  u[c,2]*u[c,1],    1+u[c,2]^2,       u[c,2]*u[c,3]  ],
441 [   -G[d]*u[d,3],  u[d,3]*u[d,1],    u[d,3]*u[d,2],    1+u[d,3]^2     ]
442 ),
443 invert_R_from_RealRefraction: invert(R_from_RealRefraction),
444 ratsimp (invert_R_from_RealRefraction . R_from_RealRefraction - ident (4)));
445matrix ([0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0]);
446
447/* invert_by_adjoint exists */
448
449(kill (foo, foo_inv),
450 foo : matrix ([1, 7, -20], [-1, 4, -2], [3, -2, 7]),
451 foo_inv : invert_by_adjoint (foo),
452 foo . foo_inv);
453matrix ([1, 0, 0], [0, 1, 0], [0, 0, 1]);
454
455/* invert_by_gausselim exists */
456
457(kill (foo, foo_inv),
458 foo : matrix ([1, 7, -20], [-1, 4, -2], [3, -2, 7]),
459 foo_inv : invert_by_gausselim (foo),
460 foo . foo_inv);
461matrix ([1, 0, 0], [0, 1, 0], [0, 0, 1]);
462
463/* end additional tests for invert */
464
465/* float inf / SIGN1 mischief -- see also rtest_extra */
466
467block ([I, M],
468  I:87^611,
469  M:91^211,
470  sign(sqrt(I)-M));
471pos;
472
473/* verify that verbify is not applied to array name
474 * SF bug #2865: "Use exp as a list"
475 */
476exp : [1, 2, 3];
477[1, 2, 3];
478
479exp[1] : 123;
480123;
481
482exp;
483[123, 2, 3];
484
485exp[1];
486123;
487
488remvalue (exp);
489[exp];
490
491/* further tests */
492
493/* verify that $array does not apply $verbify to array name */
494
495(kill (foo),
496 apply (array, [nounify ('foo), 10]),
497 member (nounify ('foo), arrays));
498true;
499
500/* verify that arrstore does not apply $verbify to array name */
501
502((nounify ('foo)) [0] :: 123,
503 first (listarray (nounify ('foo))));
504123;
505
506'foo[1];
507'foo[1];
508
509ev ('foo[1], nouns);
510foo[1];
511
512/* verify that mdefine does not apply $verbify to array function name */
513
514(kill (foo, x, a), define (funmake (arraymake (nounify ('foo), [x]), [a]), x^a));
515'foo[x](a) := x^a;
516
517member (nounify (foo), arrays);
518true;
519
520/* verify that mdefine does not apply $verbify to function name */
521
522(kill (bar), define (funmake (nounify ('bar), ['x]), 2*'x), 0);
5230;
524
525(nounify ('bar))(u);
526'bar(u);
527
528ev ('bar(u), nouns);
529bar(u);
530
531/* verify that consfundef does not apply $verbify to function name */
532
533apply (fundef, [nounify ('bar)]);
534'bar(x) := 2*x;
535
536/* verify that consfundef does not apply $verbify to array function name */
537
538apply (fundef, [nounify ('foo)]);
539'foo[x](a) := x^a;
540
541/* function and array names which are strings are verbified, however */
542
543"bar"(x) := 1 + x;
544bar(x) := 1 + x;
545
546bar(u);
5471 + u;
548
549"bar"(z);
5501 + z;
551
552fundef (bar);
553bar(x) := 1 + x;
554
555fundef ("bar");
556bar(x) := 1 + x;
557
558member ('(bar(x)), functions);
559true;
560
561"baz"[x] := x^2;
562baz[x] := x^2;
563
564baz[10];
565100;
566
567"baz"[11];
568121;
569
570"quux"[u](v) := u - v;
571quux[u](v) := u - v;
572
573member (baz, arrays);
574true;
575
576member (quux, arrays);
577true;
578
579/* Bug #481: ('m)[1] (meval) */
580('m)[1];
581m[1];
582
583('m)(1);
584m(1);
585
586kill (functions, arrays);
587done;
588
589/* Verify that we catch malformed lambda expressions when they are simplified.
590 * These used to only be checked for if/when lambda expressions were applied
591 * to arguments.
592 */
593
594/* no parameter list */
595errcatch (lambda ())$
596[];
597
598/* empty body */
599errcatch (lambda ([x]))$
600[];
601
602/* non-symbol in parameter list */
603errcatch (lambda ([42], 'foo))$
604[];
605
606/* misplaced list parameter (for optional arguments) */
607errcatch (lambda ([[l], x], 'foo))$
608[];
609
610/* invalid list parameter (for optional arguments) */
611errcatch (lambda ([[l1, l2]], 'foo))$
612[];
613
614/* attempting to bind a constant is allowed */
615block ([c],
616  local (c),
617  declare (c, constant),
618  errcatch (lambda ([c], c)))$
619[lambda([c], c)];
620
621/* Verify that the parameter/variable lists of functions, lambda expressions,
622 * macros and blocks cannot contain duplicate variables.  Lots of cases...
623 */
624
625errcatch (foo (x, x) := x)$
626[];
627
628errcatch (foo (x, 'x) := x)$
629[];
630
631errcatch (foo (x, [x]) := x)$
632[];
633
634errcatch (foo (x, ['x]) := x)$
635[];
636
637errcatch (foo [x, x] := x)$
638[];
639
640errcatch (foo [x] (y, y) := y)$
641[];
642
643errcatch (foo [x] (y, 'y) := y)$
644[];
645
646errcatch (foo [x] (y, [y]) := y)$
647[];
648
649errcatch (foo [x] (y, ['y]) := y)$
650[];
651
652errcatch (foo [x, x] (y) := y)$
653[];
654
655errcatch (lambda ([x, x], x))$
656[];
657
658errcatch (lambda ([x, 'x], x))$
659[];
660
661errcatch (lambda ([x, [x]], x))$
662[];
663
664errcatch (lambda ([x, ['x]], x))$
665[];
666
667errcatch (foo (x, x) ::= x)$
668[];
669
670errcatch (foo (x, [x]) ::= x)$
671[];
672
673errcatch (block ([x, x], x))$
674[];
675
676errcatch (block ([x, x:foo], x))$
677[];
678
679/* try to verify that save(...) handles arrays correctly */
680
681(kill (all),
682 myundeclared[1234, foo] : sin(bar),
683 myundeclared[bar, baz + 1] : cos(quux),
684 myundeclared[1729, u^2] : tan(bar),
685 myundeclared[mumble, 7] : 9876,
686 array (mydeclared, 1, 1, 1),
687 mydeclared[0, 0, 0] : bar + 0,
688 mydeclared[0, 0, 1] : bar + 1,
689 mydeclared[0, 1, 0] : bar + 2,
690 mydeclared[0, 1, 1] : bar + 3,
691 mydeclared[1, 0, 0] : bar + 4,
692 mydeclared[1, 0, 1] : bar + 5,
693 mydeclared[1, 1, 0] : bar + 6,
694 mydeclared[1, 1, 1] : bar + 7,
695 myvalue : make_array (fixnum, 2, 2, 2, 2),
696 myvalue[0, 0, 0, 0] : %pi + 1,
697 myvalue[1, 1, 1, 1] : %e - 1,
698 use_fast_arrays : true,
699 myfast[foo, bar, baz] : blurf,
700 myfast["mumble", "abc", "xy", "Z"] : 2*blurf,
701 myfast[sin(foo), 1 - baz] : 3*blurf,
702 reset (use_fast_arrays),
703 save (sconcat (maxima_tempdir, "/tmpsavearrays.lisp"), arrays, values),
704 0);
7050;
706
707[arrays, values];
708[[myundeclared, mydeclared], [myvalue, myfast]];
709
710arrayinfo (myundeclared);
711[hashed, 2, [1234,foo],[1729,u^2],[bar,baz+1],[mumble,7]];
712
713listarray (myundeclared);
714[sin(bar), tan(bar), cos(quux), 9876];
715
716arrayinfo (mydeclared);
717[declared, 3, [1, 1, 1]];
718
719listarray (mydeclared);
720[bar, bar + 1, bar + 2, bar + 3, bar + 4, bar + 5, bar + 6, bar + 7];
721
722arrayinfo (myvalue);
723[declared, 4, [1, 1, 1, 1]]; /* "declared" seems wrong here ... oh well */
724
725listarray (myvalue);
726[%pi + 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, %e - 1];
727
728(arrayinfo (myfast), [%%[1], %%[2], sort (rest (%%, 2))]);
729[hash_table, true, [["mumble", "abc", "xy", "Z"], [foo, bar, baz], [sin(foo), 1 - baz]]];
730
731sort (listarray (myfast));
732[blurf, 2*blurf, 3*blurf];
733
734kill (myundeclared, mydeclared, myvalue, myfast);
735done;
736
737[arrays, values];
738[[], []];
739
740[errcatch (arrayinfo (myundeclared)),
741 errcatch (arrayinfo (mydeclared)),
742 errcatch (arrayinfo (myvalue)),
743 errcatch (arrayinfo (myfast))];
744[[], [], [], []];
745
746[errcatch (listarray (myundeclared)),
747 errcatch (listarray (mydeclared)),
748 errcatch (listarray (myvalue)),
749 errcatch (listarray (myfast))];
750[[], [], [], []];
751
752(load (sconcat (maxima_tempdir, "/tmpsavearrays.lisp")), 0);
7530;
754
755[arrays, values];
756[[myundeclared, mydeclared], [myvalue, myfast]];
757
758arrayinfo (myundeclared);
759[hashed, 2, [1234,foo],[1729,u^2],[bar,baz+1],[mumble,7]];
760
761listarray (myundeclared);
762[sin(bar), tan(bar), cos(quux), 9876];
763
764arrayinfo (mydeclared);
765[declared, 3, [1, 1, 1]];
766
767listarray (mydeclared);
768[bar, bar + 1, bar + 2, bar + 3, bar + 4, bar + 5, bar + 6, bar + 7];
769
770arrayinfo (myvalue);
771[declared, 4, [1, 1, 1, 1]]; /* "declared" seems wrong here ... oh well */
772
773listarray (myvalue);
774[%pi + 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, %e - 1];
775
776(arrayinfo (myfast), [%%[1], %%[2], sort (rest (%%, 2))]);
777[hash_table, true, [["mumble", "abc", "xy", "Z"], [foo, bar, baz], [sin(foo), 1 - baz]]];
778
779sort (listarray (myfast));
780[blurf, 2*blurf, 3*blurf];
781
782/* firstn, lastn */
783
784(kill (a, b, z), mylist : [a, 3, %pi, b, -2, %e, z, 0, %phi], 0);
7850;
786
787firstn (mylist, 0);
788[];
789
790firstn (mylist, 1);
791[a];
792
793firstn (mylist, 4);
794[a, 3, %pi, b];
795
796firstn (mylist, 8);
797[a, 3, %pi, b, -2, %e, z, 0];
798
799firstn (mylist, 9);
800[a, 3, %pi, b, -2, %e, z, 0, %phi];
801
802firstn (mylist, 10);
803[a, 3, %pi, b, -2, %e, z, 0, %phi];
804
805firstn (mylist, 100);
806[a, 3, %pi, b, -2, %e, z, 0, %phi];
807
808errcatch (firstn (mylist, -4));
809[];
810
811lastn (mylist, 0);
812[];
813
814lastn (mylist, 1);
815[%phi];
816
817lastn (mylist, 4);
818[%e, z, 0, %phi];
819
820lastn (mylist, 8);
821[3, %pi, b, -2, %e, z, 0, %phi];
822
823lastn (mylist, 9);
824[a, 3, %pi, b, -2, %e, z, 0, %phi];
825
826lastn (mylist, 10);
827[a, 3, %pi, b, -2, %e, z, 0, %phi];
828
829lastn (mylist, 100);
830[a, 3, %pi, b, -2, %e, z, 0, %phi];
831
832errcatch (lastn (mylist, -4));
833[];
834
835(kill (foo), fooexpr : funmake (foo, mylist), 0);
8360;
837
838firstn (fooexpr, 0);
839foo();
840
841firstn (fooexpr, 1);
842foo(a);
843
844firstn (fooexpr, 4);
845foo(a, 3, %pi, b);
846
847firstn (fooexpr, 8);
848foo(a, 3, %pi, b, -2, %e, z, 0);
849
850firstn (fooexpr, 9);
851foo(a, 3, %pi, b, -2, %e, z, 0, %phi);
852
853firstn (fooexpr, 10);
854foo(a, 3, %pi, b, -2, %e, z, 0, %phi);
855
856firstn (fooexpr, 100);
857foo(a, 3, %pi, b, -2, %e, z, 0, %phi);
858
859errcatch (firstn (fooexpr, -4));
860[];
861
862lastn (fooexpr, 0);
863foo();
864
865lastn (fooexpr, 1);
866foo(%phi);
867
868lastn (fooexpr, 4);
869foo(%e, z, 0, %phi);
870
871lastn (fooexpr, 8);
872foo(3, %pi, b, -2, %e, z, 0, %phi);
873
874lastn (fooexpr, 9);
875foo(a, 3, %pi, b, -2, %e, z, 0, %phi);
876
877lastn (fooexpr, 10);
878foo(a, 3, %pi, b, -2, %e, z, 0, %phi);
879
880lastn (fooexpr, 100);
881foo(a, 3, %pi, b, -2, %e, z, 0, %phi);
882
883errcatch (lastn (fooexpr, -4));
884[];
885
886(kill (bar), nary ("bar"));
887"bar";
888
889barexpr : apply ("bar", mylist);
890a bar 3 bar %pi bar b bar (-2) bar %e bar z bar 0 bar %phi;
891
892firstn (barexpr, 0);
893"bar"();
894
895firstn (barexpr, 1);
896"bar"(a);
897
898firstn (barexpr, 4);
899a bar 3 bar %pi bar b;
900
901firstn (barexpr, 8);
902a bar 3 bar %pi bar b bar (-2) bar %e bar z bar 0;
903
904firstn (barexpr, 9);
905a bar 3 bar %pi bar b bar (-2) bar %e bar z bar 0 bar %phi;
906
907firstn (barexpr, 10);
908a bar 3 bar %pi bar b bar (-2) bar %e bar z bar 0 bar %phi;
909
910firstn (barexpr, 100);
911a bar 3 bar %pi bar b bar (-2) bar %e bar z bar 0 bar %phi;
912
913errcatch (firstn (barexpr, -4));
914[];
915
916lastn (barexpr, 0);
917"bar"();
918
919lastn (barexpr, 1);
920"bar"(%phi);
921
922lastn (barexpr, 4);
923%e bar z bar 0 bar %phi;
924
925lastn (barexpr, 8);
9263 bar %pi bar b bar (-2) bar %e bar z bar 0 bar %phi;
927
928lastn (barexpr, 9);
929a bar 3 bar %pi bar b bar (-2) bar %e bar z bar 0 bar %phi;
930
931lastn (barexpr, 10);
932a bar 3 bar %pi bar b bar (-2) bar %e bar z bar 0 bar %phi;
933
934lastn (barexpr, 100);
935a bar 3 bar %pi bar b bar (-2) bar %e bar z bar 0 bar %phi;
936
937errcatch (lastn (barexpr, -4));
938[];
939
940(kill ("bar", baz, k), bazexpr : funmake (baz[k], mylist), 0);
9410;
942
943firstn (bazexpr, 0);
944baz[k]();
945
946firstn (bazexpr, 1);
947baz[k](a);
948
949firstn (bazexpr, 4);
950baz[k](a, 3, %pi, b);
951
952firstn (bazexpr, 8);
953baz[k](a, 3, %pi, b, -2, %e, z, 0);
954
955firstn (bazexpr, 9);
956baz[k](a, 3, %pi, b, -2, %e, z, 0, %phi);
957
958firstn (bazexpr, 10);
959baz[k](a, 3, %pi, b, -2, %e, z, 0, %phi);
960
961firstn (bazexpr, 100);
962baz[k](a, 3, %pi, b, -2, %e, z, 0, %phi);
963
964errcatch (firstn (bazexpr, -4));
965[];
966
967lastn (bazexpr, 0);
968baz[k]();
969
970lastn (bazexpr, 1);
971baz[k](%phi);
972
973lastn (bazexpr, 4);
974baz[k](%e, z, 0, %phi);
975
976lastn (bazexpr, 8);
977baz[k](3, %pi, b, -2, %e, z, 0, %phi);
978
979lastn (bazexpr, 9);
980baz[k](a, 3, %pi, b, -2, %e, z, 0, %phi);
981
982lastn (bazexpr, 10);
983baz[k](a, 3, %pi, b, -2, %e, z, 0, %phi);
984
985lastn (bazexpr, 100);
986baz[k](a, 3, %pi, b, -2, %e, z, 0, %phi);
987
988errcatch (lastn (bazexpr, -4));
989[];
990
991/* bug reported to mailing list 2017-09-05: "compiling a function including sorting" */
992
993(kill(test), test(L):=sort(L,lambda([s,t],s[2]>t[2])), 0);
9940;
995
996(L:[[b,2],[a,1],[d,4],[c,3]], test (L));
997[[d,4],[c,3],[b,2],[a,1]];
998
999(compile (test), test (L));
1000[[d,4],[c,3],[b,2],[a,1]];
1001
1002(mycmp (s, t) := s[2]>t[2], 0);
10030;
1004
1005sort (L, mycmp);
1006[[d,4],[c,3],[b,2],[a,1]];
1007
1008(compile (mycmp), is (?fboundp (mycmp) # false));
1009true;
1010
1011sort (L, mycmp);
1012[[d,4],[c,3],[b,2],[a,1]];
1013
1014
1015/* bug reported to mailing list 2017-12-21: "array() with use_fast_arrays: true" */
1016
1017(use_fast_arrays : true,
1018 kill (a, b, a0, a1),
1019 array ([a, b], 1),
1020 map (?arrayp, %%));
1021[true, true];
1022
1023map (?length, [a, b]);
1024[2, 2];
1025
1026(a[0] : a0,
1027 a[1] : a1,
1028 listarray (a));
1029[a0, a1];
1030
1031listarray (b);
1032[false, false];
1033
1034(reset (), 0);
10350;
1036
1037/* mailing list 2017-12-28: "op('(do 1))" */
1038
1039(e : '(for i:111 thru 114 do push(i, L)),
1040 [op(e), args(e)]);
1041["do", [i,111,false,false,114,false,push(i,L)]];
1042
1043e1 : apply (funmake, ["do", '[i,111,false,false,114,false,push(i,L)]]);
1044for i:111 thru 114 do push(i, L);
1045
1046block ([L:[]], ev(e1), L);
1047[114, 113, 112, 111];
1048
1049(e : '(for x in [100, 200, 300] do push(x - 1, M)),
1050 [op(e), args(e)]);
1051["do_in", [x,[100,200,300],false,false,false,false,push(x-1,M)]];
1052
1053e1 : apply (funmake, ["do_in", '[x,[100,200,300],false,false,false,false,push(x-1,M)]]);
1054for x in [100, 200, 300] do push(x - 1, M);
1055
1056block ([M:[]], ev(e1), M);
1057[299, 199, 99];
1058
1059(reset(use_fast_arrays),1);
10601;
1061
1062/* ensure arraymake cannot create an invalid expression
1063 * follow-on work to commit 3140a35
1064 */
1065
1066(kill (foo, bar, baz, quux),
1067 arraymake (foo, [bar]));
1068foo[bar];
1069
1070arraymake (foo[bar], [baz]);
1071foo[bar][baz];
1072
1073arraymake (foo[bar][baz], [quux]);
1074foo[bar][baz][quux];
1075
1076arraymake (foo[bar, baz], [111, 222]);
1077foo[bar, baz][111, 222];
1078
1079errcatch (arraymake (1, [bar]));
1080[];
1081
1082errcatch (arraymake ("foo", [bar]));
1083[];
1084
1085errcatch (arraymake (lambda([x], x), [bar]));
1086[];
1087
1088/* slightly unusual constructs which do yield valid expressions, and are therefore permissable */
1089
1090arraymake (foo + bar, [baz]);
1091(foo + bar)[baz];
1092
1093arraymake (foo(bar), [baz]);
1094foo(bar)[baz];
1095
1096