/*  Make an n element Maxima list [a,a,...a].
*/

duplicate(a,n) := block([k],
  makelist(a,k,1,n));

/* Compute f(m), where m is a diagonalizable matrix and
   f : C -> C.
*/

matrix_function(f, m) := block([z, s, v, k],
  load("eigen"),
  z : UNITEIGENVECTORS(m),
  if (NONDIAGONALIZABLE) then (
      false
  ) else (
      s : z[1],              /* s = the spectrum of m */
      s : apply(append, fullmapl(duplicate, s[1], s[2])),
      v : rest(z),          /* v = normalized eigenvectors of m */
      sum(f(s[k]) * transpose(v[k]) . CONJUGATE(v[k]),k,1,length(v))));

/* Compute exp(m), where m is a diagonalizable matrix.
*/

matrix_exp(m) := matrix_function(exp,m);