function [Y,Z]=schurlog(A,p,q,v)
% Y=SCHURLOG(A,p) computes A^(1/p) using a binary powering technique.
% Y=SCHURLOG(A,p,q) computes A^(q/p)=(A^(1/p))^q using a 
%   binary powering technique. 
% Y=SCHURLOG(A,p,q,v) computes A^(q/p) where the determinations of
%   the p-th roots are chosen according to the vector v(1),...,v(s)
%   where s is the number of diagonal blocks in the real Schur form
%   of A and v(i) is the determination of the p-th root of one  
%   eigenvalue on the i-th block obtained by schur(A), in the complex
%   conjugate case two conjugate determinations are chosen in order
%   to get real roots and fractional powers
%   v(i)=0 for the principal determination
%   v(i)=1,...,p-1, for the other determinations.
%
% Reference.
% B. Iannazzo and C. Manasse, "A Schur logarithmic algorithm for
% fractional powers of a matrix", SIAM J. Matrix Anal. Appl, 34-2,
% pp. 794-813, 2013.
%
% Warning! This is just a straightforwad implementation. 
% The code is not optimized.
% works for real matrices

n=size(A,1);
tol=1e-15;

if nargin<3
  q=1;
end

% compute the Schur form of A and the block structure of T
[Q,T]=schur(A);
[m,fb,fe,s]=block_structure(T);

if nargin<4
  v=zeros(m,1);
else
  if length(v)~=m
    disp('Error: length of v different from the number of blocks of schur(A)');
  end
end

% compute the binary expansion of p
[b,c]=binary(p);
tt=length(b); % number of binary digits of p or floor(log2(p))+1
ts=length(c); % number of nonzero binary digits of p

Y=zeros(n);

V=zeros(n,n,tt+1);
W=zeros(n,n,ts+1);
U=zeros(n,n,ts+1);

V(:,:,1)=eye(n);
W(:,:,1)=eye(n);
W(:,:,ts+1)=T;
U(:,:,1)=T;
U(:,:,ts+1)=eye(n); 

