|
14 | 14 | #include <Python.h> |
15 | 15 | #include <stdint.h> |
16 | 16 | #include <math.h> |
| 17 | +#include <fenv.h> |
17 | 18 | #include <errno.h> |
18 | 19 | #include <string.h> |
19 | 20 | #include <stdlib.h> |
@@ -398,6 +399,175 @@ static inline int tp_mod_f64(double a, double b, double *out) |
398 | 399 | return 0; |
399 | 400 | } |
400 | 401 |
|
| 402 | +/* ------------------------------------------------------------------------------------------ |
| 403 | + * numpy scalar arithmetic (#220, step 2). An operation on a value derived from a numpy element is |
| 404 | + * numpy's, not CPython's: no ZeroDivisionError, inf / nan results, int64 division by zero, and a |
| 405 | + * RuntimeWarning per floating-point exception. tp_np_f64 returns the exceptions that numpy reports |
| 406 | + * (TP_NP_* bits). A pure function deopts on any bit (the interpreter then warns); an impure one calls |
| 407 | + * tp_np_warn. The algorithms are numpy's (npymath: npy_divmod / npy_floor_divide / npy_remainder). |
| 408 | + * ------------------------------------------------------------------------------------------ */ |
| 409 | + |
| 410 | +#define TP_NP_DIVZERO 1 |
| 411 | +#define TP_NP_OVERFLOW 2 |
| 412 | +#define TP_NP_INVALID 4 |
| 413 | +enum { TP_NP_ADD, TP_NP_SUB, TP_NP_MUL, TP_NP_DIV, TP_NP_FLOORDIV, TP_NP_MOD }; |
| 414 | + |
| 415 | +#if defined(__GNUC__) || defined(__clang__) |
| 416 | +#define TP_NOINLINE __attribute__((noinline, unused)) |
| 417 | +#else |
| 418 | +#define TP_NOINLINE |
| 419 | +#endif |
| 420 | + |
| 421 | +static inline double tp_np_divmod(double a, double b, double *modulus) |
| 422 | +{ |
| 423 | + double div, mod, floordiv; |
| 424 | + mod = fmod(a, b); |
| 425 | + if (!b) { |
| 426 | + *modulus = mod; |
| 427 | + return a / b; |
| 428 | + } |
| 429 | + div = (a - mod) / b; |
| 430 | + if (mod) { |
| 431 | + if ((b < 0) != (mod < 0)) { |
| 432 | + mod += b; |
| 433 | + div -= 1.0; |
| 434 | + } |
| 435 | + } |
| 436 | + else { |
| 437 | + mod = copysign(0.0, b); |
| 438 | + } |
| 439 | + if (div) { |
| 440 | + floordiv = floor(div); |
| 441 | + if (div - floordiv > 0.5) { |
| 442 | + floordiv += 1.0; |
| 443 | + } |
| 444 | + } |
| 445 | + else { |
| 446 | + floordiv = copysign(0.0, a / b); |
| 447 | + } |
| 448 | + *modulus = mod; |
| 449 | + return floordiv; |
| 450 | +} |
| 451 | + |
| 452 | +/* The volatile operands and result keep the compiler from moving, folding or sharing the operation |
| 453 | + * across the flag clear and test (the function is also noinline). */ |
| 454 | +static TP_NOINLINE int tp_np_f64(int op, double a, double b, double *out) |
| 455 | +{ |
| 456 | + volatile double va = a, vb = b, vr = 0.0, vd = 0.0; |
| 457 | + double mod; |
| 458 | + int flags, force = 0, clear = 0; |
| 459 | + feclearexcept(FE_ALL_EXCEPT); |
| 460 | + switch (op) { |
| 461 | + case TP_NP_ADD: vr = va + vb; break; |
| 462 | + case TP_NP_SUB: vr = va - vb; break; |
| 463 | + case TP_NP_MUL: vr = va * vb; break; |
| 464 | + case TP_NP_DIV: vr = va / vb; break; |
| 465 | + case TP_NP_FLOORDIV: |
| 466 | + if (!vb) { |
| 467 | + feclearexcept(FE_ALL_EXCEPT); |
| 468 | + vr = va / vb; |
| 469 | + if (!va) { |
| 470 | + force = TP_NP_INVALID; |
| 471 | + } |
| 472 | + } |
| 473 | + else { |
| 474 | + vr = tp_np_divmod(va, vb, &mod); |
| 475 | + } |
| 476 | + break; |
| 477 | + default: |
| 478 | + if (!vb) { |
| 479 | + vr = fmod(va, vb); |
| 480 | + clear = 1; |
| 481 | + } |
| 482 | + else { |
| 483 | + vd = tp_np_divmod(va, vb, &mod); /* numpy computes the quotient too: its flags count */ |
| 484 | + (void)vd; |
| 485 | + vr = mod; |
| 486 | + } |
| 487 | + break; |
| 488 | + } |
| 489 | + flags = 0; |
| 490 | + if (!clear) { |
| 491 | + int f = fetestexcept(FE_DIVBYZERO | FE_OVERFLOW | FE_INVALID); |
| 492 | + flags = ((f & FE_DIVBYZERO) ? TP_NP_DIVZERO : 0) | ((f & FE_OVERFLOW) ? TP_NP_OVERFLOW : 0) |
| 493 | + | ((f & FE_INVALID) ? TP_NP_INVALID : 0); |
| 494 | + } |
| 495 | + *out = vr; |
| 496 | + return flags | force; |
| 497 | +} |
| 498 | + |
| 499 | +/* numpy's RuntimeWarning for each bit of `flags`, in numpy's order, attributed to `file`:`line` (the |
| 500 | + * operation's own source position, as the interpreter's frame would give). 0, or -1 with the |
| 501 | + * exception set (a warning turned into an error). `op` is numpy's scalar name ("add", "divide"...). */ |
| 502 | +static TP_NOINLINE int tp_np_warn(PyObject *tp_module, int flags, const char *op, PyObject *file, int line) |
| 503 | +{ |
| 504 | + PyObject *dict = PyModule_GetDict(tp_module); /* borrowed */ |
| 505 | + PyObject *reg = NULL, *modname = NULL; |
| 506 | + static const struct { int bit; const char *text; } kinds[3] = { |
| 507 | + {TP_NP_DIVZERO, "divide by zero"}, {TP_NP_OVERFLOW, "overflow"}, {TP_NP_INVALID, "invalid value"}}; |
| 508 | + int k, rc = 0; |
| 509 | + if (PyDict_GetItemStringRef(dict, "__warningregistry__", ®) < 0) { |
| 510 | + return -1; |
| 511 | + } |
| 512 | + if (reg == NULL) { |
| 513 | + PyObject *fresh = PyDict_New(), *key = PyUnicode_FromString("__warningregistry__"), *got = NULL; |
| 514 | + int rr = -1; |
| 515 | + if (fresh != NULL && key != NULL) { |
| 516 | + rr = PyDict_SetDefaultRef(dict, key, fresh, &got); |
| 517 | + } |
| 518 | + Py_XDECREF(fresh); |
| 519 | + Py_XDECREF(key); |
| 520 | + if (rr < 0) { |
| 521 | + return -1; |
| 522 | + } |
| 523 | + reg = got; |
| 524 | + } |
| 525 | + if (PyDict_GetItemStringRef(dict, "__name__", &modname) < 0) { |
| 526 | + Py_DECREF(reg); |
| 527 | + return -1; |
| 528 | + } |
| 529 | + for (k = 0; k < 3 && rc == 0; k++) { |
| 530 | + if (flags & kinds[k].bit) { |
| 531 | + PyObject *msg = PyUnicode_FromFormat("%s encountered in scalar %s", kinds[k].text, op); |
| 532 | + if (msg == NULL) { |
| 533 | + rc = -1; |
| 534 | + break; |
| 535 | + } |
| 536 | + rc = PyErr_WarnExplicitObject(PyExc_RuntimeWarning, msg, file, line, modname, reg); |
| 537 | + Py_DECREF(msg); |
| 538 | + } |
| 539 | + } |
| 540 | + Py_XDECREF(modname); |
| 541 | + Py_DECREF(reg); |
| 542 | + return rc; |
| 543 | +} |
| 544 | + |
| 545 | +/* int64 on the numpy path (pure functions): 0 ok, 1 deopt (the interpreter warns). */ |
| 546 | +static inline int tp_np_truediv_i64(int64_t a, int64_t b, double *out) |
| 547 | +{ |
| 548 | + if (b == 0) { |
| 549 | + return 1; |
| 550 | + } |
| 551 | + *out = (double)a / (double)b; /* numpy converts both to double; CPython rounds exactly */ |
| 552 | + return 0; |
| 553 | +} |
| 554 | + |
| 555 | +static inline int tp_np_floordiv_i64(int64_t a, int64_t b, int64_t *out) |
| 556 | +{ |
| 557 | + if (b == 0) { |
| 558 | + return 1; |
| 559 | + } |
| 560 | + return tp_floordiv_i64(a, b, out); /* INT64_MIN // -1 is 1 (numpy's overflow warning) */ |
| 561 | +} |
| 562 | + |
| 563 | +static inline int tp_np_mod_i64(int64_t a, int64_t b, int64_t *out) |
| 564 | +{ |
| 565 | + if (b == 0) { |
| 566 | + return 1; |
| 567 | + } |
| 568 | + return tp_mod_i64(a, b, out); |
| 569 | +} |
| 570 | + |
401 | 571 | /* ------------------------------------------------------------------------------------------ |
402 | 572 | * Comparisons: exact int/float comparison (float_richcompare), never via (double)int. |
403 | 573 | * ------------------------------------------------------------------------------------------ */ |
|
0 commit comments