function [N,PHI0] = phi1m(A)
%
% function N = phi1m(A)
%

% Scale A by power of 2 so that its norm is < 1/2 .
[f,e] = log2(norm(A,1));
s = min(max(0,e+1),1023);
A = A/2^s;
% Pade approximation for phik(A)
ID = eye(size(A));
n = [1,1/30,1/30,1/936,1/4680,1/171600,1/3603600,1/259459200];
d = [1,-7/15,1/10,-1/78,1/936,-1/17160,1/514800,-1/32432400];
q = length(n);
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 = full(D\N);
% Undo scaling by repeated squaring
PHI0 = A*N+ID;
for i = 1:s
  N = (PHI0+ID)*N/2;
  PHI0 = PHI0*PHI0;
end
