function GF_demo

	% caratteristica del campo GF(q) (potenza di un numero primo: p^m)
	p = 2;
	m = 8;
	GF = GF_init(p,m);
	% caratteristica del campo GF(q)
	q = p^m;
	GFx = GFx_init(GF);
	t = 8;

	% polinomio generatore di un codice BCH
	g_BCH = BCH_poly(GFx,t);
	fprintf(['generatore di un codice BCH con rate (%d,%d)'...
		' e capacità correttiva %d (%d bit): %s\n'],...
		GFx.q-numel(g_BCH),GFx.q-1,t,t*log2(p),poly2str(g_BCH));

	% polinomio generatore di un codice Reed-Solomon
	g_RS = RS_poly(GFx,t);
	fprintf(['generatore di un codice Reed-Solomon con rate %d x (%d,%d)'...
		' e capacità correttiva %d (tra %d e %d bit): %s\n'],...
		m,GFx.q-numel(g_RS),GFx.q-1,t,t,t*log2(q),xpoly2str(g_RS));

	subplot(2,1,1)
	G = GFx_dft(GFx,poly2xpoly(GFx,g_BCH));
	stem(-1:GFx.q-2,log2(G+2))
	set(gca,'XLim',[-1,GFx.q-2])
	ylabel('DFT')

	subplot(2,1,2)
	G = GFx_dft(GFx,g_RS);
	stem(-1:GFx.q-2,log2(G+2))
	set(gca,'XLim',[-1,GFx.q-2])
	ylabel('DFT')

% 	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% GF(p) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 	disp(['============== polinomi di grado ',num2str(m-1),...
% 		' su GF(',num2str(p),') ==============']);
% 	for num_a = 0:q-1
% 		for num_b = 0:q-1
% 			a = num2poly(GF,num_a);
% 			b = num2poly(GF,num_b);
% 			a_plus_b=poly_add(GF,a,b);
% 			a_minus_b=poly_sub(GF,a,b);
% 			a_times_b=poly_mul(GF,a,b);
% 			[quot,rem]=poly_div(GF,a,b);
% 			disp(['a=(',poly2str(a),') b=(',poly2str(b),...
% 				') a+b=(',poly2str(a_plus_b),...
% 				') a-b=(',poly2str(a_minus_b),...
% 				') a*b=(',poly2str(a_times_b),...
% 				') a/b=(',poly2str(quot),') (rem: ',poly2str(rem),')']);
% 		end
% 	end
% 	
% 	% stampa i primi polinomi (in ordine di valore numerico convenzionale)
% 	n_poly = p^(m+1);
% 	str = ['i primi ',num2str(n_poly),' polinomi:'];
% 	for i = 0:n_poly-1
% 		str = [str,' (',poly2str(num2poly(GF,i)),')'];
% 	end
% 	disp(str);
% 	
% 	% trova e stampa i polinomi irriducibili di grado fino ad m
% 	irr_poly_num = list_irr_poly(GF,m);
% 	str = ['i polinomi a coefficienti in GF(',num2str(p),...
% 		'), irriducibili e con grado fino a ',num2str(m),':'];
% 	for i = 1:numel(irr_poly_num)
% 		poly = num2poly(GF,irr_poly_num(i));
% 		str = [str,' (',poly2str(poly),')[',...
% 			num2str(order_of(GF,poly,m)),']'];
% 	end
% 	disp(str);
% 	
% 	% trova e stampa il primo polinomio primitivo (ovvero rispetto al quale
% 	% gli elementi del campo diversi da 0 ed 1 hanno ordine ordine p^m-1)
% 	prim_poly_num = first_prim_poly(GF,m);
% 	str = ['il primo polinomio a coefficienti in GF(',num2str(p),...
% 		'), primitivo e di grado ',num2str(m),':'];
% 	str = [str,' (',poly2str(num2poly(GF,prim_poly_num)),')'];
% 	disp(str);
% 	
% 	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% GF(p^m) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 	disp(['============== campo finito GF(',...
% 		num2str(p),'^',num2str(m),') ==============']);
% 	disp(['primitive polynomial: ',poly2str(GF.prim_poly)]);
% 	for num_a = 0:q-1
% 		for num_b = 0:q-1
% 			a = num2poly(GF,num_a);
% 			b = num2poly(GF,num_b);
% 			a_plus_b=GF_add(GF,a,b);
% 			a_minus_b=GF_sub(GF,a,b);
% 			a_times_b=GF_mul(GF,a,b);
% 			a_over_b=GF_div(GF,a,b);
% 			disp(['a=(',poly2str(a),') b=(',poly2str(b),...
% 				') a+b=(',poly2str(a_plus_b),...
% 				') a-b=(',poly2str(a_minus_b),...
% 				') a*b=(',poly2str(a_times_b),...
% 				') a/b=(',poly2str(a_over_b),')']);
% 		end
% 	end
% 	
% 	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% GF(q^r) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 	disp(['============== polinomi di grado ',num2str(r-1),...
% 		' su GF(',num2str(q),') ==============']);
% 	for num_a = 0:q^r-1
% 		for num_b = 0:q^r-1
% 			a = num2xpoly(GFx,num_a);
% 			b = num2xpoly(GFx,num_b);
% 			a_plus_b=xpoly_add(GFx,a,b);
% 			a_minus_b=xpoly_sub(GFx,a,b);
% 			a_times_b=xpoly_mul(GFx,a,b);
% 			[quot,rem]=xpoly_div(GFx,a,b);
% 			disp(['a=(',xpoly2str(a),') b=(',xpoly2str(b),...
% 				') a+b=(',xpoly2str(a_plus_b),...
% 				') a-b=(',xpoly2str(a_minus_b),...
% 				') a*b=(',xpoly2str(a_times_b),...
% 				') a/b=(',xpoly2str(quot),') (rem: ',xpoly2str(rem),')']);
% 		end
% 	end
end