% work a block column at a time, j is the column index
for j=1:m

	indj=fb(j):fe(j);
	% first compute the diagonal block (j,j) of V(:,:,k) and W(:,:,k)
	V(indj,indj,2)=rootblock(T(indj,indj),p,1,v(j));
	for k=3:tt+1
		V(indj,indj,k)=V(indj,indj,k-1)*V(indj,indj,k-1);
	end
	W(indj,indj,2)=V(indj,indj,tt+1);
	for k=3:ts
		W(indj,indj,k)=W(indj,indj,k-1)*V(indj,indj,c(k-1)+2);
	end
	% Just the diagonal of U is defined
	for k=ts:-1:1
		U(indj,indj,k)=U(indj,indj,k+1)*V(indj,indj,c(k)+2);
	end
	
	% work on the block (i,j) from the diagonal to the top of column j
	for i=j-1:-1:1
		indi=fb(i):fe(i);
		% compute B^{(h)}
		clear B,
		for h=2:tt		
			B(:,:,h)...
			=V(indi,fb(i+1):fe(j-1),h)*V(fb(i+1):fe(j-1),indj,h);
		end		

		% compute C^{(h)}
		clear C
		C(:,:,1)=zeros(s(i),s(j));
		for h=2:ts
			C(:,:,h)...
			=W(indi,fb(i+1):fe(j-1),h)*V(fb(i+1):fe(j-1),indj,c(h)+2);
		end
				
		% compute the right hand side (rhs) of the system to be solved 

		% rhs = T_{ij}
		rhs=reshape(T(indi,indj),s(i)*s(j),1);
		rhs1=rhs;
		
		% rhs = rhs - terms containing B
		clear sigma
		Ii=eye(s(i));Ij=eye(s(j));
		for h=1:tt-1
			sigma(:,:,h)...
			=kron(V(indj,indj,h+1).',Ii)+kron(Ij,V(indi,indi,h+1));
		end
		clear emme
		tmp=1;
		emme(:,:,tt)=kron(U(indj,indj,2).',Ii);
		for h=tt-1:-1:1
			emme(:,:,h)=emme(:,:,h+1)*sigma(:,:,h);
			if (b(tt-h+1)==1)
				tmp=tmp+1;
				emme(:,:,h)=emme(:,:,h)...
				+kron(U(indj,indj,tmp+1).',W(indi,indi,tmp));
			end
		end
		
		for h=2:tt			
			rhs=rhs-emme(:,:,h)*reshape(B(:,:,h),s(i)*s(j),1);
		end 
		
		% rhs = rhs - terms containing C
		tmp=eye(s(j));
		tmp2=C(:,:,ts);
		for h=ts-1:-1:2			
			tmp2=tmp2+C(:,:,h)*U(indj,indj,h+1);
		end
		rhs=rhs-reshape(tmp2,s(i)*s(j),1);
		
		% solve the linear system and unvec the solution
		V(indi,indj,2)=reshape(emme(:,:,1)\rhs,s(i),s(j));
		
		% now compute the blocks of V(:,:,k) and W(:,:,k)
		for k=3:tt+1
			V(indi,indj,k)=V(indi,indi,k-1)*V(indi,indj,k-1)...
			+V(indi,indj,k-1)*V(indj,indj,k-1)+B(:,:,k-1);
		end
		
		W(indi,indj,2)=V(indi,indj,tt+1);
		for k=3:ts
			W(indi,indj,k)=W(indi,indi,k-1)*V(indi,indj,c(k-1)+2)...
			+W(indi,indj,k-1)*V(indj,indj,c(k-1)+2)+C(:,:,k-1);
		end
		
	end
end

% Z is the solution of the quasi triangular equation
Z=V(:,:,2);
 
% if q is defined then compute Z^q
% Warning: the implementation of this part has computational cost 
% O(log^2 p), this choice is for simplicity, a O(log p) implementation
% is possible
if (nargin==3 || q>1)
	[bq,cq]=binary(q);
	ts=length(cq);
	Z=eye(n);
	tmp=0;
	for h=1:ts
		tmp=tmp+2^cq(h);
		for j=1:m
			indj=fb(j):fe(j);
			ZN(indj,indj)=rootblock(T(indj,indj),p,tmp);
			for i=j-1:-1:1
				indi=fb(i):fe(i);
				ZN(indi,indj)=zeros(s(i),s(j));
				ZN(indi,indj)=ZN(indi,indj)...
				+Z(indi,fb(i):fe(j))*V(fb(i):fe(j),indj,cq(h)+2);
			end
		end
		Z=ZN;
	end
end

% now get back Y from Z which is quasi triangular
Y=Q*Z*Q';

end		
	
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [m,fb,fe,s]=block_structure(T)
% [m,fb,fe,s]=block_structure(T) computes the block structure of 
%		of the quasi-upper triangular matrix T
%		m is the number of diagonal blocks
%		fb is an array containing the begin of each block
%               fe is an array containing the end of each block
%		s is an array containing the sizes of the diagonal blocks
	
	n=length(T);
	tol=1e-15;

	i=1;j=1;
	
	while(i<n)
		v(j)=i;
		if(abs(T(i+1,i))<tol)
			i=i+1;s(j)=1;
		else
		i=i+2;s(j)=2;
		end
		j=j+1;
	end
	if(i==n)
		v(j)=n;s(j)=1;
	end

	m=length(v);
	for i=1:m
		fb(i)=v(i);
		fe(i)=v(i)+s(i)-1;
	end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [b,c]=binary(p)
% [b,c]=binary(p) computes the binary structure of p
%		b is the array of the binary digits of p
%		c is the array of the position of ones in the binary expansion of p
 
	flag=1;
	k=1;
	h=1;
	while (flag==1)
		if (mod(p,2)==1)
			b(k)=1;
			m(h)=k;
			h=h+1;
		else 
			b(k)=0;
		end
		p=floor(p/2);
		k=k+1;
		if (p==0)
			flag=0;
		end
	end

	s=length(b);
	for k=1:s
		d(s-k+1)=b(k);
	end
	b=d;

	s=length(m);
	for k=1:s
		c(s-k+1)=m(k)-1;
	end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function R=rootblock(R,p,s,k);
% R=rootblock(R,p,s,k) computes R^(p/s) by a direct formula
%   if R is 1x1 or 2x2 and k is the determination of the root
				
  
	if length(R)==1 
		R=R^(1/p)*exp(2*pi*i*k/p);
		R=R^s;
	else 
		theta=1/2*(R(1,1)+R(2,2));
		mu=1/2*sqrt((R(1,1)-R(2,2))^2-4*R(1,2)*R(2,1));
		eigen=theta+i*mu;
	
		eigenroot=(eigen^(1/p)*exp(2*pi*i*k/p))^s;
		alfa=real(eigenroot);
		beta=imag(eigenroot);
	
		R=alfa*eye(2)+(beta/mu)*(R-theta*eye(2));
	end
end