## maxima-commits

 ```Update of /cvsroot/maxima/maxima/share/linearalgebra In directory sc8-pr-cvs1.sourceforge.net:/tmp/cvs-serv3795 Added Files: announcement.txt conjugate.lisp linalg-utilities.lisp linalg.demo linalg.mac linalg.usage lu.lisp matrixexp.lisp matrixexp.usage polynomialp.lisp test-conjugate.mac test-linalg.mac test-lu.mac test-matrixexp.mac test-polynomialp.mac Log Message: Linear algebra functions by Barton Willis. Copied verbatim from http://www.unk.edu/uploadedFiles/facstaff/profiles/willisb/linearalgebra(3).zip . --- NEW FILE: announcement.txt --- The new version of linalg can be downloaded from http://www.unk.edu/facstaff/profiles/willisb/. A tour of new functions: (1) The file linearalgebra/conjugate.lisp has an alternative to the /eigen/conjugate function. The /eigen/conjugate function only does the substitution %i --> -%i. The conjugate function in linearalgebra is (I hope) somewhat smarter. (%o1) c:/maxima/linearalgebra/linalg.mac (%i2) declare(z,complex)\$ (%i3) conjugate([z,z^2,sqrt(z)]); (%o3) [CONJUGATE(z),CONJUGATE(z)^2,CONJUGATE(sqrt(z))] Off the branch cut of sqrt, conjugate and sqrt commute. Telling Maxima that z isn't on branch cut allows it to transform conjugate(sqrt(z)) --> sqrt(conjugate(z)). (%i4) assume(imagpart(z) > 0)\$ (%i5) conjugate(sqrt(z)); (%o5) sqrt(CONJUGATE(z)) (2) New functions for generating matrices: (%i6) hilbert_matrix(2); (%o6) matrix([1,1/2],[1/2,1/3]) (%i7) vandermonde_matrix([a,b]); (%o7) matrix([1,a],[1,b]) (3) A matrix kronnecker_product (%i9) kronnecker_product(hilbert_matrix(2), vandermonde_matrix([5,6])); (%o9) matrix([1,5,1/2,5/2],[1,6,1/2,3],[1/2,5/2,1/3,5/3],[1/2,3,1/3,2]) (4) A function that converts block matrices to unblocked matrices (%i12) matrix([matrix([4,5]), matrix([6,7])]); (%o12) matrix([matrix([4,5]),matrix([6,7])]) (%i13) mat_unblocker(%); (%o13) matrix([4,5,6,7]) (5) LU factorization There is code for LU factorization. It uses partial pivoting when the matrix is numeric. (%i1) load("linalg"); Warning - you are redefining the MACSYMA function RANK (%o1) c:/maxima/linearalgebra/linalg.mac (%i2) lu_factor(matrix([a,b],[c,d])); (%o2) [matrix([a,b],[c/a,d-(b*c)/a]),[1,2]] Do the LU factorization using CRE expressions. Try this using the general representation -- it is very slow. (%i4) lu_factor(vandermonde_matrix([a,b,c,d,e]),crefield)\$ (%i5) m : first(%)\$ (%i6) factor(product(m[i,i],i,1,5)); (%o6) (b-a)*(c-a)*(c-b)*(d-a)*(d-b)*(d-c)*(e-a)*(e-b)*(e-c)*(e-d) LU factor using CL complex numbers (%i14) w[i,j] := ?random(1.0) + %i * ?random(1.0)\$ Evaluation took 0.00 seconds (0.00 elapsed) (%i15) m : genmatrix(w,100,100)\$ Evaluation took 0.45 seconds (0.45 elapsed) (%i16) lu_factor(m,complexfloatfield)\$ Evaluation took 1.73 seconds (1.73 elapsed) (This code will not win any speed awards.) There is also code for solving linear systems using the LU factorization (%i19) m : hilbert_matrix(5)\$ (%i20) m : lu_factor(m)\$ (%i21) lu_backsub(m, matrix([1],[1],[1],[1],[1])); (%o21) matrix([5],[-120],[630],[-1120],[630]) (%i22) hilbert_matrix(5) . %; (%o22) matrix([1],[1],[1],[1],[1]) (6) My spectral representation code is now included in linalg. As always, let me know how I can improve this code. Barton Willis --- NEW FILE: conjugate.lisp --- ;; Copyright 2005 by Barton Willis ;; This is free software; you can redistribute it and/or ;; modify it under the terms of the GNU General Public License, ;; http://www.gnu.org/copyleft/gpl.html. ;; This software has NO WARRANTY, not even the implied warranty of ;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. (\$put '\$conjugate 1 '\$version) (defprop \$conjugate tex-postfix tex) (defprop \$conjugate ("^\\star") texsym) (defprop \$conjugate 160. tex-lbp) ;; A fix feature(x,complex) --> true for all symbols bug. (defmfun \$featurep (e ind) (cond ((not (symbolp ind)) (merror "~M is not a symbolic atom - `featurep'." ind)) ((eq ind '\$integer) (maxima-integerp e)) ((eq ind '\$noninteger) (nonintegerp e)) ((eq ind '\$even) (mevenp e)) ((eq ind '\$odd) (moddp e)) ((eq ind '\$real) (if (atom e) (or (numberp e) (kindp e '\$real) (numberp (numer e))) (free (\$rectform e) '\$%i))) ((symbolp e) (kindp e ind)))) (defun op-equalp (e &rest op) (and (consp e) (consp (car e)) (some #'(lambda (s) (equal (caar e) s)) op))) (defprop \$conjugate simp-conjugate operators) (setf (get 'mabs 'realvaluedp) t) (setf (get '%realpart 'realvaluedp) t) (setf (get '%imagpart 'realvaluedp) t) ;; When a function commutes with the conjugate, give the function the ;; commutes-with-conjugate property. The log function commutes with ;; the conjugate on all of C except on the negative real axis. Thus ;; log does not get the commutes-with-conjugate property. Instead, ;; log gets the conjugate-function property. (setf (get 'mplus 'commutes-with-conjugate) t) (setf (get 'mtimes 'commutes-with-conjugate) t) (setf (get '%cosh 'commutes-with-conjugate) t) (setf (get '%sinh 'commutes-with-conjugate) t) (setf (get '%tanh 'commutes-with-conjugate) t) (setf (get '%sech 'commutes-with-conjugate) t) (setf (get '%csch 'commutes-with-conjugate) t) (setf (get '%coth 'commutes-with-conjugate) t) (setf (get '%cos 'commutes-with-conjugate) t) (setf (get '%sin 'commutes-with-conjugate) t) (setf (get '%tan 'commutes-with-conjugate) t) (setf (get '%sec 'commutes-with-conjugate) t) (setf (get '%csc 'commutes-with-conjugate) t) (setf (get '%cot 'commutes-with-conjugate) t) (setf (get '\$matrix 'commutes-with-conjugate) t) (setf (get 'mlist 'commutes-with-conjugate) t) (setf (get 'mequal 'commutes-with-conjugate) t) ;; When a function has the conjugate-function property, ;; use a non-generic function to conjugate it. (setf (get '%log 'conjugate-function) 'conjugate-log) (setf (get 'mexpt 'conjugate-function) 'conjugate-mexpt) ;; Return true iff Maxima can prove that z is not on the ;; negative real axis. (defun off-negative-real-axisp (z) (setq z (\$rectform z)) (let ((x (\$realpart z)) (y (\$imagpart z)) (\$prederror nil)) (or (eq t (mevalp `((mgreaterp) ,y 0))) (eq t (mevalp `((mgreaterp) 0 ,y))) (eq t (mevalp `((mgreaterp) ,x 0)))))) ;; Return conjugate(log(x)). Actually, x is a lisp list (x). (defun conjugate-log (x) (setq x (car x)) (if (off-negative-real-axisp x) (simplifya `((%log) ,(mfuncall '\$conjugate x)) nil) `((\$conjugate simp) ((%log) ,x)))) ;; Return conjugate(x^p), where e = (x, p) (defun conjugate-mexpt (e) (let ((\$prederror nil) (x (first e)) (p (second e))) (cond ((or (integerp p) (and (\$ratnump p) (off-negative-real-axisp x))) (power (mfuncall '\$conjugate x) p)) ((eq t (mevalp `((mgreaterp) ,x 0))) (power x (mfuncall '\$conjugate p))) (t `((\$conjugate simp) ,(power x p)))))) (defun simp-conjugate (e y z) (declare (ignore y)) (oneargcheck e) (let ((f)) (cond ((not (memq 'simp (car e))) (setq e (simplifya (nth 1 e) z)) (cond ((complexp e) (conjugate e)) ;; never happens, but might someday. ((\$mapatom e) (cond ((like e '\$%i) (mult -1 '\$%i)) ((\$featurep e '\$complex) `((\$conjugate simp) ,e)) (t e))) ((get (mop e) 'realvaluedp) e) ((op-equalp e '\$conjugate) (car (margs e))) ((get (mop e) 'commutes-with-conjugate) (mfuncall '\$map '\$conjugate e)) ((setq f (get (mop e) 'conjugate-function)) (funcall f (margs e))) (t `((\$conjugate simp) ,e)))) (t e)))) --- NEW FILE: linalg-utilities.lisp --- ;; Copyright 2005 by Barton Willis ;; This is free software; you can redistribute it and/or ;; modify it under the terms of the GNU General Public License, ;; http://www.gnu.org/copyleft/gpl.html. ;; This software has NO WARRANTY, not even the implied warranty of ;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. (\$put '\$linalgutilities 1 '\$version) (defmacro while (cond &rest body) `(do () ((not ,cond)) ,@body)) ;; Return true iff e is an expression with operator op1, op2,...,or opn. (defun op-equalp (e &rest op) (and (consp e) (consp (car e)) (some #'(lambda (s) (equal (caar e) s)) op))) ;; Return true iff the operator of e is a Maxima relation operator. (defun mrelationp (a) (op-equalp a 'mlessp 'mleqp 'mequal 'mgeqp 'mgreaterp)) (defun \$listp (e &optional (f nil)) (and (op-equalp e 'mlist) (or (eq f nil) (every #'(lambda (s) (eq t (mfuncall f s))) (margs e))))) (defun \$matrixp (e &optional (f nil)) (and (op-equalp e '\$matrix) (every #'(lambda (s) (\$listp s f)) (margs e)))) (defun \$require_nonempty_matrix (m pos fun) (if (not (and (\$matrixp m) (> (\$length m) 0) (> (\$length (\$first m)) 0))) (merror "The ~:M argument of the function ~:M must be a nonempty matrix" pos fun))) (defun \$require_matrix (m pos fun) (if (not (\$matrixp m)) (merror "The ~:M argument of the function ~:M must be a matrix" pos fun))) (defun \$require_square_matrix (m pos fun) (if (not (and (\$matrixp m) (= (\$length m) (\$length (\$first m))))) (merror "The ~:M argument of the function ~:M must be a square matrix" pos fun))) (defun \$matrix_size(m) (\$require_matrix m "\$first" "\$matrix_size") `((mlist) ,(\$length m) ,(\$length (\$first m)))) (defun \$require_list (lst pos fun) (if (not (\$listp lst)) (merror "The ~:M argument of the function ~:M must be a list" pos fun))) (defun \$require_posinteger(i pos fun) (if (not (and (integerp i) (> i 0))) (merror "The ~:M argument of the function ~:M must be a positive integer" pos fun))) (defun \$ctranspose (m) (let ((\$matrix_element_transpose `((lambda simp) ((mlist) x) ((mcond) ((\$matrixp) x) ((\$ctranspose) x) t ((\$conjugate) x))))) (mfuncall '\$transpose m))) --- NEW FILE: linalg.demo --- load("linalg"); m : matrix([1,2],[1,2]); nullspace(m); columnspace(m); ptriangularize(m-z*ident(2),z); m : matrix([1,2,3],[4,5,6],[7,8,9]) - z * ident(3); mm : ptriangularize(m,z); algebraic : true; tellrat(mm[3,3]); mm : ratsimp(mm); nullspace(mm); m : matrix([1,2,3,4],[5,6,7,8],[9,10,11,12],[13,14,15,16]); column_space(m); apply('orthogonal_complment, args(nullspace(transpose(m)))); --- NEW FILE: linalg.mac --- /* Author Barton Willis University of Nebraska at Kearney Copyright (C) 2005, Barton Willis Brief Description: Maxima code for finding the nullspace and column space of a matrix, along with code that triangularizes matrices using the method described in ``Eigenvalues by row operations,'' by Barton Willis. This is free software you can redistribute it and/or modify it under the terms of the GNU General Public License, http://www.gnu.org/copyleft/gpl.html. This software has NO WARRANTY, not even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. */ put('linalg,1,'version); /* The translator generates code that doesn't work when translate_fast_arrays is true. So set translate_fast_arrays to false. */ eval_when([translate, compile], translate_fast_arrays : false); eval_when(translate, declare_translated(require_matrix,matrix_size,ptriangularize_with_proviso,locate_matrix_entry, nullspace,require_integer,columnspace,rowswap,rowop,require_symbol, column_reduce,good_pivot,hipow_gzero,mat_norm,require_list))\$ eval_when([batch,load,translate,compile], (if get('nset, 'version) = false then load("nset"), if get('conjugate, 'version) = false then load("conjugate"), if get('polynomialp, 'version) = false then load("polynomialp"), if get('lu, 'version) = false then load("lu"), if get('matrixexp, 'version) = false then load("matrixexp"), if get('linalgutilities, 'version) = false then load("linalg-utilities"))); require_integer(i, pos, fn) := if integerp(i) then true else error("The", pos, "argument to" ,fn, "must be an integer"); require_list(l, pos, fn) := if listp(l) then true else error("The", pos, "argument to", fn, "must be a list"); require_symbol(x, pos, fn) := if symbolp(x) or subvarp(x) then true else error("The", pos, "argument to", fn, "must be a symbol"); request_rational_matrix(m, pos, fn) := if every('identity, map(lambda([s], every('ratnump,s)), args(m))) then true else print("Some entries in the matrix are not rational numbers. The result might be wrong."); matrix_size(m) := ( require_matrix(m,"first", "matrix_size"), [length(m), length(first(m))]); dotproduct(a,b) := block([scalarmatrixp : true], require_matrix(a,"first", "dotproduct"), require_matrix(b,"second", "dotproduct"), if matrix_size(a) # matrix_size(b) or second(matrix_size(a)) # 1 then error("Arguments to dotproduct must be vectors with the same size"), ctranspose(a) . b); /* I use algsys instead of linsolve because it seems that linsolve doesn't pay attention to tellrat information. */ nullspace(m) := block([n, x, e, vo, %r, proviso, nc, nr, pos, scalarmatrixp : false, domxmxops : true], require_matrix(m, "first", "nullspace"), nc : length(first(m)), nr : length(m), if nc = 0 or nr = 0 then error("The argument to 'nullspace' must be a matrix with one or more rows and columns"), m : ptriangularize_with_proviso(m, ?gensym()), proviso : second(m), m : first(m), for i : 1 thru nr do ( if listp(pos : locate_matrix_entry(m, i, 1, i, nc, lambda([x],x # 0), 'bool)) then ( if not(constantp(x : part(part(m, first(pos)),second(pos)))) then proviso : adjoin(x # 0,proviso))), if not(emptyp(proviso)) then print("Proviso: ",proviso), x : apply('matrix, makelist([?gensym()],i,1,nc)), e : apply('append, args(m.x)), e : map('rat,e), %r : setify(%rnum_list), e : first(algsys(e,(listofvars(x)))), %r : listify(setdifference(setify(%rnum_list), %r)), vo : if emptyp(%r) then [ ] else map("=",%r, makelist(0,i,1,length(%r))), x : subst(e, x), x : sort(makelist(subst(vo, subst(vi = 1, x)), vi, %r)), funmake('span, x)); nullity(m) := length(nullspace(m)); orthogonal_complement([v]) := block([sz], require_matrix(first(v)," ","orthogonal_complement"), sz : matrix_size(first(v)), if second(sz) # 1 then error("Each argument must be a column vector"), for vi in v do ( require_matrix(vi," ","orthogonal_complement"), if matrix_size(vi) # sz then error("Each vector must have the same dimension")), nullspace(apply('matrix, map('first, map('transpose, v))))); locate_matrix_entry(m, r1, c1, r2, c2, fn, rel) := block([im, cm, mf, e, nr, nc, ok : true, frel], require_matrix(m, "first", "locate_matrix_entry"), require_integer(r1, "second", "locate_matrix_entry"), require_integer(r2, "second", "locate_matrix_entry"), require_integer(c1, "third", "locate_matrix_entry"), require_integer(r2, "fourth", "locate_matrix_entry"), nr : length(m), nc : length(first(m)), if (c1 > c2) or (r1 > r2) or (r1 < 1) or (c1 < 1) or (r2 < 1) or (c2 < 1) or (c1 > nc) or (c2 > nc) or (r1 > nr) or (r2 > nr) then error("Bogus submatrix"), mf : if rel = 'min then 'inf else if rel = 'max then 'minf else 0, frel : if rel = 'max then ">" else if rel = 'min then "<" else if rel = 'bool then lambda([x,y],x) else error("Last argument must be 'max' or 'min'"), im : 0, cm : 0, for i : r1 thru r2 while ok do ( for j : c1 thru c2 while ok do ( e : apply(fn, [m[i,j]]), if is(apply(frel, [e, mf])) then ( im : i, cm : j, if rel = 'bool then ok : false, mf : e))), if im > 0 and cm > 0 then [im, cm] else false); columnspace(a) := block([m, nc, nr, cs : [], pos, x, proviso : [], algebraic : true], require_matrix(a, "first","columnspace"), /* request_rational_matrix(a, "first", "columnspace"), */ nc : length(first(a)), nr : length(a), if nc = 0 or nr = 0 then error("The argument to 'columnspace' must be a matrix with one or more rows and columns"), m : ptriangularize_with_proviso(a, ?gensym()), proviso : second(m), m : first(m), for i : 1 thru nr do ( if listp(pos : locate_matrix_entry(m, i, 1, i, nc, lambda([x],x # 0), 'bool)) then ( if not(constantp(x : part(part(m, first(pos)),second(pos)))) then proviso : adjoin(x # 0,proviso), cs : cons(col(a,second(pos)),cs))), cs : sort(cs), if not(emptyp(proviso)) then print("Proviso: ",proviso), funmake('span,cs)); /* Maxima's rank function doesn't complain with symbolic matrices. I don't approve of rank(matrix([a,b],[c,d])) --> 2. */ rank(m) := length(columnspace(m)); rowswap(m,i,j) := block([n, p, r], require_matrix(m, "first", "rowswap"), require_integer(i, "second", "rowswap"), require_integer(j, "third", "rowswap"), n : length(m), if (i < 1) or (i > n) or (j > n) then error("Array index out of bounds"), p : copymatrix(m), r : p[i], p[i] : p[j], p[j] : r, p); columnswap(m,i,j) := transpose(rowswap(transpose(m),i,j)); /* row(m,i) <-- row(m,i) - theta * row(m,i) */ rowop(m,i,j,theta) := block([p : copymatrix(m)], p[i] : p[i] - theta * p[j], p); /* col(m,i) <-- col(m,i) - theta * col(m,i) */ columnop(m,i,j,theta) := transpose(rowop(transpose(m),i,j,theta)); hipow_gzero(e,x) := block([n : hipow(e,x)], if n > 0 then n else 'inf); good_pivot(e,x) := freeof(x,e) and e # 0; ptriangularize(m,v) := block([p : ptriangularize_with_proviso(m,v)], if not(emptyp(p[2])) then print("Proviso: ",p[2]), p[1]); ptriangularize_with_proviso(m,v) := block([nr, nc, proviso : [], mp], require_matrix(m, "first", "ptriangularize"), require_symbol(v, "second", "ptriangularize"), nr : length(m), nc : length(first(m)), for i : 1 thru min(nr, nc) do ( mp : column_reduce(m,i,v), m : first(mp), proviso : append(proviso, second(mp))), proviso : setify(proviso), [expand(m), proviso]); column_reduce(m,i,x) := block([nc, nr, pos, q, proviso : []], /* require_matrix(m, "first", "column_reduce"), */ nc : length(first(m)), nr : length(m), /* require_integer(i, "second", "column_reduce"), if (i < 1) or (i > (length(first(m)))) then error("Bogus matrix column"), */ while not(good_pivot(m[i,i],x)) and (i+1 <= nr) and listp(pos : locate_matrix_entry(m,i+1,i,nr,i, lambda([e],e # 0), 'bool)) do ( pos : locate_matrix_entry(m,i,i,nr,i,lambda([e], good_pivot(e,x)), 'bool), if listp(pos) then ( m : rowswap(m,i,pos[1])) else ( m : expand(m), pos : locate_matrix_entry(m, i, i, nr, i, lambda([e], hipow(e,x)), 'max), m : rowswap(m,i,pos[1]), pos : locate_matrix_entry(m, i+1, i, nr, i, lambda([e], hipow_gzero(e,x)), 'min), if listp(pos) then ( q : divide(m[i,i], m[pos[1],i],x), m : rowop(m,i,pos[1],first(q))))), pos : locate_matrix_entry(m, i, i, nr, i, lambda([e], good_pivot(e,x)), 'bool), if listp(pos) then ( m : rowswap(m,i,first(pos)), for j : i+1 thru nr do ( if not(constantp(ratnumer(m[i,i]))) then ( proviso : cons(ratnumer(m[i,i]) # 0, proviso)), if not(constantp(ratdenom(m[i,i]))) then ( proviso : cons(ratdenom(m[i,i]) # 0, proviso)), m : rowop(m,j,i,m[j,i]/m[i,i]))) else ( pos : locate_matrix_entry(m, i, i, nr, i, lambda([e], e # 0), 'bool), if listp(pos) then m : rowswap(m,i, first(pos))), [m, proviso]); zeromatrixp(m) := if matrixp(m) then every(lambda([e],is(equal(rectform(e),0))), apply('append,args(m))) else error("Argument to 'zeromatrixp' must be a matrix"); mat_norm(m, p) := block([listarith : true, d], require_matrix(m,"first","mat_norm"), if p = 1 then mat_norm(transpose(m), 'inf) else if p = 'inf then ( d : 0, for ri in m do ( d : max(d, apply("+", map('cabs, ri)))), d) else if p = 'frobenius then apply("+", apply('append, args(cabs(m)^2))) else error("Not able to compute the ",p," norm")); mat_unblocker(m):=block( require_matrix(m, "first", "mat_unblocker"), if every('identity, apply('append, args(matrixmap('matrixp, m)))) then ( m : apply('addrow, args( map(lambda([x], apply('addcol, x)), m)))) else m); kronnecker_product(a,b) := ( require_matrix(a,"first", "kronnecker_product"), require_matrix(b,"second", "kronnecker_product"), mat_unblocker(outermap('matrix_element_mult, a,b))); hilbert_matrix(n) := ( require_posinteger(n,"first","hilbert_matrix"), funmake('matrix, makelist(makelist(1/(i + j - 1),i,1, n),j, 1, n))); /* To avoid a "too many arguments" error, use use funmake('matrix, ....) instead of apply('matrix, ...). */ vandermonde_matrix(vars) := block([n], require_list(vars,"first","vandermonde_matrix"), n : length(vars), funmake('matrix, makelist(makelist(vars[i]^j,j,0,n-1),i,1,n))); polytocompanion(p,x) := block([n], if not POLYNOMIALP(p,[x], lambda([e], freeof(x,e))) then error("First argument to 'polytocompanion' must be a polynomial"), p : expand(p), n : hipow(p,x), p : multthru(p / coeff(p,x,n)), genmatrix(lambda([i,j], if i-j=1 then 1 else if j = n then -coeff(p,x,i-1) else 0),n)); --- NEW FILE: linalg.usage --- Installation To use linalg, you'll need to have nset version 1.203. If you are using Maxima 5.9.1cvs or better, you already have the required version of nset. Otherwise you'll need to download and install nset. You may download nset from http://www.unk.edu/facstaff/profiles/willisb/. You will need to append the path to linearalgebra to file_search_maxima and file_search_lisp. If the path to the directory linearalgebra is "c:/maxima/linearalgebra", you will use the Maxima commands (%i1) file_search_maxima : cons("c:/maxima/linearalgebra/###.{mac}",file_search_maxima)\$ (%i2) file_search_lisp : cons("c:/maxima/linearalgebra/###.{lisp}",file_search_lisp)\$ Now you should be able to load and use linalg. (%i3) load("linalg"); Warning - you are redefining the MACSYMA function EIGENVALUES Warning - you are redefining the MACSYMA function EIGENVECTORS Warning - you are redefining the MACSYMA function RANK (%o3) c:/maxima/linearalgebra/linalg.mac (%i4) m : matrix([1-z,2],[5,8-z]); (%o4) matrix([1-z,2],[5,8-z]) (%i5) ptriangularize(m,z); (%o5) matrix([5,8-z],[0,-z^2/5+(9*z)/5+2/5]) columnop(m,i,j,theta) If m is a matrix, return the matrix that results from doing the column operation Ci <- Ci - theta * Cj. If m doesn't have a row i or j, signal an error. columnswap(m,i,j) If m is a matrix, swap columns i and j. If m doesn't have a column i or j, signal an error. columnspace(m) If m is a matrix, return span(v1,v2,...,vn), where the set {v1,v2,...,vn} is a basis for the column space of m. The span of the empty set is {0}. Thus, when the column space has only one member, return span(). ctranspose(m) Return the conjugate of the transpose of the matrix m. If m is a block matrix, transpose and conjugate each block. dotproduct(u,v) Return the dotproduct of vectors u and v. This is the same as conj(transpose(u)) . v. The arguments u and v must be n x 1 matrices. get_lu_factors(x) When x = lu_factor(A), then get_lu_factors returns a list of the form [P, L, U], where P is a permutation matrix, L is lower triangular with ones on the digonal, and U is upper triangular. And A = P L U. hilbert_matrix(n) Return the n by n Hilbert matrix. When n isn't a nonnegative integer, signal an error. kronnecker_product(a,b) Return the Kronnecker product of the matrices a and b. listp (e, {pred}) Given an optional argument pred, return true if e is a Maxima list and pred evaluates to true for every list element. When listp is not given the optional argument, return true if e is a Maxima list. In all other cases, return false. locate_matrix_entry(m, r1, c1, r2, c2, fn, rel) The first argument must be a matrix; the arguments r1 through c2 determine a sub-matrix of m that consists of rows r1 through r2 and columns c1 through c2. Find a entry in the sub-matrix m that satisfies some property. Three cases: (1) rel = bool and fn a predicate: Scan the sub-matrix from left to right then top to bottom, and return the index of the first entry that satisfies the predicate fn. If no matrix entry satisfies fn, return false. (2) rel = 'max and fn real-valued: Scan the sub-matrix looking for an entry that maximizes fn. Return the index of a maximizing entry. (3) rel = 'min and fn real-valued: Scan the sub-matrix looking for an entry that minimizes fn. Return the index of a minimizing entry. lu_backsub(m, b fld) When m = lu_factor(A), then lu_backsub(m,b, fld) solves the linear systme A x = b. lu_factor(m, generalfield) Return a list of the form [LU, perm, {cnd}], where (1) The matrix LU contains the factorization of m in a packed form. Packed form means three things: First, the rows of LU are permuted according to the list perm. If, for example, perm is the list [3,2,1], the actual first row of the LU factorization is the third row of the matrix LU. Second, the lower triangular factor of m is the lower triangular part of LU with the diagonal entries replaced by all ones. Third, the upper triangular factor of m is the upper triangular part of LU. (2) When the field is either 'floatfield' or 'complexfloatfield,' the number 'cnd' is an upper bound for the infinity norm condition number of m. For all other fields, the condition number isn't estimated. For these fields, lu_factor returns a two item list. The argument m must be a square matrix. The argument 'fld' must be a symbol that determines a field. The pre-defined fields are: (a) generalfield -- the field of Maxima expressions, (b) floatfield -- the field of floating point numbers of the type double, (c) complexfloatfield -- the field of complex floating point numbers of the type double, (d) crefield -- the field of Maxima CRE expressions, (e) rationalfiled -- the field of rational numbers. When the field is floatfield or complexfloatfield, the algorithm uses partial pivoting; when the field is generalfield, rows are switched only when needed to avoid a zero pivot. Floating point addition arithmetic isn't associative, so the meaning of 'field' differs from the mathematical definition. There is no user-interface for defining a new field. A user that is familiar with Common Lisp should be able to define a new field. To compute the factorization, the first task is to convert each matrix entry to a member of the indicated field. When conversion isn't possible, the factorization halts with an error message. Members of the field needn't be Maxima expressions. Members of the complexfloatfield, for example, are Common Lisp complex numbers. Thus after computing the factorization, the matrix entries must be converted to Maxima expressions. See also get_lu_factors. Examples (%i1) w[i,j] := ?random(1.0) + %i * ?random(1.0)\$ Evaluation took 0.00 seconds (0.00 elapsed) (%i2) m : genmatrix(w,100,100)\$ Evaluation took 00.30 seconds (00.30 elapsed) (%i3) lu_factor(m,complexfloatfield)\$ Evaluation took 1.45 seconds (1.45 elapsed) (%i4) lu_factor(m,generalfield)\$ Evaluation took 12.10 seconds (12.10 elapsed) (%i1) m : matrix([1-z,3],[3,8-z]); (%o1) matrix([1-z,3],[3,8-z]) (%i2) lu_factor(m,generalfield); (%o2) [matrix([1-z,3],[3/(1-z),-z-9/(1-z)+8]),[1,2]] (%i3) get_lu_factors(%); (%o3) [matrix([1,0],[0,1]),matrix([1,0],[3/(1-z),1]),matrix([1-z,3],[0,-z-9/(1-z)+8])] (%i4) %[1] . %[2] . %[3]; (%o4) matrix([1-z,3],[3,8-z]) mat_norm(m, p) Return the matrix p-norm of the matrix m. The allowed values for p are 1, inf, and frobenius (the Frobenius matrix norm). The matrix m should be an unblocked matrix. matrixp (e, {pred}) Given an optional argument pred, return true if e is a matrix and pred evaluates to true for every matrix element. When matrixp is not given an optional argument, return true if e is a matrix. In all other cases, return false. nonnegintegerp (n) Return true if and only if n >= 0 and n is an integer. nullspace(m) If m is a matrix, return span(v1,v2,...,vn), where the set {v1,v2,...,vn} is a basis for the nullspace of m. The span of the empty set is {0}. Thus, when the nullspace has only one member, return span(). nullity(m) If m is a matrix, return the dimension of the nullspace of m. orthogonal_complement(v1,v2,...,vn) Return span(u1,u2,...,um), where the set {u1,u2,...,um} is a basis for the orthogonal complement of the set(v1,v2,...,vn). Each vector v1 through vn must be a n x 1 matrix. polynomialp(p, l, {(coeffp 'constantp), (exponp 'nonnegintegerp)}) Return true if p is a polynomial in the variables in the Maxmia list l, The predicate coeffp (default constantp) must evaluate to true for each coefficient, and the predicate exponp must evaluate to true for all exponents of the variables in the list l. If you want to use a non-default value for exponp, you must supply coeffp with a value even if you want to use the default for coeffp. The polynomial needn't be expanded: (%i1) polynomialp((x+1)*(x+2),[x]); (%o1) true (%i2) polynomialp((x+1)*(x+2)^a,[x]); (%o2) false An example using non-default values for coeffp and exponp: (%i3) polynomialp((x+1)*(x+2)^(3/2),[x],numberp,numberp); (%o3) true (%i4) polynomialp((x^(1/2)+1)*(x+2)^(3/2),[x],numberp,numberp); (%o4) true Polynomials with two variables: (%i5) polynomialp(x^2 + 5*x*y + y^2,[x]); (%o5) false (%i6) polynomialp(x^2 + 5*x*y + y^2,[x,y]); (%o6) true polytocompanion(p,x) If p is a polynomial in x, return the companion matrix of p. For a monic polynomial p of degree n, we have p = (-1)^n charpoly(polytocompanion(p, x)). When p isn't a polynomial in x, signal an error. ptriangularize(m,v) If m is a matrix with each entry a polynomial in v, return a matrix m' such that (1) m' is upper triangular, (2) m' = E_n ... E_1 m, where E1 through En are elementary matrices whose entries are polynomials in v, (3) |det(m)| = |det(m')|, Note: This function doesn't check that every entry is a polynomial in v. rowop(m,i,j,theta) If m is a matrix, return the matrix that results from doing the row operation Ri <- Ri - theta * Rj. If m doesn't have a row i or j, signal an error. rowswap(m,i,j) If m is a matrix, swap rows i and j. If m doesn't have a row i or j, signal an error. mat_unblocker(m) If m is a block matrix, return its unblocked form. If m is a matrix, return m; otherwise, signal an error. If you use block matrices, most likely you'll want to set matrix_element_mult : "." and matrix_element_transpose : 'transpose. Example (%i1) a : matrix([1,2],[3,4])\$ (%i2) b : matrix([7,8],[9,10])\$ (%i3) matrix([a,b]); (%o3) matrix([matrix([1,2],[3,4]),matrix([7,8],[9,10])]) (%i4) mat_unblocker(%); (%o4) matrix([1,2,7,8],[3,4,9,10]) vandermonde_matrix([x1,x2,...,xn]) Return a n by n matrix whose i-th row is [1, xi, xi^2, ... xi^(n-1)]. zeromatrixp(m) If is(equal(e,0)) is true for each entry of the matrix m, return true; otherwise, return false. When m is not a matrix, signal an error. Short demo: (%i1) batch("linalg.demo"); batching #pc:/maxima/linearalgebra/linalg.demo (%i2) LOAD(linalg) (%o2) c:/maxima/linearalgebra/linalg.mac (%i3) m:matrix([1,2],[1,2]) (%o3) matrix([1,2],[1,2]) (%i4) nullspace(m) (%o4) span(matrix([1],[-1/2])) (%i5) columnspace(m) (%o5) span(matrix([1],[1])) (%i6) ptriangularize(m-z*IDENT(2),z) (%o6) matrix([1,2-z],[0,3*z-z^2]) (%i7) m:matrix([1,2,3],[4,5,6],[7,8,9])-z*IDENT(3) (%o7) matrix([1-z,2,3],[4,5-z,6],[7,8,9-z]) (%i8) mm:ptriangularize(m,z) (%o8) matrix([4,5-z,6],[0,66/49,-z^2/7+(102*z)/49+132/49],[0,0,(49*z^3)/264-(245*z^2)/88-(147*z)/44]) (%i9) ALGEBRAIC:true (%o9) true (%i10) TELLRAT(mm[3,3]) (%o10) [z^3-15*z^2-18*z] (%i11) mm:RATSIMP(mm) (%o11) matrix([4,5-z,6],[0,66/49,-(7*z^2-102*z-132)/49],[0,0,0]) (%i12) nullspace(mm) (%o12) span(matrix([1],[(z^2-14*z-16)/8],[-(z^2-18*z-12)/12])) (%i13) m:matrix([1,2,3,4],[5,6,7,8],[9,10,11,12],[13,14,15,16]) (%o13) matrix([1,2,3,4],[5,6,7,8],[9,10,11,12],[13,14,15,16]) (%i14) column_space(m) (%o14) column_space(matrix([1,2,3,4],[5,6,7,8],[9,10,11,12],[13,14,15,16])) (%i15) APPLY('orthogonal_complment,ARGS(nullspace(TRANSPOSE(m)))) (%o15) orthogonal_complment(matrix([0],[1],[-2],[1]),matrix([1],[0],[-3],[2])) --- NEW FILE: lu.lisp --- ;; Copyright 2005 by Barton Willis ;; This is free software; you can redistribute it and/or ;; modify it under the terms of the GNU General Public License, ;; http://www.gnu.org/copyleft/gpl.html. ;; This software has NO WARRANTY, not even the implied warranty of ;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. (\$put '\$lu 1 '\$version) (defun \$mytest (fn) (let ((*collect-errors* nil)) (setq fn (\$file_search1 fn '((mlist) \$file_search_maxima))) (test-batch fn nil :show-all t :show-expected t))) (defstruct mfield coerce-to-lisp-float abs great add div mult sub add-id mult-id max fzerop maxima-to-mfield mfield-to-maxima) (defun \$require_field (f pos fun) (if (not (mfield-p f)) (merror "The ~:M argument of the function ~:M must be an ordered field" pos fun))) (defun list-swap (p i j) (let ((x (nth i p))) (setf (nth i p) (nth j p)) (setf (nth j p) x))) (defparameter \$floatfield (make-mfield :coerce-to-lisp-float #'(lambda (s) s) :abs #'abs :great #'> :add #'+ :div #'/ :mult #'* :sub #'- :add-id 0.0 :mult-id 1.0 :max #'max :fzerop #'(lambda (s) (= s 0.0)) :mfield-to-maxima #'(lambda (s) s) :maxima-to-mfield #'(lambda (s) (setq s (\$float s)) (if (floatp s) s (merror "Unable to convert matrix entry to a float"))))) (defparameter \$complexfloatfield (make-mfield :coerce-to-lisp-float #'(lambda (s) s) :abs #'abs :great #'> :add #'+ :div #'/ :mult #'* :sub #'- :add-id 0.0 :mult-id 1.0 :max #'max :fzerop #'(lambda (s) (= s 0.0)) :mfield-to-maxima #'(lambda (s) (add (realpart s) (mul '\$%i (imagpart s)))) :maxima-to-mfield #'(lambda (s) (setq s (\$rectform s)) (if (and (floatp (\$float (\$realpart s))) (floatp (\$float (\$imagpart s)))) (complex (\$float (\$realpart s)) (\$float (\$imagpart s))) (merror "Unable to convert matrix entry to a complex float"))))) (defparameter \$rationalfield (make-mfield :coerce-to-lisp-float #'(lambda (s) (\$float s)) :abs #'abs :great #'> :add #'+ :div #'/ :mult #'* :sub #'- :add-id 0 :mult-id 1 :max #'max :fzerop #'(lambda (s) (= s 0)) :mfield-to-maxima #'(lambda (s) (simplify `((rat) ,(numerator s) ,(denominator s)))) :maxima-to-mfield #'(lambda (s) (if (\$ratnump s) (if (integerp s) s (/ (\$num s) (\$denom s))) (merror "Unable to convert matrix entry to a rational number"))))) (defparameter \$crefield (make-mfield :coerce-to-lisp-float nil :abs #'(lambda (s) (simplify (mfuncall '\$cabs s))) :great #'(lambda (a b) (like b 0)) :add #'add :div #'div :mult #'mult :sub #'sub :add-id 0 :mult-id 1 :max #'(lambda (a b) (maximum (list a b))) :fzerop #'(lambda (s) (like s 0)) :mfield-to-maxima #'(lambda (s) s) :maxima-to-mfield #'(lambda (s) (\$rat s)))) (defparameter \$generalfield (make-mfield :coerce-to-lisp-float nil :abs #'(lambda (s) (simplify (mfuncall '\$cabs s))) :great #'(lambda (a b) (like b 0)) :add #'(lambda (a b) (\$rectform (add a b))) :div #'(lambda (a b) (\$rectform (div a b))) :mult #'(lambda (a b) (\$rectform (mult a b))) :sub #'(lambda (a b) (\$rectform (sub a b))) :add-id 0 :mult-id 1 :max #'(lambda (a b) (maximum (list a b))) :fzerop #'(lambda (s) (like s 0)) :mfield-to-maxima #'(lambda (s) s) :maxima-to-mfield #'(lambda (s) s))) (defparameter \$bigfloatfield (make-mfield :coerce-to-lisp-float #'(lambda (s) (setq s (\$rectform (\$float s))) (complex (\$realpart s) (\$imagpart s))) :abs #'(lambda (s) (simplify (mfuncall '\$cabs s))) :great #'mgrp :add #'(lambda (a b) (\$rectform (add a b))) :div #'(lambda (a b) (\$rectform (div a b))) :mult #'(lambda (a b) (\$rectform (mult a b))) :sub #'(lambda (a b) (\$rectform (sub a b))) :add-id 0 :mult-id 1 :max #'(lambda (a b) (maximum (list a b))) :fzerop #'(lambda (s) (like s bigfloatzero)) :mfield-to-maxima #'(lambda (s) s) :maxima-to-mfield #'(lambda (s) (setq s (\$bfloat s)) (if (\$bfloatp s) s (merror "Unable to convert matrix entry to a big float"))))) ;; Map the lisp function fn over the r by c Maxima matrix m. (defun matrix-map (m r c fn) (loop for i from 1 to r do (loop for j from 1 to c do (setmatelem m (funcall fn (nth j (nth i m))) i j)))) ;; Return the i,j entry of the Maxima matrix m. The rows of m have been permuted according ;; to the Maxima list p. (defun m-elem (m p i j) (nth j (nth (nth i p) m))) ;; Set the i j entry of the Maxima matrix m to x. (defun setmatelem (m x i j) (setf (nth j (nth i m)) x)) ;; Return m[perm[i], k] - sum(m[perm[i],s] * m[perm[s],k],s,1,n) (defun partial-matrix-prod (m p i k n fadd fsub fmult add-id) (loop for s from 1 to n do (setq add-id (funcall fadd add-id (funcall fmult (m-elem m p i s) (m-elem m p s k))))) (setmatelem m (funcall fsub (m-elem m p i k) add-id) (nth i p) k)) (defun \$lu_factor (mm &optional (fld \$generalfield)) (\$require_square_matrix mm "\$first" "\$lu_factor") (\$require_nonempty_matrix mm "\$first" "\$lu_factor") (\$require_field fld "\$second" "\$lu_factor") (let ((m (copy-tree mm)) (c (\$length mm)) (perm)) (matrix-map m c c (mfield-maxima-to-mfield fld)) (loop for i from 1 to c do (push i perm)) (setq perm (reverse perm)) (push '(mlist) perm) (lu-factor m perm c fld))) (defun \$get_lu_factors (x) (let ((mat (\$first x)) (mp) (p (\$second x)) (perm) (r) (c) (id) (lower) (upper)) (setq r (\$matrix_size mat)) (setq c (\$second r)) (setq r (\$first r)) (setq id (\$args (\$ident c))) (loop for i from 1 to c do (push (nth (nth i p) id) perm) (push (nth (nth i p) mat) mp)) (setq perm (reverse perm)) (setq mp (reverse mp)) (push '(\$matrix) perm) (push '(\$matrix) mp) (setq lower (copy-tree mp)) (setq upper (copy-tree mp)) (loop for i from 1 to r do (loop for j from 1 to c do (if (= i j) (setmatelem lower 1 i j)) (if (> i j) (setmatelem upper 0 i j)) (if (< i j) (setmatelem lower 0 i j)))) `((mlist) ,(\$transpose perm) ,lower ,upper))) (defun lu-factor (m perm c fld) (let ((pos) (kp1) (mx) (cnd) (fadd (mfield-add fld)) (fsub (mfield-sub fld)) (fmult (mfield-mult fld)) (fdiv (mfield-div fld)) (fabs (mfield-abs fld)) (fgreat (mfield-great fld)) (fzerop (mfield-fzerop fld)) (add-id (mfield-add-id fld))) (loop for k from 1 to c do (loop for i from k to c do (partial-matrix-prod m perm i k (- k 1) fadd fsub fmult add-id)) (setq mx (funcall fabs (m-elem m perm k k))) (setq pos k) (loop for s from k to c do (if (funcall fgreat (funcall fabs (m-elem m perm s k)) mx) (setq pos s mx (funcall fabs (m-elem m perm s k))))) (list-swap perm k pos) (setq kp1 (+ 1 k)) (loop for i from kp1 to c do (if (funcall fzerop (m-elem m perm k k)) (merror "Unable to compute the LU factorization")) (setmatelem m (funcall fdiv (m-elem m perm i k) (m-elem m perm k k)) (nth i perm) k) (partial-matrix-prod m perm k i (- k 1) fadd fsub fmult add-id))) (cond ((not (eq nil (mfield-coerce-to-lisp-float fld))) (setq cnd (mat-cond-by-lu m perm c (mfield-coerce-to-lisp-float fld))) (matrix-map m c c (mfield-mfield-to-maxima fld)) `((mlist) ,m ,perm ,cnd)) (t `((mlist) ,m ,perm))))) (defun row-sum (mat i start stop fn) (let ((acc 0.0)) (setq mat (nth i mat)) (loop for i from start to stop do (setq acc (+ acc (abs (funcall fn (nth i mat)))))) acc)) (defun mat-cond-by-lu (mat perm n fn) (let ((l-norm) (l-cond) (u-norm) (uu-norm) (u-cond) (d)) (setq l-norm 0.0) (loop for i from 2 to n do (setq l-norm (max l-norm (row-sum mat (nth i perm) 1 (- i 1) fn)))) (setq l-cond (* (+ 1 l-norm) (/ (- 1 (expt l-norm n)) (- 1 l-norm)))) (setq u-norm 0.0) (setq uu-norm 0.0) (loop for i from 1 to n do (setq u-norm (max u-norm (row-sum mat (nth i perm) i n fn))) (setq d (abs (funcall fn (m-elem mat perm i i)))) (cond ((= d 0.0) (setq uu-norm '\$inf) (return nil))) (setq uu-norm (max uu-norm (/ (row-sum mat (nth i perm) (+ i 1) n fn) d)))) (cond ((eq uu-norm '\$inf) '\$inf) (t (setq u-cond (* uu-norm (/ (- 1 (expt uu-norm n)) (- 1 uu-norm)))) (* l-cond u-cond))))) (defun \$lu_backsub(m b1 &optional (fld \$generalfield)) (\$require_list m "\$first" "\$lu_backsub") (\$require_field fld "\$second" "\$lu_backsub") (let ((mat) (n) (r) (c) (bb) (acc) (perm) (id-perm) (b) (fadd (mfield-add fld)) (fsub (mfield-sub fld)) (fmult (mfield-mult fld)) (fdiv (mfield-div fld)) (add-id (mfield-add-id fld))) (setq mat (copy-tree (\$first m))) (setq perm (\$second m)) (setq n (\$matrix_size mat)) (setq r (\$first n)) (setq c (\$second n)) (matrix-map mat r c (mfield-maxima-to-mfield fld)) (setq b (copy-tree b1)) (matrix-map b r 1 (mfield-maxima-to-mfield fld)) (setq bb (copy-tree b)) (loop for i from 1 to r do (setmatelem bb (m-elem b perm i 1) i 1)) (setq id-perm (mfuncall '\$makelist 'i 'i 1 r)) (loop for i from 2 to c do (setq acc add-id) (loop for j from 1 to (- i 1) do (setq acc (funcall fadd acc (funcall fmult (m-elem mat perm i j) (m-elem bb id-perm j 1))))) (setmatelem bb (funcall fsub (m-elem bb id-perm i 1) acc) i 1)) (loop for i from r downto 1 do (setq acc (m-elem bb id-perm i 1)) (loop for j from (+ 1 i) to c do (setq acc (funcall fsub acc (funcall fmult (m-elem mat perm i j) (m-elem bb id-perm j 1))))) (setmatelem bb (funcall fdiv acc (m-elem mat perm i i)) i 1)) (matrix-map bb r 1 (mfield-mfield-to-maxima fld)) bb)) --- NEW FILE: matrixexp.lisp --- ;; Copyright 2004 by Barton Willis ;; This is free software; you can redistribute it and/or ;; modify it under the terms of the GNU General Public License, ;; http://www.gnu.org/copyleft/gpl.html. ;; This software has NO WARRANTY, not even the implied warranty of ;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. ;; Innocent looking problems, such as matrixexp(matrix([x,1,0],[1,1,1],[0,1,1])), ;; generate huge (and most likely worthless) expressions. These huge ;; expressions strain Maxima's rational function code; to avoid errors ;; such as "..quotient by polynomial of higher degree" I found it necessary to ;; ;; (1) set gcd == spmod and algebraic == true, ;; (2) always give ratsimp and fullratsimp a value (even if nil) for its ;; optional second argument, ;; (3) set ratvars to an empty list at the start of most functions. ;; I didn't try to find the real cause for these bugs. The function ;; spectral_rep does check its output. Although these checks are slow, ;; I recommend that they stay. ;; Map the CL function f over the elements of a Maxima matrix mat and ;; return a Maxima matrix. The first argument should be a CL function. (\$put '\$matrixexp 1 '\$version) (defun my-matrix-map (f mat) (setq mat (mapcar 'cdr (cdr mat))) (setq mat (mapcar #'(lambda (s) (mapcar f s)) mat)) (setq mat (mapcar #'(lambda (s) (cons `(mlist simp) s)) mat)) `((\$matrix simp) ,@mat)) ;; Signal an error if 'mat' is not a square matrix; otherwise, return true. (defun require-square-matrix (mat pos fun-name) (if (not (\$matrixp mat)) (merror "The ~:M argument to ~:M must be a matrix" pos fun-name)) (if (not (= (\$length (\$args mat)) (\$length (\$first (\$args mat))))) (merror "The ~:M argument to ~:M must be a square matrix" pos fun-name))) ;; When mat is a square matrix, return exp(mat * x). The second ;; argument is optional and it defaults to 1. (defun \$matrixexp (mat &optional (x 1)) (let ((sp) (d) (p) (id) (n (\$length (\$args mat))) (f)) (\$ratvars) (require-square-matrix mat "\$first" "\$matrixexp") (setq mat (\$spectral_rep mat)) (setq sp (\$first mat)) (setq p (\$second mat)) (setq sp (cons '(mlist simp) (mapcar #'(lambda (s) (simplify (\$exp (mult s x)))) (cdr sp)))) (setq d (mult x (\$third mat))) (setq id (\$ident n)) (setq f id) (setq n (+ n 1)) (dotimes (i n) (setq f (add id (div (ncmul2 d f) (- n i))))) (\$fullratsimp (ncmul2 (ncmul2 sp p) f) nil))) ;; Let f(var) = expr. This function returns f(mat), where 'mat' is a ;; square matrix. Here expr is an expression---it isn't a function! (defun require-lambda (e n pos fun-name) (let ((var)) (if (and (consp e) (consp (car e)) (eq 'lambda (mop e)) (= 3 (length e)) (\$listp (nth 1 e)) (setq var (cdr (nth 1 e))) (= n (length var)) (every #'(lambda (s) (or (symbolp s) (\$subvarp s))) var)) (list var (nth 2 e)) (merror "The ~:M argument to `~:M' must be a lambda form with ~:M variable(s)" pos fun-name n)))) (defun \$matrixfun (lamexpr mat) (let ((z (gensym)) (expr) (var) (sp) (d) (p) (di) (n (\$length (\$args mat))) (f 0)) (require-square-matrix mat "\$second" "\$matrixexp") (setq expr (require-lambda lamexpr 1 "\$first" "\$matrixfun")) (setq var (nth 0 (nth 0 expr))) (setq expr (nth 1 expr)) (\$ratvars) (setq expr (\$substitute z var expr)) (setq mat (\$spectral_rep mat)) (setq sp (\$first mat)) (setq p (\$second mat)) (setq d (\$third mat)) (setq di (\$ident n)) (setq sp (cdr sp)) (dotimes (i (+ n 1)) (setq f (add f (ncmul2 di (ncmul2 (cons '(mlist simp) (mapcar #'(lambda (s) (\$substitute s z expr)) sp)) p)))) (setq di (ncmul2 di d)) (setq expr (div (\$diff expr z) (factorial (+ i 1))))) (\$fullratsimp (simplify f) nil))) ;; Return the residue of the rational expression e with respect to the ;; variable var at the point pt. Assumptions: ;; (1) the denominator of e divides ker, ;; (2) e is a rational expression, ;; (3) ker is a polynomial, ;; (4) pt is z zero of ker and ord is its order. (defun rational-residue (e var pt ker ord) (let ((\$gcd '\$spmod) (\$algebraic t) (\$ratfac nil) (p) (q) (f (sub var pt))) (\$ratvars) (setq e (\$fullratsimp e var)) (setq p (\$num e)) (setq q (\$denom e)) (setq p (mult p (\$quotient ker q var))) (setq e (\$fullratsimp (div p (\$quotient ker (power f ord) var)) var)) (\$fullratsimp (\$substitute pt var (div (\$diff e var (- ord 1)) (factorial (- ord 1)))) nil))) (defun \$spectral_rep (mat) (require-square-matrix mat "\$first" "\$spectral_rep") (let ((\$gcd '\$spmod) (\$algebraic t) (\$resultant '\$subres) (ord) (zi) (\$ratfac nil) (z (gensym)) (res) (m) (n (\$length (\$args mat))) (p) (p1) (p2) (sp) (proj)) (\$ratvars) (setq p (\$newdet (sub mat (mult z (\$ident n))))) (if (oddp n) (setq p (mult -1 p))) ;; p1 = (z - z1)(z - z2) ... (z - zk), where z1 thru zk are ;; the distinct zeros of p. (setq p1 (\$first (\$divide p (\$gcd p (\$diff p z) z) z))) (setq p2 (\$resultant p1 (\$diff p1 z) z)) (cond ((and (not (\$constantp p2)) (not (like 0 p2))) (\$ratvars) (setq p2 (\$sqfr p2)) (if (mminusp p2) (setq p2 (mult -1 p2))) (setq p2 (if (mtimesp p2) (margs p2) (list p2))) (setq p2 (mapcar #'(lambda (s) (if (mexptp s) (nth 1 s) s)) p2)) (setq p2 `((mtimes simp) ,@p2)) (mtell "Proviso: assuming ~:M" `((mnotequal simp) ,p2 0)))) (setq sp (\$solve p z)) (setq sp (mapcar '\$rhs (cdr sp))) (cond ((not (eq n (apply #'+ (cdr \$multiplicities)))) (print `(ratvars = ,\$ratvars gcd = '\$gcd algebraic = ,\$algebraic)) (print `(ratfac = ,\$ratfac)) (merror "Unable to find the spectrum"))) (setq res (\$fullratsimp (ncpower (sub (mult z (\$ident n)) mat) -1) z)) (setq m (length sp)) (dotimes (i m) (setq zi (nth i sp)) (setq ord (nth (+ i 1) \$multiplicities)) (push (my-matrix-map #'(lambda (e) (rational-residue e z zi p ord)) res) proj)) (setq proj (nreverse proj)) (setq m (length proj)) (dotimes (i m) (setq mat (sub mat (mult (nth i sp) (nth i proj))))) (cond ((check-spectral-rep proj n) (push `(mlist simp) proj) (push `(mlist simp) sp) `((mlist simp) ,sp ,proj, (\$fullratsimp mat nil))) (t (merror "Unable to find the spectral representation"))))) (defun check-spectral-rep (proj n) (let* ((m (length proj)) (okay) (zip (\$zeromatrix n n)) (qi) (\$gcd '\$spmod) (\$algebraic t) (\$ratfac nil)) (\$ratvars) (setq proj (mapcar #'(lambda (s) (\$fullratsimp s nil)) proj)) (setq okay (like (\$ident n) (\$fullratsimp (apply 'add proj) nil))) (dotimes (i m) (setq qi (nth i proj)) (setq qi (\$fullratsimp (sub qi (ncmul2 qi qi)) nil)) (setq okay (and okay (like zip qi)))) (dotimes (i m) (setq qi (nth i proj)) (dotimes (j i) (setq okay (and okay (like zip (\$fullratsimp (ncmul2 qi (nth j proj)) nil)))) (setq okay (and okay (like zip (\$fullratsimp (ncmul2 (nth j proj) qi) nil)))))) okay)) --- NEW FILE: matrixexp.usage --- Every square matrix can be expressed as a linear combination of commuting projection matrices plus a nilpotent matrix. Specifically, given a square matrix M, there are scalars z1, z2, ... zn and matrices P1, P2, ... , Pn, D, such that (1) M = z1 P1 + z2 P2 + ... + zn Pn + D, (2) spectrum(M) = {z1,z2, ..., zn}, (3) Pi Pj = Pj Pi = kron_delta(i,j) Pi, for 1 <= i,j <= n, (4) D is nilpotent, (5) D Pi = Pi D, for 1 <= i <= n, (6) P1 + P2 + ... + Pn = I. We call (1) the spectral representation of M. Recall that a matrix is D is nilpotent means that there is a positive integer k such that D^k = 0. The spectral representation of a matrix is not a continuous function of the matrix arguments. Here is an example (C1) spectral_rep(matrix([1,0],[0,x])); Proviso: assuming x-1 # 0 [ 0 0 ] [ 1 0 ] [ 0 0 ] (D1) [[x, 1], [[ ], [ ]], [ ]] [ 0 1 ] [ 0 0 ] [ 0 0 ] Maxima does give a warning that it has assumed that x does not equal 1. Unfortunately, the assumption x-1 # 0 isn't available to the user. When x = 1, the spectral representation is (C2) spectral_rep(matrix([1,0],[0,1])); [ 1 0 ] [ 0 0 ] (D2) [[1], [[ ]], [ ]] [ 0 1 ] [ 0 0 ] (C3) Functions: spectral_rep(mat) When mat is a square matrix, return a list of the form [[z1,z2, ... zn], [P1, P2, ... ,Pn], D], where z1 P1 + z2 P2 + ... + zn Pn + D is the spectral representation of M. matrixexp(mat, t) When mat is a square matrix, return exp(mat * t). The second argument is optional and it defaults to 1. matrixfun(lambda-expr, mat) When mat is a square matrix and lambda-expr is a Maxima lambda expression, extend the function lambda-expr to matrices and return lambda-expr(mat). Thus matrixfun(lambda([x],x^2), mat) == mat . mat and matrixfun(lambda([x], 1/x), mat) == mat^^-1. If the lambda-expr involves an unknown function, you may need to load pdiff; otherwise, you'll get an error. Here is an example (C1) load("matexp.lisp")\$ (C2) m : matrix([1,2],[0,1])\$ (C3) matrixfun(lambda([x], f(x)),m); Attempt to differentiate with respect to a number: 1 -- an error. Quitting. To debug this try DEBUGMODE(TRUE);) (C4) load("pdiff")\$ (C5) matrixfun(lambda([x], f(x)),m); [ f(1) 2 f (1) ] (D5) [ (1) ] [ ] [ 0 f(1) ] (C6) Barton Willis wrote matexp. This code is covered by the GNU public license. Please send bug reports to the Maxima mailing list. --- NEW FILE: polynomialp.lisp --- ;; Author Barton Willis ;; University of Nebraska at Kearney ;; Copyright (C) 2004, 2005, Barton Willis ;; Brief Description: polynomial predicate function. ;; This program is free software; you can redistribute it and/or modify ;; it under the terms of the GNU General Public License as published by ;; the Free Software Foundation; either version 2 of the License, or ;; (at your option) any later version. ;; This program is distributed in the hope that it will be useful, ;; but WITHOUT ANY WARRANTY; without even the implied warranty of ;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the ;; GNU General Public License for more details. ;; You should have received a copy of the GNU General Public License ;; along with this program; if not, write to the Free Software ;; Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. (eval-when (:compile-toplevel :load-toplevel :execute) (if (not (functionp 'op-equalp)) (\$load "linalg-utilities")) (if (not (functionp 'require-list-or-set)) (\$load "nset"))) (\$put '\$polynomialp 1 '\$version) ;; Return true iff n is an integer and n >= 0. (defun \$nonnegintegerp (n) (and (integerp n) (>= n 0))) (defun \$polynomialp (p vars &optional (coeffp '\$constantp) (exponp '\$nonnegintegerp)) (setq vars (require-list-or-set vars "\$polynomialp")) (setq vars (mapcar '\$ratdisrep vars)) (if (every #'(lambda (s) (or (\$symbolp s) (\$subvarp s))) vars) (polynomialp (\$ratdisrep p) vars coeffp exponp) (merror "The second argument to polynomialp must be a list of symbols"))) (defun polynomialp (p vars coeffp exponp) (or (mfuncall coeffp p) (if (member p vars :test #'eq) t nil) (and (op-equalp p 'mtimes 'mplus) (every #'(lambda (s) (polynomialp s vars coeffp exponp)) (margs p))) (and (op-equalp p 'mexpt) (polynomialp (car (margs p)) vars coeffp exponp) (mfuncall exponp (cadr (margs p)))))) --- NEW FILE: test-conjugate.mac --- conjugate(-1); -1\$ conjugate(6); 6\$ conjugate(-2/3); -2/3\$ conjugate(sqrt(5)); sqrt(5)\$ conjugate(sqrt(3) + sqrt(7)); sqrt(3) + sqrt(7)\$ conjugate(5.6); 5.6\$ conjugate(-5.6); -5.6\$ conjugate(5.6 * %i); -%i * 5.6\$ conjugate(8.23b3); 8.23b3\$ conjugate(1.4b0 + %i * 5.7b0); 1.4b0 - %i * 5.7b0\$ conjugate(sqrt(7 + %i * 8)); sqrt(7 - %i * 8); conjugate(log(1 - %i)); log(1 + %i)\$ conjugate(%e); %e\$ conjugate(-%e); -%e\$ conjugate(%pi); %pi\$ conjugate(%phi); %phi\$ conjugate(x); x\$ conjugate(-x); -x\$ conjugate(1+x*(1+x)); 1+x*(1+x)\$ conjugate(sqrt(1+x^2)); sqrt(1+x^2)\$ conjugate(x+abs(x)); x+abs(x)\$ conjugate(cos(x)); cos(x)\$ conjugate(sin(x)); sin(x)\$ conjugate(tan(x)); tan(x)\$ conjugate(sec(x)); sec(x)\$ conjugate(cot(x)); cot(x)\$ conjugate(csc(x)); csc(x)\$ conjugate(cosh(x)); cosh(x)\$ conjugate(sinh(x)); sinh(x)\$ conjugate(tanh(x)); tanh(x)\$ conjugate(sech(x)); sech(x)\$ conjugate(coth(x)); coth(x)\$ conjugate(csch(x)); csch(x)\$ conjugate(exp(x)); exp(x)\$ conjugate(sqrt(x)); conjugate(sqrt(x))\$ conjugate(abs(x))\$ abs(x)\$ (declare(z,complex),0); 0\$ conjugate(conjugate(z)); z\$ conjugate(conjugate(conjugate(z))); conjugate(z)\$ conjugate(5*z); 5 * conjugate(z); conjugate(5 + z); 5 + conjugate(z); conjugate(5 + 7 * z); 5 + 7 * conjugate(z); conjugate(%i * z); -%i * conjugate(z)\$ conjugate(x + z); x + conjugate(z)\$ conjugate(x * z); x * conjugate(z)\$ conjugate(cos(z)); cos(conjugate(z))\$ conjugate(sin(z)); sin(conjugate(z))\$ conjugate(tan(z)); tan(conjugate(z))\$ conjugate(sec(z)); sec(conjugate(z))\$ conjugate(cot(z)); cot(conjugate(z))\$ conjugate(csc(z)); csc(conjugate(z))\$ conjugate(cosh(z)); cosh(conjugate(z))\$ conjugate(sinh(z)); sinh(conjugate(z))\$ conjugate(tanh(z)); tanh(conjugate(z))\$ conjugate(sech(z)); sech(conjugate(z))\$ conjugate(coth(z)); coth(conjugate(z))\$ conjugate(csch(z)); csch(conjugate(z))\$ conjugate(exp(z)); exp(conjugate(z))\$ conjugate(realpart(z)); realpart(z)\$ conjugate(imagpart(z)); imagpart(z)\$ conjugate(x+z); x+conjugate(z)\$ conjugate(x*z); x*conjugate(z)\$ conjugate(%i * z); -%i * conjugate(z)\$ conjugate(cos(x*z)); cos(x*conjugate(z))\$ (assume(imagpart(z) > 0),0); 0\$ conjugate(log(z)); log(conjugate(z))\$ conjugate(sqrt(z)); sqrt(conjugate(z))\$ conjugate(log(5 + z)); log(5 + conjugate(z))\$ (declare(w,complex),0); 0\$ (assume(realpart(w) > 0),0); 0\$ conjugate(log(w)); log(conjugate(w))\$ conjugate(sqrt(w)); sqrt(conjugate(w))\$ conjugate(log(7 + w)); log(7 + conjugate(w))\$ (declare(q,complex),0); 0\$ conjugate(q^2); conjugate(q)^2\$ conjugate(1/q); 1/conjugate(q)\$ conjugate(1 + q + q^2); 1 + conjugate(q) + conjugate(q)^2\$ conjugate(q*(1+q)); conjugate(q) * (1 + conjugate(q))\$ conjugate((1+q)/(1-q)); (1 + conjugate(q))/(1 - conjugate(q))\$ --- NEW FILE: test-linalg.mac --- /*--- ctranspose -----*/ (m : matrix([%i,1],[2,%i],[0,2]),0); 0\$ ctranspose(m); matrix([-%i,2,0],[1,-%i,2])\$ zeromatrixp(m - ctranspose(ctranspose(m))); true\$ (m : matrix([m,m]),0); 0\$ zeromatrixp(ctranspose(mat_unblocker(m)) - mat_unblocker(ctranspose(m))); true\$ ctranspose(matrix([%i])); matrix([-%i])\$ /*--- matrix norm ----*/ (m : matrix([2,3],[6,7]), 0); 0\$ mat_norm(m,1); 10\$ mat_norm(m,'inf); 13\$ mat_norm(m,'frobenius); 2^2 + 3^2 + 6^2 + 7^2\$ (m : matrix([2*%i,3],[6,7]), 0); 0\$ mat_norm(m,1); 10\$ mat_norm(m,'inf); 13\$ mat_norm(m,'frobenius); 2^2 + 3^2 + 6^2 + 7^2\$ /*--- matrix size ----*/ matrix_size(matrix([1])); [1,1]\$ matrix_size(matrix([false],[false])); [2,1]\$ matrix_size(matrix([false,false])); [1,2]\$ /*--- matrix rank and nullity -----*/ rank(matrix([0])); 0\$ nullity(matrix([0])); 1\$ rank(matrix([1])); 1\$ nullity(matrix([1])); 0\$ rank(matrix([0,0])); 0\$ nullity(matrix([0,0])); 2\$ rank(matrix([0,1])); 1\$ nullity(matrix([0,1])); 1\$ rank(matrix([0],[0])); 0\$ nullity(matrix([0],[0])); 1\$ rank(matrix([1],[0])); 1\$ nullity(matrix([1],[0])); 0\$ rank(matrix([0],[1])); 1\$ nullity(matrix([0],[1])); 0\$ rank(matrix([1],[1])); 1\$ nullity(matrix([1],[1])); 0\$ rank(matrix([0],[%i])); 1\$ nullity(matrix([0],[%i])); 0\$ rank(matrix([1],[1])); 1\$ nullity(matrix([1],[1])); 0\$ rank(matrix([1],[%i])); 1\$ nullity(matrix([1],[%i])); 0\$ rank(matrix([1,2],[0,0])); 1\$ nullity(matrix([1,2],[0,0])); 1\$ rank(matrix([1,2],[0,1])); 2\$ nullity(matrix([1,2],[0,1])); 0\$ /*--- dotproducts ------*/ dotproduct(matrix([1]),matrix([1])); 1\$ dotproduct(matrix([%i]),matrix([%i])); 1\$ dotproduct(matrix([1],[2]),matrix([-2],[1])); 0\$ dotproduct(matrix([1],[%i]),matrix([1],[%i])); 2\$ /*--- matrix nullity ----*/ nullity(matrix([0])); 1\$ nullity(matrix([1])); 0\$ nullity(matrix([0,0])); 2\$ nullity(matrix([0,1])); 1\$ nullity(matrix([1,0])); 1\$ nullity(matrix([1,2],[0,0])); 1\$ nullity(matrix([1,2],[0,1])); 0\$ /*---orthogonal complements ---*/ (m : matrix([1,2,3]),0); 0\$ (s1 : args(columnspace(m)),0); 0\$ (s2 : args(nullspace(ctranspose(m))),0); 0\$ every(lambda([e], e=0), map(lambda([e], apply('dotproduct,e)), cartesian_product(setify(s1),setify(s2)))); true\$ /*----*/ (m : matrix([0,0,3]),0); 0\$ (s1 : args(columnspace(m)),0); 0\$ (s2 : args(nullspace(ctranspose(m))),0); 0\$ every(lambda([e], e=0), map(lambda([e], apply('dotproduct,e)), cartesian_product(setify(s1),setify(s2)))); true\$ /*----*/ (m : matrix([0,0,0]),0); 0\$ (s1 : args(columnspace(m)),0); 0\$ (s2 : args(nullspace(ctranspose(m))),0); 0\$ every(lambda([e], e=0), map(lambda([e], apply('dotproduct,e)), cartesian_product(setify(s1),setify(s2)))); true\$ /*----*/ (m : matrix([0,0,3],[0,0,1]),0); 0\$ (s1 : args(columnspace(m)),0); 0\$ (s2 : args(nullspace(ctranspose(m))),0); 0\$ every(lambda([e], e=0), map(lambda([e], apply('dotproduct,e)), cartesian_product(setify(s1),setify(s2)))); true\$ /*----*/ (m : matrix([1,2,3],[%i,2,1]),0); 0\$ (s1 : args(columnspace(m)),0); 0\$ (s2 : args(nullspace(ctranspose(m))),0); 0\$ every(lambda([e], e=0), map(lambda([e], apply('dotproduct,e)), cartesian_product(setify(s1),setify(s2)))); true\$ /*----*/ (m : matrix([1,2,3],[%i,2,1],[1,2,3],[1,2,3]),0); 0\$ (s1 : args(columnspace(m)),0); 0\$ (s2 : args(nullspace(ctranspose(m))),0); 0\$ every(lambda([e], e=0), map(lambda([e], apply('dotproduct,e)), cartesian_product(setify(s1),setify(s2)))); true\$ /*----*/ (m : matrix([1,2,3],[4,5,6],[7,8,9]),0); 0\$ (s1 : args(columnspace(m)),0); 0\$ (s2 : args(nullspace(ctranspose(m))),0); 0\$ every(lambda([e], e=0), map(lambda([e], apply('dotproduct,e)), cartesian_product(setify(s1),setify(s2)))); true\$ /*-----*/ ptriangularize(matrix([1]),z); matrix([1])\$ ptriangularize(matrix([0]),z); matrix([0])\$ ptriangularize(matrix([z]),z); matrix([z])\$ (m : matrix([1,2],[3,4]) - z*ident(2),0); 0\$ abs(ratsimp(determinant(m) / determinant(ptriangularize(m,z)))); 1\$ (m : matrix([1,2,3],[4,5,6], [7,8,9]) - z*ident(3),0); 0\$ abs(ratsimp(determinant(m) / determinant(ptriangularize(m,z)))); 1\$ (m : matrix([1-z,2,3-z],[4,5-z,6], [7-z,8,9+z]),0); 0\$ abs(ratsimp(determinant(m) / determinant(ptriangularize(m,z)))); 1\$ locate_matrix_entry(matrix([1]),1,1,1,1,lambda([e],e # 0), 'bool); [1,1]\$ locate_matrix_entry(matrix([0]),1,1,1,1,lambda([e],e # 0), 'bool); false\$ locate_matrix_entry(matrix([1]),1,1,1,1,identity, 'min); [1,1]\$ locate_matrix_entry(matrix([1]),1,1,1,1,identity, 'max); [1,1]\$ locate_matrix_entry(matrix([1,1]),1,1,1,2,identity, 'min); [1,1]\$ locate_matrix_entry(matrix([1,1]),1,1,1,2,identity, 'max); [1,1]\$ locate_matrix_entry(matrix([1,2]),1,1,1,2,identity, 'max); [1,2]\$ locate_matrix_entry(matrix([1,2]),1,1,1,2,identity, 'min); [1,1]\$ (m : matrix([1,2,3],[4,5,6],[7,8,9]),0); 0\$ locate_matrix_entry(m,1,1,3,3,lambda([e],e+%pi), min); [1,1]\$ locate_matrix_entry(m,1,1,3,3,lambda([e],e+%pi), max); [3,3]\$ locate_matrix_entry(m,1,2,3,3,lambda([e],e+%pi), min); [1,2]\$ locate_matrix_entry(m,1,2,3,3,lambda([e],e+%pi), max); [3,3]\$ columnop(matrix([2,1]),1,2,-2); matrix([4,1])\$ columnop(matrix([2,1]),2,1,x); matrix([2,1-2*x])\$ columnswap(matrix([a,b]),1,2); matrix([b,a])\$ columnswap(matrix([a,b,c]),2,3); matrix([a,c,b])\$ columnswap(columnswap(matrix([a,b],[c,d]),1,2),1,2); matrix([a,b],[c,d])\$ ctranspose(matrix([5+%i*6])); matrix([5-%i*6])\$ ctranspose(matrix([1],[%i])); matrix([1,-%i])\$ rowop(matrix([a],[b]),1,2,x); matrix([a-b*x],[b]); rowswap(rowswap(matrix([a,b],[c,d]),1,2),1,2); matrix([a,b],[c,d])\$ (m : matrix([0,1,1],[0,0,0],[0,0,0]),0); 0\$ nullspace(m); funmake('span, sort([matrix([1],[0],[0]),matrix([0],[1],[-1])]))\$ columnspace(m); span(matrix([1],[0],[0]))\$ zeromatrixp(ctranspose(apply(addcol,args(columnspace(ctranspose(m))))) . apply(addcol,args(nullspace(m)))); true\$ (m : matrix([0,1,2,3],[0,0,1,0],[0,0,0,0]),0); 0\$ zeromatrixp(ctranspose(apply(addcol,args(columnspace(ctranspose(m))))) . apply(addcol,args(nullspace(m)))); true\$ (m : matrix([0,0,2,3],[0,0,1,0],[0,0,0,0]),0); 0\$ zeromatrixp(ctranspose(apply(addcol,args(columnspace(ctranspose(m))))) . apply(addcol,args(nullspace(m)))); true\$ (m : matrix([0,0,%i,3],[0,0,%i,0],[0,0,0,0]),0); 0\$ zeromatrixp(ctranspose(apply(addcol,args(columnspace(ctranspose(m))))) . apply(addcol,args(nullspace(m)))); true\$ (remvalue(s1,s2,m),0); 0\$ hilbert_matrix(2); matrix([1,1/2],[1/2,1/3])\$ kronnecker_product(matrix([a]),matrix([b])); matrix([a*b]); kronnecker_product(matrix([1,2],[3,4]),matrix([a,b],[c,d])); mat_unblocker(matrix([matrix([a,b],[c,d]), 2*matrix([a,b],[c,d])],[3*matrix([a,b],[c,d]),4*matrix([a,b],[c,d])]))\$ listp([-1,-2,3],lambda([e], is(e < 0))); false\$ listp([-1,-2,3],lambda([e], ratnump(e))); true\$ matrixp(matrix([1.3b0,7.8b0],[5.7b0, 9.1b0]), bfloatp); true\$ matrixp(matrix([1.3b0,7.8b0],[5.7b0, 7/8]), bfloatp); false\$ polytocompanion(a*x^3+b*x^2+c*x+d,x); matrix([0,0,-d/a],[1,0,-c/a],[0,1,-b/a])\$ vandermonde_matrix([a]); matrix([1])\$ vandermonde_matrix([a,b]); matrix([1,a],[1,b])\$ --- NEW FILE: test-lu.mac --- (kill(values),0); 0\$ (lu_test(m) := block([mp], mp : get_lu_factors(lu_factor(m,generalfield)), zeromatrixp(m - mp[1] . mp[2] . mp[3])),0); 0\$ lu_test(matrix([1])); true\$ lu_test(matrix([0,1],[1,0])); true\$ lu_test(matrix([0,a],[b,0])); true\$ lu_test(matrix([a,b],[c,d])); true\$ lu_test(matrix([a,b],[b,d])); true\$ lu_test(matrix([1/2,b],[b,2/3])); true\$ lu_test(matrix([1/2 + %i, %i - 1],[9 - %i, 19 + %i])); true\$ lu_test(matrix([0,1,0],[0,0,1],[1,2,4])); true\$ lu_test(matrix([0,a,0],[b,0,0],[0,0,c])); true\$ lu_test(matrix([a,b,c],[d,e,f],[g,h,i])); true\$ lu_test(matrix([0,0,1],[1,0,0],[0,1,0])); true\$ lu_test(matrix([0,0,1],[1,0,0],[0,1,x])); true\$ lu_test(hilbert_matrix(3)); true\$ lu_test(hilbert_matrix(4)); true\$ lu_test(hilbert_matrix(5)); true\$ lu_test(vandermonde_matrix([1,2,%i])); true\$ lu_test(vandermonde_matrix([1 + %i,1 - %i, 8, 32])); true\$ (lu_test(m, tol) := block([mp], mp : get_lu_factors(lu_factor(m,bigfloatfield)), matrixp(m - mp[1] . mp[2] . mp[3], lambda([s], is(cabs(s) < tol)))),0); 0\$ (fpprec : 20,0); 0\$ lu_test(vandermonde_matrix([1,1/2,1/3,1/4]), 1.0e-18); true\$ (fpprec : 40,0); 0\$ lu_test(vandermonde_matrix([1,1/2,1/3,1/4]), 1.0e-38); true\$ lu_test(hilbert_matrix(13), 1.0e-38); true\$ (lu_test(m, tol) := block([mp], mp : get_lu_factors(lu_factor(m,complexfloatfield)), matrixp(m - mp[1] . mp[2] . mp[3], lambda([s], is(cabs(s) < tol)))),0); 0\$ lu_test(hilbert_matrix(7) + %i * vandermonde_matrix([1,2,3,4,5,6,7]), 1.0e-10); true\$ (m : matrix([1,2,3],[4,5,6],[7,8,10]),0); 0\$ (b : m . matrix([x],[y],[z]),0); 0\$ lu_backsub(lu_factor(m,crefield), b, crefield); matrix([x],[y],[z])\$ (m : matrix([0,2,3],[0,5,6],[7,8,10]),0); 0\$ (b : m . matrix([x],[y],[z]),0); 0\$ lu_backsub(lu_factor(m, crefield), b, crefield); matrix([x],[y],[z])\$ (m : matrix([0,2,3],[0,5,6],[7,8,10]),0); 0\$ (b : m . matrix([x],[y],[z]),0); 0\$ lu_backsub(lu_factor(m, crefield), b, crefield); matrix([x],[y],[z])\$ (m : matrix([0,1,2],[1,0,1],[0,0,8]),0); 0\$ (b : m . matrix([x],[y],[z]),0); 0\$ lu_backsub(lu_factor(m, crefield), b, crefield); matrix([x],[y],[z])\$ --- NEW FILE: test-matrixexp.mac --- (load("matrixexp.lisp"),0); 0\$ (remvalue(a,b,c,x),0); 0\$ (xequal(a,b) := block([gcd : 'spmod, algebraic : true], ratvars(), is(fullratsimp(a) = fullratsimp(b)) or is(fullratsimp(rectform(a)) = fullratsimp(rectform(b)))),0); 0\$ scalarmatrixp : false; false\$ (check(mat,n) := block([s, okay, z : ?gensym(), i : 0, m], s : spectral_rep(mat), okay : xequal(mat - s[1] . s[2] - s[3] , zeromatrix(n,n)), print(i : i + 1, okay), okay : okay and xequal(mat . s[3], s[3] . mat), print(i : i + 1, okay), okay : okay and xequal(s[3]^^n, zeromatrix(n,n)), print(i : i + 1, okay), okay : okay and xequal(matrixfun(lambda([x],1), mat), ident(n)), print(i : i + 1, okay), okay : okay and xequal(matrixfun(lambda([x], 1+x+x^2), mat), ident(n) + mat + mat . mat), print(i : i + 1, okay), okay : okay and xequal(matrixfun(lambda([x], x^2/5), mat), mat . mat / 5), print(i : i + 1, okay), okay : okay and xequal(matrixfun(lambda([x],x^8),mat), mat^^8), print(i : i + 1, okay), m : errcatch(matrixfun(lambda([z],sqrt(z)), mat)), okay : okay and if m = [ ] then member(0, s[1]) else xequal(m[1] . m[1], mat), print(i : i + 1, okay), m : errcatch(matrixfun(lambda([mat], 1/mat^2), mat)), okay : okay and if m = [ ] then member(0, s[1]) else xequal(m[1], mat^^-2), print(i : i + 1, okay), m : fullratsimp(matrixexp(mat,z) . matrixexp(mat, -z)), okay : okay and xequal(m, ident(n)), print(i : i + 1, okay), m : matrixexp(mat,z), okay : okay and xequal(diff(m,z) - m . mat, zeromatrix(n,n)), print(i : i + 1, okay), okay), 0); 0\$ check(matrix([1]),1); true\$ check(matrix([x]),1); true\$ check(matrix([1,0],[0,1]),2); true\$ check(matrix([1,0],[0,0]),2); true\$ check(matrix([0,0],[0,1]),2); true\$ check(matrix([0,1],[1,0]),2); true\$ check(matrix([0,1],[-1,0]),2); true\$ check(matrix([0,%i],[-%i,0]),2); true\$ check(matrix([0,a],[0,b]),2); true\$ check(matrix([1,2],[2,1]),2); true\$ check(matrix([1,%i],[-%i,1]),2); true\$ check(matrix([1,2],[3,1]),2); true\$ check(matrix([1,6],[0,1]),2); true\$ check(matrix([x,1],[1,1]), 2); true\$ check(matrix([sqrt(3), sqrt(2)],[sqrt(2), 9]),2); true\$ check(matrix([sqrt(6),1],[1, -9]),2); true\$ check(matrix([sqrt(3), sqrt(2)],[sqrt(2), 9]),2); true\$ check(matrix([a,b],[c,d]),2); true\$ check(matrix([1,0,0],[0,1,0],[0,0,1]),3); true\$ check(matrix([1,0,0],[0,1,0],[0,0,-1]),3); true\$ check(matrix([1,1,0],[1,1,1],[0,1,1]),3); true\$ check(matrix([1,%i,0],[-%i,1,%i],[0,-%i,1]),3); true\$ check(matrix([1,2,3],[4,5,6],[7,8,9]),3); true\$ check(matrix([1,a,b],[0,1,c],[0,0,2]),3); true\$ check(matrix([1,a,b],[0,1,c],[0,0,1]),3); true\$ check(matrix([1,a,b],[0,2,c],[0,0,1]),3); true\$ check(matrix([2,a,b],[0,2,c],[0,0,1]),3); true\$ check(matrix([1,a,b],[0,1,c],[0,0,2]),3); true\$ check(matrix([1,0,0],[a,1,0],[b,c,1]),3); true\$ check(matrix([1,0,0],[a,%i,0],[b,c,%i]),3); true\$ check(matrix([1,0,0],[a,%i,0],[b,c,1]),3); true\$ check(matrix([1,0,0],[a,%i,0],[b,c,sqrt(7)]),3); true\$ check(matrix([1,0,0],[a,%i,0],[b,c,%i]),3); true\$ /* spectral_rep checks its results--so this really does testing.*/ (spectral_rep(matrix([x,1,0],[1,1,1],[0,1,1])),0); 0\$ check(matrix([1,1,0],[1,1,1],[0,1,1]),3); true\$ (s : matrix([1,2],[2,3]), 0); 0\$ (m : s . matrix([1,0],[0,-1]) . s^^-1, 0); 0\$ xequal(s . matrix([exp(1),0],[0,exp(-1)]) . s^^-1, matrixexp(m)); true\$ /* Matrices from Gosei Furuya's Jordan form code. */ (m : matrix([2,0,0,0,0,0,0,0], [1,2,0,0,0,0,0,0], [-4,1,2,0,0,0,0,0], [2,0,0,2,0,0,0,0], [-7,2,0,0,2,0,0,0], [9,0,-2,0,1,2,0,0], [-34,7,1,-2,-1,1,2,0], [145,-17,-16,3,9,-2,0,3]), 0)\$ 0\$ xequal(matrixfun(lambda([x],x^2), m), m.m); true\$ (m : matrix([2,1,0,0,0,0], [-1,4,0,0,0,0], [-1,1,2,1,0,0], [-1,1,-1,4,0,0], [-1,1,-1,1,3,0], [-1,1,-1,1,1,2]),0)\$ 0\$ xequal(matrixfun(lambda([x],x^2), m), m.m); true\$ (m : matrix([0,1,0], [0,0,1], [1,-3,3]),0); 0\$ xequal(matrixfun(lambda([x],1+x^5), m), ident(3)+m^^5); true\$ (m : matrix([0,1,2], [0,0,0], [0,0,0]),0); 0\$ xequal(matrixfun(lambda([x],1+x^5), m), ident(3)+m^^5); true\$ (m : matrix([3,1,0,0], [-4,-1,0,0], [7,1,2,1], [-17,-6,-1,0]),0); 0\$ xequal(matrixfun(lambda([x],1+x^2), m), ident(4)+m^^2); true\$ (m : matrix([2,1,2,0], [-2,2,1,2], [-2,-1,-1,1], [3,1,2,-1]),0); 0\$ xequal(matrixfun(lambda([x],1+x^2), m), ident(4)+m^^2); true\$ (m : matrix([0,0,1,1,1], [0,0,0,1,1], [0,0,0,0,1], [0,0,0,0,0], [0,0,0,0,0]),0); 0\$ xequal(matrixfun(lambda([x],1+x^2), m), ident(5)+m^^2); true\$ (m : matrix([0,1,0], [0,0,1], [-1,-3,-3]),0); 0\$ xequal(matrixfun(lambda([x],1+x^5), m), ident(3)+m^^5); true\$ --- NEW FILE: test-polynomialp.mac --- (load("polynomialp.lisp"),0); 0\$ polynomialp(1,[x])\$ true\$ polynomialp(%pi,[x])\$ true\$ polynomialp(sqrt(23),[x])\$ true\$ polynomialp(1+x,[x])\$ true\$ polynomialp(1 + x * (sqrt(5) + x),[x])\$ true\$ polynomialp(1 + x * (sqrt(5) + x),[y])\$ false\$ polynomialp(1 + x * (sqrt(5) + y),[x,y])\$ true\$ polynomialp(1 + sqrt(x),[x], 'numberp, 'numberp); true\$ polynomialp(1 + sqrt(1 + sqrt(x)),[x], 'numberp, 'numberp); true\$ polynomialp(1 + sqrt(x + sqrt(1 + x*y)),[x,y], 'numberp, 'numberp); true\$ polynomialp(cos(x),[x]); false\$ polynomialp(cos(x),[x], 'numberp, 'numberp); false\$ polynomialp([x],[x], 'numberp, 'numberp); false\$ polynomialp((1+x)^a,[x], 'constantp, lambda([e],freeof(x,e))); true\$ polynomialp((1+x)^a,[x], 'constantp); false\$ ```