Power of a binary matrix
6 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
ENRIQUE SANCHEZ CARDOSO
am 2 Dez. 2022
Kommentiert: ENRIQUE SANCHEZ CARDOSO
am 2 Dez. 2022
Hello, I'm modeling a PRBS (Pseudo Random Bit Secuence) used in electronics.
This is defined in a for loop:
PRBS=zeros(1705,1);
pilotos=ones(11,1);
for indice=1:1705
PRBS(indice)=pilotos(11);
aux=xor(pilotos(9), pilotos(11));
pilotos(2:11)=pilotos(1:10);
pilotos(1)=aux;
end
Now I want to eliminate the for loop, so I decided to build a matrix that gives me the i-th element of PRBS vector. This matrix is shifting rows of pilotos and doing xor (matrix a does the sum of the 9th and 11th bit and later I do module 2 to get the same result of xor) of the 9th and 11th element.
a=[0 1 0 0 0 0 0 0 0 0 0;
0 0 1 0 0 0 0 0 0 0 0;
0 0 0 1 0 0 0 0 0 0 0;
0 0 0 0 1 0 0 0 0 0 0;
0 0 0 0 0 1 0 0 0 0 0;
0 0 0 0 0 0 1 0 0 0 0;
0 0 0 0 0 0 0 1 0 0 0;
0 0 0 0 0 0 0 0 1 0 0;
1 0 0 0 0 0 0 0 0 1 0;
0 0 0 0 0 0 0 0 0 0 1;
1 0 0 0 0 0 0 0 0 0 0];
PRBS(1)=pilotos(11);
aux2=pilotos*a;
PRBS(2)=aux2(11);
a2=mod(a^2, 2)
aux3=pilotos*a2;
PRBS(3)=aux3(11);
I want result=mod(a^1705, 2) to be a binary matrix, so I tried a1705=mod(a^1705, 2) but it seems that it first does a^1705 (very big numbers) and later mod 2, so I dont get the result I want beacause I lost information.
I know I can do it with a for loop, but I dont want it, I need a more efficient way.
My expected output is:
result = [1 0 1 0 1 0 0 0 0 0 1;
1 1 0 1 0 1 0 0 0 0 0;
0 1 1 0 1 0 1 0 0 0 0;
0 0 1 1 0 1 0 1 0 0 0;
0 0 0 1 1 0 1 0 1 0 0;
1 0 0 0 1 1 0 1 0 1 0;
0 1 0 0 0 1 1 0 1 0 1;
0 0 1 0 0 0 1 1 0 1 0;
0 0 0 1 0 0 0 1 1 0 1;
1 0 1 0 0 0 0 0 1 1 1;
0 1 0 1 0 0 0 0 0 1 1];
%I cant get it using:
result = a;
for k=1:1704
result = mod(result*a, 2);
end
Thank you!
3 Kommentare
Akzeptierte Antwort
John D'Errico
am 2 Dez. 2022
This is really pretty simple. Far simpler than you might think.
You want to do a powermod operation, on a matrix, so
mod(A^n,b)
where b is the modulus. Here, we have b==2 and n = 1705, but b and n could be any values.
The trick is to decompose the power n into a binary form. This leaves us with an efficient scheme to write the powermod operation. That is,
find(dec2bin(1705) == '1')
So effectively, we know that
1705 = sum(2.^[1 2 4 6 8 11])
In turn, that tells us that we can form your matrix power into a product of only 6 matrices. We can see how this works now in the code I'll write as MPowerMod.
Now test the code.
A = [0 1 0 0 0 0 0 0 0 0 0;
0 0 1 0 0 0 0 0 0 0 0;
0 0 0 1 0 0 0 0 0 0 0;
0 0 0 0 1 0 0 0 0 0 0;
0 0 0 0 0 1 0 0 0 0 0;
0 0 0 0 0 0 1 0 0 0 0;
0 0 0 0 0 0 0 1 0 0 0;
0 0 0 0 0 0 0 0 1 0 0;
1 0 0 0 0 0 0 0 0 1 0;
0 0 0 0 0 0 0 0 0 0 1;
1 0 0 0 0 0 0 0 0 0 0];
n = 1705;
b = 2;
res1 = MPowerMod(A,n,b)
Verification:
res2 = mod(sym(A)^n,b)
res1-res2
Of course, MPowerMod works for any modulus and power. It should also be pretty fast.
timeit(@() MPowerMod(A,n,b))
As yuo can see, it is efficient since only very few matrix multiplies are done.
function res = MPowerMod(A,n,b)
% compute mod(A^n,b), using a powermod scheme to avoid overflows
%
% arguments: (INPUT)
% A - square (INTEGER) matrix
% n - non-negative integer, no larger than 2^53-1
% b - positive integer >= 2
%
[r,c] = size(A);
if r ~= c
error('A must be a square matrix')
end
res = eye(r);
if n == 0
return
end
A = mod(A,b);
% binary expansion of n
nbin = flip(dec2bin(n)) - '0';
for i = 1:numel(nbin)
% does this power appear in the binary expansion of n?
if nbin(i)
res = mod(res*A,b);
end
% square A only if we are not done
if i < numel(nbin)
A = mod(A*A,b);
end
end
end
Weitere Antworten (3)
Jan
am 2 Dez. 2022
Bearbeitet: Jan
am 2 Dez. 2022
[EDITED] This is not the multiplication wanted by the OP - I leave it for educational reasons.
Assuming, that you mean the usual definition of multiplying binary matrices:
A = randi([0,1], 11, 11); % Test data, use your "a"
tic % Simple loop approach:
B = A;
for k = 1:1704
B = MMBool(B, A);
end
toc
tic % Smart collection of squaring:
C = MPowerBool(A, 1705);
toc
isequal(B, C)
function C = MMBool(A, B)
% C = A * B for binary matrices
% Input:
% A, B: Matrices with matching sizes, logial or numerical, but with
% values 0 and 1 only.
% Output:
% C = A * B
%
% (C) 2022 Jan, Heidelberg, CC BY-SA 3.0
mA = height(A);
nA = width(A);
nB = width(B);
C = false(mA, nB);
for i = 1:mA
for j = 1:nB
for k = 1:nA
if A(i,k) && B(k,j)
C(i, j) = true;
break;
end
end
end
end
end
function Y = MPowerBool(X, p)
% Y = X^p for binary matrices
% Input:
% X: Square binary matrix, logical or numerical, but with values 0 or 1 only
% p: Integer > 0
% Output:
% Y: X^p, same class as input.
%
% (C) 2022 Jan, Heidelberg, CC BY-SA 3.0
% See: James Tursa's https://www.mathworks.com/matlabcentral/fileexchange/25782
nBit = floor(log2(p)) + 1;
b = rem(floor(p ./ pow2(0:nBit - 1)), 2);
Y = [];
for k = 1:nBit
if b(k)
if isempty(Y)
Y = X;
else
Y = MMBool(Y, X);
end
end
if k ~= nBit
X = MMBool(X, X);
end
end
end
This can be accelerated, if the inputs are numerical matrices by using this instead of MMBool:
C = double((A * B) == 1);
0 Kommentare
Image Analyst
am 2 Dez. 2022
What is your expected output? Do you want element by element multiplication or matrix multiplication? Here is element by element multiplication
a=[0 1 0 0 0 0 0 0 0 0 0;
0 0 1 0 0 0 0 0 0 0 0;
0 0 0 1 0 0 0 0 0 0 0;
0 0 0 0 1 0 0 0 0 0 0;
0 0 0 0 0 1 0 0 0 0 0;
0 0 0 0 0 0 1 0 0 0 0;
0 0 0 0 0 0 0 1 0 0 0;
0 0 0 0 0 0 0 0 1 0 0;
1 0 0 0 0 0 0 0 0 1 0;
0 0 0 0 0 0 0 0 0 0 1;
1 0 0 0 0 0 0 0 0 0 0];
b = logical(a .^ 1704)
whos a
whos b
Here is matrix multiplication:
a=[0 1 0 0 0 0 0 0 0 0 0;
0 0 1 0 0 0 0 0 0 0 0;
0 0 0 1 0 0 0 0 0 0 0;
0 0 0 0 1 0 0 0 0 0 0;
0 0 0 0 0 1 0 0 0 0 0;
0 0 0 0 0 0 1 0 0 0 0;
0 0 0 0 0 0 0 1 0 0 0;
0 0 0 0 0 0 0 0 1 0 0;
1 0 0 0 0 0 0 0 0 1 0;
0 0 0 0 0 0 0 0 0 0 1;
1 0 0 0 0 0 0 0 0 0 0];
b = logical(a ^ 1704)
whos a
whos b
4 Kommentare
Image Analyst
am 2 Dez. 2022
But where is the exponentiation taking place? You said "I want result=a^1705", in other words a raised to the 1705'th power. I don't see that in your loop. How does your loop with mod accomplish exponentiation?
Jan
am 2 Dez. 2022
Bearbeitet: Jan
am 2 Dez. 2022
If you really define the matrix multiplication with mod(A*B, 2):
A = randi([0,1], 11, 11); % Test data, use your "a"
tic % Simple loop approach:
B = A;
for k = 1:1704 % 1705 - 1 !
B = mod(B * A, 2);
end
toc
tic % Smart collection of squaring:
C = MPower_mod2(A, 1705);
toc
isequal(B, C)
function Y = MPower_mod2(X, p)
% Y = X^p with mod(M, 2) in each multiplication
% Input:
% X: Square binary matrix, numerical with values 0 or 1 only
% p: Integer > 0
% Output:
% Y: X^p, same class as input.
%
% (C) 2022 Jan, Heidelberg, CC BY-SA 3.0
% See: James Tursa's https://www.mathworks.com/matlabcentral/fileexchange/25782
nBit = floor(log2(p)) + 1;
b = rem(floor(p ./ pow2(0:nBit - 1)), 2); % p as binary vector
Y = [];
for k = 1:nBit
if b(k)
if isempty(Y)
Y = X;
else
Y = mod(Y * X, 2);
end
end
if k ~= nBit
X = mod(X * X, 2);
end
end
end
Siehe auch
Kategorien
Mehr zu Number Theory finden Sie in Help Center und File Exchange
Produkte
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!