Feat: numpy array parameters, step 2: warnings and division (#220)

docs/SPEC.md

@@ -689,7 +689,7 @@ All in `PM/commonTest`, so they run wherever the interpreter loads.
   | Array parameters only indexed, stored, `len`-ed; no aliasing (entry guard) | not compiled; diagnostic |
   | Integer ops proven safe by the verifier's intervals (constants, `len`, `range` variables, guards such as `if x < e` / `while lo < hi`, joined local assignments, `//` and `%` by constants) | emitted unchecked (a shift/mask for `//`/`%` by a power of two of a non-negative value); the sanitizer CI run builds with `TP_CHECK_PROVEN` and turns a wrong proof into a SystemError (#148) |
   | `array.array[float]` / `array.array[int]` parameters (exact `array.array`, typecode `d` or `q`/`l`, 8-byte items) are read and written in place: the buffer export is held for the call and released on every exit, and every element access is one relaxed atomic 8-byte load or store, so each access is atomic like the interpreter's `a[i]`; no copy-in and no limit on the number of buffers. Open (#149): while the export is held, another thread's resize of that array raises BufferError. numpy arrays: see the next row | anything else: object path |
-  | numpy `npt.NDArray[np.float64]` / `NDArray[np.int64]` parameters (exact `numpy.ndarray`, 1-D, C-contiguous, native byte order, writable when stored to) take the same buffer path (#220, step 1) when the function accesses the array purely by element (`x[i]`, `len(x)`); a function that uses it as a whole object (`np.sum(x)`, `x.T`, passing `x` on, a slice) compiles as before with the parameter an opaque object. Other dtypes, including bool, are opaque objects. Results follow the interpreter: a value derived from an element is a numpy scalar there, so it is boxed as `np.float64` / `np.int64` / `np.bool` when it escapes (return, `Box`), decided per value at run time by a taint flag; a function declared `-> float` / `int` / `bool` that holds a numpy array returns the interpreter's object. `float(v)` of such a value is a plain float. An out-of-range index raises numpy's `IndexError("index i is out of bounds for axis 0 with size n")`. The numpy RuntimeWarnings (overflow, invalid value) and numpy's division are not emulated in this step: a pure function deopts on a non-finite float `+ - *` result or an unproved int overflow, so the interpreter emits them; `/`, `//`, `%` on a value derived from an element, and `+ - *` on one in an impure function, are not compiled | anything else: object path |
+  | numpy `npt.NDArray[np.float64]` / `NDArray[np.int64]` parameters (exact `numpy.ndarray`, 1-D, C-contiguous, native byte order, writable when stored to) take the same buffer path (#220, step 1) when the function accesses the array purely by element (`x[i]`, `len(x)`); a function that uses it as a whole object (`np.sum(x)`, `x.T`, passing `x` on, a slice) compiles as before with the parameter an opaque object. Other dtypes, including bool, are opaque objects. Results follow the interpreter: a value derived from an element is a numpy scalar there, so it is boxed as `np.float64` / `np.int64` / `np.bool` when it escapes (return, `Box`), decided per value at run time by a taint flag; a function declared `-> float` / `int` / `bool` that holds a numpy array returns the interpreter's object. `float(v)` of such a value is a plain float. An out-of-range index raises numpy's `IndexError("index i is out of bounds for axis 0 with size n")`. Arithmetic on a value derived from an element follows numpy's scalar semantics (#220, step 2): float `/`, `//`, `%` never raise ZeroDivisionError (inf / nan, with numpy's result bits and sign rules, `//` and `%` as numpy's `divmod`), int64 `//` and `%` by zero give 0, int64 `/` converts both operands to double. numpy's RuntimeWarning (category and text exactly as numpy's, e.g. `overflow encountered in scalar add`, `divide by zero encountered in scalar divide`, `invalid value encountered in scalar floor_divide`) is emitted once per flagged operation, with the source file and line of the operation. A pure function does not emit it itself: on any floating-point exception (overflow, invalid, divide by zero) or int64 wrap or division by zero it deopts, so the interpreter emits it; an impure function cannot deopt after an effect, so the compiled code emits the warning (`PyErr_WarnExplicitObject`, honouring `warnings` filters and the module's `__warningregistry__`). Float `+ - * / // %` on such values is compiled in both kinds of function, including a store of the result into the array (`x[i] -= d * m`); int64 `/ // %` is compiled in a pure function only. Int64 `+ - *` and `// %` of an element in an impure function (a native INT would have to be a numpy scalar) stay on the object path | anything else: object path |
   | A function holding an array parameter runs no user code: closed, and no slot read or `C(args)` whose slow path runs user code, in its whole call tree (#203) | not compiled; diagnostic |
   | OBJ references follow one ownership rule; no IR form for free, borrowed alias or NULL read | not compiled; diagnostic |
   | Integer overflow checked and promoted (deopt); an op marked `proven` must re-prove its interval fits i64 | not compiled; diagnostic `verify/proven-overflow` |

python-multiplatform-ksp/src/main/python/typedpython/cgen.py

@@ -202,6 +202,15 @@ class CGenError(Exception):
            ir.BinOpKind.MOD: "tp_mod_i64", ir.BinOpKind.TRUEDIV: "tp_truediv_i64"}
 _F64_OP = {ir.BinOpKind.TRUEDIV: "tp_truediv_f64", ir.BinOpKind.FLOORDIV: "tp_floordiv_f64",
            ir.BinOpKind.MOD: "tp_mod_f64"}
+_NP_I64_OP = {ir.BinOpKind.TRUEDIV: "tp_np_truediv_i64", ir.BinOpKind.FLOORDIV: "tp_np_floordiv_i64",
+              ir.BinOpKind.MOD: "tp_np_mod_i64"}
+_NP_F64_CODE = {k: f"TP_NP_{k.name}" for k in (
+    ir.BinOpKind.ADD, ir.BinOpKind.SUB, ir.BinOpKind.MUL, ir.BinOpKind.TRUEDIV,
+    ir.BinOpKind.FLOORDIV, ir.BinOpKind.MOD)}
+_NP_F64_CODE[ir.BinOpKind.TRUEDIV] = "TP_NP_DIV"
+_NP_NAME = {ir.BinOpKind.ADD: "add", ir.BinOpKind.SUB: "subtract", ir.BinOpKind.MUL: "multiply",
+            ir.BinOpKind.TRUEDIV: "divide", ir.BinOpKind.FLOORDIV: "floor_divide",
+            ir.BinOpKind.MOD: "remainder"}
 _F64_INLINE = {ir.BinOpKind.ADD: "+", ir.BinOpKind.SUB: "-", ir.BinOpKind.MUL: "*"}
 _OBJ_OP = {k: f"TP_BINOP_{k.name}" for k in ir.BinOpKind}
 # SPEC N-13: the INT operations (tp_runtime.h "Native ints"); TRUEDIV of ints is an OBJ operation
@@ -482,6 +491,7 @@ def visit(name: str) -> int | None:
 class _ModGen:
     def __init__(self, module: ir.Module):
         self.module = module
+        self.display_path = ""
         self.functions = {f.name: f for f in module.functions}
         if len(self.functions) != len(module.functions):
             raise CGenError("duplicate compiled function names")
@@ -631,6 +641,7 @@ def generate(module: ir.Module, source_path: Path, display_path: str | None = No
     display_path = source_path.name if display_path is None else display_path
     source = source_path.read_text(encoding="utf-8")
     mg = _ModGen(module)
+    mg.display_path = display_path       # the file name of a compiled numpy RuntimeWarning (#220)
     for f in module.functions:
         _check_function(f, mg.functions)
     for f in module.functions:
@@ -1042,16 +1053,28 @@ def ev_BinOp(self, e: ir.BinOp) -> str:
                     self.emit(f"{t} = {a} & {_i64_lit((1 << k) - 1)};")
             elif e.proven:
                 self.emit(f"tp_proven({e.op.name.lower()}, {a}, {b}, {t});")
+            elif e.np and e.op in _NP_I64_OP:
+                # #220: numpy's int64 / // % (division by zero, INT64_MIN // -1): a deopt, the
+                # interpreter then warns (the verifier allows it in a pure function only)
+                self.check(f"{_NP_I64_OP[e.op]}({a}, {b}, &{t})")
             else:
                 self.check(f"{_I64_OP[e.op]}({a}, {b}, &{t})")
             return t
+        if lt is Type.F64 and e.np:
+            t = self.tmp("double")
+            self.emit(f"tp_s = tp_np_f64({_NP_F64_CODE[e.op]}, {a}, {b}, &{t});")
+            if self.f.pure:
+                # #220: numpy would warn: the interpreter redoes the call and emits it
+                self.emit("if (tp_s != 0) { tp_rc = 1; goto tp_exit; }")
+            else:
+                fi = self.mg.name_index(self.mg.display_path)
+                self.emit(f"if (tp_s != 0 && tp_np_warn(tp_module, tp_s, \"{_NP_NAME[e.op]}\", "
+                          f"tp_st->names[{fi}], {e.line}) != 0) {{ tp_rc = -1; goto tp_exit; }}")
+            return t
         if lt is Type.F64:
             t = self.tmp("double")
             if e.op in _F64_INLINE:
                 self.emit(f"{t} = {a} {_F64_INLINE[e.op]} {b};")
-                if e.np:
-                    # #220: numpy would warn (overflow / invalid value); the interpreter redoes the call
-                    self.emit(f"if (!isfinite({t})) {{ tp_rc = 1; goto tp_exit; }}")
             else:
                 self.check(f"{_F64_OP[e.op]}({a}, {b}, &{t})")
             return t

python-multiplatform-ksp/src/main/python/typedpython/frontend.py

@@ -1825,6 +1825,13 @@ def to_int(self, e: ir.Expr) -> ir.Expr:
 
     def arith(self, op: ir.BinOpKind, left: ir.Expr, right: ir.Expr, node: ast.AST,
               inplace: bool = False) -> ir.Expr:
+        e = self._arith(op, left, right, node, inplace)
+        if isinstance(e, ir.BinOp) and e.type in ir.SCALARS and self.numpy:
+            e = dataclasses.replace(e, line=node.lineno)      # the warning site (#220)
+        return e
+
+    def _arith(self, op: ir.BinOpKind, left: ir.Expr, right: ir.Expr, node: ast.AST,
+               inplace: bool = False) -> ir.Expr:
         lt, rt = left.type, right.type
         if Type.BOOL in (lt, rt) or Type.NONE in (lt, rt):
             raise Skip(node.lineno, f"arithmetic on {TYPE_NAME[Type.BOOL if Type.BOOL in (lt, rt) else Type.NONE]} "

python-multiplatform-ksp/src/main/python/typedpython/ir.py

@@ -160,7 +160,13 @@ class BinOp(Expr):
     # F64 ADD/SUB/MUL only (#220): an operand may be a numpy scalar in the interpreter, which warns
     # (RuntimeWarning: overflow / invalid value) where a double does not. A non-finite result →
     # deopt, so the interpreter redoes the call and emits the warning; legal only in a pure function.
+    # Step 2 (#220): `np` also covers F64 TRUEDIV / FLOORDIV / MOD (numpy's results: inf / nan, no
+    # ZeroDivisionError) and unproven I64 TRUEDIV / FLOORDIV / MOD (numpy's int64: division by zero
+    # gives 0 / inf / nan; both are a deopt, and legal only in a pure function). Any floating-point
+    # exception (overflow, invalid, divide by zero) is a deopt in a pure function; in an impure F64
+    # operation the code emits numpy's RuntimeWarning itself, at `line` of the cgen display path.
     np: bool = False
+    line: int = 0      # source line of the operation (the warning site of an impure `np` operation)
 
 
 class UnaryOpKind(Enum):

python-multiplatform-ksp/src/main/python/typedpython/runtime/tp_runtime.h

@@ -14,6 +14,7 @@
 #include <Python.h>
 #include <stdint.h>
 #include <math.h>
+#include <fenv.h>
 #include <errno.h>
 #include <string.h>
 #include <stdlib.h>
@@ -398,6 +399,175 @@ static inline int tp_mod_f64(double a, double b, double *out)
     return 0;
 }
 
+/* ------------------------------------------------------------------------------------------
+ * numpy scalar arithmetic (#220, step 2). An operation on a value derived from a numpy element is
+ * numpy's, not CPython's: no ZeroDivisionError, inf / nan results, int64 division by zero, and a
+ * RuntimeWarning per floating-point exception. tp_np_f64 returns the exceptions that numpy reports
+ * (TP_NP_* bits). A pure function deopts on any bit (the interpreter then warns); an impure one calls
+ * tp_np_warn. The algorithms are numpy's (npymath: npy_divmod / npy_floor_divide / npy_remainder).
+ * ------------------------------------------------------------------------------------------ */
+
+#define TP_NP_DIVZERO 1
+#define TP_NP_OVERFLOW 2
+#define TP_NP_INVALID 4
+enum { TP_NP_ADD, TP_NP_SUB, TP_NP_MUL, TP_NP_DIV, TP_NP_FLOORDIV, TP_NP_MOD };
+
+#if defined(__GNUC__) || defined(__clang__)
+#define TP_NOINLINE __attribute__((noinline, unused))
+#else
+#define TP_NOINLINE
+#endif
+
+static inline double tp_np_divmod(double a, double b, double *modulus)
+{
+    double div, mod, floordiv;
+    mod = fmod(a, b);
+    if (!b) {
+        *modulus = mod;
+        return a / b;
+    }
+    div = (a - mod) / b;
+    if (mod) {
+        if ((b < 0) != (mod < 0)) {
+            mod += b;
+            div -= 1.0;
+        }
+    }
+    else {
+        mod = copysign(0.0, b);
+    }
+    if (div) {
+        floordiv = floor(div);
+        if (div - floordiv > 0.5) {
+            floordiv += 1.0;
+        }
+    }
+    else {
+        floordiv = copysign(0.0, a / b);
+    }
+    *modulus = mod;
+    return floordiv;
+}
+
+/* The volatile operands and result keep the compiler from moving, folding or sharing the operation
+ * across the flag clear and test (the function is also noinline). */
+static TP_NOINLINE int tp_np_f64(int op, double a, double b, double *out)
+{
+    volatile double va = a, vb = b, vr = 0.0, vd = 0.0;
+    double mod;
+    int flags, force = 0, clear = 0;
+    feclearexcept(FE_ALL_EXCEPT);
+    switch (op) {
+    case TP_NP_ADD: vr = va + vb; break;
+    case TP_NP_SUB: vr = va - vb; break;
+    case TP_NP_MUL: vr = va * vb; break;
+    case TP_NP_DIV: vr = va / vb; break;
+    case TP_NP_FLOORDIV:
+        if (!vb) {
+            feclearexcept(FE_ALL_EXCEPT);
+            vr = va / vb;
+            if (!va) {
+                force = TP_NP_INVALID;
+            }
+        }
+        else {
+            vr = tp_np_divmod(va, vb, &mod);
+        }
+        break;
+    default:
+        if (!vb) {
+            vr = fmod(va, vb);
+            clear = 1;
+        }
+        else {
+            vd = tp_np_divmod(va, vb, &mod);     /* numpy computes the quotient too: its flags count */
+            (void)vd;
+            vr = mod;
+        }
+        break;
+    }
+    flags = 0;
+    if (!clear) {
+        int f = fetestexcept(FE_DIVBYZERO | FE_OVERFLOW | FE_INVALID);
+        flags = ((f & FE_DIVBYZERO) ? TP_NP_DIVZERO : 0) | ((f & FE_OVERFLOW) ? TP_NP_OVERFLOW : 0)
+                | ((f & FE_INVALID) ? TP_NP_INVALID : 0);
+    }
+    *out = vr;
+    return flags | force;
+}
+
+/* numpy's RuntimeWarning for each bit of `flags`, in numpy's order, attributed to `file`:`line` (the
+ * operation's own source position, as the interpreter's frame would give). 0, or -1 with the
+ * exception set (a warning turned into an error). `op` is numpy's scalar name ("add", "divide"...). */
+static TP_NOINLINE int tp_np_warn(PyObject *tp_module, int flags, const char *op, PyObject *file, int line)
+{
+    PyObject *dict = PyModule_GetDict(tp_module);          /* borrowed */
+    PyObject *reg = NULL, *modname = NULL;
+    static const struct { int bit; const char *text; } kinds[3] = {
+        {TP_NP_DIVZERO, "divide by zero"}, {TP_NP_OVERFLOW, "overflow"}, {TP_NP_INVALID, "invalid value"}};
+    int k, rc = 0;
+    if (PyDict_GetItemStringRef(dict, "__warningregistry__", &reg) < 0) {
+        return -1;
+    }
+    if (reg == NULL) {
+        PyObject *fresh = PyDict_New(), *key = PyUnicode_FromString("__warningregistry__"), *got = NULL;
+        int rr = -1;
+        if (fresh != NULL && key != NULL) {
+            rr = PyDict_SetDefaultRef(dict, key, fresh, &got);
+        }
+        Py_XDECREF(fresh);
+        Py_XDECREF(key);
+        if (rr < 0) {
+            return -1;
+        }
+        reg = got;
+    }
+    if (PyDict_GetItemStringRef(dict, "__name__", &modname) < 0) {
+        Py_DECREF(reg);
+        return -1;
+    }
+    for (k = 0; k < 3 && rc == 0; k++) {
+        if (flags & kinds[k].bit) {
+            PyObject *msg = PyUnicode_FromFormat("%s encountered in scalar %s", kinds[k].text, op);
+            if (msg == NULL) {
+                rc = -1;
+                break;
+            }
+            rc = PyErr_WarnExplicitObject(PyExc_RuntimeWarning, msg, file, line, modname, reg);
+            Py_DECREF(msg);
+        }
+    }
+    Py_XDECREF(modname);
+    Py_DECREF(reg);
+    return rc;
+}
+
+/* int64 on the numpy path (pure functions): 0 ok, 1 deopt (the interpreter warns). */
+static inline int tp_np_truediv_i64(int64_t a, int64_t b, double *out)
+{
+    if (b == 0) {
+        return 1;
+    }
+    *out = (double)a / (double)b;          /* numpy converts both to double; CPython rounds exactly */
+    return 0;
+}
+
+static inline int tp_np_floordiv_i64(int64_t a, int64_t b, int64_t *out)
+{
+    if (b == 0) {
+        return 1;
+    }
+    return tp_floordiv_i64(a, b, out);     /* INT64_MIN // -1 is 1 (numpy's overflow warning) */
+}
+
+static inline int tp_np_mod_i64(int64_t a, int64_t b, int64_t *out)
+{
+    if (b == 0) {
+        return 1;
+    }
+    return tp_mod_i64(a, b, out);
+}
+
 /* ------------------------------------------------------------------------------------------
  * Comparisons: exact int/float comparison (float_richcompare), never via (double)int.
  * ------------------------------------------------------------------------------------------ */