function [f,g,H]=logist_X_fgH(w, X,Y,alpha)
  % w is k(C-1) by 1
  % X is n by k
  % Y is n by 1, taking values from 1..C
  % alpha is k(C-1) by 1
  
  % constants
  eps_inc=1E-5;
  
  % get inputs
  C=max(Y);
  [n,k]=size(X);
  kC=length(w);
  k2=[0:C-1]*k;
  y=double(repmat(Y(:),1,C)==repmat([1:C],n,1));

  % compute u and mu
  u=zeros(n,C-1);
  mu=zeros(n,C);
  for i=1:C-1
    u(:,i)=X*w(k2(i)+1:k2(i+1));
  end
  expu=exp(u);
  mu=zeros(n,C);
  mu(:,1)=1./(1+sum(expu,2));
  mu(:,2:C)=expu./repmat(1+sum(expu,2),1,C-1);
  assert(norm(sum(mu,2)-1)<1E-3);
  
  % get f g and H
  if nargout>=1
    f=-sum(sum(log(mu+eps_inc).*y))+1/2*sum(w.^2.*alpha);
  end
  if nargout>=2
    g=zeros(kC,1);
    for i=1:C-1
      g(k2(i)+1:k2(i+1))=X'*(mu(:,i+1)-y(:,i+1));
    end
    g=g+w.*alpha;
  end
  if nargout>=3
    H=zeros(kC,kC);
    for i=1:C-1
      for j=1:C-1
        delta=double(i==j);
        t=mu(:,i+1).*(delta-mu(:,j+1));
        Hblock=X'*(repmat(t,1,k).*X);
        H(k2(i)+1:k2(i+1),k2(j)+1:k2(j+1))=Hblock;
      end
    end
    H=H+diag(alpha);
  end
