Power of a binary matrix
Show older comments
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 Comments
Jan
on 2 Dec 2022
If would be useful, if you post the code with the loop to clarify, what you want to achieve.
Is modulo 2 really wanted? Usually the multiplication of binary matrices is defined as: C = double((A * B) == 1), not C = mod(A * B, 2). Please explain this explicitly.
Jan
on 2 Dec 2022
@ENRIQUE SANCHEZ CARDOSO: Your code calculates a^1706, not 1705.
ENRIQUE SANCHEZ CARDOSO
on 2 Dec 2022
Accepted Answer
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
1 Comment
ENRIQUE SANCHEZ CARDOSO
on 2 Dec 2022
More Answers (3)
[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);
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 Comments
ENRIQUE SANCHEZ CARDOSO
on 2 Dec 2022
Edited: Jan
on 2 Dec 2022
Torsten
on 2 Dec 2022
for indice3=1:1704
instead of
for indice3=1:1705
Otherwise you have to start with a^0 = eye(11)
Image Analyst
on 2 Dec 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?
ENRIQUE SANCHEZ CARDOSO
on 2 Dec 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
1 Comment
ENRIQUE SANCHEZ CARDOSO
on 2 Dec 2022
Categories
Find more on Number Theory in Help Center and File Exchange
Products
on 2 Dec 2022
on 2 Dec 2022
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
