funm_randomized

function F = funm_randomized(A,fun)
%FUNM_RANDOMIZED  Evaluate general matrix function using
% randomized approximate diagonalization method of Davies (2007).
tol = 8*eps(A(1,1));  % Tolerance for single or double precision A.
E = randn(size(A)); 
[V,D] = eig(A + (tol*norm(A,'fro')/norm(E,'fro'))*E);
F = V*diag(fun(diag(D)))/V;

Thanks, Pavel - that's an excellent point. E should indeed be of the same precision as A and your change achieves that very neatly.