function N = phi1m(A)
%
% function N = phi1m(A)
%

% Scale A by power of 2 so that its norm is < 1/2 .
[f,e] = log2(norm(A,'inf'));
s = min(max(0,e),1023);
A = A/2^s;
% Pade approximation for phik(A)
ID = eye(size(A));
X = A;
q = 6;
n = [1,1/26,5/156,1/858,1/5720,1/205920,1/8648640];
d = [1,-6/13,5/52,-5/429,1/1144,-1/25740,1/1235520];
N = n(1)*ID;
D = d(1)*ID;
X = ID;
for i = 2:q
   X = A*X;
   N = N+n(i)*X;
   D = D+d(i)*X;
end
N = D\N;
% Undo scaling by repeated squaring
PHI0 = A*N+ID;
for i = 1:s
  N = (PHI0+ID)*N/2;
  PHI0 = PHI0*PHI0;
end
