Numpy gcd function

Question

Does numpy have a gcd function somewhere in its structure of modules?

I'm aware of fractions.gcd but thought a numpy equivalent maybe potentially quicker and work better with numpy datatypes.

I have been unable to uncover anything on google other than this link which seems out of date and I don't know how I would access the _gcd function it suggests exists.

Naively trying:

np.gcd
np.euclid

hasn't worked for me...

I think the `_gcd` function you're talking about is at `numpy.core._internal._gcd`, but it's in pure Python (and so not too speedy) and doesn't handle numpy arrays in any case. — DSM, Mar 22 '13 at 11:43
is the de facto approach to use fractions.gcd? I assume that will be a faster C implementation — bph, Mar 22 '13 at 11:45

score 12 · Accepted Answer · answered Mar 22 '13 at 12:46

12

You can write it yourself:

def numpy_gcd(a, b):
    a, b = np.broadcast_arrays(a, b)
    a = a.copy()
    b = b.copy()
    pos = np.nonzero(b)[0]
    while len(pos) > 0:
        b2 = b[pos]
        a[pos], b[pos] = b2, a[pos] % b2
        pos = pos[b[pos]!=0]
    return a

Here is the code to test the result and speed:

In [181]:
n = 2000
a = np.random.randint(100, 1000, n)
b = np.random.randint(1, 100, n)
al = a.tolist()
bl = b.tolist()
cl = zip(al, bl)
from fractions import gcd
g1 = numpy_gcd(a, b)
g2 = [gcd(x, y) for x, y in cl]
print np.all(g1 == g2)

True

In [182]:
%timeit numpy_gcd(a, b)

1000 loops, best of 3: 721 us per loop

In [183]:
%timeit [gcd(x, y) for x, y in cl]

1000 loops, best of 3: 1.64 ms per loop

answered Mar 22 '13 at 12:46

HYRY

94,853
25
187
187

two and a bit times quicker for your numpy version – bph Mar 22 '13 at 13:07
the function you wrote doesn't work for numpy_gcd(10, 5) or am I doing something wrong? I get 0-d arrays can't be indexed – evan54 Apr 15 '14 at 21:14
Its because the function expects numpy arrays (and probably numpy scalars). You should be able to do numpy_gcd[10], [5]) or modify the function to work directly with Python scalars. – coderforlife Apr 22 '15 at 05:07

score 12 · Answer 2 · answered Oct 25 '16 at 23:31

12

Public service announcement for anyone using Python 3.5

from math import gcd
gcd(2, 4)

And if you want to write it yourself in a one-liner:

def gcd(a: int, b: int): return gcd(b, a % b) if b else a

answered Oct 25 '16 at 23:31

lababidi

2,654
1
22
14

score 9 · Answer 3 · answered Mar 22 '13 at 11:50

9

It seems there is no gcd function yet in numpy. However, there is a gcd function in fractions module. If you need to perform gcd on numpy arrays, you could build a ufunc using it:

gcd = numpy.frompyfunc(fractions.gcd, 2, 1)

answered Mar 22 '13 at 11:50

Charles Brunet

21,797
24
83
124

8

I'm surprised numpy doesn't contain a gcd function – bph Mar 22 '13 at 11:59
Nice. To take the gcd of an entire list, you can use `numpy.ufunc.reduce(gcd, m)` where `m` is the list and `gcd` is as defined in this answer. – Ian Hincks Jun 06 '17 at 13:27

score 6 · Answer 4 · edited Oct 26 '21 at 11:27

6

The functions gcd (Greatest Common Divisor) and lcm (Lowest Common Multiple) have been added to numpy in version 1.15.

You can use them both "as is" on a pair of scalars

import numpy as np

np.gcd(-5, 10)  # yields '5'

or on a list or array using .reduce:

np.gcd.reduce(np.array([-5, 10, 0, 5]))  # yields '5'

edited Oct 26 '21 at 11:27

smarie

4,568
24
39

answered Dec 06 '18 at 11:50

johannes_r

61
1
1

score 1 · Answer 5 · answered Nov 02 '17 at 20:12

In case the desired result is not an element-wise gcd but rather the gcd of all numbers in the array, you may use the code below.

import numpy as np
from math import gcd as mathgcd

def numpy_set_gcd(a):
    a = np.unique(a)
    if not a.dtype == np.int or a[0] <= 0:
        raise ValueError("Argument must be an array of positive " +
                         "integers.")

    gcd = a[0]
    for i in a[1:]:
        gcd = mathgcd(i, gcd)
        if gcd == 1:
            return 1 

    return gcd

Depending on the use case, it can be faster to omit the sorting step a = np.unique(a).

An alternative (maybe more elegant but slower) implementation using ufuncs is

import numpy as np
from math import gcd as mathgcd
npmathgcd = np.frompyfunc(mathgcd, 2, 1)

def numpy_set_gcd2(a):
    a = np.unique(a)
    if not a.dtype == np.int or a[0] <= 0:
        raise ValueError("Argument must be an array of positive " +
                         "integers.")
    npmathgcd.at(a[1:], np.arange(a.size-1), a[:-1])
    return a[-1]

Numpy gcd function

5 Answers5

Linked