classdef Elem
properties (SetAccess = public),
field;
prvek; %prvek = element (in Czech)
end
methods
function result = Elem(field,prvek)
%Elem - constructor Elem(field,prvek)
%element of a field
result.field = field;
if ~isvector(prvek)
error('Elem:badInput','input must be polynomial (array) or scalar');
end
[~,r] = deconv(prvek,field.pol);
if r ~= prvek %if the polynomial is greater than irreducible pol.
error('Elem:badInput','the element is greater then irreducible polynomial');
end
if sum(prvek >= field.p) > 0
result.prvek = mod(prvek,field.p);
warning('Elem:badInput','koef. of element must be smaller than "p". It was repaired by mod p');
end
result.prvek = prvek;
end
function result = sameField(elem1,elem2)
%sameField - are the elements from the same field?
%works even for n-dimensional matrices
if ~isa(elem1,'Elem')|| ~isa(elem2,'Elem')
error('not instances of Elem');
end
if isequal(elem1.field,elem2.field)
result = true;
else
error('elements from different fields');
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% OVERLOADING OPERATORS %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function result = plus(a,b)
%plus - overloading of plus operator
%works for all dimensions
if prod(double(sameField(a,b))) %from the same field
if numel(a)==1 && numel(b)==1 %both scalars
result = a;
result.prvek = Field.polAdd(a.prvek,b.prvek);
result.prvek = mod(result.prvek,result.field.p);
else
result = Field.distrib(@plus,a,b);
end
end
end
function result = minus(a,b)
%minus - minus overloading
%works for all dimensions
result = plus(a,-b);
end
function result = uminus(a)
%uminus - unary minus
%works for all dimensions
result = a;
for i=1:numel(result)
result(i).prvek = mod(-result(i).prvek,result(1).field.p);
end
end
function result = uplus(a)
%uplus - unary plus overloading
result = a;
end
function result = times(a,b)
%times - times overloading
%works for all dimensions
if numel(a)==1 && numel(b)==1 %if scalars
if sameField(a,b) %from the same field
result = Elem(a.field,a.prvek);
[~,result.prvek] = deconv(conv(a.prvek,b.prvek),result.field.pol);
result.prvek = mod(result.prvek,result.field.p);
result.prvek = Field.cutPol(result.prvek);
end
else
result = Field.distrib(@times,a,b);
end
end
function result = mtimes(a,b)
%mtimes - mtimes overloading
%for matrix and vectors multipling, 2-D maximum
siz1 = size(a); siz2 = size(b);
if length(siz1)>2 || length(siz1)>2
error('mtimes:toMuchDimension','at least one variable has more than 2 dimensions');
end
if siz1(2) == siz2(1)
result = zeros(a(1).field,[siz1(1),siz2(2)]);
for i=1:siz1(1)
for j=1:siz2(2)
pom = a(i,:).*(b(:,j))';
result(i,j)=sum(pom);
end
end
else
error('mtimes:badDimensions','bad lenght of matrix of vector');
end
end
function result = rdivide(a,b)
%rdivide - rdivide overloading
%works for arrays of the same lenght, one or both can be scalars
if numel(a)==1 && numel(b)==1 %if scalars
if sameField(a,b)
if nonZeroElem(b)
result = a*inverz(b);
else
error('mrdivide:zerodividing','b is zero element');
end
end
else
result = Field.distrib(@rdivide,a,b);
end
end
function result = ldivide(a,b)
%ldivide - ldivide overloding
%works for arrays of the same lenght, one or both can be scalars
result = rdivide(b,a);
end
function result = mrdivide(a,b)
%mrdivide -
%works for arrays of the same lenght, one or both can be scalars
if numel(b)==1
result = rdivide(a,b);
return
end
if size(a)==size(b)
result = a*b^(-1);
else
error('mrdivide:badDimensions','bad matrix size');
end
end
function result = mldivide(a,b)
%mldivide - dlen matic zleva
%funguje pro pole stejnch rozmr, nebo jmenovatel je skalr
result = mrdivide(b,a);
end
function result = power(a,mocnina)
%power - peten opertoru mocnn jednotlivch prvk skalrem
%a me bt pole libovoln dimenze, mocnina mus bt cel slo nebo
%matice celch sel stejnho rozmru
if numel(mocnina) == 1
mocnina = ones(size(a))*mocnina;
else
if size(mocnina) ~= size(a)
error('mpower:badPowerMatDimension','mocnina mus bt matice stejnch rozmr, nebo skalr');
end
end
if int8(mocnina) ~= pom
error('mpower:powerIsNotInteger','mocnina mus bt cel slo');
end
mocnina = int8(mocnina);
siz = size(a);
result=ones(a(1).field,[1,numel(a)]);
for j=1:numel(a) %pouito linern indexovn obecn matice
if mocnina(j) == 0
continue
else
result(j) = a(j);
end
for i=2:abs(mocnina(j))
result(j) = result(j).*a(j);
end
if mocnina < 0
result(j) = inverz(result(j));
end
end
result = reshape(result,siz);
end
function result = mpower(a,mocnina)
%mpower - peten opertoru mocnn matice
%a mus bt tvercov matice
%mocnina mus bt cel slo
if ~isscalar(mocnina)
error('mpower:powerIsNotScalar','mocnina mus bt skalr');
end
if int8(mocnina)~=mocnina
warning('mpower:powerIsNotInt','mocnina byla pevedena na cel slo');
end
mocnina = int8(mocnina);
if mocnina == 0
result = ones(a(1).field, size(a));
return;
end
result=a;
for i=1:mocnina-1
result = result*a;
end
if mocnina < 0
result = inv(result);
end
end
%relacni operatory
function result = ne(a,b)
%ne - peten opertoru nonekvivalence
%vstupy mus mt stejny rozmr
if eq(a,b) == true;
result = false;
else
result = true;
end
end
function result = eq(a,b)
%eq - peten opertoru ekvivalence
%vstupy mus mt stejny rozmr
if numel(a)==1 && numel(b)==1 %jsou-li to prvky a ne matice
if sameField(a,b) %ze stejnho tlesa
if a.prvek == b.prvek
result = true;
else
result = false;
end
end
else %je-li vice prvku
result = Field.distrib(@eq,a,b);
end
end
function result = colon(a,b,c)
%colon - peten opertoru dvojteka
%vstupem jsou poten a konen hodnota vektoru, ppadn i
%prostedn vstup, kter znamen inkrementan krok, vyjde vektor dlky
%maximln charakteristika tlesa
if nargin == 3
elemd = b; %inkrementan krok
b = c; %horn mez
sameField(a,elemd);
sameField(a,b);
end
if nargin == 2
if sameField(a,b)
elemd = ones(a.field,[1,1]);
end
end
if nargin>3 || nargin<2
error('colon:badArguments','mus zadat 2 nebo 3 argumenty');
end
result(1) = a;
i = 1;
while result(i)~=b && i<=a.field.p
i = i+1;
result(i) = result(i-1)+ elemd;
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% OVERLOADING OPERATORS END %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function display(E)
% display - petypovn funkce pro vpis instance
siz = size(E);
disp(' ');
for i=1:numel(E)
[I,J] = ind2sub(siz,i);
I = num2str(I);
J = num2str(J);
disp(['Elem: ',inputname(1),'(',I,',',J,') = [',num2str(E(i).prvek),']']);
end
if ~isempty(E)
display2(E(1).field);
end
disp(' ');
end
function result = sum(elem)
%sum - peten funkce sum
%funguje pro libovoln dimenze
elem = reshape(elem,[],1);
result = elem(1);
for i=2:numel(elem)
result = result + elem(i);
end
end
function result = prod(elem)
%prod - peten funkce prod
%funguje pro libovoln dimenze
elem = reshape(elem,[],1);
result = ones(elem(1).field,[1,1]);
for i=1:numel(elem)
result = result .* elem(i);
end
end
function result = nonZeroElem(elem)
%nonZeroElem - zjist zda se jedn o nulov prvek. Kdy ano, vrt false.
pomocny = zeros(elem(1).field, size(elem));
if elem ~= pomocny
result = true;
else
result = false;
end
end
function result = inverz(a)
%inverzni prvek ke skalaru
if ~isscalar(a.prvek)
[~,u] = deconv(a.prvek,a.field.pol); v = a.field.pol;
x1 = 1; x2 = 0;
while ~(u(end) == 1 && sum(u)==1) %Eukliduv algoritmus pro polynomy
u = Field.cutPol(u);
v = Field.cutPol(v);
[pom1,pom2] = deconv(v,u);
v = u;u = pom2;
pom3 = Field.polAdd(x2,(-1)*conv(pom1,x1));
x2 = x1; x1 = pom3;
u = mod(u,a.field.p);
end
result = Elem(a.field,mod(x1,a.field.p));
for i=1:length(result.prvek)
if mod(result.prvek(i),1)~=0
celaCast = floor(result.prvek(i));
zlomek = result.prvek(i)- celaCast;
pomocnyEl = Elem(Field(result.field.p,1),round(1/zlomek));
pomocnyEl = inverz(pomocnyEl);
if celaCast == 0
celaCast = 1;
end
result.prvek(i) = mod(pomocnyEl.prvek*celaCast,result.field.p);
end
end
else
u = mod(a.prvek,a.field.p); v = a.field.p;
x1 = 1; x2 = 0;
while u ~= 1 %Eukliduv algoritmus pro skalry
pom1 = floor(v/u);
pom2 = v - pom1*u;
v = u;u = pom2;
pom3 = x2 - pom1*x1;
x2 = x1; x1 = pom3;
end
result = Elem(a.field,mod(x1,a.field.p));
end
end
function A = rref(A)
%rref - reduced row echelon form, provede Gauss-Jordanovu eliminaci
zeroEl = zeros(A(1).field,[1,1]);
oneEl = ones(A(1).field,[1,1]);
A = Gauss(A);
s = size(A);
a = min(s);
for i=1:a
A(i,[1:i-1 i+1:end])= A(i,[1:i-1 i+1:end])/A(i,i);
A(i,i)=oneEl;
ss=A(i,i+1:end);
for j=[1:i-1 i+1:a]
if ~isempty(A(j,i+1:end))
A(j,i+1:end)= A(j,i+1:end)-A(j,i).*ss;
end
end
A([1:i-1 i+1:end],i)=zeroEl;
end
end
function [A,det] = Gauss(A)
%Gauss - provede Gaussovu eliminaci na horn trojhelnkovou matici
zeroEl = zeros(A(1).field,[1,1]);
oneEl = ones(A(1).field,[1,1]);
s = 0;
det = oneEl;
[m,~]=size(A);
for j=1:m-1
if A(j,j)== zeroEl
k=j;
for k=k+1:m
if A(k,j)== zeroEl
continue
end
break
end
if A(k,k)==zeroEl
break
end
det = -det;
B=A(j,:);
A(j,:)=A(k,:);
A(k,:)=B;
end
for i=1+s:m-1
L=A(i+1,j)/A(j,j);
A(i+1,:)=A(i+1,:)-L*A(j,:);
end
s=s+1;
end
for i=1:m
det = det*A(i,i);
end
end
function result = det(A)
%determinant matice prvk konench tles
%vyuv funkce rref
s=size(A);
if s(1) ~= s(2)
error('det:notSquare','matice mus bt tvercov');
end
[~,result] = Gauss(A);
end
function result = rank(A)
%hodnost matice
%vyuv funkce rref
zeroEl = zeros(A(1).field,[1,1]);
A = Gauss(A);
s=size(A);
a=min(s);
result = a;
for i=1:a
if A(i,i) == zeroEl
result = i-1;
break;
end
end
end
function A = inv(A)
%inv - inverzn matice. funguje jen pro regulrn (tvercov) matice,
%jinak vrac chybu
zeroEl = zeros(A(1).field,[1,1]);
s = size(A);
if s(1) ~= s(2)
error('inv:notSquare','matice mus bt tvercov');
end
if det(A) ~= zeroEl
A = rref([A,eye(A(1).field,s(1),s(2))]);
A = A(:,s(2)+1:2*s(2));
else
warning('inv:notRegular','matice mus bt regulrn');
end
end
function result = skalarniSoucin(a,b)
%skalarniSoucin - pot skalrn souin dvou stejn dlouhch vektor
if ~isvector(a) || ~isvector(b)
error('hammingNorm:notVector','jeden ze vstup nen vektor');
end
if size(a)~=size(b)
error('hammingNorm:badSize','vektory nejsou stejn velk');
end
result = a*b;
end
function result = linNezavislost(varargin)
%linNezavislost - pov zda je libovoln poet vektor lin. nezvisl
prvni = varargin{1};
s = size(prvni);
if s(1) == 1
A(:,1) = prvni';
for i=2:nargin
if ~isvector(varargin{i})
error('linNezavislost:notVector','jeden ze vstup nen vektor');
end
if length(varargin{i}) ~= length(prvni)
error('linNezavislost:badSize','vektory nejsou stejn velk');
end
A(:,i)= varargin{i}';
end
else
A(:,1) = prvni;
for i=2:nargin
if ~isvector(varargin{i})
error('linNezavislost:notVector','jeden ze vstup nen vektor');
end
if length(varargin{i}) ~= length(prvni)
error('linNezavislost:badSize','vektory nejsou stejn velk');
end
A(:,i)= varargin{i};
end
end
if rank(A)==nargin
result = true;
else
result = false;
end
end
function result = hammingDist(a,b)
%hammingDist - Hammingova vzdlenost dvou vektor. k v kolika prvcch se
% li dva stejn dlouh vektory (kdov slova)
if ~isvector(a) || ~isvector(b)
error('hammingNorm:notVector','jeden ze vstup nen vektor');
end
if size(a)~=size(b)
error('hammingNorm:badSize','vektory nejsou stejn velk');
end
L = sum(double(a==b));
result = length(a)-L;
end
function X = pinv(A)
%pseudoinverz - Moore-Penrosova inverze Grevillovym algoritmem
oneEl = ones(A(1).field,[1,1]);
[~,n] = size(A);
zeroEl = zeros(A(1).field,[1,1]);
c = A(:,1);
if c == zeroEl
X = c';
else
X = c'/(c'*c);
end
for j = 2:n
d = X*A(:,j);
c = A(:,j) - A(:,1:(j-1))*d;
if c == zeroEl
bT = (d'*X)/(oneEl+d'*d);
else
bT = c'/(c'*c);
end
X = [X - d*bT; bT];
end
end
function U = triu(A)
%triu - horn trojhelnkov st
U = zeros(A(1).field,size(A));
[m,n]=size(A);
for i=1:m
for j=i+1:n
U(i,j) = A(i,j);
end
end
end
function L = tril(A)
%triu - horn trojhelnkov st
L = zeros(A(1).field,size(A));
m = size(A,1);
for i=1:m
for j=1:i-1
L(i,j) = A(i,j);
end
end
end
function D = diag(A)
%triu - horn trojhelnkov st
D = zeros(A(1).field,size(A));
[m,~]=size(A);
for i=1:m
D(i,i) = A(i,i);
end
end
function [L,U] = lu(A)
%lu - provede LU-dekonpozici tvercov matice
zeroEl = zeros(A(1).field,[1,1]);
[m,n]=size(A);
if m~=n
error('lu:notSquare','matice nen tvercov');
end
for q = 1:m-1
for i = q+1:m
if A(q,q) == zeroEl
error('lu:nullDividing','na diagonle je nulov prvek');
else
A(i,q) = A(i,q)/A(q,q);
for j = q+1:m
A(i,j) = A(i,j) - A(i,q)*A(q,j);
end
end
end
end
L = eye(A(1).field,m,m);
L = L + tril(A);
U = zeros(A(1).field,[m,m]) + triu(A) + diag(A);
end
function [Particular, NullSpace] = Frobenius(A,B)
%Frobenius - vypoet soustavy linernch rovnic. "b" me bt i matice,
%potom vsledn partikulrn een bude tak matice
[m,n] = size(A);
if nargin == 1
B = zeros(A(1).field,[m,1]);
end
if m ~= size(B,1)
error('Frobenius:badInput','A a B nem stejn poet dku');
end
hodnost = rank(A);
volnost = n-hodnost;
if hodnost ~= rank([A,B])
error('Frobenius:overDetermined','soustava nem een. jeden nebo vce sloupc z B nevyhovj');
end
Aplus = pinv(A);
Particular = Aplus*B;
NullSpace = (eye(A(1).field,n,n)-Aplus*A);
NullSpace = NullSpace(:,1:volnost);
end
end
end