% inizializzazione delle variabili che descrivono il campo finito GF(p^m)
%
% il parametro p deve essere un numero primo, altrimenti la costruzione
% della tabella di divisione fallisce (e GF(p) non è un campo)
%
% il parametro prim_poly definisce (nella notazione numerica convenzionale
% di num2poly/poly2num) il polinomio primitivo da utilizzare per la
% costruzione del campo GF(p^m) (in particolare per ottenere la
% moltiplicazione invertibile); se non viene passato un valore allora viene
% usato il primo (sempre rispetto a num2poly) polinomio primitivo di grado
% pari ad m
%
function GF = GF_init(p,m,prim_poly_num)

	% Un numero primo
	GF.p = p;

	% Il numero di componenti dei vettori del campo.
	GF.m = m;
	% se riguardati come polinomi, il grado massimo dei polinomi è m - 1 ed
	% il polinomio primitivo che serve a definire la moltiplicazione
	% invertibile tra elementi deve avere grado m

	% Numero di elementi del campo
	q = p^m;
	GF.q = q;

	% Nella rappresentazione polinomiale del campo occorre sapere
	% effettuare le operazioni aritmetiche tra gli elementi di GF(p) e le
	% operazioni polinomiali tra polinomi a coefficienti in GF(p).
	%
	% Le operazioni aritmetiche +, - e * si effettuano in modulo p, mentre
	% per la divisione occorre una ricerca. Per efficienza, tabelliamo i
	% valori degli inversi degli elementi non nulli di GF(p)
	inv = zeros(1,p);
	for i = 0:p-1
		for j = 0:p-1
			if mod(i*j,p) == 1
				inv(i) = j;
				inv(j) = i;
			end
		end
	end
	GF.inv = inv;

	% se il polinomio primitivo non viene passato come parametro, trovalo
	if nargin < 3
		if false
		elseif q == 9; prim_poly_num = 14; % X^2 + X + 2
		else
			prim_poly_num = first_prim_poly(GF,m);
			disp(['default polynomial [',num2str(prim_poly_num),']: ',...
				poly2str(num2poly(GF,prim_poly_num))]);
		end
	end
	GF.prim_poly = num2poly(GF,prim_poly_num);

	% 	% tabella dei numeri convenzionali degli inversi dei polinomi non
	% 	% nulli del campo GF(p^m)
	% 	inv_poly = zeros(1,q-1);
	% 	for num_a = 1:q-1
	% 		if inv_poly(num_a) == 0
	% 			a = num2poly(GF,num_a);
	% 			for num_b = 1:q-1
	% 				b = num2poly(GF,num_b);
	% 				a_times_b = GF_mul(GF,a,b);
	% 				if numel(a_times_b) == 1 && a_times_b == 1
	% 					inv_poly(num_a) = num_b;
	% 					inv_poly(num_b) = num_a;
	% 				end
	% 			end
	% 		end
	% 	end
	% 	GF.inv_poly = inv_poly;
end

