function y=ln(x)
% 10^-4 <= x <= 10^4
x1=x(:);
assert(all(1.0e-04 <= x1));
assert(all(x1 <= 1.0E+04));
% find k such that e^k <= x < e^k+1
% note: exp(-10) ~ 4.54E-05
%       exp(10)  ~ 2.20E+04
powers  = (-10:10)';
edges   = exp(powers);
k       = discretize(x1,edges)
x2      = x1./edges(k);
% now x2 is in the fundamental domain [1,e]
% build the approximating polynomial based on N Chebyshev nodes
% from Exersise 3.3.7 we know the worst-case error formula:
% error <= (e-1)^N/N/2^(2N-1)
% Thus N=25 guarantees error < 10^(-10)
N  = 25;
Xb = chebyshev(exp(0),exp(1),N);
Yb = log(Xb);
c  = newtdd(Xb,Yb,N);
y  = arrayfun(@(z)nest(N-1,c,z,Xb),x2)+powers(k);
y  = reshape(y,size(x));


function bins = discretize(x, edges)
%  discretize  Group numeric data into bins or categories.
%     BINS = discretize(X,EDGES) returns the indices of the bins that the 
%     elements of X fall into.  EDGES is a numeric vector that contains bin 
%     edges in monotonically increasing order. An element X(i) falls into 
%     the j-th bin if EDGES(j) <= X(i) < EDGES(j+1), for 1 <= j < N where 
%     N is the number of bins and length(EDGES) = N+1. The last bin includes 
%     the right edge such that it contains EDGES(N) <= X(i) <= EDGES(N+1). For 
%     out-of-range values where X(i) < EDGES(1) or X(i) > EDGES(N+1) or 
%     isnan(X(i)), BINS(i) returns NaN.     
% 
%
% I have this function in my Matlab installation
% here is a "poor man's" implementation
%
% It would be much more efficient to do a binary search in a
% sorted array of edges instead of using find
%
bins = NaN*ones(size(x));
for i=1:length(x(:)),
    J = find(edges <= x(i));
    if isempty(J)
        continue;
    elseif J(end)==length(edges) && edges(J(end))<x(i),
        continue;
    else
        bins(i)=J(end);
    end
end