function [T,H,p,q] = THsplit(A)
% returns matrices T, H and numbers p, q, such that 
%
%   A = T + H + p*P + q*Q
% 
%  where 
%  P = matrix of all ones = ones(N) 
%  Q = ones with alternating signs
% 
%  The input matrix A is assumed to be a square matrix A that is the sum of
%  a Toeplitz matrix and a Hankel matrix.
% 

pass_test = TplusHtest(A);
if ~pass_test
    error('Input matrix must be Toeplitz plus Hankel.')
end

N = size(A,1); 

%% get the T + H splitting
% There are two matrices in the intersection 
% of Toeplitz matrices and Hankel matrices: P and Q.
% So we find the projections of our matrix onto the intersection.
% Get projection p onto all 1s: 
P = ones(N);
p = sum(sum(A.*P))/N^2; % of course this is not efficient 
% Get projection q onto alternating +/-1s :
Q = (-1).^(1:N); Q = Q.'*Q;
q = sum(sum(A.*Q))/N^2;

% The remainder of the matrix is strictly Toeplitz plus strictly Hankel
TH = A - p*P - q*Q;

%% Find the Hankel part
%  Differences along diagonals are due to the Hankel part 
d = diag(TH,0);  % main diagonal
d = diff(d); 
d = [0
     d];
d = cumsum(d);
e = diag(TH,-1); % sub diagonal
e = diff(e);
e = [0 
     e];
e = cumsum(e);

a = zeros(2*N,1);
a(1:2:2*N) = d;
a(2:2:2*N-1) = e;
r = a(1:N); 
c = a(N+1:end-1); 
c = [r(N)
     c];
r = flipud(r);
H = toeplitz(c,r);
% Turn Toeplitz into Hankel by flipping:
H = fliplr(H);

% The Toeplitz part is what is left after we take out the Hankel part:
T = TH - H;