MATLAB method for determining the highest common factor of a matrix

Example:
Data Types:
|
|
|
|
|
|
|
|
|
Output Arguments
collapse all
— Greatest common divisor
real, nonnegative integer values
Greatest common divisor, returned as an array of real nonnegative
integer values. Thus, effectively the Euclidean algorithm, applied to a vector of numbers…
d = min(A);
while
true
r = mod(A,d);
if
~any(r)
break
end
% find the smallest non-zero element of r
r(r == 0) = inf;
d = min(r);
end
When the loop terminates, which it must do as long as all of the elements of A are positive integers, then d is the GCD of the group.


Solution 1:

If loops are not your preference, would you be interested in using recursive functions instead?

iif = @(varargin) varargin{2 * find([varargin{1:2:end}], 1, 'first')}();
gcdrec=@(v,gcdr) iif(length(v)==1,v, ...
                     v(1)==1,1, ...
                     length(v)==2,@()gcd(v(1),v(2)), ...
                     true,@()gcdr([gcd(v(1),v(2)),v(3:end)],gcdr));
mygcd=@(v)gcdrec(v(:)',gcdrec);
A=[0,0,0; 2,4,2;-2,0,8];
divisor=mygcd(A);
A=A/divisor;

The function defined by

iif

enables the creation of a recursive function,

gcdrec

, which can be used to compute the greatest common divisor of an array. To accomplish this,

iif

implements an inline conditional construct that tests the first argument to determine if it is

true

. If it is, the function returns the second argument. If not, it progresses to the third argument and repeats the process until it finds an argument that satisfies

true

. It is important to use

@()

to safeguard recursive functions and other quantities used within them to prevent errors.

By implementing

iif

, the recursive function

gcdrec

operates in the following manner:

  • It returns the input vector as it is, if it is a scalar.
  • In case the initial element of the vector is 1, recovery becomes impossible and the function returns 1 quickly, especially for bigger matrices.
  • If the length of the input vector is two, the function will calculate the greatest common divisor using

    gcd

    .
  • Otherwise, it invokes itself using a truncated vector where the initial two elements are replaced by their highest common factor.

The purpose of

mygcd

is to serve as a user-friendly interface.

The speed should be satisfactory, with the only potential issue being the recursion depth for exceptionally large problems. To compare with @Adriaan’s looping version, I conducted a brief timing evaluation using

A=randi(100,N,N)-50

,

N=100

,

N=1000

,

N=5000

,

tic

, and

toc

.


  1. N=100

    includes the following:
  2. looping 0.008 seconds
  3. recursive 0.002 seconds

  4. N=1000

    includes the following items:
  5. looping 0.46 seconds
  6. recursive 0.04 seconds

  7. N=5000

    includes the following items:
  8. looping 11.8 seconds
  9. recursive 0.6 seconds

The update reveals that the recursion limit, which is set at 500 by default, was not exceeded due to the absence of a common divisor in the data. However, doubling a random matrix and applying

N=100

would cause the recursion limit to be reached quickly, making it unsuitable for large matrices. On the other hand, @Adriaan’s solution works well for small matrices.

In an attempt to address the recursion limit issue, I experimented with rewriting the input vector to a smaller size at each recursive step. However, this approach, although effective, resulted in a significant decrease in speed, taking 2 seconds for

N=100

and a whopping 261 seconds for

N=1000

. I am still searching for a solution that strikes a balance between a reasonably large matrix size and an acceptable runtime.


Solution 2:

 A = [0,0,0; 2,4,2;-2,0,8];
 B = 1;
 kk = max(abs(A(:))); % start at the end
 while B~=0 && kk>=0
     tmp = mod(A,kk);
     B = sum(tmp(:));
     kk = kk - 1;
 end
 kk = kk+1;


To compute the sum of all elements in your matrix after applying the modulus operation,

B

initializes a counter. The

kk

is used to run through integers while

mod(A,kk)

computes the modulus after division for each element in

A

. If all elements in the matrix are divisible by 2,

mod(A,kk)

will return 0 for each element. The resulting modulo-matrix is transformed into a single column by

sum(tmp(:))

and summed to obtain a number. If this number is 0, it indicates that all elements in

A

are divisible by

kk

. The loop stops as soon as this happens, and the number in

kk

is the common divisor. Since

kk

is decreased every count, it is one value too low, thus one is added. Even though this might not be the quickest approach, it will suffice for now.

The loop has been modified to run in reverse as you are searching for the largest CD and not the smallest one. In case you have a matrix similar to

[4,8;16,8]

, the loop will halt at

2

instead of

4

. I apologize for the mistake, but this solution is now functional, despite the fact that the other two solutions presented here are much quicker.

One way to perform matrix division is shown below.

divided_matrix = A/kk;


Solution 3:


I concur, loops are not favorable to me either. We should eliminate them.

unqA = unique(abs(A(A~=0))).';             %//'
col_extent = [2:max(unqA)]';               %//'
B = repmat(col_extent,1,numel(unqA));
B(bsxfun(@gt,col_extent,unqA)) = 0;
divisor = find(all(bsxfun(@times,bsxfun(@rem,unqA,B)==0,B),2),1,'first');
if isempty(divisor)
    out = A;
else
    out = A/divisor;
end


Sample runs

Case #1:

A =
     0     0     0
     2     4     2
    -2     0     8
divisor =
     2
out =
     0     0     0
     1     2     1
    -1     0     4

Case #2:

A =
     0     3     0
     5     7     6
    -5     0    21
divisor =
     1
out =
     0     3     0
     5     7     6
    -5     0    21

Frequently Asked Questions