% operazioni aritmetiche su polinomi a coefficienti in GF(p)
%
% NOTE: In questa implementazione i polinomi a coefficienti in GF(p) sono
%       rappresentati con i vettori dei coefficienti in ordine di grado
%       Il vettore 'a' rappresenta il polinomio:
%
%          a(1) + a(2) X + a(3) X^2 + ... a(end) X^(numel(a)-1)
%
%       Nelle funzioni si assume che il coefficiente di grado più alto,
%       a(end), sia diverso da zero; il polinomio nullo è rappresentato dal
%       vettore vuoto, []
%
function c = poly_add(GF,a,b)
	na = numel(a);
	nb = numel(b);
	if na > nb
		c = a;
		if nb > 0; c(1:nb) = mod(c(1:nb) + b, GF.p); end
	elseif na < nb
		c = b;
		if na > 0; c(1:na) = mod(c(1:na) + a, GF.p); end
	elseif na == 0 % && nb == 0
		c = [];
	else
		c = mod(a + b, GF.p);
		temp = numel(c);
		while temp > 0 && c(temp) == 0; temp = temp - 1; end
		c(temp+1:end) = [];
	end
end
% function c = poly_sub(GF,a,b)
% 	na = numel(a);
% 	nb = numel(b);
% 	if na > nb
% 		c = a;
% 		if nb > 0; c(1:nb) = mod(c(1:nb) - b, GF.p); end
% 	elseif na < nb
% 		c = mod(-b,GF.p);
% 		if na > 0; c(1:na) = mod(c(1:na) + a, GF.p); end
% 	elseif na == 0 % && nb == 0
% 		c = [];
% 	else
% 		c = mod(a - b, GF.p);
% 		temp = numel(c);
% 		while temp > 0 && c(temp) == 0; temp = temp - 1; end
% 		c(temp+1:end) = [];
% 	end
% end
function c = poly_mul(GF,a,b)
	na = numel(a);
	nb = numel(b);
	if na == 0 || nb == 0
		c = [];
	elseif na == 1 && nb == 1
		c = mod(a * b, GF.p);
	else
		c = mod(conv(a,b), GF.p);
	end
end
function [q,r] = poly_div(GF,a,b)
	na = numel(a);
	nb = numel(b);
	if nb == 0
		display('*** division by zero ***');
		q = [];
		r = a;
	elseif na < nb || na == 0
		q = [];
		r = a;
	elseif nb == 1
		q = mod(a * GF.inv(b),GF.p);
		r = [];
	else
		q = zeros(1,na-nb+1);
		r = a;
		shift = na - nb;
		while shift >= 0
			if r(nb+shift)
				q(1+shift) = mod(r(nb+shift) * GF.inv(b(nb)), GF.p);
				% r = mod(r(1:nb+shift)-[zeros(1,shift),q(1+shift)*b], GF.p);
				for i = shift+1:shift+nb
					r(i) = mod(r(i)-q(1+shift)*b(i-shift), GF.p);
				end
			end
			shift = shift - 1;
		end
		temp = numel(r);
		while temp > 0 && r(temp) == 0; temp = temp - 1; end
		r(temp+1:end) = [];
	end
end

% ritorna una stringa che rappresenta il polinomio in notazione "umana"
%
function str = poly2str(a)
	if numel(a) == 0
		str = '0';
		return
	end
	str=[];
	first = true;
	for i = numel(a):-1:1
		if a(i) > 0
			% separa ogni monomio dopo il primo con un '+'
			if ~first
				str=[str,' + '];
			end
			% stampa il monomio di grado (i-1)
			if a(i) > 1 || i == 1
				str=[str,num2str(a(i))];
			end
			if i > 2
				str=[str,'X^',num2str(i-1)];
			elseif i == 2
				str=[str,'X'];
			end
			first = false;
		end
	end
end

% conversione da numero a polinomio a coefficienti in GF(p) (dec2base)
% e viceversa (base2dec)
%
function a = num2poly(GF,n)
	if n == 0
		a = [];
	else
		% na = ceil(log(n)/log(GF.p));
		x = 1;
		na = 0;
		while x <= n
			x = x * GF.p;
			na = na + 1;
		end
		% potenze successive di GF.p
		x = zeros(1,na);
		x(1) = 1;
		for i = 2:na
			x(i) = x(i-1) * GF.p;
		end
		% conversione
		a = zeros(1,na);
		for i = na:-1:1
			a(i) = floor(n / x(i));
			n = n - x(i) * a(i);
		end
	end
