Here is a function that beats those above by a mile:
function S = DiscSampVec2(x,p,ns);
[~,idx] = histc(rand(1,ns),[0,cumsum(p)]);
S = x(idx);
This is a slight modification of the function on page 47, Kroese, Taimre, and Botev, Handbook of Monte Carlo Methods, Wiley, 2011.
>> ns=10^6; n=10^3; x=1:n;
>> p = rand(1,n); p = p/sum(p);
>> S = DiscSampVec2(x,p,ns);
@Derek, I am aware and glad that there are other formulations of this task. 'N' and 'M' are clearly defined in the description above. There is no overt relationship between P and N & M. P is a distribution, in the mathematical sense of the word, and because it describes the probability of picking an index of P, the sum of those values must equal one, i.e. it must be normalized.