diff --git a/.github/CODEOWNERS b/.github/CODEOWNERS index 269770261b8..d072b133376 100644 --- a/.github/CODEOWNERS +++ b/.github/CODEOWNERS @@ -687,6 +687,7 @@ peps/pep-0808.rst @FFY00 peps/pep-0809.rst @zooba peps/pep-0810.rst @pablogsal @DinoV @Yhg1s peps/pep-0811.rst @sethmlarson @gpshead +peps/pep-0812.rst @serhiy-storchaka peps/pep-0814.rst @vstinner @corona10 # ... peps/pep-2026.rst @hugovk diff --git a/peps/pep-0812.rst b/peps/pep-0812.rst new file mode 100644 index 00000000000..dab1cee57d8 --- /dev/null +++ b/peps/pep-0812.rst @@ -0,0 +1,701 @@ +PEP: 812 +Title: Imaginary type and IEC 60559-compatible complex arithmetic +Author: Sergey B Kirpichev +Sponsor: Serhiy Storchaka +Discussions-To: Pending +Status: Draft +Type: Standards Track +Created: 31-Oct-2025 +Python-Version: 3.15 +Post-History: `10-Sep-2023 `__, + `14-Jan-2025 `__ + + +Abstract +======== + +This PEP alters the meaning of the *imaginary literal* in Python. From now it will +denote a newly added :external+py3.14:class:`complex` *subtype* (called therein +imaginary or pure-imaginary complex number) with a zero real component. This +component of imaginary numbers will be *ignored* in arithmetic operations when +the other operand is either a :external+py3.14:class:`complex` (but not imaginary) +or a :external+py3.14:class:`float`. + + +Motivation +========== + +Complex numbers were part of the Python language from very early releases: +since v1.2 in the stdlib and since v1.3 as the +:external+py3.14:class:`complex` builtin. But until recently their +implementation was simplistic, based on the approach (dating back to Fortran +66) that just adds a complex data type as an ordered pair of +:external+py3.14:class:`float` numbers, doing arithmetic on them by the usual +formulae. When operands are of mixed types, this assumes [1]_ we first do type +coercions, i.e. conversion of integers and floating-point numbers to complex +numbers with a zero (0.0) imaginary component. The imaginary unit (``1j`` in +Python) in such a system is a complex number with zero real component. + +Unfortunately, this doesn't work well if components of the complex number are +floating-point numbers, used as a model of the *extended* real line. Such as +ones, specified by the IEC 60559 standard, which beyond normal (finite and +nonzero) numbers have special values like signed zero, infinities and NaNs. +Lets take simple examples with multiplication in Python-like pseudocode (where +``~>`` denotes type coercion and ``~=`` is substitution of ``complex(0, y)`` +for ``yj``): + +.. code:: python + + 2.0 * (inf+3j) ~> complex(2, 0) * complex(inf, 3) + == complex(2.0*inf - 0.0*3.0, 2.0*3.0 + 0.0*inf) + == complex(inf, nan) + +.. code:: python + + 2j * (inf+3j) ~= complex(0, 2) * complex(inf, 3) + == complex(0.0*inf - 2.0*3.0, 2.0*inf + 0.0*3.0) + == complex(nan, inf) + +.. code:: python + + 2j * complex(-0.0, 3) ~= complex(0, 2) * complex(-0.0, 3) + == complex(-0.0*0.0 - 2.0*3.0, -2.0*0.0 + 0.0*3.0) + == complex(-6.0, 0.0) + + +Users with some mathematical background instead would expect, that ``x+yj`` +notation in Python is equal to ``complex(x, y)`` and is just the usual +rectangular form for a complex number with components ``x`` (real) and ``y`` +(imaginary) and that the above examples will work, following the rules of +elementary algebra (with the assumption that ``1j**2`` is ``-1``): + +.. code:: python + + 2.0 * (inf+3j) == complex(2.0*inf, 2.0*3.0) == complex(inf, 6) + 2j * (inf+3j) == complex(-2.0*3.0, 2.0*inf) == complex(-6, inf) + 2j * complex(-0.0, 3) == complex(-2.0*3.0, -2.0*0.0) == complex(-6, -0.0) + + +The same affects addition (or subtraction): + +.. code:: python + + 2.0 - 0j ~= 2.0 - complex(0, 0) ~> complex(2.0, 0) - complex(0, 0) + == complex(2.0 - 0.0, 0.0 - 0.0) == complex(2, 0) + +.. code:: python + + -0.0 + 2j ~= -0.0 + complex(0, 2) ~> complex(-0.0, 0) + complex(0, 2) + == complex(-0.0 + 0.0, 0.0 + 2.0) == complex(0, 2) + + +The simplistic approach for complex arithmetic is the underlying reason for numerous +and recurring issues in the CPython bug tracker (here is an incomplete list: +:cpython-issue:`61538`, :cpython-issue:`66738`, :cpython-issue:`67418`, +:cpython-issue:`69639`, :cpython-issue:`70026`, :cpython-issue:`71550`, +:cpython-issue:`84450`, :cpython-issue:`85657`, :cpython-issue:`105027`, +:cpython-issue:`107854`, :cpython-issue:`112176`, :cpython-issue:`122615`) and +mathematical libraries in the Python ecosystem (e.g. `numpy/numpy#26310 +`_, `mpmath/mpmath#473 +`_, `mpmath/mpmath#774 +`_). Among others, the broken +``repr(eval(repr(x))) == repr(x)`` invariant for complex numbers: + +.. code:: pycon + + >>> -0.0+1j + 1j + >>> complex(-0.0, 1) # note funny signed integer zero + (-0+1j) + + +To work around the described behavior, it's required to check the operands of +arithmetic expressions and handle special numbers separately. For example, you can't +just blindly take analytic identity from the textbook and use it to implement +some mathematical function. Let's see how this already happens in the stdlib +on `cmath.asin() +`_ +example. Here is a pure-Python version of the same code: + +.. code:: python + + def asin(z): + # asin(z) = -i asinh(iz) + z = complex(-z.imag, z.real) # z -> iz + z = cmath.asinh(z) + return complex(z.imag, -z.real) # z -> -iz + +Note that here we are essentially doing component-wise computations, complex +arithmetic is not used at all. In other words, being simple to implement --- +it's less useful to end users. + +A more modern approach [2]_, reflecting advances in the IEC standard for real +floating-point arithmetic, instead avoids coercion of reals to complexes +(``~>``) and use a separate data type (imaginary) to represent the imaginary +unit, *ignoring its real component in arithmetic* (i.e. no implicit cast +(``~=``) to a complex number with zero real part). The ``cmath.asin()`` would +be implemented with this approach simply by: + +.. code:: python + + def asin(z): + return -1j*cmath.asinh(1j*z) + + +It's pioneered by the C99 standard ([3]_, [4]_). This is also how complex +arithmetic is implemented in the Ada language, see the `Ada 2022 Reference +Manual, Annex G `_. Some +mathematical libraries (like the GNU GSL, see the `GSL Reference Manual, §5.5 +`_) or the GNU MPC, see +the `MPC Reference Manual, §5.7 +`_) have special +routines to implement mixed-mode arithmetic for complex numbers, i.e., when one +operand either pure-real or pure-imaginary. As a side effect, this also +introduces some performance boost for operations with mixed types (e.g., +multiplication of complex and real numbers costs only two real multiplications, +not four). Although it's more important that in the IEC floating-point +arithmetic results here are uniquely determined by the usual mathematical +formulae. + +For a first step, :cpython-pr:`124829` added in the CPython 3.14 mixed-mode +rules for complex arithmetic, combining real and complex operands. So, some +examples from above now are working correctly: + +.. code:: pycon + + >>> from cmath import inf + >>> 2*(inf+3j) + (inf+6j) + >>> 2-0j + (2-0j) + + +Unfortunately, this is only a half-way solution. To fix the rest of examples +we need a separate type for pure-imaginary complex numbers. + + +Rationale +========= + +Let's collect here some arguments (with possible answers) against adoption of +a newer approach to complex arithmetic in Python. + + +Special cases are rare +---------------------- + +That's not true, as special numbers come not just from input, but also +*during* computations (e.g. with underflow or overflow). Robust software must +account for them and currently this usually requires reinventing complex +arithmetic in an application, i.e. use the :external+py3.14:class:`complex` type +just as a bag for its components. The only known cure for this is presented +by the PEP: + + Generally, mixed-mode arithmetic combining real and complex variables + should be performed directly, not by first coercing the real to complex, + lest the sign of zero be rendered uninformative; the same goes for + combinations of pure imaginary quantities with complex variables. + + -- Kahan, W: `Branch cuts for complex elementary functions + `_. + + +The complex facility should be simple +------------------------------------- + +Simplicity is a goal, but most importantly for the end user, not implementers +of the complex arithmetic. + + +Upcoming C2y standard abandons the imaginary type +------------------------------------------------- + +That [5]_ might be viewed as a failure of the new approach: no compiler from major +players had a correct implementation of the C99 Annex G. + +On the other hand, this might be also viewed as an indication of poor adoption of +Annex G itself. Notably, MSVC lacks its support. So, neither the CPython, nor +any other Python implementation (per author's knowledge), uses native complex +arithmetic from the C language, and things will hardly be changed soon. + +It's also important to note that removal documents from the C language +committee don't discuss mathematical arguments for the imaginary types at all +[6]_. + + +Specification +============= + +The :external+py3.14:ref:`imaginary literals ` create instances of +imaginary type: + +.. class:: imaginary(x=0.0) + + Create an imaginary number from a number or a string. + + This is equivalent of ``float(x)*1j``. + + +Strictly speaking, we need only one such object, the imaginary unit ``1j``, +with a property: + +.. code:: + + 1j*1j == -1.0 + + +An arbitrary complex value is a direct sum of a pure-real +(:external+py3.14:class:`float` number) part and a pure-imaginary complex +number and the following identities hold (assuming ``x`` and ``y`` are +:external+py3.14:class:`float`\ s): + +.. code:: python + + complex(x, y) == x + y*1j + repr(complex(x, y)) == repr(x + y*1j) + + +Imaginary and complex numbers will have disinct string representations: + +.. code:: python + + repr(complex(x, y)) = ("(" + + format(x, ".0f" if x and x.is_integer() else "") + + ("+" if math.copysign(1, y) == 1 else "") + + repr(y*1j) + ")") + repr(x*1j) = (repr(x) + "j").replace(".0j", "j") + +Parsing strings with the integer "negative zero" in real part (i.e. ``"-0+1j"`` +or ``"(-0+1j)"``) will be deprecated in the :external+py3.14:class:`complex` +constructor. + +The :mod:`marshal` module will be adjusted to support new type. + +The :func:`ast.literal_eval` and :func:`compile` will reject non-imaginary +complex values for :class:`ast.Constant` node type, e.g. +``ast.Constant(value=1+2j)`` (now permitted). + + +Complex Arithmetic +------------------ + +The tables below define unary operations, additive operators (binary ``+`` and +``-``) and multiplicative operators (``*`` and ``/``) for all possible +combinations of types (integer operand values will be implicitly converted to +:external+py3.14:class:`float`'s). In all cases the result approximates the +real and imaginary parts, respectively, of the mathematical formula to be +computed. + +.. table:: Unary operations + :align: left + + +--------+--------+--------------+---------------+-------------+ + | z | +z | -z | z.conjugate() | abs(z) | + +--------+--------+--------------+---------------+-------------+ + | x | x | -x | x | abs(x) | + +--------+--------+--------------+---------------+-------------+ + | yj | yj | (-y)j | (-y)j | abs(y) | + +--------+--------+--------------+---------------+-------------+ + | x + yj | x + yj | (-x) + (-y)j | x + (-y)j | hypot(x, y) | + +--------+--------+--------------+---------------+-------------+ + + +.. table:: Addition and subtraction + :align: left + + +----------+------------+--------------+--------------------+ + | ± | u | vj | u + vj | + +----------+------------+--------------+--------------------+ + | x | x ± u | x + (±v)j | (x ± u) + (±v)j | + +----------+------------+--------------+--------------------+ + | yj | ±u + yj | (y ± v)j | ±u + (y ± v)j | + +----------+------------+--------------+--------------------+ + | x + yj | x ± u + yj | x + (y ± v)j | (x ± u) + (y ± v)j | + +----------+------------+--------------+--------------------+ + + +If both operands have imaginary type, then the result has imaginary type. If +one operand has real type and the other operand has imaginary type, or if +either operand has complex type, then the result has complex type. + + +.. table:: Multiplication + :align: left + + +----------+----------------+-----------------+----------------------------+ + | ``*`` | u | vj | u + vj | + +----------+----------------+-----------------+----------------------------+ + | x | x*u | (x*v)j | (x*u) + (x*v)j | + +----------+----------------+-----------------+----------------------------+ + | yj | (y*u)j | -y*v | (-y*v) + (y*u)j | + +----------+----------------+-----------------+----------------------------+ + | x + yj | (x*u) + (y*u)j | (-y*v) + (x*v)j | (x*u - y*v) + (y*u + x*v)j | + +----------+----------------+-----------------+----------------------------+ + + +.. table:: Division (assuming ``w = u**2 + v**2``) + :align: left + + +----------+----------------+---------------+----------------------------------+ + | / | u | vj | u + vj | + +----------+----------------+---------------+----------------------------------+ + | x | x/u | (-x/v)j | x*u/w + (-x*v/w)j | + +----------+----------------+---------------+----------------------------------+ + | yj | (y/u)j | y/v | y*v/w + (y*u/w)j | + +----------+----------------+---------------+----------------------------------+ + | x + yj | (x/u) + (y/u)j | y/v + (-x/v)j | (x*u + y*v)/w + ((y*u - x*v)/w)j | + +----------+----------------+---------------+----------------------------------+ + + +If one operand has real type and the other operand has imaginary type, then the +result has imaginary type. If both operands have imaginary type, then the +result has real type. If either operand has complex type, then the result has +complex type. + +Results of all mixed-mode operations, except for division when the right +operand is a complex number, are unambigously defined by floating-point +arithmetic. This specification does not indicate how exactly the results are +to be evaluated [7]_ for complex multiplication (when *both* operands are +complex numbers) and for division (when the *right* operand is a complex +number). Though, the following holds. + +For complex multiplication [8]_: + +- + .. code:: + + (z * w).conjugate() == z.conjugate() * w.conjugate() + z * w = w * z + (-z) * w = -(z * w) + +- if one operand is an infinity [9]_ and the other operand is a nonzero finite + number or an infinity, then the result is an infinity. + +For division: + +- with a nonzero right operand + + .. code:: + + (z / w).conjugate() == z.conjugate() / w.conjugate() + (-z) / w = -(z / w) + z / (-w) = -(z / w) + +- if the first operand is an infinity and the second operand is a finite + number, then the result is an infinity; + +- if the first operand is a finite number and the second operand is an + infinity, then the result is a zero; + +- if the first operand is a nonzero finite number or an infinity and the second + operand is a zero, then the result is an infinity. + + +New C API +--------- + +.. c:type:: PyImaginaryObject + + This subtype of :c:type:`PyComplexObject` represents a Python imaginary + number object. + + +.. c:var:: PyTypeObject PyImaginary_Type + + This instance of :c:type:`PyTypeObject` represents purely imaginary numbers, + the Python complex number type *without* real component. + + +.. c:function:: int PyImaginary_Check(PyObject *p) + + Return true if its argument is a :c:type:`PyImaginaryObject` or a subtype of + :c:type:`PyImaginaryObject`. This function always succeeds. + + +.. c:function:: int PyImaginary_CheckExact(PyObject *p) + + Return true if its argument is a :c:type:`PyImaginaryObject`, but not a + subtype of :c:type:`PyImaginaryObject`. This function always succeeds. + + +.. c:function:: PyObject* PyImaginary_FromDouble(double imag) + + Return a new :c:type:`PyImaginaryObject` object with *imag* imaginary + component. Return ``NULL`` with an exception set on error. + + The imaginary component value of a new object could be taken with + :c:func:`PyComplex_ImagAsDouble`. + + +In conformance with the recent C-API group `decision +`__ we don't offer API +to do arithmetic on low-level representation of complex numbers in CPython. +Instead, it's expected that C-API users either will use `PyNumber_* +`__ API or will export numbers +from Python objects and do arithmetic with some external library (like the GNU +GSL), then import back. + + +It's not Magic +============== + +New arithmetic rules correct some more examples, where using known analytic +identities produced wrong results. Here is an example with :func:`~cmath.atan` +near branch cut: + +.. code:: pycon + + >>> import cmath + >>> z = 2j - 0 # or complex(-0.0, 2) + >>> cmath.atan(z) + (-1.5707963267948966+0.5493061443340549j) + >>> atan = lambda z: 1j*(cmath.log(1 - 1j*z) - cmath.log(1 + 1j*z))/2 + >>> atan(z) # was "(1.5707963267948966+0.5493061443340549j)" + (-1.5707963267948966+0.5493061443340549j) + + +Though, we should mention that floating-point arithmetic is not a replacement +for ``limit()`` facilities of computer algebra systems. Using the same identity +near a real line will show wrong results: + +.. code:: pycon + + >>> z = 2+0j + >>> cmath.atan(z) + (1.1071487177940904+0j) + >>> atan(z) + (1.1071487177940904+0j) + >>> z = 2-0j + >>> cmath.atan(z) + (1.1071487177940904-0j) + >>> atan(z) + (1.1071487177940904+0j) + + +Of course, similar happens already for real floating-point arithmetic: + +.. code:: pycon + + >>> f = lambda x: (1 + x)/(1 - x) - 1 + >>> f(1e-15) + 2.220446049250313e-15 + >>> f(0.0) + 0.0 + >>> f(-0.0) + 0.0 + >>> f(-1e-15) + -2.1094237467877974e-15 + +Applications must carefully choose expressions from equivalent forms. + + +Backwards Compatibility +======================= + +In one sense, this PEP should have relatively low impact for end users. + +Indeed, no new syntax is introduced. The results for complex arithmetic will be +different if computation triggers some corner cases, where before either +meaningless values were obtained (``nan``\ s) or wrong zero signs. In the +latter case, the results will be indistinguishable from equality (``==``) testing. + +Major differences impose the new rule ``1j*1j -> float(-1)`` (was +``complex(-1, 0.0)``). Though, there is again no difference for equality. + +Here we list variants of backward incompatible behavior: + +.. code:: pycon + + >>> type(1j) # was "" + + >>> type(-123j) # was "" + + >>> -123j # was "(-0-123j)" + -123j + >>> complex(+0.0, 1) # was "1j" + (0.0+1j) + >>> complex(-0.0, 1) # was "(-0+1j)" + (-0.0+1j) + >>> complex('1j') # was "1j" + (0.0+1j) + >>> format(1j, "f") # was '0.000000+1.000000j' + '1.000000j' + >>> format(-1j, "f") # was '-0.000000-1.000000j' + '-1.000000j' + >>> +0.0+1j # was "1j" + (0.0+1j) + >>> -0.0+1j # was "1j" + (-0.0+1j) + >>> float('inf')*1j # was "(nan+infj)" + infj + >>> float('nan')*1j # was "(nan+nanj)" + nanj + >>> -0.0*1j # was "(-0-0j)" + -0j + + +Working with the initial implementation shows that most test failures in the +CPython test suite come from cases where imaginary literals are used just as +"some complex numbers" to produce exceptions. Running the `mpmath +`_ (development version) test suite shows +only two test failures. `NumPy `_ (v2.3.4) +has around 16 broken tests. + + +How to Teach This +================= + +While internally complex arithmetic will be more complicated (but not too much, +see `Reference Implementation `_), its semantics +will be closer to the usual mathematical notation in textbooks on complex +analysis, with much less room for confusion for newcomers. Roughly speaking, it +will be the floating-point arithmetic, augmented by the special algebraic +symbol ``1j``, whose square is ``-1.0``. + + +Reference Implementation +======================== + +A draft implementation is available at +https://github.com/skirpichev/cpython/pull/1. + + +Rejected Ideas +============== + +We might try to implement complex arithmetic that will specially treat complex +numbers of the form ``complex(0.0, y)`` --- just like instances of imaginary type +proposed. But such a proposal alters the arithmetic rules on the set of +complex numbers itself, in particular ``complex(a, b) + complex(c, d)`` will +not be exactly equal to ``complex(a + c, b + d)`` anymore. + +The set of imaginary numbers with special treatment in complex arithmetic +might be implemented differently, as a distinct form for the complex type +constructor, say ``complex(imag=y, pure=True)``. However, experiments show +that such an implementation is more complicated internally and harder to +explain then a dedicated concept of imaginary numbers as a subtype of complex. + + +Open Issues +=========== + +The PEP expose a new subtype as the :class:`imaginary` builtin, but maybe we +shouldn't? This looks redundant, as all imaginary values could be obtained by +scaling the imaginary unit, i.e. ``imaginary(x) == float(x)*1j``. + + +Acknowledgements +================ + +Thanks to Mark Dickinson for pointing to the right solution and helpful +discussion on various earlier versions of this idea. + + +Footnotes +========= + +.. [1] `The Fortran 2023 standard + `_ (ISO/IEC 1539:2023) + §10.1.5.2.1 says: + + Except for a value of type real or complex raised to an integer + power, if the operands have different types or kind type + parameters, the effect is as if each operand that differs in type + or kind type parameter from those of the result is converted to the + type and kind type parameter of the result before the operation is + performed. + + + Although it's not specified how exactly operations are implemented: + + Quite apart from the fact of the exclusion, the Fortran standard + itself contains no specification or requirement on the algorithm + used to calculate complex multiplication. + + As was pointed out in email, there are algorithms for complex + multiply other than the "traditional" one. One such algorithm + omits parts of the traditional calculation when the real or + imaginary part of one of the operands is known to be zero. + + Furthermore, as the standard contains no specification or + requirement, it thus contains no requirement that the same + algorithm be used at all times. Thus anything + "processor-dependent" can depend on "the phase of the moon" or + indeed anything else. + + -- https://j3-fortran.org/doc/year/24/24-179.txt + +.. [2] W. Kahan and J. W. Thomas. `Augmenting a Programming Language with + Complex Arithmetic + `_. + Technical Report UCB/CSD 91/667, Univ. of Calif. at Berkeley, December, + 1991. + +.. [3] ISO/IEC 9899:1999, Annex G. See `N1256 (final draft) + `_. + https://open-std.org/, WG14. 2007. + +.. [4] See `Rationale for C99 + `_. + https://open-std.org/, WG14, 2003 and `Issues Regarding Imaginary Types for + C and C++ + `_, + by Jim Thomas and Jerome T. Coonen, The Journal of C Language Translation, + Volume 5, Number 3, March 1994. + +.. [5] `N3274: Remove imaginary types + `_. + https://open-std.org/, WG14. June 14, 2024. + +.. [6] See `N3206: The future of imaginary types + `_, WG14. 2023. + +.. [7] Another alternative for the multiplication of complex numbers + (with only three multiplies), see e.g. "Handbook of Floating-Point + Arithmetic" by Muller at al, 2010, Algorithm 4.8: + + .. code:: python + + def karatsuba_mul(z, w): + x, y = z.real, z.imag + u, v = w.real, w.imag + p1 = (x + y)*(u + v) + p2 = x*u + p3 = y*v + return complex(p2 - p3, p1 - p2 - p3) + + Other variants include using a fused multiply-add (FMA) instruction. + +.. [8] Its easy to smash such properties, e.g. the commutativity, as shows the + `following algorithm `_: + + .. code:: python + + def mul3(z, w): + x, y = z.real, z.imag + u, v = w.real, w.imag + t = x*(u + v) + return complex(t - v*(x + y), t + u*(y - x)) + + Example on IEC 60559-conformant system: + + .. code:: pycon + + >>> z = 1e308+1.6e308j + >>> w = 0.1+1j + >>> mul3(z, w) + (-inf+1.1600000000000002e+308j) + >>> mul3(w, z) + (inf+infj) + +.. [9] A complex value with at least one infinite part is regarded as an + infinity (even if its other part is a NaN). A complex value is a finite + number if each of its parts is a finite number (neither infinite nor NaN). + A complex value is a zero if each of its parts is a zero. + + +Copyright +========= + +This document is placed in the public domain or under the +CC0-1.0-Universal license, whichever is more permissive.