end
function n = poly2num(GF,a)
	na = numel(a);
	n = 0;
	for i = na:-1:1
		n = n * GF.p + a(i);
	end
end

% genera casualmente un polinomio di grado n a coefficienti in GF(p)
%
% function a = poly_rand(GF,n)
% 	na = 1+n;
% 	a = floor(GF.p * rand(1,na));
% 	while a(end) == 0
% 		a = floor(GF.p * rand(1,na));
% 	end
% end

% ricerca dei polinomi monici irriducibili (divisibili solo per se stessi e
% per 1) fino ad un grado massimo assegnato
%
% i primi polinomi sono:
% '0'             non irriducibile in quanto divisibile con resto nullo da
%                 qualsiasi polinomio
% '1'             irriducibile per definizione
% 'X + a'         irriducibile
% 'X^2 + a X + b': rid. se esistono c,d in GF(p) con c+d=a e c*d=b
% 'X^3 + a X^2 + b X + c': rid. se esistono d,e,f con d+f=a, e+d*f=b, e*f=c
%
% algoritmo: crivello di Eratostene: una tabella di variabili booleane che
%            vengono azzerate se il polinomio ad esse corrispondente è
%            multiplo di un polinomio irriducibile
function irr_poly_num = list_irr_poly(GF,n_max,trace)
	if nargin < 3
		trace = [];
	end
	irr_poly_num = [];
	[ip,tab] = first_irr_poly(GF,n_max,trace);
	while ip
		irr_poly_num = [irr_poly_num,ip];
		[ip,tab] = next_irr_poly(GF,n_max,trace,ip,tab);
	end
end
function [ip,tab] = first_irr_poly(GF,n_max,trace)
	if nargin < 3
		trace = [];
	end

	tab = ones(1,GF.p^(n_max+1));

	% il polinomio 1 è un caso particolare di polinomio irriducibile
	tab(1) = 0;

	% trova il prossimo (ovvero il primo)
	[ip,tab] = next_irr_poly(GF,n_max,trace,1,tab);
end
function [ip,tab] = next_irr_poly(GF,n_max,trace,last_ip,tab)
	ip = 0;
	found = false;
	i = last_ip + 1;
	while ~found && i < GF.p^(n_max+1)
		a = num2poly(GF,i);
		if a(end) == 1 % il polinomio è monico
			if tab(i) == 1 % il polinomio è irriducibile

				% se il grado del polinomio è quello desiderato, aggiungilo
				% alla lista da ritornare
				if (numel(trace)==0 || a(end-1) == trace)
					ip = i;

					% esci dal loop
					found = true;
				end

				% marca come non irriducibili tutti i multipli di grado non
				% maggiore di n_max
				na = numel(a);
				for j = 2:GF.p^(n_max-na+2)-1
					tab(poly2num(GF,poly_mul(GF,a,num2poly(GF,j)))) = 0;
				end
			end
		end
		i = i + 1;
	end
end

% calcola l'ordine moltiplicativo del polinomio X
%
function count = order_of(GF,poly)
	% il polinomio X
	X = [0 1];
	% il polinomio a = X^0
	a = 1;

	% in GF.p^n - 1 passi non si devono ottenere valori ripetuti
	tab = zeros(1,GF.p^(numel(poly)-1));
	count = 0;
	while true
		val = poly2num(GF,a);
		if val==0 || tab(val)
			break
		end

		count = count + 1;
		tab(val) = 1;

		% a <- a*X mod p(X)
		[ignored,a] = poly_div(GF,poly_mul(GF,a,X),poly);
	end
end

% ricerca di un polinomio monico primitivo di grado e traccia assegnata
%
% un polinomio p(X) è primitivo quando X^k MOD p(X) al variare di k genera
% tutti gli elementi non nulli del campo (ovvero ha periodo massimo, pari a
% GF.p^n-1)
%
% function prim_poly = list_prim_poly(GF,n,trace)
% 	if nargin < 3
% 		trace = [];
% 	end
% 
% 	% la ricerca è effettuata nell'insieme dei polinomi irriducibili
% 	prim_poly = [];
% 	[pp,tab_irr] = first_prim_poly(GF,n,trace);
% 	while pp
% 		prim_poly = [prim_poly,pp];
% 		[pp,tab_irr] = next_prim_poly(GF,n,trace,pp,tab_irr);
% 	end
% end
function [pp,tab_irr] = first_prim_poly(GF,n,trace)
	if nargin < 3
		trace = [];
	end

	tab_irr = ones(1,GF.p^(n+1));

	% il polinomio 1 è un caso particolare di polinomio irriducibile
	tab_irr(1) = 0;

	% trova il prossimo (ovvero il primo)
	[pp,tab_irr] = next_prim_poly(GF,n,trace,1,tab_irr);
