Reciprocal.impl computes np.float32(1.0) / x: float32 precision for float64 scalars under NEP 50 (Python linker) #2430

Summary

scalar.basic.Reciprocal.impl is np.float32(1.0) / x. Under NumPy 2's promotion rules (NEP 50) a NumPy float32 scalar divided by a Python float gives a float32, so the Python implementation of the reciprocal of a float64 scalar carries float32 precision. The C implementation and np.reciprocal (used for arrays through nfunc_spec) are exact, so this shows only where the scalar Python impl runs — the Python linker, and anywhere scalar constants are folded through impl.

Reproduce

pytensor 3.3.0, numpy 2.4.6, Python 3.12, macOS arm64.

import numpy as np, pytensor, pytensor.tensor as pt
a = pt.dscalar("a")
f = pytensor.function([a], [1 / a, pt.reciprocal(a), a ** -1.0, 2.0 / a], mode="FAST_RUN")  # with PYTENSOR_FLAGS=linker=py
for name, v in zip(["1 / a", "reciprocal(a)", "a ** -1.0", "2.0 / a"], f(0.657)):
    print(name, repr(float(v)), abs(float(v) - eval(name.replace("reciprocal(a)", "1 / a").replace("a", "0.657"))) / (1 / 0.657))

1 / a          1.522070050239563  2.3e-08
reciprocal(a)  1.522070050239563  2.3e-08
a ** -1.0      1.522070050239563  2.3e-08
2.0 / a        3.0441400304414    0.0

The first three are rewritten to Reciprocal (pytensor.dprint shows the single node) and come back 2.3e-8 off; 2.0 / a is a true_div and exact. Directly in NumPy:

>>> np.float32(1.0) / 0.657
np.float32(1.52207)
>>> np.float32(1.0) / np.float64(0.657)
np.float64(1.5220700152207)

Fix

return 1.0 / x (or np.reciprocal(x)) in Reciprocal.impl. np.float32(1.0) was presumably there to keep float32 inputs float32 under the old value-based casting; with NEP 50 it now downcasts Python floats instead.

Where it was noticed

pymc-marketing's BinomialAdstock (1 / alpha - 1) and the half-life parameterisation of its GeometricAdstock (0.5 ** (1.0 / halflife)), and any InverseGamma prior (beta / x), evaluated under the Python linker: the log density of a whole MMM is off by 6e-7 and its gradient by 4e-8 relative, against a float64 reference that NumPy and an independent implementation agree on.