Repository: numpy/numpy
Pull request: API,BUG: Fix copyto (and ufunc) handling of scalar cast safety
This defines scalar cast safety to be "safe", "equiv", or "unsafe" (for Python int, float, complex).
Although, equiv is a bit of an odd one.  We define it by using the common DType with the desired one, and assuming that this cast is at least safe.  With that, we can create an appropriate array for the scalar (and the rest works fine).

This means that any python float can be safely cast to e.g. float32 and any Python integer to any NumPy integer even though the actual assignment may overflow (floats) or fail (integers). `can_cast` itself is not defined here.  It could be, but has the disadvantage that it won't error out on assignment later for integers.

Closes gh-26381, gh-27075

---

~Still needs release notes~, but ready for review.  Sorry it ended up as one big commit, but besides moving some helper out of ufunc.c looking at it per-file is just as well.

Problem statement:
BUG: NEP 50 inconsistency in `copyto`
As @ev-br noticed:  Copyto is currently inconsistent in NEP 50, which will be noticed relatively rarely, because the default is `same-kind` casting, where it just doesn't matter.

It turns out, `copyto` is the only (meaningful) user of `PyArray_AssignRawScalar` when the input is 0-D, and that function is the only *direct* user of the old scalar logic.

Now, `copyto` runs into the same difficulty as ufuncs (even more niche), about how you define "safe" casting for Python int/float/complex.  In principle, we could:
1. Just designate that the `src` is an array, and ignore any subtleties (right now it uses old logic, which also kicks in for 0-d arrays).
   * Means that `np.copyto(uint8_arr, 1000)` *wraps* the integer (which is what happens right now anyway).
2. Make a decision on how to define "safe" for Python scalars:
   * For integers, safe could be "fits the result" (based on the value).
   * For floats, safe is tricky, they always lose precision.

I don't care too much which way we go here, but I think we should punt this to NumPy 2.1 and then need to make a decision though.  We were bound to have such issues, that it is actually using old logic is annoying, but not much worse than similar ones we must still have in Python.

The advantage of defining things as 2. is that it can make `copyto` behave like a ufunc: `np.positive(src, out=dst, casting=...)`  (positive is just effectively an identity ufunc).
That part has the subtlety about mixing `casting=` with Python scalars (but not the issue about 0-D arrays).

---

I would be much obliged for input, there is a reason this is still there: Deciding what to do with scalars and cast-safety is not trivial and probably more of an "expectation" decision than trying to reason what is "right".

Hints from issue discussion:
> The advantage of defining things as 2. is that it can make `copyto` behave like a ufunc: `np.positive(src, out=dst, casting=...)`

Agreed with this - they should behave the same as much as possible, so I think we should go for 2.

But your question about `casting="safe"` for the python case seems relevant. In principle, I guess one could have it mean the same as `np.copyto(a, a.dtype.type(python-int))`. Though perhaps it is also not unreasonable to just document that option, since if we were to go that route, it would seem not inappropriate to do that everywhere, i.e., pick the minimum size dtype that actually is able to hold the value and then get the final loop from there (though typically the other dtypes may then lead to a preference for unsigned, if possible).
+1 for `np.copyto(dst, src, casting)` being `np.positive(src, out=dst, casting=casting).

If this means that short ints may raise, OK, this is self-consistent. An alternative is to basically set `casting='unsafe'` whenever a scalar is involved. Then all operations succeed but short/unsigned ints may wrap. This is also consistent. (and this is what pytorch does I believe).
But then it's indeed best to make the change across the board indeed. Currently,

```
In [16]: np.asarray(2001, dtype=np.uint8)
---------------------------------------------------------------------------
OverflowError                             Traceback (most recent call last)
Cell In[16], line 1
----> 1 np.asarray(2001, dtype=np.uint8)

OverflowError: Python integer 2001 out of bounds for uint8
```
Right, but two question:
* do you guys care about the upcoming release?  Because I still lean towards waiting because it is rather niche.  I.e. I doubt it matters to `cupy` (we either need to ignore it or need to work around one way or another).  But changing it might break already released downstream.
* I agree with `np.positive` being the right behavior, but unfortunately that doesn't lead to a definition for `casting=` behavior for Python scalars there.  (I have some tendencies for what to do, but they are not strong.)
  * This is an opinionated decisions, IMO.  I have not been able to reason clearly one way or another for Python scalars.  Beyond that the default behavior should be what it is now (raise an error).  (Additionally, we probably won't yet have the API to get what we land on.)
IMO, copyto(python scalar) can indeed wait until numpy 2.1 or something. It is annoying yes, but can be worked around.

Why does `np.positive` not lead to a definition of casting for scalars? This seems consistent:

```
In [9]: In [8]: for src in [101, 1001]:
   ...:    ...:   for casting in ['no', 'equiv', 'safe', 'same_kind', 'unsafe']:
   ...:    ...:     dst = np.ones(1, dtype=np.uint8)
   ...:    ...:     try:
   ...:    ...:         np.positive(src, out=dst, casting=casting)
   ...:    ...:         result = dst.copy()
   ...:    ...:     except Exception as e:
   ...:    ...:         result = e
   ...:    ...:     print(src, casting, result)
   ...: 
