function logf = logmultigaussian(x,mu,S)
  % Returns the logarithm of a multivariate Gaussian density
  % mu must be a 1xd row vector, S a dxd square covariance matrix
  % x a nxd matrix, with one point per row
  % returns a nx1 column vector of log-densities, calculated at these points
  [n,d] = size(x);
  xmu = op(x,'-',mu);
  logf = (-d/2)*log(2*pi)-sum((xmu*inv(S)).*xmu,2)/2-log(det(S))/2;