end
function [pp,tab_irr] = next_prim_poly(GF,n,trace,last_pp,tab_irr)
	pp = 0;

	% la ricerca è effettuata nell'insieme dei polinomi irriducibili
	[ip,tab_irr] = next_irr_poly(GF,n,trace,last_pp,tab_irr);
	found = false;
	while ~found && ip
		temp = num2poly(GF,ip);
		if numel(temp) == n+1

			% calcola l'ordine moltiplicativo degli elementi di GF(p^m)
			% rispetto al polinomio pp
			if order_of(GF,temp) == GF.p^n-1
				pp = ip;

				% esci dal loop
				found = true;
			end
		end
		if ~found
			% passa al prossimo polinomio irriducibile
			[ip,tab_irr] = next_irr_poly(GF,n,trace,ip,tab_irr);
		end
	end
end

% operazioni aritmetiche nel campo finito GF(p^m)
%
function c = GF_add(GF,a,b)
	c = poly_add(GF,a,b);
end
% function c = GF_sub(GF,a,b)
% 	c = poly_sub(GF,a,b);
% end
function c = GF_mul(GF,a,b)
	[ignored,c] = poly_div(GF,poly_mul(GF,a,b),GF.prim_poly);
end
% function c = GF_div(GF,a,b)
% 	if numel(b) == 0
% 		display('*** division by zero ***');
% 		c = [];
% 	else
% 		b_inv = num2poly(GF,GF.inv_poly(poly2num(GF,b)));
% 		[ignored,c] = poly_div(GF,poly_mul(GF,a,b_inv),GF.prim_poly);
% 	end
% end

% un elemento casuale di GF(p^m)
%
% function a = GF_rand(GF)
% 	a = num2poly(GF,floor(GF.q * rand));
% end

% inizializzazione delle variabili che descrivono il campo finito GF(q)
%
% il valore di q è il numero di elementi del campo finito GF(p^m) a cui
% appartengono i coefficienti
%
% La notazione scelta per rappresentare gli elementi di GF(p^m) è quella
% esponenziale/logaritmica, in cui tutti gli elementi non nulli del campo
% sono rappresentati come potenze di un elemento primitivo (il polinomio X)
%
% L'elemento nullo del campo (elemento neutro rispetto alla somma) è
% rappresentato da alpha^-inf e convenzionalmente dall'esponente -1, gli
% altri elementi del campo hanno esponenti compresi tra 0 (l'elemento
% neutro rispetto alla somma) e q-2 (sono q in totale)
%
function GFx = GFx_init(GF)

	% numero di elementi del campo
	q = GF.q;
	GFx.q = q;

	% viene anche memorizzata la tabella dei logaritmi degli elementi di
	% GF(p) in base alpha
	p = GF.p;
	GFx.p = p;
	GFx.log = zeros(p,1);
	GFx.log(1+0) = -1;
	GFx.exp = -1*ones(q,1);
	GFx.exp(2+-1) = 0;

	% costruzione della tabella dei logaritmi di Zech:
	%            alpha^z_n = 1 + alpha^n
	% z_n è memorizzato in z(2+n) con n = -1, 0, 1, ... q-2
	%
	% i logaritmi di Zech sono utili ad accelerare la somma, in quanto
	%          a^n + a^m = a^n (1 + a^(m-n)) = a^(n+z(m-n))
	%
	% per costruire la tabella, costruisco le due tabelle dei valori
	% numerici convenzionali delle potenze alpha^n (potenze successive del
	% polinomio X) e di 1+alpha^n (le stesse sommate al polinomio 1)
	%
	num_alpha  = zeros(q,1); % num_alpha(2+n)  = poly2num(alpha^n)
	num_alphaz = zeros(q,1); % num_alphaz(2+n) = poly2num(1+alpha^n == alpha^z_n)

	% casi particolari del polinomio nullo
	%
	num_alpha(-1+2) = 0;  % num(alpha^-inf) == 0
	num_alphaz(-1+2) = 1; % num(1+alpha^-inf) == 1

	a = 1;     % alpha^0
	X = [0 1]; % alpha^1
	for n = 0:q-2
		% alpha^n, alpha^z(n)
		num_alpha(n+2) = poly2num(GF,a);
		num_alphaz(n+2) = poly2num(GF,GF_add(GF,a,1));

		% salva i logaritmi degli elementi di GF(p)
		if num_alpha(n+2) < p
			GFx.log(1+num_alpha(n+2)) = n;
			GFx.exp(2+n) = num_alpha(n+2);
		end

		a = GF_mul(GF,a,X);
	end

	% la tabella dei logaritmi di Zech si ottiene trovando le
	% corrispondenze tra i due vettori costruiti
	z = zeros(q,1);
	for n = -1:q-2
		temp = num_alphaz(n+2);
		for j = -1:q-2
			if temp == num_alpha(j+2)
				z(n+2) = j;
				break
			end
		end
	end
	GFx.z = z;

	% ricerca dei polinomi minimali (ovvero identificazione di quale dei
	% polinomi irriducibili si annulla in quale elemento del campo)
	min_pol = zeros(q,1);
	irr_poly = list_irr_poly(GF,GF.m);
	for ip_num = 1:numel(irr_poly)
		ip = poly2xpoly(GFx,num2poly(GF,irr_poly(ip_num)));
		for i = -1:q-2
			if min_pol(2+i) == 0
				if xpoly_eval(GFx,ip,i) == -1
					min_pol(2+i) = ip_num;
					if i>0
						j=mod(i*p,q-1);
						while min_pol(2+j) == 0
							min_pol(2+j) = ip_num;
							j=mod(j*p,q-1);
						end
					end
				end
			end
		end
	end
	GFx.irr_poly = irr_poly;
	GFx.min_pol = min_pol;

	% ---------------------------------------------------------------------
	% altri campi della struttura GF che occorrono per potere passare GFx
	% alle funzioni che usano GF
	GFx.prim_poly = GF.prim_poly;
	GFx.inv = GF.inv;
