Power of a binary matrix

8 ビュー (過去 30 日間)
ENRIQUE SANCHEZ CARDOSO
ENRIQUE SANCHEZ CARDOSO 2022 年 12 月 2 日
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 件のコメント
ENRIQUE SANCHEZ CARDOSO
ENRIQUE SANCHEZ CARDOSO 2022 年 12 月 2 日
@Jan I've just edited my question to show you what I am trying to do.

サインインしてコメントする。

採用された回答

John D'Errico
John D'Errico 2022 年 12 月 2 日
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')
ans = 1×6
1 2 4 6 8 11
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)
res1 = 11×11
0 1 0 1 0 0 0 0 0 1 1 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 1 0 0 0 0 0 1 1 1 0
Verification:
res2 = mod(sym(A)^n,b)
res2 = 
res1-res2
ans = 
Of course, MPowerMod works for any modulus and power. It should also be pretty fast.
timeit(@() MPowerMod(A,n,b))
ans = 7.2756e-05
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
  1 件のコメント
ENRIQUE SANCHEZ CARDOSO
ENRIQUE SANCHEZ CARDOSO 2022 年 12 月 2 日
Thank you!

サインインしてコメントする。

その他の回答 (3 件)

Jan
Jan 2022 年 12 月 2 日
編集済み: Jan 2022 年 12 月 2 日
[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
Elapsed time is 0.039653 seconds.
tic % Smart collection of squaring:
C = MPowerBool(A, 1705);
toc
Elapsed time is 0.007874 seconds.
isequal(B, C)
ans = logical
1
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);

Image Analyst
Image Analyst 2022 年 12 月 2 日
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)
b = 11×11 logical array
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
whos a
Name Size Bytes Class Attributes a 11x11 968 double
whos b
Name Size Bytes Class Attributes b 11x11 121 logical
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)
b = 11×11 logical array
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
whos a
Name Size Bytes Class Attributes a 11x11 968 double
whos b
Name Size Bytes Class Attributes b 11x11 121 logical
  4 件のコメント
ENRIQUE SANCHEZ CARDOSO
ENRIQUE SANCHEZ CARDOSO 2022 年 12 月 2 日
@Image Analyst I've just edited my question to let you know what Im trying to do. I want result=mod(a^1705, 2) like the for loop does. Ty.

サインインしてコメントする。


Jan
Jan 2022 年 12 月 2 日
編集済み: Jan 2022 年 12 月 2 日
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
Elapsed time is 0.011007 seconds.
tic % Smart collection of squaring:
C = MPower_mod2(A, 1705);
toc
Elapsed time is 0.003123 seconds.
isequal(B, C)
ans = logical
1
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
  1 件のコメント
ENRIQUE SANCHEZ CARDOSO
ENRIQUE SANCHEZ CARDOSO 2022 年 12 月 2 日
Thanks! @Jan

サインインしてコメントする。

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by