101 no Cannot cast ufunc 'positive' output from dtype('int64') to dtype('uint8') with casting rule 'no'
101 equiv Cannot cast ufunc 'positive' output from dtype('int64') to dtype('uint8') with casting rule 'equiv'
101 safe Cannot cast ufunc 'positive' output from dtype('int64') to dtype('uint8') with casting rule 'safe'
101 same_kind Cannot cast ufunc 'positive' output from dtype('int64') to dtype('uint8') with casting rule 'same_kind'
101 unsafe [101]
1001 no Cannot cast ufunc 'positive' output from dtype('int64') to dtype('uint8') with casting rule 'no'
1001 equiv Cannot cast ufunc 'positive' output from dtype('int64') to dtype('uint8') with casting rule 'equiv'
1001 safe Cannot cast ufunc 'positive' output from dtype('int64') to dtype('uint8') with casting rule 'safe'
1001 same_kind Cannot cast ufunc 'positive' output from dtype('int64') to dtype('uint8') with casting rule 'same_kind'
1001 unsafe [233]
```

---------------------------------------------

What looks like a bug to me is value-based logic in 0D arrays. As discussed offline: 

```
In [7]: In [8]: for src in [np.asarray(101), np.asarray(1001)]:
   ...:    ...:   for casting in ['no', 'equiv', 'safe', 'same_kind', 'unsafe']:
   ...:    ...:     dst = np.ones(1, dtype=np.uint8)
   ...:    ...:     try:
   ...:    ...:         np.copyto(dst, src, casting)
   ...:    ...:         result = dst.copy()
   ...:    ...:     except Exception as e:
   ...:    ...:         result = e
   ...:    ...:     print(src, casting, result)
   ...: 