end

% conversione da polinomio a coefficienti in GF(p) a polinomio a
% coefficienti in GF(q)
%
function xpoly = poly2xpoly(GFx,poly)
	xpoly = zeros(1,numel(poly));
	for i = 1:numel(poly)
		xpoly(i) = GFx.log(1+poly(i));
	end
end

% valutazione di un polinomio a coefficienti in GF(q)
%
function val = xpoly_eval(GFx,poly,x)
	q = GFx.q;
	val = poly(end);
	for i = numel(poly)-1:-1:1
		% % val = val * x + poly(i);
		% val = xscalar_add(GFx,xscalar_mul(GFx,x,val),poly(i));
		if x < 0 || val < 0
			x_times_val = -1;
		else
			x_times_val = mod(x + val, q-1);
		end
		if x_times_val < 0
			val = poly(i);
		elseif poly(i) < 0
			val = x_times_val;
		else
			t = GFx.z(2+abs(poly(i) - x_times_val));
			if t < 0
				val = -1;
			else
				val = mod(min(x_times_val,poly(i)) + t, q-1);
			end
		end
	end
end

% operazioni artimetiche tra scalari in GF(p^m)
%
function c = xscalar_add(GFx,a,b)
	if a < 0
		c = b;
	elseif b < 0
		c = a;
	else
		t = GFx.z(2+abs(b-a));
		if t < 0
			c = -1;
		else
			c = mod(min(a,b) + t,GFx.q-1);
		end
	end
end
% function c = xscalar_sub(GFx,a,b)
% 	% somma l'elemento opposto
% 	if b == -1
% 		minus_b = b;
% 	else
% 		minus_b = mod(b + GFx.log(end),GFx.q-1);
% 	end
% 	if a < 0
% 		c = minus_b;
% 	elseif minus_b < 0
% 		c = a;
% 	else
% 		t = GFx.z(2+abs(minus_b-a));
% 		if t < 0
% 			c = -1;
% 		else
% 			c = mod(min(a,minus_b) + t,GFx.q-1);
% 		end
% 	end
% end
function c = xscalar_mul(GFx,a,b)
	if a < 0 || b < 0
		c = -1;
	else
		c = mod(a + b, GFx.q-1);
	end
end
% function c = xscalar_div(GFx,a,b)
% 	if a < 0
% 		c = -1;
% 	elseif b < 0
% 		display('*** division by zero ***');
% 		c = -1;
% 	else
% 		c = mod(a - b, GFx.q-1);
% 	end
% end

