BUG: np.polyval returns float64 for float32 coefficients

[ev-br]

Describe the issue:

Under NEP 50, I'd expect that np.polyval with single precision coefficients and a scalar evaluation value preserves single precision, but it converts to float64 instead.

Reproduce the code example:

In [6]: a1 = np.array([1., 2., 3., 4., 5.], dtype=np.float32)
   ...: a2 = 5.0
   ...: np.polyval(a1, a2).dtype
Out[6]: dtype('float64')

In [7]: np.__version__
Out[7]: '2.0.0rc2'

Error message:

No response

Python and NumPy Versions:

2.0.0rc2
3.11.0 | packaged by conda-forge | (main, Jan 14 2023, 12:27:40) [GCC 11.3.0]

Runtime Environment:

No response

Context for the issue:

No response

[luxedo]

This function is quite simple. The error occurs at NX.asanyarray and NX.zeros_like that casts python's float to np.float64 instead of np.float16.

    p = NX.asarray(p)
    if isinstance(x, poly1d):
        y = 0
    else:
        x = NX.asanyarray(x)
        y = NX.zeros_like(x)
    for pv in p:
        y = y * x + pv
    return y

One possible solution would be to test if x is a numpy array and pass the dtype argument.

    p = NX.asarray(p)
    if isinstance(x, poly1d):
        y = 0
    elif isinstance(x, ndarray):  # Preserve dtypes for NEP 50
        x = NX.asanyarray(x, dtype=x.dtype)
        y = NX.zeros_like(x, dtype=x.dtype)
    else:
        x = NX.asanyarray(x)
        y = NX.zeros_like(x)
    for pv in p:
        y = y * x + pv
    return y

Or maybe asanyarray and zeros_like should return the type correctly. With some guidance I'm can to build a PR to address this.

[ngoldbaum]

Maybe something like:

x = NX.asanyarray(x, dtype=p.dtype)

Note how I'm getting the dtype from p, not x. You also don't need to explicitly specify the dtype in zeros_like, since zeros_like gets the dtype from the input array-like if none is explicitly specified.

[seberg]

The best thing is probably to try to not convert things, the _array_converter was designed to help with that, but otherwise also if type(y) not in (int, float, complex): y = np.asarray(y) may be good (with some related fixes likely).

using dtype=p.dtype is unclear to me, since x could be an array with higher precision.

[eendebakpt]

As seberg already noted, using

x = NX.asanyarray(x, dtype=p.dtype)

will give incorrect results for something like polyval([1, 2, 3], 2.5]) (the 2.5 is converted to 2).

This works on the cases I tested so far:

    p = NX.asarray(p)
    if isinstance(x, poly1d):
        y = 0
    elif isinstance(x, (float, int, complex)):
        x = x
        y = NX.zeros( (1,), dtype=p.dtype)
    else:
        x = NX.asanyarray(x)
        y = NX.zeros_like(x)
    for pv in p:
        y = y * x + pv
    return y

@luxedo Would you like to make a PR see whether this works?

[luxedo]

That's right! I was testing opposite condition.

I tested a bit and the return value had the correct dtype but returned an array, so I just expanded the condition with the new types (float, int, complex) to keep the return consistent. I'm sending the PR right now.

    p = NX.asarray(p)
    if isinstance(x, (poly1d, float, int, complex)):
        y = 0
    else:
        x = NX.asanyarray(x)
        y = NX.zeros_like(x)
    for pv in p:
        y = y * x + pv
    return y

>>> np.polyval(np.array([5.], dtype=np.float32), 5)  # Suggestion
array([5.], dtype=float32)  # Returns an array<float32>

>>> np.polyval(np.array([5.], dtype=np.float32), 5)  # PR
np.float32(5.0)  # Returns np.float32

[eendebakpt]

The same issues occurs with the np.polynomial.Polynomial class. A minimal example:

import numpy as np
from numpy.polynomial import polyutils as pu

p=np.polynomial.Polynomial(np.array([1, 2], dtype=np.float32))
type(p(2)) # np.float64

# the reason is already in pu.mapparms
off, scl = pu.mapparms(p.domain, p.window)
type(off) # float64

p.domain.dtype, p.window.dtype # int32, int32

@seberg As the expert on NEP 50, if possible this should also be fixed right?

[eendebakpt]

There might be some more cases. For example

import numpy as np
arr = np.polydiv(1, np.float32(1))
print(arr)
print(arr[0].dtype) # float64

[seberg]

Yeah, sounds right, although not sure we should backport them. Not sure I am understanding the first one, int32 sounds weird there (unless that is a 32bit build?).

[eendebakpt]

I checked a bit more:

import numpy as np
arr = np.polydiv(1, np.float32(1))
print(arr)
print(arr[0].dtype) # float64

This is easy to fix: inside polydiv replace

    u = atleast_1d(u) + 0.0
    v = atleast_1d(v) + 0.0

with

    u = atleast_1d(u + 0.0)
    v = atleast_1d(v + 0.0)

(which, as a nice bonus, is also faster)

Another case:

np.polyadd(1, np.array([2], dtype=np.float16)).dtype # float64

To fix that, we need to have some special cases to check whether one of the input arguments is a scalar. One could also argue that polyadd expects an array input (and a list[float] or tuple[int] is also fine), and inputting a scalar is not a supported use case.