Search code examples
matlaboctavegreatest-common-divisor

How to compute the GCD of a vector in GNU Octave / Matlab


gcd (A1, A2, ...) computes the GCD of elements A1(1), A2(1), .... Being the elements stored in a vector A, how to compute gcd (A)?
(I mean, gcd (4, 2, 8) = 2, gcd ([4, 2, 8] will raise an error in GNU Octave 4.0.0).


Solution

  • With cell array expansion

    Here is a one-liner, valid only in octave (thanks to nirvana-msu for pointing out matlab's limitation):

    A = [10 25 15];
    gcd(num2cell(A){:})
    # ans =  5
    

    This use cell array expansion, which is a bit hidden there :

    Accessing multiple elements of a cell array with the ‘{’ and ‘}’ operators will result in a comma separated list of all the requested elements

    so here A{:} is interpreted as A(1), A(2), A(3), and thus gcd(A{:}) as gcd(A(1), A(2), A(3))


    Performance

    Still under octave

    A = 3:259;
    tic; gcd(num2cell(A){:}); toc
    
    Elapsed time is 0.000228882 seconds.
    

    while with the gcd_vect in @nirvana_msu answer,

    tic; gcd_vect(A); toc
    
    Elapsed time is 0.0184669 seconds.
    

    This is because using recursion implies a high performance penalty (at least under octave). And actually for more than 256 elements in A, recursion limit is exhausted.

    tic; gcd_vect(1:257); toc
    
    <... snipped bunch of errors as ...>
    error: evaluating argument list element number 2
    error: called from
    gcd_vect at line 8 column 13
    

    This could be improved a lot by using a Divide and conquer algorithm

    While the cell array expansion (octave only) scales well:

    A = 127:100000;
    tic; gcd(num2cell(A){:}); toc
    Elapsed time is 0.0537438 seconds.
    

    Divide and conquer algorithm (best)

    This one should work under matlab too (not tested though. Feedback welcome).

    It uses recursion too, like in other answers, but with Divide and conquer

    function g = gcd_array(A)
      N = numel(A);
    
      if (mod(N, 2) == 0)
        % even number of elements
        % separate in two parts of equal length
        idx_cut = N / 2;
        part1 = A(1:idx_cut);
        part2 = A(idx_cut+1:end);
        % use standard gcd to compute gcd of pairs
        g = gcd(part1(:), part2(:));
        if ~ isscalar(g)
           % the result was an array, compute its gcd
           g = gcd_array(g);
        endif
      else
        % odd number of elements
        % separate in one scalar and an array with even number of elements
        g = gcd(A(1), gcd_array(A(2:end)));
      endif
    endfunction
    

    timings:

    A = 127:100000;
    tic; gcd_array(A); toc
    Elapsed time is 0.0184278 seconds.
    

    So this seems even better than cell array expansion.