#5708 matrixfun gives wrong results for matrices with %i (final simplification runs at domain: real)

[Ilia Gradina]

Following up on #2596, as suggested there, one ticket per remaining bug. Commit 48a0fb fixes matrixexp, but matrixfun on a matrix containing %i is still wrong on master (6a3bd21, SBCL):

(%i1) load(linearalgebra)$
(%i2) m : %i*matrix([1,2,3],[3,2,0],[1,2,1])$
(%i3) d : float(rectform(matrixfun(lambda([x], exp(x)), m) - matrixexp(m)))$
(%i4) lmax(map(lambda([r], lmax(map(cabs, r))), args(d)));
(%o4)                          31.1083788527916

matrixexp(m) itself is right: it agrees with the power series to 6e-40 at fpprec 40. matrixfun(sin) and matrixfun(cos) on the same m are off by 33.8 and 32.4 against sin(A) = (exp(%iA) - exp(-%iA))/(2%i) and cos(A) = (exp(%iA) + exp(-%iA))/2 built from matrixexp. For the real matrix, %imatrixfun(sinh) agrees with the same formula to 2e-37, so the spectral data is fine.

The cause: $matrixfun binds domain: complex in its inner let, but calls (final-simp ff) after that let has closed, so the final ratsimp/factor runs at domain: real. matrixexp doesn't have this problem because matrixexp-simp binds domain: complex itself.

Binding domain: complex around final-simp unconditionally fixes it, but the symbolic result of matrixfun(lambda([x], x+1/x), matrix([x,1,0],[1,1,1],[0,1,1])) grows from 41 to 36,236 characters. So the patch below binds it only when the matrix contains %i, and adds two tests to rtest_matrixexp.mac:

diff --git a/share/linearalgebra/matrixexp.lisp b/share/linearalgebra/matrixexp.lisp
index f8e4fdd..dd676e9 100644
--- a/share/linearalgebra/matrixexp.lisp
+++ b/share/linearalgebra/matrixexp.lisp
@@ -178,7 +178,7 @@
       (merror "The ~:M argument to `~:M' must be a lambda form with ~:M variable(s)" pos fun-name n))))

 (defmfun $matrixfun (lamexpr mat)

-  (let ((ff))
+  (let ((ff) (cplx (not ($freeof '$%i mat))))
   (let (($gcd '$spmod) ($ratmx nil) (z (gensym)) (expr) (var) (sp) (d) (p) (di) 
          (n ($length ($args mat))) (f 0)  ($matrix_element_mult "*") ($matrix_element_add "+")
         ($domain '$complex))
@@ -201,7 +201,9 @@
       (setq di (ncmul2 di d))
       (setq expr (div ($diff expr z) (factorial (+ i 1)))))
     (setq ff f))
-    (final-simp ff)))
+    (if cplx
+        (let (($domain '$complex)) (final-simp ff))
+        (final-simp ff))))

 ;; Return the residue of the rational expression e with respect to the
 ;; variable var at the point pt.  Assumptions:
diff --git a/share/linearalgebra/rtest_matrixexp.mac b/share/linearalgebra/rtest_matrixexp.mac
index 1514834..99f9dba 100644
--- a/share/linearalgebra/rtest_matrixexp.mac
+++ b/share/linearalgebra/rtest_matrixexp.mac
@@ -329,3 +329,15 @@ matrix([0,0,0],[0,0,0],[0,0,0])$
 (kill(s,m,mat,xxx,zzz,rrr,xequal,check), 0);
 0$

+
+/* matrixfun on a matrix with %i: the final simplification must run with domain : complex */
+block([m : %i*matrix([1,2,3],[3,2,0],[1,2,1]), d],

+  d : float(rectform(matrixfun(lambda([x], exp(x)), m) - matrixexp(m))),
+  is(lmax(map(lambda([r], lmax(map(cabs, r))), args(d))) < 1.0e-10));
+true$
+
+block([m : %i*matrix([1,2,3],[3,2,0],[1,2,1]), d],
+  d : float(rectform(matrixfun(lambda([x], sin(x)), m)
+                     - (matrixexp(%i*m) - matrixexp(-%i*m))/(2*%i))),
+  is(lmax(map(lambda([r], lmax(map(cabs, r))), args(d))) < 1.0e-10));
+true$

With it, matrixfun(exp) agrees with matrixexp exactly, sin and cos agree with the exp formula to 2e-37, and the symbolic example stays at 41 characters. rtest_matrixexp gives 83/83 with the patch; without the fix the two new tests fail (81/83).