Difference between numpy power and ** for certain values

Question

I have a numpy array where the entries in f**2 differ from f[i]**2, but only for some specific value.

import numpy as np
np.set_printoptions(precision = 16)

f = np.array([ -40709.6555510835, -40708.6555510835, -33467.081758611654, -27653.379955714125])

f2 = f**2
# f2 = np.power(f,2)

print("outside loop", np.abs(f[1]**2 - f2[1]), np.abs(1.0 - f2[1] / f[1]**2), f.dtype, f[1].dtype, f2.dtype, f2[1].dtype)

for i, val in enumerate(f):        
    print("inside loop", i, np.abs(val**2 - f2[i]), np.abs(1.0 - f2[i] / val**2), val.dtype, f2.dtype, f2[i].dtype)

Produces output:

outside loop 2.384185791015625e-07 2.220446049250313e-16 float64 float64 float64 float64
inside loop 0 0.0 0.0 float64 float64 float64
inside loop 1 2.384185791015625e-07 2.220446049250313e-16 float64 float64 float64
inside loop 2 0.0 0.0 float64 float64 float64
inside loop 3 0.0 0.0 float64 float64 float64

I do note that this is a relative error on the order of epsilon.

This issue goes away when using np.power instead of ** in the definition of f2. Even so, why is f[i]**2 not the same as the ith value of f**2 (even if only for certain values in f).

I'm using python 3.10.6 and the latest numpy 1.26.4.

Edit:

The fundamental issue is captured in:

import numpy as np
f = np.array([-40708.6555510835])
print((f[0])**2 - (f**2)[0])

which displays a value of

-2.384185791015625e-07

I would like to know why that specific number has this specific issue. If you'd like confirmation, or to try different values for f, see this demo.

Answer 1

The results are different because f**2 calls numpy.square, while f[0]**2 and numpy.power(f, 2) call numpy.power.

---

numpy.ndarray.__pow__ is written in C. It looks like this:

static PyObject *
array_power(PyObject *a1, PyObject *o2, PyObject *modulo)
{
    PyObject *value = NULL;

    if (modulo != Py_None) {
        /* modular exponentiation is not implemented (gh-8804) */
        Py_INCREF(Py_NotImplemented);
        return Py_NotImplemented;
    }

    BINOP_GIVE_UP_IF_NEEDED(a1, o2, nb_power, array_power);
    if (fast_scalar_power(a1, o2, 0, &value) != 0) {
        value = PyArray_GenericBinaryFunction(a1, o2, n_ops.power);
    }
    return value;
}

The value = PyArray_GenericBinaryFunction(a1, o2, n_ops.power); is a Python function call to the numpy.power ufunc object, but first, it tries a fast_scalar_power function. This function tries to optimize exponentiation with common scalar powers, such as 2.

For the f**2 operation, fast_scalar_power detects the exponent of 2, and delegates the operation to numpy.square:

else if (exponent ==  2.0) {
    fastop = n_ops.square;
}

For numpy.power(f, 2), this is of course a direct call to numpy.power. numpy.power doesn't go through fast_scalar_power, and doesn't have any special handling for an exponent of 2. (Depending on what underlying power implementation it hits, that implementation might still have special handling for 2, though.)

For scalars, I believe numpy.float64.__pow__ actually just calls array_power:

static PyObject *
gentype_power(PyObject *m1, PyObject *m2, PyObject *modulo)
{
    if (modulo != Py_None) {
        /* modular exponentiation is not implemented (gh-8804) */
        Py_INCREF(Py_NotImplemented);
        return Py_NotImplemented;
    }

    BINOP_GIVE_UP_IF_NEEDED(m1, m2, nb_power, gentype_power);
    return PyArray_Type.tp_as_number->nb_power(m1, m2, Py_None);
}

so it hits fast_scalar_power, but one of the first checks in fast_scalar_power is

if (PyArray_Check(o1) &&

An instance of numpy.float64 does not pass PyArray_Check, which checks for objects whose type is exactly numpy.ndarray. Thus, the scalar goes through the general numpy.power code path.