function [F,n_terms] = fun_atom(T,fun,tol,prnt)
%FUNM_ATOM  Function of triangular matrix with nearly constant diagonal.
%           [F, N_TERMS] = FUNM_ATOM(T, FUN, TOL, PRNT) evaluates function
%           FUN at the upper triangular matrix T, where T has nearly constant
%           diagonal.  A Taylor series is used.
%           FUN(X,K) must return the K'th derivative of
%           the function represented by FUN evaluated at the vector X.
%           TOL is a convergence tolerance for the Taylor series,
%           defaulting to EPS.
%           If PRNT ~= 0 trace information is printed.
%           N_TERMS is the number of terms taken in the Taylor series.
%           N_TERMS  = -1 signals lack of convergence.

if nargin < 3 | isempty(tol), tol = eps; end
if nargin < 4, prnt = 0; end

if isequal(fun,@fun_log)   % LOG is special case.
   [F,n_terms]  = logm_isst(T,prnt);
   return
end

itmax = 500;

n = length(T);
if n == 1, F = feval(fun,T,0); n_terms = 1; return, end

lambda = sum(diag(T))/n;
F = eye(n)*feval(fun,lambda,0);
f_deriv_max = zeros(itmax+n-1,1);
N = T - lambda*eye(n);
mu = norm( (eye(n)-abs(triu(T,1)))\ones(n,1),inf );

P = N;
max_d = 1;

for k = 1:itmax

    f = feval(fun,lambda,k);
    F_old = F;
    F = F + P*f;
    rel_diff = norm(F - F_old,inf)/(tol+norm(F_old,inf));
    if prnt
        fprintf('%3.0f: coef = %5.0e', k, abs(f)/factorial(k));
        fprintf('  N^k/k! = %7.1e', norm(P,inf));
        fprintf('  rel_d = %5.0e',rel_diff);
        fprintf('  abs_d = %5.0e',norm(F - F_old,inf));
    end
    P = P*N/(k+1);

    if rel_diff <= tol

      % Approximate the maximum of derivatives in convex set containing
      % eigenvalues by maximum of derivatives at eigenvalues.
      for j = max_d:k+n-1
          f_deriv_max(j) = norm(feval(fun,diag(T),j),inf);
      end
      max_d = k+n;
      omega = 0;
      for j = 0:n-1
          omega = max(omega,f_deriv_max(k+j)/factorial(j));
      end

      trunc = norm(P,inf)*mu*omega;  % norm(F) moved to RHS to avoid / 0.
      if prnt
          fprintf('  [trunc,test] = [%5.0e %5.0e]', ...
                   trunc, tol*norm(F,inf))
      end
      if prnt == 5, trunc = 0; end % Force simple stopping test.
      if trunc <= tol*norm(F,inf)
         n_terms = k;
         if prnt, fprintf('\n'), end, return
      end
    end

    if prnt, fprintf('\n'), end

end
n_terms = -1;