René's URL Explorer Experiment

{"@context":"https://schema.org","@type":"DiscussionForumPosting","headline":"Invalid corner cases (resulting in nan+nanj) in _Py_c_prod()","articleBody":"# Bug report\r\n\r\n### Bug description:\r\n\r\nReproducer:\r\n\r\n```pycon\r\n\u003e\u003e\u003e z = 1e300+1j\r\n\u003e\u003e\u003e z*complex('(inf+infj)')  # should be (nan+infj)\r\n(nan+nanj)\r\n```\r\n\r\n\u003cdetails\u003e\r\n\r\n\u003csummary\u003ec.f. C code\u003c/summary\u003e\r\n\r\n```c\r\n#include \u003cstdio.h\u003e\r\n#include \u003ccomplex.h\u003e\r\n#include \u003cmath.h\u003e\r\nint main (void)\r\n{\r\n    complex double z = CMPLX(1e300, 1);\r\n    z = CMPLX(INFINITY, INFINITY)*z;\r\n    printf(\"%e %e\\n\", creal(z), cimag(z));\r\n    return 0;\r\n}\r\n```\r\n```\r\n$ cc a.c \u0026\u0026 ./a.out\r\n-nan inf\r\n```  \r\n\r\n\u003c/details\u003e\r\n\r\nThat quickly leads to problems in optimized algorithm for exponentiation (by squaring), that claims to be \"more accurate\":\r\n```pycon\r\n\u003e\u003e\u003e z**5\r\n(nan+nanj)\r\n\u003e\u003e\u003e z**5.0000001  # generic algorithm\r\nTraceback (most recent call last):\r\n  File \"\u003cpython-input-5\u003e\", line 1, in \u003cmodule\u003e\r\n    z**5.0000001\r\n    ~^^~~~~~~~~~\r\nOverflowError: complex exponentiation\r\n\u003e\u003e\u003e _testcapi._py_c_pow(z, 5)  # for integer exponent it also signals overflow\r\n((inf+infj), 34)\r\n```\r\n\r\nFollowing patch (mostly a literal translation of the ``_Cmultd()`` routine from the C11 Annex G.5.2) should solve the issue.\r\n\u003cdetails\u003e\r\n\r\n\u003csummary\u003ea patch\u003c/summary\u003e\r\n\r\n```diff\r\ndiff --git a/Objects/complexobject.c b/Objects/complexobject.c\r\nindex 59c84f1359..1d9707895c 100644\r\n--- a/Objects/complexobject.c\r\n+++ b/Objects/complexobject.c\r\n@@ -53,11 +53,48 @@ _Py_c_neg(Py_complex a)\r\n }\r\n\r\n Py_complex\r\n-_Py_c_prod(Py_complex a, Py_complex b)\r\n+_Py_c_prod(Py_complex z, Py_complex w)\r\n {\r\n     Py_complex r;\r\n-    r.real = a.real*b.real - a.imag*b.imag;\r\n-    r.imag = a.real*b.imag + a.imag*b.real;\r\n+    double a = z.real, b = z.imag, c = w.real, d = w.imag;\r\n+    double ac = a*c, bd = b*d, ad = a*d, bc = b*c;\r\n+\r\n+    r.real = ac - bd;\r\n+    r.imag = ad + bc;\r\n+\r\n+    if (isnan(r.real) \u0026\u0026 isnan(r.imag)) {\r\n+        /* Recover infinities that computed as (nan+nanj) */\r\n+        int recalc = 0;\r\n+        if (isinf(a) || isinf(b)) {  /* z is infinite */\r\n+            /* \"Box\" the infinity and change nans in the other factor to 0 */\r\n+            a = copysign(isinf(a) ? 1.0 : 0.0, a);\r\n+            b = copysign(isinf(b) ? 1.0 : 0.0, b);\r\n+            if (isnan(c)) c = copysign(0.0, c);\r\n+            if (isnan(d)) d = copysign(0.0, d);\r\n+            recalc = 1;\r\n+        }\r\n+        if (isinf(c) || isinf(d)) {  /* w is infinite */\r\n+            /* \"Box\" the infinity and change nans in the other factor to 0 */\r\n+            c = copysign(isinf(c) ? 1.0 : 0.0, c);\r\n+            d = copysign(isinf(d) ? 1.0 : 0.0, d);\r\n+            if (isnan(a)) a = copysign(0.0, a);\r\n+            if (isnan(b)) b = copysign(0.0, b);\r\n+            recalc = 1;\r\n+        }\r\n+        if (!recalc \u0026\u0026 (isinf(ac) || isinf(bd) || isinf(ad) || isinf(bc))) {\r\n+            /* Recover infinities from overflow by changing nans to 0 */\r\n+            if (isnan(a)) a = copysign(0.0, a);\r\n+            if (isnan(b)) b = copysign(0.0, b);\r\n+            if (isnan(c)) c = copysign(0.0, c);\r\n+            if (isnan(d)) d = copysign(0.0, d);\r\n+            recalc = 1;\r\n+        }\r\n+        if (recalc) {\r\n+            r.real = Py_INFINITY*(a*c - b*d);\r\n+            r.imag = Py_INFINITY*(a*d + b*c);\r\n+        }\r\n+    }\r\n+\r\n     return r;\r\n }\r\n\r\n```\r\n\r\nc.f. clang's version:\r\nhttps://github.com/llvm/llvm-project/blob/4973ad47181710d2a69292018cad7bc6f95a6c1a/libcxx/include/complex#L712-L792\r\n\u003c/details\u003e\r\n\r\nAlso, maybe ``_Py_c_prod()`` code should set errno on overflows, just as ``_Py_c_abs()``.  E.g. with above version we have:\r\n```pycon\r\n\u003e\u003e\u003e z**5\r\nTraceback (most recent call last):\r\n  File \"\u003cpython-input-1\u003e\", line 1, in \u003cmodule\u003e\r\n    z**5\r\n    ~^^~\r\nOverflowError: complex exponentiation\r\n\u003e\u003e\u003e z*z*z*z*z\r\n(-inf+infj)\r\n```\r\n\r\nIf this issue does make sense, I'll provide a patch.\r\n\r\n### CPython versions tested on:\r\n\r\nCPython main branch\r\n\r\n### Operating systems tested on:\r\n\r\n_No response_\n\n\u003c!-- gh-linked-prs --\u003e\n### Linked PRs\n* gh-120287\n\u003c!-- /gh-linked-prs --\u003e\n","author":{"url":"https://github.com/skirpichev","@type":"Person","name":"skirpichev"},"datePublished":"2024-06-04T03:08:28.000Z","interactionStatistic":{"@type":"InteractionCounter","interactionType":"https://schema.org/CommentAction","userInteractionCount":4},"url":"https://github.com/120010/cpython/issues/120010"}