function LT_demo
	
	NAIF = false;
	
	% non usare K grandi con l'implementazione Naif, che ha
	% complessità O(K^2)
	if NAIF
		K = 64;
		N_base = 96;
		N_overhead = 96;
	else
		K = 1000;
		N_base = 1100;
		N_overhead = 400;
% 		K = 10000;
% 		N_base = 11000;
% 		N_overhead = 2000;
	end
	p = 1 - K/N_base;
	N = N_base+N_overhead;
	fprintf(['Digital fountain LT code demo: Erasure rate = %g '...
		'-- K/N = %d/%d (%.1f%% of BEC Capacity)\n'],p,K,N,K/N/(1-p)*100);
	
	% vettore con gli estremi destri degli intervalli associati alle
	% probabilità (valori della funzione di distribuzione cumulativa)
% 	R = soliton(K);
	R = robust_soliton(K);
	
	% attendi questo numero di eventi di errore
	failures_threshold = 30;
	
	% attendi al più questo tempo
	seconds_budget = 30;
	tic;
	
	trials = 0;
	failures = 0;
	undelivered = 0;
	while failures < failures_threshold && toc < seconds_budget
		trials = trials + 1;
		
		if NAIF
			
			G = [];
			for n = 1:N
				
				% trova il grado della colonna
				L = degree_from_distribution(R);
				
				% genera le posizioni di L uno a caso
				g = [ones(L,1);zeros(K-L,1)];
				for i = 1:K-1
					dif = floor((K-i+1)*rand);
					if dif
						g([i,i+dif]) = g([i+dif,i]);
					end
				end
				
				% canale BEC
				if rand > p
					G = [G,g];
				end
				
			end
			
			% decodifica a sostituzione 'semplice'
			while numel(G)
				
				found = false;
				for n = 1:size(G,2)
					if sum(G(:,n)) == 1
						found = true;
						% trova la riga in cui sta l'uno
						for k = 1:size(G,1)
							if G(k,n) == 1
								break
							end
						end
						break
					end
				end
				
				if found
					% elimina riga k-esima e colonna n-esima
					G(:,n) = [];
					G(k,:) = [];
				else
					failures = failures + 1;
					undelivered = undelivered + size(G,1);
					fprintf(['Overhead: %g%% -- Pr{decoding failure} = %g (%d/%d)'...
						' -- Undelivered nodes per failure: %g\n'],...
						N_overhead/N*100,failures/trials,failures,trials,...
						undelivered/failures);
					break
				end
				
			end
			
		else
			
			% quante parole vengono ricevute delle N trasmesse?
			N_rx = N;
			for n = 1:N
				if rand < p
					N_rx = N_rx - 1;
				end
			end
			
			% costruzione della matrice G (vengono memorizzate
			% solo le posizioni degli '1' nelle righe e nelle
			% colonne)
			G_col = cell(N_rx,1);
			deg_col = zeros(1,N_rx);
			G_row = cell(K,1);
			deg_row = zeros(1,K);
			used = zeros(1,K);
			for n = 1:N_rx
				
				% grado della colonna
				L = degree_from_distribution(R);
				deg_col(n) = L;
				
				% L indici distinti tra 1 e K
				ind = 1:L;
				used(1:L) = 1;
				for i = 1:L
					ind_in = ind(i);
					ind_out = ind_in+floor((K-i+1)*rand);
					if ~used(ind_out)
						used(ind_in) = 0;
						deg_row(ind_in) = deg_row(ind_in) - 1;
						used(ind_out) = 1;
						deg_row(ind_out) = deg_row(ind_out) + 1;
						ind(i) = ind_out;
					end
				end
				G_col{n} = sort(ind);
				used(ind) = 0;
			end
			for i = 1:K
				G_row{i} = zeros(1,deg_row(i));
				deg_row(i) = 0;
			end
			for n = 1:N_rx
				for i = 1:deg_col(n)
					k = G_col{n}(i);
					deg_row(k) = deg_row(k) + 1;
					G_row{k}(deg_row(k)) = n;
				end
			end
			
			% decodifica 'semplice' (risolvi le colonne che
			% coinvolgono una sola incognita, e abbassa di grado
			% le colonne che la contengono)
			used = cell(N_rx,1);
			ripple = zeros(1,N_rx);
			ripple_k = zeros(1,N_rx);
			NR = 0;
			for n = 1:N_rx
				used{n} = zeros(1,deg_col(n));
				if deg_col(n) == 1
					NR = NR + 1;
					ripple(NR) = n;
					ripple_k(NR) = G_col{n};
				end
			end
			while NR > 0
				next_ripple = ripple;
				next_ripple_k = ripple_k;
				next_NR = 0;
				for r = 1:NR
					n = ripple(r);
					k = ripple_k(r);
					% questo può capitare perché eliminiamo tutte
					% le colonne nel ripple a ciascuna passata,
					% ed una colonna potrebbe apparire più di una
					% volta tra i vicini di colonne nel ripple
					if deg_col(n) < 1
						continue
					end
					for j = 1:deg_row(k)
						% elimina l'incognita k dalla colonna i
						i = G_row{k}(j);
						% ricerca binaria di k nella colonna i
						M = numel(G_col{i});
						m = ceil(M/2);
						d = ceil(m/2);
						while G_col{i}(m) ~= k
							if G_col{i}(m) > k
								m = max(m-d,1);
							else
								m = min(m+d,M);
							end
							d = ceil(d/2);
						end
						used{i}(m) = 1;
						deg_col(i) = deg_col(i) - 1;
						if deg_col(i) == 1
							next_NR = next_NR + 1;
							next_ripple(next_NR) = i;
							for temp_k = 1:numel(used{i})
								if used{i}(temp_k) == 0
									break
								end
							end
							next_ripple_k(next_NR) = G_col{i}(temp_k);
						end
					end
					deg_row(k) = 0;
				end
				NR = next_NR;
				ripple = next_ripple;
				ripple_k = next_ripple_k;
			end
			if sum(deg_row) ~= 0
				failures = failures + 1;
				for i = 1:numel(G_row)
					undelivered = undelivered ...
						+ (numel(G_row{i}) > 0);
				end
				fprintf(['Overhead: %g%% -- Pr{decoding failure} = %g (%d/%d)'...
					' -- Undelivered nodes per failure: %g\n'],...
					N_overhead/N*100,failures/trials,failures,trials,...
					undelivered/failures);
			end
		end
	end
	fprintf(['Failure Rate: %.1f%% -- %g encodings-decodings/s\n'],...
		failures/trials*100,trials/toc);



