function p = dmp(N,H)
	
	if H > log2(N)
		fprintf('H <= %g with N = %g\n',log2(N),N);
		p = ones(1,N)/N;
		return
	end
	
	% la densita' di massa di probabilita' ritornata e' una
	% permutazione casuale di una distribuzione
	%
	% a^(0:N-1)*(1-a)/(1-a^N)
	%
	% con il parametro 'a' determinato (col metodo delle secanti) in
	% modo da ottenere l'entropia richiesta
	
	x0 = 0.75;
	y0 = entropia(x0.^(0:N-1)*(1-x0)/(1-x0^N));
	x1 = 0.5;
	y1 = entropia(x1.^(0:N-1)*(1-x1)/(1-x1^N));
	while abs(y1-y0) > 1e-10
		x = x0+(H-y0)/(y1-y0)*(x1-x0);
		if x < 0
			x = x0/10;
		end
		p = x.^(0:N-1)*(1-x)/(1-x^N);
		y = entropia(p);
		x0 = x1;
		y0 = y1;
		x1 = x;
		y1 = y;
	end
	
	for i = 1:N
		j = i+floor((N-i+1)*rand);
		if i ~= j
			p([i,j]) = p([j,i]);
		end
	end
end

% calcola l'entropia del vettore di distribuzione di massa di
% probabilita' p (evitando gli zeri)
%
function val = entropia(p)
    val = 0;
	for i = 1:numel(p)
		if p(i) > 0
			val = val - p(i)*log2(p(i));
		end
	end
end
