function d = simplex_proj(c,tau)

n = max(size(c));
p = -c;
pmax = max(p);
sm = sum(p);
if sm >= n*pmax - tau
  lambda = (tau+sm)/n;
  d = max(0, c+lambda*ones(n,1));
  clear p;
  return;
end 

p = sort(p);
phi = zeros(n,1);

sm = 0;
for i = 1:n
  sm = sm + p(i);
  phi(i) = i*p(i) - sm;
end

lid = 1; rid = n;
while (rid - lid > 1)
  mid = ceil((lid+rid)/2);
  if phi(mid) < tau
    lid = mid;
  else
    rid = mid;
  end
end
k = lid;
delta = (tau - phi(k))/k;
lambda = p(k) + delta;
d = max(0, c+lambda*ones(n,1));
clear p phi;
