function dis = distillable_coherence(rho,tol)
% distillable_coherence(rho,tol) computes the one-shot distillable coherence under MIO/DIO
% rho: dxd density matrix or normalised state vector
% tol: error tolerance 0<=tol<=1

dim = size(rho);
d = max(dim);
if dim(1)==1; rho=rho'*rho; 
elseif dim(2)==1; rho=rho*rho'; end
rho = (rho+rho')/2; % to avoid numerical problems

cvx_begin sdp quiet

    variable G(d,d) hermitian
    variable k nonnegative
    0 <= G <= eye(d)
    trace(rho*G) >= 1-tol;
    diag(G) == k*ones(d,1)

    minimise k

cvx_end

dis = log2(floor(1/cvx_optval));

end