101 no [101]
101 equiv [101]
101 safe [101]
101 same_kind [101]
101 unsafe [101]
1001 no Cannot cast scalar from dtype('int64') to dtype('uint8') according to the rule 'no'
1001 equiv Cannot cast scalar from dtype('int64') to dtype('uint8') according to the rule 'equiv'
1001 safe Cannot cast scalar from dtype('int64') to dtype('uint8') according to the rule 'safe'
1001 same_kind [233]
1001 unsafe [233]
```

OTOH, if it's more convenient to fix this together with python scalars, downstream can adapt.
Argg.  Partially, because `np.positive` is more strict there than we may want it to be.

We allow: `np.add(101, 101, dtype=np.uint8)` because not doing so would also not allow `np.add(101, np.uint8(0))` also note try using `int8` instead of `uint8` since `int8` and `int64` are same-kind.  Making that distinction for Python scalars doesn't seem right.

Now, the closest to just saying "all fine", to me seems like saying: OK, if casting is passed explicitly, we just use the default dtype (even for uint that isn't same-kind to Python int).  But if `casting` is _not_ passed the cast-safety is actually "same-kind or try convert py-scalar".
TBH and FWIW, what you write sounds too tricky. I've no idea how to predict what's going on or how to explain what's going on (the "how to teach this" PEP rubric). For `np.copyto`, if `np.positive(src, out=dst, casting=casting)` does not cut it, I think there are two reasonable options:
- convert the scalar to the default dtype --- equivalent to `np.asarray(scalar)`? --- and proceed from there with casting as specified, explicitly or via a default kwarg
- convert the src scalar to `dst.dtype` --- equivalent to `np.asarray(scalar, dtype=dst.dtype)`? --- and proceed from there with the casting as specified, explicitly or via a default kwarg 
Well, it's what we have:  I have struggled with this and brought up the fact that it is tricky many times and over a very long time...
* You probably want `101 + uint8(1)` to work and return a uint8
* And most people seem to want `np.add(101, uint8(1))` to be identical to the first.

So you have no choice but to do the promotion by default by basically ignoring the value and saying a Python integer could be any integer really.  And that works fine, except that there is a line when code calls `np.asarray(python_int)`, which is a problem because after that you cannot also behave like a ufunc anymore.

But back to the problem: This is mostly fine for *promotion*, but if you look at casting things get tricky because while promotion used to be implemented via `can_cast` it isn't the same thing and we don't have a definition for it for Python scalars that makes much sense.

Unless you are suggesting to not do *any* NEP 50 style weak promotion at all, which of course is the only truly consistent thing without pitfalls, but...
FWIW, my understanding is that in the last meeting it seemed like there was a preference for just considering these casts safe (in particular `pyfloat -> float32`).
@charris I made some progress on this today.  It is a bit round-about, but fixes the general issue of cast-safety in ufuncs first (this part works, I need to use it also for copyto; PR upcoming tomorrow latest).

It would be neat if we can manage to get it into the rc.  I'll have a PR ready tomorrow (probably not tonight).

To sketch the idea briefly (and yes, it seems very round-about maybe, but for our types pretty much every step is cut short):
* `np.copyto(uint8_arr, 3.0)`, needs to find a `dtype` to use for the 3.0
* We use `common = common_dtype(UInt8, PyFloatDType)` (which is float64 here, or Float32 for `float32_arr, 3.0`).  I.e. this step is pretty trivial in practice.
* We use `common` to create the temporary array for the 3.0.  We assumes that the common DType operation can only reasonably give us something that must be considered a `"safe"` cast.
* Work with it as normal.

A slightly awkward thing, is whether to bother to keep allowing `copyto(float64_arr, 3.0, casting="equiv")` (equivalent casting).  To do that, the easiest is to do `np.asarray()` first, and then see if the above gave the same dtype or not.
(That is fine, and easier right now anyway.  It is just a bit round-about in the longer term.)

Commit message: Apply suggestions from code review

Co-authored-by: Nathan Goldbaum <nathan.goldbaum@gmail.com>

Patch to review:
diff --git a/doc/release/upcoming_changes/27091.change.rst b/doc/release/upcoming_changes/27091.change.rst
new file mode 100644
index 000000000000..5b71692efabd
--- /dev/null
+++ b/doc/release/upcoming_changes/27091.change.rst
@@ -0,0 +1,24 @@
+Cast-safety fixes in ``copyto`` and ``full``
+--------------------------------------------
+``copyto`` now uses NEP 50 correctly and applies this to its cast safety.
+Python integer to NumPy integer casts and Python float to NumPy float casts
+are now considered "safe" even if assignment may fail or precision may be lost.
+This means the following examples change slightly:
+
+* ``np.copyto(int8_arr, 1000)`` previously performed an unsafe/same-kind cast
+   of the Python integer.  It will now always raise, to achieve an unsafe cast
+   you must pass an array or NumPy scalar.
+* ``np.copyto(uint8_arr, 1000, casting="safe")`` will raise an OverflowError
+  rather than a TypeError due to same-kind casting.
+* ``np.copyto(float32_arr, 1e300, casting="safe")`` will overflow to ``inf``
+  (float32 cannot hold ``1e300``) rather raising a TypeError.
+
+Further, only the dtype is used when assigning NumPy scalars (or 0-d arrays),
+meaning that the following behaves differently:
+
+* ``np.copyto(float32_arr, np.float64(3.0), casting="safe")`` raises.
+* ``np.coptyo(int8_arr, np.int64(100), casting="safe")`` raises.
+  Previously, NumPy checked whether the 100 fits the ``int8_arr``.
+
+This aligns ``copyto``, ``full``, and ``full_like`` with the correct NumPy 2
+behavior.
\ No newline at end of file
diff --git a/numpy/_core/src/multiarray/abstractdtypes.c b/numpy/_core/src/multiarray/abstractdtypes.c
index 8d00084f0efe..214833737792 100644
--- a/numpy/_core/src/multiarray/abstractdtypes.c
+++ b/numpy/_core/src/multiarray/abstractdtypes.c
@@ -378,3 +378,117 @@ NPY_NO_EXPORT PyArray_DTypeMeta PyArray_PyComplexDType = {{{
     .dt_slots = &pycomplexdtype_slots,
     .scalar_type = NULL,  /* set in initialize_and_map_pytypes_to_dtypes */
 };
+
+
+/*
+ * Additional functions to deal with Python literal int, float, complex
+ */
+/*
+ * This function takes an existing array operand and if the new descr does
+ * not match, replaces it with a new array that has the correct descriptor
+ * and holds exactly the scalar value.
+ */
+NPY_NO_EXPORT int
+npy_update_operand_for_scalar(
+    PyArrayObject **operand, PyObject *scalar, PyArray_Descr *descr,
+    NPY_CASTING casting)
+{
+    if (PyArray_EquivTypes(PyArray_DESCR(*operand), descr)) {
+        /*
+        * TODO: This is an unfortunate work-around for legacy type resolvers
+        *       (see `convert_ufunc_arguments` in `ufunc_object.c`), that
+        *       currently forces us to replace the array.
+        */
+        if (!(PyArray_FLAGS(*operand) & NPY_ARRAY_WAS_PYTHON_INT)) {
+            return 0;
+        }
+    }
+    else if (NPY_UNLIKELY(casting == NPY_EQUIV_CASTING) &&
+             descr->type_num != NPY_OBJECT) {
+        /*
+         * increadibly niche, but users could pass equiv casting and we
+         * actually need to cast.  Let object pass (technically correct) but
+         * in all other cases, we don't technically consider equivalent.
+         * NOTE(seberg): I don't think we should be beholden to this logic.
+         */
+        PyErr_Format(PyExc_TypeError,
+            "cannot cast Python %s to %S under the casting rule 'equiv'",
+            Py_TYPE(scalar)->tp_name, descr);
+        return -1;
+    }
+
+    Py_INCREF(descr);
+    PyArrayObject *new = (PyArrayObject *)PyArray_NewFromDescr(
+            &PyArray_Type, descr, 0, NULL, NULL, NULL, 0, NULL);
+    Py_SETREF(*operand, new);
+    if (*operand == NULL) {
+        return -1;
+    }
+    if (scalar == NULL) {
+        /* The ufunc.resolve_dtypes paths can go here.  Anything should go. */
+        return 0;
+    }
+    return PyArray_SETITEM(new, PyArray_BYTES(*operand), scalar);
+}
+
+
+/*
+ * When a user passed a Python literal (int, float, complex), special promotion
+ * rules mean that we don't know the exact descriptor that should be used.
+ *
+ * Typically, this just doesn't really matter.  Unfortunately, there are two
+ * exceptions:
+ * 1. The user might have passed `signature=` which may not be compatible.
+ *    In that case, we cannot really assume "safe" casting.
+ * 2. It is at least fathomable that a DType doesn't deal with this directly.
+ *    or that using the original int64/object is wrong in the type resolution.
+ *
+ * The solution is to assume that we can use the common DType of the signature
+ * and the Python scalar DType (`in_DT`) as a safe intermediate.
+ */
+NPY_NO_EXPORT PyArray_Descr *
+npy_find_descr_for_scalar(
+    PyObject *scalar, PyArray_Descr *original_descr,
+    PyArray_DTypeMeta *in_DT, PyArray_DTypeMeta *op_DT)
+{
+    PyArray_Descr *res;
+    /* There is a good chance, descriptors already match... */
+    if (NPY_DTYPE(original_descr) == op_DT) {
+        Py_INCREF(original_descr);
+        return original_descr;
+    }
+
+    PyArray_DTypeMeta *common = PyArray_CommonDType(in_DT, op_DT);
+    if (common == NULL) {
+        PyErr_Clear();
+        /* This is fine.  We simply assume the original descr is viable. */
+        Py_INCREF(original_descr);
+        return original_descr;
+    }
+    /* A very likely case is that there is nothing to do: */
+    if (NPY_DTYPE(original_descr) == common) {
+        Py_DECREF(common);
+        Py_INCREF(original_descr);
+        return original_descr;
+    }
+    if (!NPY_DT_is_parametric(common) ||
+            /* In some paths we only have a scalar type, can't discover */
+            scalar == NULL ||
+            /* If the DType doesn't know the scalar type, guess at default. */
+            !NPY_DT_CALL_is_known_scalar_type(common, Py_TYPE(scalar))) {
+        if (common->singleton != NULL) {
+            Py_INCREF(common->singleton);
+            res = common->singleton;
+            Py_INCREF(res);
+        }
+        else {
+            res = NPY_DT_CALL_default_descr(common);
+        }
+    }
+    else {
+        res = NPY_DT_CALL_discover_descr_from_pyobject(common, scalar);
+    }
+
+    Py_DECREF(common);
+    return res;
+}
diff --git a/numpy/_core/src/multiarray/abstractdtypes.h b/numpy/_core/src/multiarray/abstractdtypes.h
index 7bf8191e6917..3c96ffe8e0ef 100644
--- a/numpy/_core/src/multiarray/abstractdtypes.h
+++ b/numpy/_core/src/multiarray/abstractdtypes.h
@@ -1,6 +1,7 @@
 #ifndef NUMPY_CORE_SRC_MULTIARRAY_ABSTRACTDTYPES_H_
 #define NUMPY_CORE_SRC_MULTIARRAY_ABSTRACTDTYPES_H_
 
+#include "numpy/ndarraytypes.h"
 #include "arrayobject.h"
 #include "dtypemeta.h"
 
@@ -68,6 +69,19 @@ npy_mark_tmp_array_if_pyscalar(
     return 0;
 }
 
+
+NPY_NO_EXPORT int
+npy_update_operand_for_scalar(
+    PyArrayObject **operand, PyObject *scalar, PyArray_Descr *descr,
+    NPY_CASTING casting);
+
+
+NPY_NO_EXPORT PyArray_Descr *
+npy_find_descr_for_scalar(
+    PyObject *scalar, PyArray_Descr *original_descr,
+    PyArray_DTypeMeta *in_DT, PyArray_DTypeMeta *op_DT);
+
+
 #ifdef __cplusplus
 }
 #endif
diff --git a/numpy/_core/src/multiarray/array_assign_scalar.c b/numpy/_core/src/multiarray/array_assign_scalar.c
index 6818c1aa2a1b..0199ba969eb9 100644
--- a/numpy/_core/src/multiarray/array_assign_scalar.c
+++ b/numpy/_core/src/multiarray/array_assign_scalar.c
@@ -243,8 +243,7 @@ PyArray_AssignRawScalar(PyArrayObject *dst,
     }
 
     /* Check the casting rule */
-    if (!can_cast_scalar_to(src_dtype, src_data,
-                            PyArray_DESCR(dst), casting)) {
+    if (!PyArray_CanCastTypeTo(src_dtype, PyArray_DESCR(dst), casting)) {
         npy_set_invalid_cast_error(
                 src_dtype, PyArray_DESCR(dst), casting, NPY_TRUE);
         return -1;
diff --git a/numpy/_core/src/multiarray/multiarraymodule.c b/numpy/_core/src/multiarray/multiarraymodule.c
index 944898ceecf7..83fc505a97e6 100644
--- a/numpy/_core/src/multiarray/multiarraymodule.c
+++ b/numpy/_core/src/multiarray/multiarraymodule.c
@@ -1929,21 +1929,55 @@ array_asfortranarray(PyObject *NPY_UNUSED(ignored),
 
 
 static PyObject *
-array_copyto(PyObject *NPY_UNUSED(ignored), PyObject *args, PyObject *kwds)
+array_copyto(PyObject *NPY_UNUSED(ignored),
+        PyObject *const *args, Py_ssize_t len_args, PyObject *kwnames)
 {
-    static char *kwlist[] = {"dst", "src", "casting", "where", NULL};
-    PyObject *wheremask_in = NULL;
-    PyArrayObject *dst = NULL, *src = NULL, *wheremask = NULL;
+    PyObject *dst_obj, *src_obj, *wheremask_in = NULL;
+    PyArrayObject *src = NULL, *wheremask = NULL;
     NPY_CASTING casting = NPY_SAME_KIND_CASTING;
+    NPY_PREPARE_ARGPARSER;
 
-    if (!PyArg_ParseTupleAndKeywords(args, kwds, "O!O&|O&O:copyto", kwlist,
-                &PyArray_Type, &dst,
-                &PyArray_Converter, &src,
-                &PyArray_CastingConverter, &casting,
-                &wheremask_in)) {
+    if (npy_parse_arguments("copyto", args, len_args, kwnames,
+            "dst", NULL, &dst_obj,
+            "src", NULL, &src_obj,
+            "|casting", &PyArray_CastingConverter, &casting,
+            "|where", NULL, &wheremask_in,
+            NULL, NULL, NULL) < 0) {
         goto fail;
     }
 
+    if (!PyArray_Check(dst_obj)) {
+        PyErr_Format(PyExc_TypeError,
+            "copyto() argument 1 must be a numpy.ndarray, not %s",
+            Py_TYPE(dst_obj)->tp_name);
+    }
+    PyArrayObject *dst = (PyArrayObject *)dst_obj;
+
+    src = (PyArrayObject *)PyArray_FromAny(src_obj, NULL, 0, 0, 0, NULL);
+    if (src == NULL) {
+        goto fail;
+    }
+    PyArray_DTypeMeta *dtype = NPY_DTYPE(PyArray_DESCR(src));
+    Py_INCREF(dtype);
+    if (npy_mark_tmp_array_if_pyscalar(src_obj, src, &dtype)) {
+        /* The user passed a Python scalar */
+        PyArray_Descr *descr = npy_find_descr_for_scalar(
+                src_obj, PyArray_DESCR(src), dtype,
+                NPY_DTYPE(PyArray_DESCR(dst)));
+        Py_DECREF(dtype);
+        if (descr == NULL) {
+            goto fail;
+        }
+        int res = npy_update_operand_for_scalar(&src, src_obj, descr, casting);
+        Py_DECREF(descr);
+        if (res < 0) {
+            goto fail;
+        }
+    }
+    else {
+        Py_DECREF(dtype);
+    }
+
     if (wheremask_in != NULL) {
         /* Get the boolean where mask */
         PyArray_Descr *dtype = PyArray_DescrFromType(NPY_BOOL);
@@ -4431,7 +4465,7 @@ static struct PyMethodDef array_module_methods[] = {
         METH_FASTCALL | METH_KEYWORDS, NULL},
     {"copyto",
         (PyCFunction)array_copyto,
-        METH_VARARGS|METH_KEYWORDS, NULL},
+        METH_FASTCALL | METH_KEYWORDS, NULL},
     {"nested_iters",
         (PyCFunction)NpyIter_NestedIters,
         METH_VARARGS|METH_KEYWORDS, NULL},
@@ -5129,7 +5163,7 @@ PyMODINIT_FUNC PyInit__multiarray_umath(void) {
 
     // initialize static reference to a zero-like array
     npy_static_pydata.zero_pyint_like_arr = PyArray_ZEROS(
-            0, NULL, NPY_LONG, NPY_FALSE);
+            0, NULL, NPY_DEFAULT_INT, NPY_FALSE);
     if (npy_static_pydata.zero_pyint_like_arr == NULL) {
         goto err;
     }
diff --git a/numpy/_core/src/umath/ufunc_object.c b/numpy/_core/src/umath/ufunc_object.c
index 92bc7793f2ad..9aac8d8cf188 100644
--- a/numpy/_core/src/umath/ufunc_object.c
+++ b/numpy/_core/src/umath/ufunc_object.c
@@ -101,8 +101,8 @@ static int
 resolve_descriptors(int nop,
         PyUFuncObject *ufunc, PyArrayMethodObject *ufuncimpl,
         PyArrayObject *operands[], PyArray_Descr *dtypes[],
-        PyArray_DTypeMeta *signature[], PyObject *inputs_tup,
-        NPY_CASTING casting);
+        PyArray_DTypeMeta *signature[], PyArray_DTypeMeta *original_DTypes[],
+        PyObject *inputs_tup, NPY_CASTING casting);
 
 
 /*UFUNC_API*/
@@ -2356,11 +2356,12 @@ reducelike_promote_and_resolve(PyUFuncObject *ufunc,
     if (evil_ndim_mutating_hack) {
         ((PyArrayObject_fields *)out)->nd = 0;
     }
-    /* DTypes may currently get filled in fallbacks and XDECREF for error: */
-    Py_XDECREF(operation_DTypes[0]);
-    Py_XDECREF(operation_DTypes[1]);
-    Py_XDECREF(operation_DTypes[2]);
+    // TODO: Clean up multiple cleanup!
     if (ufuncimpl == NULL) {
+        /* DTypes may currently get filled in fallbacks and XDECREF for error: */
+        Py_XDECREF(operation_DTypes[0]);
+        Py_XDECREF(operation_DTypes[1]);
+        Py_XDECREF(operation_DTypes[2]);
         return NULL;
     }
 
@@ -2372,9 +2373,17 @@ reducelike_promote_and_resolve(PyUFuncObject *ufunc,
      * (although this should possibly happen through a deprecation)
      */
     if (resolve_descriptors(3, ufunc, ufuncimpl,
-            ops, out_descrs, signature, NULL, casting) < 0) {
+            ops, out_descrs, signature, operation_DTypes, NULL, casting) < 0) {
+        /* DTypes may currently get filled in fallbacks and XDECREF for error: */
+        Py_XDECREF(operation_DTypes[0]);
+        Py_XDECREF(operation_DTypes[1]);
+        Py_XDECREF(operation_DTypes[2]);
         return NULL;
     }
+    /* DTypes may currently get filled in fallbacks and XDECREF for error: */
+    Py_XDECREF(operation_DTypes[0]);
+    Py_XDECREF(operation_DTypes[1]);
+    Py_XDECREF(operation_DTypes[2]);
 
     /*
      * The first operand and output should be the same array, so they should
@@ -4023,12 +4032,12 @@ static int
 resolve_descriptors(int nop,
         PyUFuncObject *ufunc, PyArrayMethodObject *ufuncimpl,
         PyArrayObject *operands[], PyArray_Descr *dtypes[],
-        PyArray_DTypeMeta *signature[], PyObject *inputs_tup,
-        NPY_CASTING casting)
+        PyArray_DTypeMeta *signature[], PyArray_DTypeMeta *original_DTypes[],
+        PyObject *inputs_tup, NPY_CASTING casting)
 {
     int retval = -1;
     NPY_CASTING safety;
-    PyArray_Descr *original_dtypes[NPY_MAXARGS];
+    PyArray_Descr *original_descrs[NPY_MAXARGS];
 
     NPY_UF_DBG_PRINT("Resolving the descriptors\n");
 
@@ -4043,12 +4052,12 @@ resolve_descriptors(int nop,
         PyObject *input_scalars[NPY_MAXARGS];
         for (int i = 0; i < nop; i++) {
             if (operands[i] == NULL) {
-                original_dtypes[i] = NULL;
+                original_descrs[i] = NULL;
             }
             else {
                 /* For abstract DTypes, we might want to change what this is */
-                original_dtypes[i] = PyArray_DTYPE(operands[i]);
-                Py_INCREF(original_dtypes[i]);
+                original_descrs[i] = PyArray_DTYPE(operands[i]);
+                Py_INCREF(original_descrs[i]);
             }
             /*
              * Check whether something is a scalar of the given type.
@@ -4067,27 +4076,71 @@ resolve_descriptors(int nop,
 
         npy_intp view_offset = NPY_MIN_INTP;  /* currently ignored */
         safety = ufuncimpl->resolve_descriptors_with_scalars(
-            ufuncimpl, signature, original_dtypes, input_scalars,
+            ufuncimpl, signature, original_descrs, input_scalars,
             dtypes, &view_offset
         );
+
+        /* For scalars, replace the operand if needed (scalars can't be out) */
+        for (int i = 0; i < nin; i++) {
+            if ((PyArray_FLAGS(operands[i]) & NPY_ARRAY_WAS_PYTHON_LITERAL)) {
+                /* `resolve_descriptors_with_scalars` decides the descr */
+                if (npy_update_operand_for_scalar(
+                        &operands[i], input_scalars[i], dtypes[i],
+                        /* ignore cast safety for this op (resolvers job) */
+                        NPY_SAFE_CASTING) < 0) {
+                    goto finish;
+                }
+            }
+        }
         goto check_safety;
     }
 
     for (int i = 0; i < nop; ++i) {
         if (operands[i] == NULL) {
-            original_dtypes[i] = NULL;
+            original_descrs[i] = NULL;
+            continue;
         }
-        else {
-            /*
-             * The dtype may mismatch the signature, in which case we need
-             * to make it fit before calling the resolution.
-             */
-            PyArray_Descr *descr = PyArray_DTYPE(operands[i]);
-            original_dtypes[i] = PyArray_CastDescrToDType(descr, signature[i]);
-            if (original_dtypes[i] == NULL) {
+        PyArray_Descr *descr = PyArray_DTYPE(operands[i]);
+
+        /*
+         * If we are working with Python literals/scalars, deal with them.
+         * If needed, we create new array with the right descriptor.
+         */
+        if ((PyArray_FLAGS(operands[i]) & NPY_ARRAY_WAS_PYTHON_LITERAL)) {
+            PyObject *input;
+            if (inputs_tup == NULL) {
+                input = NULL;
+            }
+            else {
+                input = PyTuple_GET_ITEM(inputs_tup, i);
+            }
+
+            PyArray_Descr *new_descr = npy_find_descr_for_scalar(
+                    input, descr, original_DTypes[i], signature[i]);
+            if (new_descr == NULL) {
                 nop = i;  /* only this much is initialized */
                 goto finish;
             }
+            int res = npy_update_operand_for_scalar(
+                &operands[i], input, new_descr, casting);
+            Py_DECREF(new_descr);
+            if (res < 0) {
+                nop = i;  /* only this much is initialized */
+                goto finish;
+            }
+
+            /* Descriptor may have been modified along the way */
+            descr = PyArray_DESCR(operands[i]);
+        }
+
+        /*
+         * The dtype may mismatch the signature, in which case we need
+         * to make it fit before calling the resolution.
+         */
+        original_descrs[i] = PyArray_CastDescrToDType(descr, signature[i]);
+        if (original_descrs[i] == NULL) {
+            nop = i;  /* only this much is initialized */
+            goto finish;
         }
     }
 
@@ -4096,7 +4149,7 @@ resolve_descriptors(int nop,
         npy_intp view_offset = NPY_MIN_INTP;  /* currently ignored */
 
         safety = ufuncimpl->resolve_descriptors(ufuncimpl,
-                signature, original_dtypes, dtypes, &view_offset);
+                signature, original_descrs, dtypes, &view_offset);
         goto check_safety;
     }
     else {
@@ -4124,7 +4177,7 @@ resolve_descriptors(int nop,
 
   finish:
     for (int i = 0; i < nop; i++) {
-        Py_XDECREF(original_dtypes[i]);
+        Py_XDECREF(original_descrs[i]);
     }
     return retval;
 }
@@ -4467,49 +4520,11 @@ ufunc_generic_fastcall(PyUFuncObject *ufunc,
 
     /* Find the correct descriptors for the operation */
     if (resolve_descriptors(nop, ufunc, ufuncimpl,
-            operands, operation_descrs, signature, full_args.in, casting) < 0) {
+            operands, operation_descrs, signature, operand_DTypes,
+            full_args.in, casting) < 0) {
         goto fail;
     }
 
-    if (promoting_pyscalars) {
-        /*
-         * Python integers need to be cast specially.  For other python
-         * scalars it does not hurt either.  It would be nice to never create
-         * the array in this case, but that is difficult until value-based
-         * promotion rules are gone.  (After that, we may get away with using
-         * dummy arrays rather than real arrays for the legacy resolvers.)
-         */
-        for (int i = 0; i < nin; i++) {
-            int orig_flags = PyArray_FLAGS(operands[i]);
-            if (!(orig_flags & NPY_ARRAY_WAS_PYTHON_LITERAL)) {
-                continue;
-            }
-            /*
-             * If descriptor matches, no need to convert, but integers may
-             * have been too large.
-             */
-            if (!(orig_flags & NPY_ARRAY_WAS_INT_AND_REPLACED)
-                    && PyArray_EquivTypes(
-                        PyArray_DESCR(operands[i]), operation_descrs[i])) {
-                continue;
-            }
-            /* Otherwise, replace the operand with a new array */
-            PyArray_Descr *descr = operation_descrs[i];
-            Py_INCREF(descr);
-            PyArrayObject *new = (PyArrayObject *)PyArray_NewFromDescr(
-                    &PyArray_Type, descr, 0, NULL, NULL, NULL, 0, NULL);
-            Py_SETREF(operands[i], new);
-            if (operands[i] == NULL) {
-                goto fail;
-            }
-
-            PyObject *value = PyTuple_GET_ITEM(full_args.in, i);
-            if (PyArray_SETITEM(new, PyArray_BYTES(operands[i]), value) < 0) {
-                goto fail;
-            }
-        }
-    }
-
     /*
      * Do the final preparations and call the inner-loop.
      */
@@ -5827,7 +5842,7 @@ ufunc_at(PyUFuncObject *ufunc, PyObject *args)
 
         /* Find the correct operation_descrs for the operation */
         int resolve_result = resolve_descriptors(nop, ufunc, ufuncimpl,
-                tmp_operands, operation_descrs, signature, NULL, NPY_UNSAFE_CASTING);
+                tmp_operands, operation_descrs, signature, operand_DTypes, NULL, NPY_UNSAFE_CASTING);
         for (int i = 0; i < 3; i++) {
             Py_XDECREF(signature[i]);
             Py_XDECREF(operand_DTypes[i]);
@@ -6152,7 +6167,7 @@ py_resolve_dtypes_generic(PyUFuncObject *ufunc, npy_bool return_context,
 
         /* Find the correct descriptors for the operation */
         if (resolve_descriptors(ufunc->nargs, ufunc, ufuncimpl,
-                dummy_arrays, operation_descrs, signature,
+                dummy_arrays, operation_descrs, signature, DTypes,
                 NULL, casting) < 0) {
             goto finish;
         }
diff --git a/numpy/_core/tests/test_api.py b/numpy/_core/tests/test_api.py
index 5b9bdb60f1b3..1ac7a49b3610 100644
--- a/numpy/_core/tests/test_api.py
+++ b/numpy/_core/tests/test_api.py
@@ -406,6 +406,37 @@ def test_copyto():
     # 'dst' must be an array
     assert_raises(TypeError, np.copyto, [1, 2, 3], [2, 3, 4])
 
+
+def test_copyto_cast_safety():
+    with pytest.raises(TypeError):
+        np.copyto(np.arange(3), 3., casting="safe")
+
+    # Can put integer and float scalars safely (and equiv):
+    np.copyto(np.arange(3), 3, casting="equiv")
+    np.copyto(np.arange(3.), 3., casting="equiv")
+    # And also with less precision safely:
+    np.copyto(np.arange(3, dtype="uint8"), 3, casting="safe")
+    np.copyto(np.arange(3., dtype="float32"), 3., casting="safe")
+
+    # But not equiv:
+    with pytest.raises(TypeError):
+        np.copyto(np.arange(3, dtype="uint8"), 3, casting="equiv")
+
+    with pytest.raises(TypeError):
+        np.copyto(np.arange(3., dtype="float32"), 3., casting="equiv")
+
+    # As a special thing, object is equiv currently:
+    np.copyto(np.arange(3, dtype=object), 3, casting="equiv")
+
+    # The following raises an overflow error/gives a warning but not
+    # type error (due to casting), though:
+    with pytest.raises(OverflowError):
+        np.copyto(np.arange(3), 2**80, casting="safe")
+
+    with pytest.warns(RuntimeWarning):
+        np.copyto(np.arange(3, dtype=np.float32), 2e300, casting="safe")
+
+
 def test_copyto_permut():
     # test explicit overflow case
     pad = 500
diff --git a/numpy/_core/tests/test_numeric.py b/numpy/_core/tests/test_numeric.py
index a9694ebee736..b665db5701a9 100644
--- a/numpy/_core/tests/test_numeric.py
+++ b/numpy/_core/tests/test_numeric.py
@@ -3472,7 +3472,11 @@ def test_empty_like(self):
     def test_filled_like(self):
         self.check_like_function(np.full_like, 0, True)
         self.check_like_function(np.full_like, 1, True)
-        self.check_like_function(np.full_like, 1000, True)
+        # Large integers may overflow, but using int64 is OK (casts)
+        # see also gh-27075
+        with pytest.raises(OverflowError):
+            np.full_like(np.ones(3, dtype=np.int8), 1000)
+        self.check_like_function(np.full_like, np.int64(1000), True)
         self.check_like_function(np.full_like, 123.456, True)
         # Inf to integer casts cause invalid-value errors: ignore them.
         with np.errstate(invalid="ignore"):
diff --git a/numpy/_core/tests/test_ufunc.py b/numpy/_core/tests/test_ufunc.py
index e777d7e07be3..a605b65079ba 100644
--- a/numpy/_core/tests/test_ufunc.py
+++ b/numpy/_core/tests/test_ufunc.py
@@ -616,6 +616,31 @@ def call_ufunc(arr, **kwargs):
         expected = call_ufunc(arr.astype(np.float64))  # upcast
         assert_array_equal(expected, res)
 
+    @pytest.mark.parametrize("ufunc", [np.add, np.equal])
+    def test_cast_safety_scalar(self, ufunc):
+        # We test add and equal, because equal has special scalar handling
+        # Note that the "equiv" casting behavior should maybe be considered
+        # a current implementation detail.
+        with pytest.raises(TypeError):
+            # this picks an integer loop, which is not safe
+            ufunc(3., 4., dtype=int, casting="safe")
+
+        with pytest.raises(TypeError):
+            # We accept python float as float64 but not float32 for equiv.
+            ufunc(3., 4., dtype="float32", casting="equiv")
+
+        # Special case for object and equal (note that equiv implies safe)
+        ufunc(3, 4, dtype=object, casting="equiv")
+        # Picks a double loop for both, first is equiv, second safe:
+        ufunc(np.array([3.]), 3., casting="equiv")
+        ufunc(np.array([3.]), 3, casting="safe")
+        ufunc(np.array([3]), 3, casting="equiv")
+
+    def test_cast_safety_scalar_special(self):
+        # We allow this (and it succeeds) via object, although the equiv
+        # part may not be important.
+        np.equal(np.array([3]), 2**300, casting="equiv")
+
     def test_true_divide(self):
         a = np.array(10)
         b = np.array(20)


Code review comments:
diff_hunk: @@ -4067,27 +4076,71 @@ resolve_descriptors(int nop,
 
         npy_intp view_offset = NPY_MIN_INTP;  /* currently ignored */
         safety = ufuncimpl->resolve_descriptors_with_scalars(
-            ufuncimpl, signature, original_dtypes, input_scalars,
+            ufuncimpl, signature, original_descrs, input_scalars,
             dtypes, &view_offset
         );
+
+        /* For scalars, replace the operand if needed (scalars can't be out) */
+        for (int i = 0; i < nin; i++) {
+            if ((PyArray_FLAGS(operands[i]) & NPY_ARRAY_WAS_PYTHON_LITERAL)) {
+                /* `resolve_descriptors_with_scalars` decides the descr */
+                if (npy_update_operand_for_scalar(
+                        &operands[i], input_scalars[i], dtypes[i],
+                        /* ignore cast safety for this op (resolvers job) */
+                        NPY_SAFE_CASTING) < 0) {
+                    goto finish;
+                }
+            }
+        }
         goto check_safety;
     }
 
     for (int i = 0; i < nop; ++i) {
         if (operands[i] == NULL) {
-            original_dtypes[i] = NULL;
+            original_descrs[i] = NULL;
+            continue;
         }
-        else {
-            /*
-             * The dtype may mismatch the signature, in which case we need
-             * to make it fit before calling the resolution.
-             */
-            PyArray_Descr *descr = PyArray_DTYPE(operands[i]);
-            original_dtypes[i] = PyArray_CastDescrToDType(descr, signature[i]);
-            if (original_dtypes[i] == NULL) {
+        PyArray_Descr *descr = PyArray_DTYPE(operands[i]);
+
+        /*
+         * If we are working with Python literals/scalars, deal with them.
+         * If needed, we create new array with the right descriptor.
+         */
+        if ((PyArray_FLAGS(operands[i]) & NPY_ARRAY_WAS_PYTHON_LITERAL)) {
+            PyObject *input;
+            if (inputs_tup == NULL) {
+                input = NULL;
+            }
+            else {
+                input = PyTuple_GET_ITEM(inputs_tup, i);
+            }
+
+            PyArray_Descr *new_descr = npy_find_descr_for_scalar(
+                    input, descr, original_DTypes[i], signature[i]);
+            if (new_descr == NULL) {
                 nop = i;  /* only this much is initialized */
                 goto finish;
             }
+            int res = npy_update_operand_for_scalar(
+                &operands[i], input, new_descr, casting);
+            Py_DECREF(new_descr);
+            if (res < 0) {
+                nop = i;  /* only this much is initialized */
path: numpy/_core/src/umath/ufunc_object.c
text: @user1:
This works but is rather ugly; can we not ensure the array is filled with NULL to start with? Or have an explicit `int ninit = 0` that gets increased as things are added?

@author:
FWIW, I added an `n_cleanup`.  I agree, it's a bit tricky with the cleanup although not sure it is a real improvement.  I don't really love defensive NULLing that much more.

@user1:
OK, I think this at least is clearer. I somewhat agree that NULLing is also not great...