Power of a binary matrix

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!

[Jan]

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]

@ENRIQUE SANCHEZ CARDOSO: Your code calculates a^1706, not 1705.

[ENRIQUE SANCHEZ CARDOSO]

@Jan I've just edited my question to show you what I am trying to do.

Answer 1

[John D'Errico]

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

[ENRIQUE SANCHEZ CARDOSO]

Thank you!

Answer 2

[Jan]

[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);

Answer 3

[Image Analyst]

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              

[ENRIQUE SANCHEZ CARDOSO]

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 indice3=1:1705
    result = mod(result*a, 2);
end

[Torsten]

for indice3=1:1704

instead of

for indice3=1:1705

Otherwise you have to start with a^0 = eye(11)

[Image Analyst]

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]

@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.

Answer 4

[Jan]

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

[ENRIQUE SANCHEZ CARDOSO]

Thanks! @Jan