Isn't this a problem with all ufuncs?

Extending your example a little, and testing for multiplication a*b:

val = 1
cases = [(' p', val)]
for dtype,name in zip([np.int64, np.uint64], ['i', 'u']):
    cases.extend([('s' + name, dtype(val)),
                  ('0' + name, np.array(val, dtype=dtype)),
                  ('1' + name, np.array([val], dtype=dtype)),
                  ('2' + name, np.array([[val]], dtype=dtype))])
names, tps = zip(*cases)
ret = [[(a*b) for n,b in cases] for n1,a in cases]
kinds = [['p' if type(k) == int else k.dtype.kind for k in r] for r in ret]
print "     " + " ".join(names)
print "\n".join(n + "    " + "  ".join(d) for n,d in zip(names, kinds))

Output: (p = python, i = int, u = uint, s = scalar, 0,1,2 = 0d,1d,2d)

      p si 0i 1i 2i su 0u 1u 2u
 p    p  i  i  i  i  f  f  u  u
si    i  i  i  i  i  f  f  u  u
0i    i  i  i  i  i  f  f  u  u
1i    i  i  i  i  i  i  i  f  f
2i    i  i  i  i  i  i  i  f  f
su    f  f  f  i  i  u  u  u  u
0u    f  f  f  i  i  u  u  u  u
1u    u  u  u  f  f  u  u  u  u
2u    u  u  u  f  f  u  u  u  u

the power operation gives the same table.

Also, this is loosely related to #8231 and the related discussions on the mailing list last year about integer powers:

First in June under subject "Integers to integer powers, let's make a decision ", then in October under the same thread plus others like "fpower".

(but to be clear, this is a somewhat different issue)

The shape inconsistencies is #10322 again

This is the result of NumPy's "type demotion" rules for mixed scalar/array operations. When a positive-dimensional array and a 0-dimensional array or array scalar participate in a binary operation, the scalar gets "type demoted" based on its value, so 3 gets demoted to uint8, and the result dtype is uint64. On the other hand, when both operands are scalar, NumPy treats 3 as having dtype np.int_, which is signed, and numpy can't find an integer dtype big enough to hold all values of uint64 and int_, so it uses float64.

This is documented in the docs for numpy.result_type. (The implementation is in convert_datatype.c.)

The Stack Overflow question also showed what looked like NumPy using a floating-point power implementation for the uint64 output case. I don't know why that was happening, but it seems to have been fixed at some point after the question was asked.

Closing as duplicate of gh-22624