% operazioni aritmetiche tra polinomi a coefficienti in GF(p^m)
%
% function c = xpoly_add(GFx,a,b)
% 	na = numel(a);
% 	nb = numel(b);
% 	if na > nb
% 		c = a;
% 		for i = 1:nb; c(i) = xscalar_add(GFx,c(i),b(i)); end
% 	elseif na < nb
% 		c = b;
% 		for i = 1:na; c(i) = xscalar_add(GFx,c(i),a(i)); end
% 	elseif na == 0 % && nb == 0
% 		c = [];
% 	else
% 		nc = na;
% 		c = zeros(1,nc); for i=1:nc; c(i) = xscalar_add(GFx,a(i),b(i)); end
% 		temp = nc;
% 		while temp > 0 && c(temp) < 0; temp = temp - 1; end
% 		c(temp+1:nc) = [];
% 	end
% end
% function c = xpoly_sub(GFx,a,b)
% 	na = numel(a);
% 	nb = numel(b);
% 	if na > nb
% 		c = a;
% 		for i = 1:nb; c(i) = xscalar_sub(GFx,c(i),b(i)); end
% 	elseif na < nb
% 		c = [a,-1*ones(1,nb-na)];
% 		for i = 1:nb; c(i) = xscalar_sub(GFx,c(i),b(i)); end
% 	elseif na == 0 % && nb == 0
% 		c = [];
% 	else
% 		nc = na;
% 		c = zeros(1,nc); for i=1:nc; c(i) = xscalar_sub(GFx,a(i),b(i)); end
% 		temp = nc;
% 		while temp > 0 && c(temp) < 0
% 			temp = temp - 1;
% 		end
% 		c(temp+1:nc) = [];
% 	end
% end
function c = xpoly_mul(GFx,a,b)
	na = numel(a);
	nb = numel(b);
	if na == 0 || nb == 0
		c = [];
	else
		% convoluzione
		c = zeros(1,na+nb-1);
		for k = 1:na+nb-1
			c(k) = -1;
			for i = 1:na
				if k-i >= 0 && 1+k-i <= nb
					c(k) = xscalar_add(GFx,c(k),xscalar_mul(GFx,a(i),b(1+k-i)));
				end
			end
		end
	end
end
% function [q,r] = xpoly_div(GFx,a,b)
% 	na = numel(a);
% 	nb = numel(b);
% 	if nb == 0
% 		display('*** division by zero ***');
% 		q = [];
% 		r = a;
% 	elseif na < nb || na == 0
% 		q = [];
% 		r = a;
% 	else
% 		q = zeros(1,na-nb+1);
% 		r = a;
% 		i = na-nb;
% 		while i >= 0
% 			if nb+i > numel(r)
% 				q(1+i) = -1;
% 			elseif r(nb+i) >= 0
% 				q(1+i) = xscalar_div(GFx,r(nb+i),b(end));
% 				r = xpoly_sub(GFx,r(1:nb+i),[-1*ones(1,i),xpoly_mul(GFx,q(1+i),b)]) ;
% 			end
% 			i = i - 1;
% 		end
% 		temp = numel(r);
% 		while temp > 0 && r(temp) < 0
% 			temp = temp - 1;
% 		end
% 		r(temp+1:numel(r)) = [];
% 	end
% end

% costruisce una stringa che rappresenta il polinomio in notazione "umana"
%
function str = xpoly2str(a)
	if numel(a) == 0
		str = '0';
		return
	end
	str=[];
	first = true;
	for i = numel(a):-1:1
		if a(i) >= 0
			% separa ogni monomio dopo il primo con un '+'
			if ~first
				str=[str,' + '];
			end
			% stampa il monomio di grado (i-1)
			if a(i) > 1
				str=[str,'alpha^',num2str(a(i))];
			elseif a(i) > 0
				str=[str,'alpha'];
			elseif i == 1
				str=[str,'1'];
			end
			if i > 1 && a(i) > 0
				str=[str,' '];
			end
			if i > 2
				str=[str,'X^',num2str(i-1)];
			elseif i == 2
				str=[str,'X'];
			end
			first = false;
		end
	end
end

