tests/rtest2.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 kill(functions,arrays,values);
  10 done$
  11 use_fast_arrays:false;
  12 false;
  13 a[n]:=n*a[n-1];
  14 a[n]:=n*a[n-1]$
  15 a[0]:1;
  16 1$
  17 a[5];
  18 120$
  19 a[n]:=n;
  20 a[n]:=n$
  21 a[6];
  22 6$
  23 a[4];
  24 24$
  25 (use_fast_arrays:true,kill(a));
  26 done;
  27 lambda([x,y,z],x^2+y^2+z^2);
  28 lambda([x,y,z],x^2+y^2+z^2)$
  29 %(1,2,a);
  30 a^2+5$
  31 1+2+a;
  32 a+3$
  33 exp:[x^2,y/3,-2];
  34 [x^2,y/3,-2]$
  35 %[1]*x;
  36 x^3$
  37 [a,exp,%];
  38 [a,[x^2,y/3,-2],x^3]$
  39 m:matrix([a,0],[b,1]);
  40 matrix([a,0],[b,1])$
  41 m^2;
  42 matrix([a^2,0],[b^2,1])$
  43 exp:m . m;
  44 matrix([a^2,0],[a*b+b,1])$
  45 m[1,1]*m;
  46 matrix([a^2,0],[a*b,a])$
  47 %-exp+1;
  48 matrix([1,1],[1-b,a])$
  49 m^^(-1);
  50 matrix([1/a,0],[-b/a,1])$
  51 [x,y] . m;
  52 matrix([b*y+a*x,y])$
  53 matrix([a,b,c],[d,e,f],[g,h,i]);
  54 matrix([a,b,c],[d,e,f],[g,h,i])$
  55 %^^2;
  56 matrix([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])$
  59 exp:x+1 = y^2;
  60 x+1 = y^2$
  61 x-1 = 2*y+1;
  62 x-1 = 2*y+1$
  63 exp+%;
  64 2*x = y^2+2*y+1$
  65 exp/y;
  66 (x+1)/y = y$
  67 1/%;
  68 y/(x+1) = 1/y$
  69 fib[n]:=if n = 1 or n = 2 then 1 else fib[n-1]+fib[n-2];
  70 fib[n]:=if n = 1 or n = 2 then 1 else fib[n-1]+fib[n-2]$
  71 fib[1]+fib[2];
  72 2$
  73 fib[3];
  74 2$
  75 fib[5];
  76 5$
  77 eta(mu,nu):=if mu = nu then mu else (if mu > nu then mu-nu else mu+nu);
  78 eta(mu,nu):=if mu = nu then mu else (if mu > nu then mu-nu else mu+nu)$
  79 eta(5,6);
  80 11$
  81 eta(eta(7,7),eta(1,2));
  82 4$
  83 if not 5 >= 2 and 6 <= 5 or 4+1 > 3 then a else b;
  84 a$
  85 kill(f);
  86 done$
  87
  88 kill(x,y,z);
  89 done$
  90 determinant(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
  93 subst(1,z,quotient(%,-54*a^2));
  94 y^3+a*x*y+x^3$
  95 f(x):=block([a,y],local(a),y:4,a[y]:x,display(a[y]));
  96 f(x):=block([a,y],local(a),y:4,a[y]:x,display(a[y]))$
  97 y:2;
  98 2$
  99 a[y+2]:0;
 100 0$
 101 f(9);
 102 done$
 103 a[y+2];
 104 0$
 105
 106 (use_fast_arrays : false, kill(a), 0);
 107 0$
 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);
 121 0;
 122
 123 /* construct a matrix of modest size */
 124 (apply (matrix, makelist ([i], i, 1, 10000)), 0);
 125 0;
 126
 127 /* verify that arguments are evaluated exactly once */
 128 block ([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 */
 132 block ([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));
 141 matrix ([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
 153 matrix([(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
 169 kill (all);
 170 done;
 171
 172 /* should trigger an error */
 173 errcatch (matrix ([1], [1, 2]));
 174 [];
 175
 176 /* should trigger an error */
 177 errcatch (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));
 188 matrix ([1, 1, 1], [1, 1, 1], [1, 1, 1]);
 189
 190 submatrix (3, 5, F, 3, 5);
 191 matrix ([1, 1, 1], [1, 1, 1], [1, 1, 1]);
 192
 193 F;
 194 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]);
 195
 196 (F : matrix ([1, 2, 3, 4], [5, 6, 7, 8], [9, 10, 11, 12]),
 197  submatrix (F, 2, 4));
 198 matrix ([1, 3], [5, 7], [9, 11]);
 199
 200 submatrix (1, 3, F);
 201 matrix ([5, 6, 7, 8]);
 202
 203 /* next one is mostly just to ensure it doesn't trigger an error */
 204 submatrix (1, 2, 3, F);
 205 matrix ();
 206
 207 /* next one is mostly just to ensure it doesn't trigger an error */
 208 submatrix (F, 1, 2, 3, 4);
 209 matrix ([], [], []);
 210
 211 F;
 212 matrix ([1, 2, 3, 4], [5, 6, 7, 8], [9, 10, 11, 12]);
 213
 214 submatrix (F);
 215 matrix ([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
 221 invert_by_lu (matrix ([v [0]]));
 222 matrix ([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);
 237 0;
 238
 239 M1 : invert (M), ratsimp;
 240 matrix([d/(a*d-b*c),-b/(a*d-b*c)],[-c/(a*d-b*c),a/(a*d-b*c)]);
 241
 242 ratsimp ([M1 . M, M . M1]);
 243 [matrix ([1, 0], [0, 1]), matrix ([1, 0], [0, 1])];
 244
 245 is (invert (M) = M^^-1);
 246 true;
 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
 252 is (invert (M) = M^^-1), detout=true, doscmxops=false, doallmxops=false;
 253 true;
 254
 255 block ([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 %%);
 257 true;
 258
 259 is (invert (M) = M^^-1), ratmx=true;
 260 true;
 261
 262 block ([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 %%);
 265 true;
 266
 267 is (invert (M) = M^^-1), ratmx=true, detout=true, doscmxops=false, doallmxops=false;
 268 true;
 269
 270 /* bigfloat elements */
 271
 272 (M : ev (M, a = 1b0, b = 2b0, c = 3b0, d = -4b0), 0);
 273 0;
 274
 275 invert (M);
 276 matrix([4.0b-1,2.0b-1],[3.0b-1,-1.0b-1]);
 277
 278 is (invert (M) = M^^-1);
 279 true;
 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
 285 is (invert (M) = M^^-1), detout=true, doscmxops=false, doallmxops=false;
 286 true;
 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);
 290 true;
 291
 292 is (invert (M) = M^^-1), ratmx=true;
 293 true;
 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 %%);
 298 true;
 299
 300 is (invert (M) = M^^-1), ratmx=true, detout=true, doscmxops=false, doallmxops=false;
 301 true;
 302
 303 /* float elements */
 304
 305 (M : float (M), 0);
 306 0;
 307
 308 invert (M);
 309 matrix([4.0e-1,2.0e-1],[3.0e-1,-1.0e-1]);
 310
 311 is (invert (M) = M^^-1);
 312 true;
 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
 318 is (invert (M) = M^^-1), detout=true, doscmxops=false, doallmxops=false;
 319 true;
 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);
 323 true;
 324
 325 is (invert (M) = M^^-1), ratmx=true;
 326 true;
 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 %%);
 331 true;
 332
 333 is (invert (M) = M^^-1), ratmx=true, detout=true, doscmxops=false, doallmxops=false;
 334 true;
 335
 336 /* handle detout=true correctly when determinant=1
 337  * reported to mailing list 2015-01-22, "Matrix inversion with detout = true?"
 338  */
 339 M : ident (4) $
 340 matrix ([1,0,0,0], [0,1,0,0], [0,0,1,0], [0,0,0,1]) $
 341
 342 M^^-1, detout=true, doscmxops=false, doallmxops=false;
 343 matrix ([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
 347 M^^-1, detout=true, doscmxops=false, doallmxops=false, invert_method='adjoint;
 348 matrix ([1,0,0,0], [0,1,0,0], [0,0,1,0], [0,0,0,1]) $
 349
 350 M^^-1, detout=true, doscmxops=false, doallmxops=false, invert_method='lu;
 351 matrix ([1,0,0,0], [0,1,0,0], [0,0,1,0], [0,0,0,1]) $
 352
 353 M^^-1, detout=true, doscmxops=false, doallmxops=false, invert_method='gausselim;
 354 matrix ([1,0,0,0], [0,1,0,0], [0,0,1,0], [0,0,0,1]) $
 355
 356 M : matrix ([17, 29], [1, 42]);
 357 matrix ([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));
 400 matrix([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));
 432 true;
 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)));
 445 matrix ([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);
 453 matrix ([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);
 461 matrix ([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
 467 block ([I, M],
 468   I:87^611,
 469   M:91^211,
 470   sign(sqrt(I)-M));
 471 pos;
 472
 473 /* verify that verbify is not applied to array name
 474  * SF bug #2865: "Use exp as a list"
 475  */
 476 exp : [1, 2, 3];
 477 [1, 2, 3];
 478
 479 exp[1] : 123;
 480 123;
 481
 482 exp;
 483 [123, 2, 3];
 484
 485 exp[1];
 486 123;
 487
 488 remvalue (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));
 498 true;
 499
 500 /* verify that arrstore does not apply $verbify to array name */
 501
 502 ((nounify ('foo)) [0] :: 123,
 503  first (listarray (nounify ('foo))));
 504 123;
 505
 506 'foo[1];
 507 'foo[1];
 508
 509 ev ('foo[1], nouns);
 510 foo[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
 517 member (nounify (foo), arrays);
 518 true;
 519
 520 /* verify that mdefine does not apply $verbify to function name */
 521
 522 (kill (bar), define (funmake (nounify ('bar), ['x]), 2*'x), 0);
 523 0;
 524
 525 (nounify ('bar))(u);
 526 'bar(u);
 527
 528 ev ('bar(u), nouns);
 529 bar(u);
 530
 531 /* verify that consfundef does not apply $verbify to function name */
 532
 533 apply (fundef, [nounify ('bar)]);
 534 'bar(x) := 2*x;
 535
 536 /* verify that consfundef does not apply $verbify to array function name */
 537
 538 apply (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;
 544 bar(x) := 1 + x;
 545
 546 bar(u);
 547 1 + u;
 548
 549 "bar"(z);
 550 1 + z;
 551
 552 fundef (bar);
 553 bar(x) := 1 + x;
 554
 555 fundef ("bar");
 556 bar(x) := 1 + x;
 557
 558 member ('(bar(x)), functions);
 559 true;
 560
 561 "baz"[x] := x^2;
 562 baz[x] := x^2;
 563
 564 baz[10];
 565 100;
 566
 567 "baz"[11];
 568 121;
 569
 570 "quux"[u](v) := u - v;
 571 quux[u](v) := u - v;
 572
 573 member (baz, arrays);
 574 true;
 575
 576 member (quux, arrays);
 577 true;
 578
 579 kill (functions, arrays);
 580 done;
 581