Essentially, math.gcd, but operating on numpy arrays with dtype= something in the int family.

>>> a = np.arange(4)
>>> b = np.arange(10,14)
>>> np.gcd(a, b)
array([10, 1, 2, 1])

Would support broadcasting:

>>> c = 2
>>> np.gcd(b, c)
array([2, 1, 2, 1])

Might want to have a version that operates on multiple values (more than 2) on an axis of a multivariable array:

>>> d = np.arange(9).reshape((3, 3))
>>> np.gcd_along_axis(d, axis=1)
array([3, 1, 1])

(This is slightly related to #4829.)