% conversione da numero a polinomio a coefficienti in GF(p^m) (dec2base)
% e viceversa (base2dec)
%
% function a = num2xpoly(GFx,n)
% 	if n == 0
% 		a = [];
% 	else
% 		% na = ceil(log(n)/log(GFx.q))
% 		x = 1;
% 		na = 0;
% 		while x <= n
% 			x = x * GFx.q;
% 			na = na + 1;
% 		end
% 		% potenze successive di GFx.q
% 		x = zeros(1,na);
% 		x(1) = 1;
% 		for i = 2:na
% 			x(i) = x(i-1) * GFx.q;
% 		end
% 		% conversione
% 		a = zeros(1,na);
% 		for i = na:-1:1
% 			a(i) = floor(n / x(i)) - 1;
% 			n = n - x(i) * (a(i)+1);
% 		end
% 	end
% end
% function n = xpoly2num(GFx,a)
% 	na = numel(a);
% 	n = 0;
% 	for i = na:-1:1
% 		n = n * GFx.q + (1+a(i));
% 	end
% end

% polinomio generatore di un codice BCH di capacità correttiva t
%
% t0 è l'esponente della prima delle 2*t potenze successive di alpha in cui
% il polinomio deve annullarsi
%
function g = BCH_poly(GFx,t,t0)

	if nargin < 3
		t0 = 1;
	end

	if 2*t+1 > GFx.q-1
		error(['la capacità correttiva richiesta eccede la capacità del'...
			' campo GF(q); usare m >= ',num2str(ceil(log(2*t+2)/log(GFx.p)))]);
	end

	% seleziona i polinomi minimali corrispondenti agli zeri prescelti
	temp = unique(sort(GFx.min_pol(2+(t0:t0+2*t-1))));

	% il polinomio generatore del codice è il prodotto dei polimoni
	% minimali
	g = 1;
	for i = 1:numel(temp)
		f = num2poly(GFx,GFx.irr_poly(temp(i)));
		g = poly_mul(GFx,g,f);
	end
end

% polinomio generatore di un codice Reed-Solomon
%
% è il polinomio equivalente a quello del codice BCH, ma a coefficienti nel
% campo esteso; nel campo esteso il polinomio minimale di alpha^i è
% semplicemente (X - alpha^i)
%
function g = RS_poly(GFx,t,t0)

	if nargin < 3
		t0 = 1;
	end

	if 2*t+1 > GFx.q
		error(['la capacità correttiva richiesta eccede la capacità del'...
			' campo GF(q); usare m >= ',num2str(ceil(log(2*t+1)/log(GFx.p)))]);
	end

	% il polinomio generatore del codice è il prodotto dei polimoni
	% minimali
	g = 0;
	for i = t0:t0+2*t-1
		g = xpoly_mul(GFx,g,[xscalar_mul(GFx,i,GFx.log(end)) 0]);
	end
end

% trasformata di Fourier discreta
%
function G = GFx_dft(GFx,g)
	q = GFx.q;
	G = zeros(1,q);
	for i = 1:q
		G(i) = xpoly_eval(GFx,g,i-2);
	end
end

% % risoluzione di un sistema lineare
% %
% function [ne,x] = GFx_solve(GFx,A,b)
% 	[n,m] = size(A);
% 
% 	% eliminazione di Gauss
% 	ne = 0;
% 	for i = 1:n
% 		% cerca un elemento non nullo nella i-esima colonna dalla i-esima
% 		% posizione in poi (è il pivot)
% 		j = i;
% 		while j <= m && A(j,i) == -1
% 			j = j + 1;
% 		end
% 		if j == m+1
% 			% non si può procedere oltre (il sistema è singolare)
% 			break
% 		else
% 			% una equazione indipendente in più nel sistema
% 			ne = ne + 1;
% 			if j ~= i
% 				temp = A(j,i:end);
% 				A(j,i:end) = A(i,i:end);
% 				A(i,i:end) = temp;
% 				temp = b(j);
% 				b(j) = b(i);
% 				b(i) = temp;
% 			end
% 			for j = i+1:n
% 				lambda = xscalar_div(GFx,A(j,i),A(i,i));
% 				for k = i:m
% 					A(j,k) = xscalar_sub(GFx,A(j,k),xscalar_mul(GFx,A(i,k),lambda));
% 				end
% 				b(j) = xscalar_sub(GFx,b(j),xscalar_mul(GFx,b(i),lambda));
% 			end
% 		end
% 	end
% 
% 	% calcola le soluzioni
% 	x = zeros(1,ne);
% 	for i = ne:-1:1
% 		temp = -1;
% 		for j = i+1:ne
% 			temp = xscalar_add(GFx,temp,xscalar_mul(GFx,A(i,j),x(j)));
% 		end
% 		x(i) = xscalar_div(GFx,xscalar_sub(GFx,b(i),temp),A(i,i));
% 	end
% end
