From db708361ea828659ac6bee7a8acec159ef3352c0 Mon Sep 17 00:00:00 2001 From: BREW <81044362+b-re-w@users.noreply.github.com> Date: Sun, 25 Oct 2026 06:04:35 +0900 Subject: [PATCH] Feat: numpy warnings, division and impure element arithmetic (#220 step 2) --- docs/SPEC.md | 2 +- .../src/main/python/typedpython/cgen.py | 29 +- .../src/main/python/typedpython/frontend.py | 7 + .../src/main/python/typedpython/ir.py | 6 + .../python/typedpython/runtime/tp_runtime.h | 170 +++++++++ .../src/main/python/typedpython/taint.py | 30 +- .../src/main/python/typedpython/verify.py | 18 +- .../src/test/python/test_buffer_numpy.py | 338 +++++++++++++++++- 8 files changed, 569 insertions(+), 31 deletions(-) diff --git a/docs/SPEC.md b/docs/SPEC.md index ab7652e5..358cc8df 100644 --- a/docs/SPEC.md +++ b/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` | diff --git a/python-multiplatform-ksp/src/main/python/typedpython/cgen.py b/python-multiplatform-ksp/src/main/python/typedpython/cgen.py index 7f3bf55c..22cd4082 100644 --- a/python-multiplatform-ksp/src/main/python/typedpython/cgen.py +++ b/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 diff --git a/python-multiplatform-ksp/src/main/python/typedpython/frontend.py b/python-multiplatform-ksp/src/main/python/typedpython/frontend.py index d34fc22b..a61478d8 100644 --- a/python-multiplatform-ksp/src/main/python/typedpython/frontend.py +++ b/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]} " diff --git a/python-multiplatform-ksp/src/main/python/typedpython/ir.py b/python-multiplatform-ksp/src/main/python/typedpython/ir.py index f337e09f..8705b8fe 100644 --- a/python-multiplatform-ksp/src/main/python/typedpython/ir.py +++ b/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): diff --git a/python-multiplatform-ksp/src/main/python/typedpython/runtime/tp_runtime.h b/python-multiplatform-ksp/src/main/python/typedpython/runtime/tp_runtime.h index 22f7d580..25aaea17 100644 --- a/python-multiplatform-ksp/src/main/python/typedpython/runtime/tp_runtime.h +++ b/python-multiplatform-ksp/src/main/python/typedpython/runtime/tp_runtime.h @@ -14,6 +14,7 @@ #include #include #include +#include #include #include #include @@ -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__", ®) < 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. * ------------------------------------------------------------------------------------------ */ diff --git a/python-multiplatform-ksp/src/main/python/typedpython/taint.py b/python-multiplatform-ksp/src/main/python/typedpython/taint.py index 581dfaa1..489c7868 100644 --- a/python-multiplatform-ksp/src/main/python/typedpython/taint.py +++ b/python-multiplatform-ksp/src/main/python/typedpython/taint.py @@ -10,9 +10,9 @@ Conservative rules (decision 2026-10-25): whatever the pass cannot reproduce exactly is refused (`Refuse`, the function stays interpreted): - * `/`, `//`, `%` with a tainted operand (numpy divides without raising and warns), - * `+ - *` on F64 with a tainted operand in an impure function (numpy warns on overflow / invalid - value; a pure function deopts on a non-finite result instead, `BinOp.np`), + (step 2: `/ // %` and F64 `+ - *` on a tainted operand are no longer refused; they are `BinOp.np` + operations with numpy's results, a deopt on a warning in a pure function and the compiled code's own + RuntimeWarning in an impure one, see ir.BinOp.) * a tainted value stored into a list or local array, passed to a compiled callee, or used as the fill of a local array (the object stored or passed would have to be the numpy scalar), * a comparison of an I64 with an F64 where one side is tainted (numpy converts the int to double, @@ -69,13 +69,6 @@ def taint(self, e: ir.Expr) -> ir.Expr | None: return None if e.explicit else self.taint(e.operand) if isinstance(e, ir.BinOp): t = _or(self.taint(e.left), self.taint(e.right)) - if t is not None and e.op in (ir.BinOpKind.TRUEDIV, ir.BinOpKind.FLOORDIV, - ir.BinOpKind.MOD): - raise Refuse(f"`{e.op.value}` of a value derived from a numpy element (numpy divides " - "without raising and warns; not reproduced yet)") - if t is not None and e.type is Type.F64 and not self.pure: - raise Refuse(f"`{e.op.value}` of a value derived from a numpy element in a function " - "with effects (numpy would warn on overflow; not reproduced yet)") return t if isinstance(e, ir.UnaryOp): return None if e.op is ir.UnaryOpKind.NOT else self.taint(e.operand) @@ -142,14 +135,17 @@ def rewrite(self, x): return stmt if isinstance(x, ir.Box) and x.operand.type in ir.SCALARS: return dataclasses.replace(x, operand=self.rewrite(x.operand), np=self.taint(x.operand)) - if isinstance(x, ir.BinOp) and x.type in ir.SCALARS: - self.taint(x) # refuses division / impure arithmetic wherever it appears - if isinstance(x, ir.BinOp) and x.type is Type.F64 and x.op in ( - ir.BinOpKind.ADD, ir.BinOpKind.SUB, ir.BinOpKind.MUL): + if isinstance(x, ir.BinOp) and x.type in ir.SCALARS and x.left.type in (Type.F64, Type.I64): tainted = _or(self.taint(x.left), self.taint(x.right)) is not None - self.deopts = self.deopts or tainted - return dataclasses.replace(x, left=self.rewrite(x.left), right=self.rewrite(x.right), - np=tainted) + division = x.op in (ir.BinOpKind.TRUEDIV, ir.BinOpKind.FLOORDIV, ir.BinOpKind.MOD) + if x.left.type is Type.F64: + np_op = tainted + else: # I64: + - * keep the checked deopt; the divisions are numpy's + np_op = tainted and division and not x.proven + if np_op: + self.deopts = self.deopts or self.pure + return dataclasses.replace(x, left=self.rewrite(x.left), right=self.rewrite(x.right), + np=True) if isinstance(x, ir.StoreIndex) and x.array not in self.numpy: self.plain(x.value, "stored into a list or local array") if isinstance(x, ir.NewArray) and x.fill is not None: diff --git a/python-multiplatform-ksp/src/main/python/typedpython/verify.py b/python-multiplatform-ksp/src/main/python/typedpython/verify.py index db8487d3..0b5f90b1 100644 --- a/python-multiplatform-ksp/src/main/python/typedpython/verify.py +++ b/python-multiplatform-ksp/src/main/python/typedpython/verify.py @@ -307,6 +307,8 @@ class Diagnostic: _FULL = (_I64_MIN, _I64_MAX) _LEN_MAX = 2 ** 62 _PROVABLE_OPS = frozenset({BinOpKind.ADD, BinOpKind.SUB, BinOpKind.MUL}) +_NP_F64_OPS = frozenset({BinOpKind.ADD, BinOpKind.SUB, BinOpKind.MUL, BinOpKind.TRUEDIV, + BinOpKind.FLOORDIV, BinOpKind.MOD}) _POW2_OPS = frozenset({BinOpKind.FLOORDIV, BinOpKind.MOD}) # `proven` means shift / mask _NEGATE = {CompareKind.EQ: CompareKind.NE, CompareKind.NE: CompareKind.EQ, CompareKind.LT: CompareKind.GE, CompareKind.GE: CompareKind.LT, @@ -774,11 +776,17 @@ def e_BinOp(self, e: BinOp, state: _State) -> ir.Expr: if type(e.np) is not bool: self.error(TYPE, f"BinOp `np` {e.np!r} is not a bool") elif e.np: - if not (e.type is Type.F64 and lt is Type.F64 and rt is Type.F64 - and e.op in (BinOpKind.ADD, BinOpKind.SUB, BinOpKind.MUL)): - self.error(TYPE, "BinOp `np` is set on something other than an f64 + - * (#220)") - else: - self.can_deopt = True # a non-finite result: the interpreter redoes the call + f64_op = e.type is Type.F64 and lt is Type.F64 and rt is Type.F64 + i64_op = (lt is Type.I64 and rt is Type.I64 and not e.proven and e.op in ( + BinOpKind.TRUEDIV, BinOpKind.FLOORDIV, BinOpKind.MOD) + and e.type is (Type.F64 if e.op is BinOpKind.TRUEDIV else Type.I64)) + if not (f64_op or i64_op): + self.error(TYPE, "BinOp `np` is set on something other than an f64 operation or an " + "unproven i64 / // % (#220)") + elif f64_op and e.op not in _NP_F64_OPS: + self.error(TYPE, f"BinOp `np` on an f64 {e.op.value}") + elif i64_op or self.fn.pure is True: + self.can_deopt = True # a warning / int64 division by zero: the interpreter redoes shift = e.op in (BinOpKind.LSHIFT, BinOpKind.RSHIFT) if lt is Type.INT: # SPEC N-13: both operands INT, result INT, never deopts (it promotes in place) diff --git a/python-multiplatform-ksp/src/test/python/test_buffer_numpy.py b/python-multiplatform-ksp/src/test/python/test_buffer_numpy.py index 30f3e681..b94a9f58 100644 --- a/python-multiplatform-ksp/src/test/python/test_buffer_numpy.py +++ b/python-multiplatform-ksp/src/test/python/test_buffer_numpy.py @@ -208,14 +208,168 @@ def nd_int_impure(x: npt.NDArray[np.int64]) -> object: def nd_scale(x: npt.NDArray[np.float64], k: float) -> None: for i in range(len(x)): x[i] = x[i] * k + + def p_div(x: npt.NDArray[np.float64]) -> object: + return x[0] / x[1] + + + def p_fdiv(x: npt.NDArray[np.float64]) -> object: + return x[0] // x[1] + + + def p_mod(x: npt.NDArray[np.float64]) -> object: + return x[0] % x[1] + + + def p_add(x: npt.NDArray[np.float64]) -> object: + return x[0] + x[1] + + + def p_sub(x: npt.NDArray[np.float64]) -> object: + return x[0] - x[1] + + + def p_mul(x: npt.NDArray[np.float64]) -> object: + return x[0] * x[1] + + + def p_pdiv(x: npt.NDArray[np.float64], k: float) -> object: + return k / x[0] + + + def p_idiv(x: npt.NDArray[np.int64]) -> object: + return x[0] / x[1] + + + def p_ifdiv(x: npt.NDArray[np.int64]) -> object: + return x[0] // x[1] + + + def p_imod(x: npt.NDArray[np.int64]) -> object: + return x[0] % x[1] + + + def p_imixed(x: npt.NDArray[np.int64]) -> object: + return x[0] // 3 + x[1] % 4 + + + def i_div(x: npt.NDArray[np.float64]) -> None: + x[2] = x[0] / x[1] + + + def i_fdiv(x: npt.NDArray[np.float64]) -> None: + x[2] = x[0] // x[1] + + + def i_mod(x: npt.NDArray[np.float64]) -> None: + x[2] = x[0] % x[1] + + + def i_add(x: npt.NDArray[np.float64]) -> None: + x[2] = x[0] + x[1] + + + def i_sub(x: npt.NDArray[np.float64]) -> None: + x[2] = x[0] - x[1] + + + def i_mul(x: npt.NDArray[np.float64]) -> None: + x[2] = x[0] * x[1] + + + def i_plain(x: npt.NDArray[np.float64], k: float) -> None: + x[2] = k / x[1] + x[3] = k * 2.0 + + + def i_multiline(x: npt.NDArray[np.float64]) -> None: + x[2] = (x[0] * + x[1]) + x[3] = (x[0] + / x[1] + + x[0] * x[1]) + + + def i_aug(x: npt.NDArray[np.float64]) -> None: + x[0] *= x[1] + x[1] /= x[2] + x[2] += x[3] + + + def i_loop(x: npt.NDArray[np.float64], y: npt.NDArray[np.float64]) -> float: + s = 0.0 + for i in range(len(x)): + y[i] = x[i] / y[i] + s += y[i] * x[i] + return s + + + def i_mixed(x: npt.NDArray[np.float64], n: int) -> None: + x[0] = n / x[1] + x[2] = x[3] * n + + + def i_ret(x: npt.NDArray[np.float64]) -> object: + x[2] = x[0] * x[1] + return x[0] / x[1] + + + def i_ifdiv(x: npt.NDArray[np.int64]) -> object: + x[2] = x[0] // x[1] + return x[2] + + + def nd_advance(dt: float, steps: int, x: npt.NDArray[np.float64], y: npt.NDArray[np.float64], + z: npt.NDArray[np.float64], vx: npt.NDArray[np.float64], + vy: npt.NDArray[np.float64], vz: npt.NDArray[np.float64], + mass: npt.NDArray[np.float64]) -> None: + for _ in range(steps): + for i in range(5): + for j in range(i + 1, 5): + dx: float = x[i] - x[j] + dy: float = y[i] - y[j] + dz: float = z[i] - z[j] + d2: float = dx * dx + dy * dy + dz * dz + mag: float = dt / (d2 * math.sqrt(d2)) + vx[i] -= dx * mass[j] * mag + vy[i] -= dy * mass[j] * mag + vz[i] -= dz * mass[j] * mag + vx[j] += dx * mass[i] * mag + vy[j] += dy * mass[i] * mag + vz[j] += dz * mass[i] * mag + for i in range(5): + x[i] += dt * vx[i] + y[i] += dt * vy[i] + z[i] += dt * vz[i] + + + def nd_energy_full(x: npt.NDArray[np.float64], y: npt.NDArray[np.float64], + z: npt.NDArray[np.float64], vx: npt.NDArray[np.float64], + vy: npt.NDArray[np.float64], vz: npt.NDArray[np.float64], + mass: npt.NDArray[np.float64]) -> float: + e: float = 0.0 + for i in range(5): + e += 0.5 * mass[i] * (vx[i] * vx[i] + vy[i] * vy[i] + vz[i] * vz[i]) + for j in range(i + 1, 5): + dx: float = x[i] - x[j] + dy: float = y[i] - y[j] + dz: float = z[i] - z[j] + e -= mass[i] * mass[j] / math.sqrt(dx * dx + dy * dy + dz * dz) + return e ''') COMPILED = { "nd_first", "nd_first_f", "nd_first_i", "nd_sum", "nd_sum_from_first", "nd_isum", "nd_gt", "nd_flag", "nd_plain", "nd_mixed", "nd_sqrt", "nd_neg", "nd_len", "nd_set", "nd_copy", "nd_iset", "nd_get", "nd_param", "nd_branch", "nd_tuple", "nd_overflow", "nd_invalid", "nd_iover", "nd_energy", "nd_copy7", + # step 2 + "nd_div", "nd_idiv", "nd_scale", "p_div", "p_fdiv", "p_mod", "p_add", "p_sub", "p_mul", "p_pdiv", + "p_idiv", "p_ifdiv", "p_imod", "p_imixed", "i_div", "i_fdiv", "i_mod", "i_add", "i_sub", "i_mul", + "i_plain", "i_multiline", "i_aug", "i_loop", "i_mixed", "i_ret", "nd_advance", "nd_energy_full", } -NOT_COMPILED = {"nd_div", "nd_idiv", "nd_scale", "nd_to_list", "nd_mixed_cmp", "nd_abs", "nd_int_impure"} # nd_abs calls user code # could warn / divide, or store a possibly warning op +IMPURE = {"i_div", "i_fdiv", "i_mod", "i_add", "i_sub", "i_mul", "i_plain", "i_multiline", "i_aug", + "i_loop", "i_mixed", "i_ret", "nd_advance", "nd_scale"} +NOT_COMPILED = {"i_ifdiv", "nd_to_list", "nd_mixed_cmp", "nd_abs", "nd_int_impure"} # nd_abs calls user code # could warn / divide, or store a possibly warning op F = [1.5, -2.25, 3.0, 4.5] I = [3, -5, 7, 11] @@ -238,7 +392,7 @@ def run(fn, args): r = ("ok", fn(*args)) except Exception as e: # noqa: BLE001 - type and text are the observable r = ("raise", type(e), str(e)) - return r, [(x.category, str(x.message), x.lineno) for x in w] + return r, [(x.category, str(x.message), x.filename, x.lineno) for x in w] def same(a, b): @@ -405,13 +559,14 @@ def test_int_overflow_warning_comes_from_the_interpreter(built): differential(m, "nd_iover", lambda: (i64(),)) -def test_division_and_impure_arithmetic_stay_on_the_object_path(built): +def test_division_matches_numpy(built): _, m, _ = built for make in (lambda: (f64([1.0, 0.0]),), lambda: (f64([-1.0, 0.0]),), lambda: (f64([0.0, 0.0]),), lambda: (f64(),)): - differential(m, "nd_div", make) - c, w = differential(m, "nd_idiv", lambda: (i64([5, 0]),)) + differential(m, "nd_div", make, expect_deopt=make is not None and make()[0][1] == 0.0) + c, w = differential(m, "nd_idiv", lambda: (i64([5, 0]),), expect_deopt=True) assert [x[0] for x in w] == [RuntimeWarning] + differential(m, "nd_idiv", lambda: (i64([5, 2]),)) differential(m, "nd_scale", lambda: (f64([BIG, 1.0]), 10.0)) differential(m, "nd_scale", lambda: (f64(), 2.0)) @@ -478,3 +633,176 @@ def _nodes(x): yield x for f in dataclasses.fields(x): yield from _nodes(getattr(x, f.name)) + + +# --- step 2: numpy division and impure arithmetic --------------------------------------------------------- + +SPECIAL = [0.0, -0.0, 1.0, -1.0, 2.5, -2.5, 1e308, -1e308, 1e-310, 5e-324, float("inf"), + float("-inf"), float("nan"), 0.3] +IMIN = -2 ** 63 +SPECIAL_I = [0, 1, -1, 7, -7, IMAX, IMIN] + + +def grid(m, name, pairs, make, pure): + """Every pair, compiled against interpreted: result, TYPE, warnings with file and line.""" + seen = set() + for a, b in pairs: + ca, ia = make(a, b), make(a, b) + before = m.__typedpython_deopts__ + c, cw = run(getattr(m, name), ca) + deopted = m.__typedpython_deopts__ - before + i, iw = run(m.__typedpython_interpreted__[name], ia) + assert c[0] == i[0], (name, a, b, c, i) + if c[0] == "ok" and c[1] is not None: + assert type(c[1]) is type(i[1]), (name, a, b) + assert same(c[1], i[1]), (name, a, b, c, i) + assert cw == iw, (name, a, b, cw, iw) + for x, y in zip(ca, ia): + if isinstance(x, np.ndarray): + assert x.dtype == y.dtype + assert all(same(p, q) and np.signbit(p) == np.signbit(q) + for p, q in zip(x.tolist(), y.tolist()) if p == p), (name, a, b, x, y) + assert all((p != p) == (q != q) for p, q in zip(x.tolist(), y.tolist())) + assert deopted == (1 if pure and cw else 0), (name, a, b, deopted, cw) + seen.update(str(w[1]) for w in cw) + return seen + + +def fpairs(): + return [(a, b) for a in SPECIAL for b in SPECIAL] + + +@pytest.mark.parametrize("name", ["p_div", "p_fdiv", "p_mod", "p_add", "p_sub", "p_mul"]) +def test_pure_float_ops_match_numpy_over_special_values(built, name): + _, m, _ = built + seen = grid(m, name, fpairs(), lambda a, b: (f64([a, b]),), pure=True) + if name == "p_div": + assert {"divide by zero encountered in scalar divide", "invalid value encountered in scalar divide", + "overflow encountered in scalar divide"} <= seen + + +@pytest.mark.parametrize("name", ["i_div", "i_fdiv", "i_mod", "i_add", "i_sub", "i_mul"]) +def test_impure_float_ops_match_numpy_and_warn_from_the_compiled_code(built, name): + _, m, _ = built + seen = grid(m, name, fpairs(), lambda a, b: (f64([a, b, 0.0]),), pure=False) + assert seen, "numpy warned for no pair: the grid tests nothing" + + +def test_impure_warning_has_the_operations_file_and_line(built): + _, m, _ = built + _, w = run(m.i_multiline, (f64([BIG, 10.0, 0.0, 0.0]),)) + lines = [x[3] for x in w] + assert len(set(lines)) >= 2, w + grid(m, "i_multiline", [(a, b) for a in (1.0, 0.0, BIG, float("inf")) for b in (0.0, 1.0, BIG)], + lambda a, b: (f64([a, b, 0.0, 0.0]),), pure=False) + grid(m, "i_aug", [(a, b) for a in (1.0, 0.0, BIG) for b in (0.0, 1.0, BIG, float("nan"))], + lambda a, b: (f64([a, b, 0.0, BIG]),), pure=False) + + +def test_impure_error_filter_raises_and_leaves_the_same_state(built): + _, m, _ = built + for fn in (m.i_div, m.__typedpython_interpreted__["i_div"]): + x = f64([1.0, 0.0, 7.0]) + with warnings.catch_warnings(): + warnings.simplefilter("error") + with pytest.raises(RuntimeWarning, match="divide by zero encountered in scalar divide"): + fn(x) + assert x.tolist() == [1.0, 0.0, 7.0] + + +def test_impure_default_filter_warns_once_per_location(built): + _, m, _ = built + out = [] + for fn in (m.i_div, m.__typedpython_interpreted__["i_div"]): + with warnings.catch_warnings(record=True) as w: + warnings.simplefilter("default") + m.__dict__.pop("__warningregistry__", None) + for _ in range(3): + fn(f64([1.0, 0.0, 0.0])) + out.append([(x.filename, x.lineno, str(x.message)) for x in w]) + assert out[0] == out[1] and len(out[0]) == 1, out + + +def test_impure_mixed_plain_and_numpy_operands(built): + _, m, _ = built + grid(m, "i_plain", [(a, b) for a in (1.0, 0.0, BIG) for b in (0.0, 2.0, -0.0)], + lambda a, b: (f64([0.0, b, 0.0, 0.0]), a), pure=False) + grid(m, "i_mixed", [(a, b) for a in (1, 0, 2 ** 62) for b in (0.0, 2.0, 1e-310)], + lambda a, b: (f64([0.0, b, 0.0, BIG]), a), pure=False) + grid(m, "p_pdiv", [(a, b) for a in (1.0, 0.0, BIG) for b in (0.0, 2.0, 1e-310)], + lambda a, b: (f64([b, 0.0]), a), pure=True) + for name in ("i_loop", "i_ret"): + for xs, ys in (([1.0, 2.0, 0.0], [0.0, 4.0, 2.0]), ([BIG, 1.0, 0.0], [1e-5, 0.0, 0.0])): + args = lambda: ((f64(xs), f64(ys)) if name == "i_loop" else (f64(xs),)) # noqa: E731 + differential(m, name, args) + + +@pytest.mark.parametrize("name", ["p_idiv", "p_ifdiv", "p_imod"]) +def test_pure_int_division_matches_numpy(built, name): + _, m, _ = built + pairs = [(a, b) for a in SPECIAL_I for b in SPECIAL_I] + grid(m, name, pairs, lambda a, b: (i64([a, b]),), pure=True) + + +def test_int_division_by_a_constant_and_mixed_ops(built): + _, m, _ = built + grid(m, "p_imixed", [(a, b) for a in SPECIAL_I for b in SPECIAL_I], + lambda a, b: (i64([a, b]),), pure=True) + differential(m, "p_idiv", lambda: (i64([2 ** 60 + 1, 3]),)) # numpy divides as doubles + + +def test_the_int_divide_by_zero_in_an_impure_function_stays_on_the_object_path(built): + proved, m, _ = built + assert "i_ifdiv" in proved.skipped + differential(m, "i_ifdiv", lambda: (i64([5, 0, 0]),)) + + +def test_a_disabled_warning_path_would_differ(built): + """The check that can fail: the compiled impure i_div must emit the warning itself.""" + _, m, _ = built + _, w = run(m.i_div, (f64([1.0, 0.0, 0.0]),)) + assert [(x[0], x[1]) for x in w] == [(RuntimeWarning, "divide by zero encountered in scalar divide")] + + +BODIES = { + "x": [0.0, 4.84143144246472090, -4.15980459285970104e-1, 1.2 , 8.34336671824457987], + "y": [0.0, -1.16032004402742839, 7.69901118419740425e-1, -1.2, 4.12479856412430479], + "z": [0.0, -1.03622044471123109e-1, -8.5e-1, 0.5, -4.0e-1], + "vx": [0.0, 1.66007664274403694e-3 * 365.24, -2.76742510726862411e-3 * 365.24, 0.1, 0.2], + "vy": [0.0, 7.69901118419740425e-3 * 365.24, 4.99852801234917238e-3 * 365.24, 0.3, 0.1], + "vz": [0.0, -6.90460016972063023e-5 * 365.24, 2.30417297573763929e-5 * 365.24, -0.2, 0.05], + "mass": [39.47841760435743, 0.03769367487038949, 0.011286326131968767, 0.0017237240570597112, 0.0020336868699246304], +} + + +def bodies(): + return [np.array(BODIES[k], dtype=np.float64) for k in ("x", "y", "z", "vx", "vy", "vz", "mass")] + + +def test_nbody_advance_and_energy_compile_and_match(built): + proved, m, _ = built + by = {f.name: f for f in proved.functions} + assert "nd_advance" in by and "nd_energy_full" in by, proved.skipped + assert by["nd_advance"].pure is False and by["nd_energy_full"].pure is True + ref = m.__typedpython_interpreted__ + a, b = bodies(), bodies() + for steps in (1, 25): + before = m.__typedpython_deopts__ + with warnings.catch_warnings(record=True) as w: + warnings.simplefilter("always") + m.nd_advance(0.01, steps, *a) + ref["nd_advance"](0.01, steps, *b) + assert m.__typedpython_deopts__ == before + assert w == [] + for p, q in zip(a, b): + assert p.tolist() == q.tolist() + e1, e2 = m.nd_energy_full(*a), ref["nd_energy_full"](*b) + assert type(e1) is type(e2) is np.float64 and e1 == e2 + # coincident bodies: the interpreter's warnings and results, from the compiled advance + a, b = bodies(), bodies() + for arrs in (a, b): + arrs[0][2] = arrs[0][1]; arrs[1][2] = arrs[1][1]; arrs[2][2] = arrs[2][1] + ca, cw = run(m.nd_advance, (0.01, 1, *a)) + ia, iw = run(ref["nd_advance"], (0.01, 1, *b)) + assert cw == iw and cw, (cw, iw) + assert all(same(p, q) for u, v in zip(a, b) for p, q in zip(u.tolist(), v.tolist()))