% questa funzione è analoga alla funzione di generazione di
% simboli da una sorgente DMS
function L = degree_from_distribution(R)
	x = rand;
	L = 1;
	while x > R(L)
		L = L + 1;
	end
	% L ha la distribuzione richiesta...



% vogliamo generare L (grado della combinazion lineare)
% secondo la distribuzione solitonica:
%  Pr{L == 1} = rho(1) = 1/K
%  Pr{L == n} = rho(n) = 1/(n*(n-1))
function R = soliton(K)
	R = zeros(1,K);
	R(1) = 1/K;
	for n = 2:K-1
		R(n) = R(n-1) + 1/(n*(n-1));
	end
	R(K) = 1;



function R = robust_soliton(K,delta,c)
	if nargin < 2
		% probabilità di fallimento tollerata (bound)
		delta = 0.5;
	end
	if nargin < 3
		% coefficiente di log(K/delta)*sqrt(K) per il numero di
		% simboli che desideriamo possano essere risolti ad ogni
		% passo
		c = 0.1;
	end	
	
	% distribuzione solitonica
	rho = zeros(1,K);
	rho(1) = 1/K;
	for k = 2:K
		rho(k) = 1/(k*(k-1));
	end
	
	% numero atteso di elementi che possono esser erisolti ad
	% ogni passo
	S = floor(c*log(K/delta) * sqrt(K));
	if S < 1
		S = 1;
	end
	
	% termine correttivo
	tau = zeros(1,K);
	for i = 1:round(K/S)-1
		tau(i) = S/(K*i);
	end
	tau(round(K/S)) = S*log(S/delta)/K;
	
	% distribuzione 'robusta'
	beta = (rho + tau)/sum(rho + tau);
	
	% estremi destri degli intervalli
	R = zeros(1,K);
	R(1) = beta(1);
	for k = 2:K-1
		R(k) = R(k-1) + beta(k);
	end
	R(K) = 1;


