gh-117999: drop special case for small integer exponents in complex_pow()

[skirpichev]

• Issue: Small integer powers of real±0j are invalid #117999

---

This has hegative impact on performance, see #118000 (comment). On another hand, it fixes issue and makes CPython pow() behave like cpow() from the C99 standard. Glibc, for example, has no special code to handle small integer exponents: https://sourceware.org/git/?p=glibc.git;a=blob;f=math/s_cpow_template.c;h=92e3e986639be118ec36fb7f991ba47d8115816e;hb=HEAD Ditto GSL: https://git.savannah.gnu.org/cgit/gsl.git/tree/complex/math.c#n398

[skirpichev]

CC @picnixz

See also #60200

[picnixz]

I had some thoughts on that matter. For small integer exponents $b = \Re(b) + \Im(b)\cdot\jmath$ with $\Im(b) = 0$ and $\Re(b) \in \mathbb{Z}$, we would now call _Py_c_pow in https://github.com/python/cpython/blob/main/Objects/complexobject.c:

vabs = hypot(a.real,a.imag);
len = pow(vabs,b.real);
at = atan2(a.imag, a.real);
phase = at*b.real;
/* if (b.imag != 0.0) {
    len /= exp(at*b.imag);
    phase += b.imag*log(vabs);
} */ // ignored since b.imag == 0
r.real = len*cos(phase);
r.imag = len*sin(phase);

Let $a = r(\cos \phi + \jmath \sin \phi)$ and $b = n\in\mathbb{Z}$. De Moivre's formula tells me that $a^{n} = r^{n}(\cos n\phi + \jmath \sin n\phi)$. So, we could technically have an alternative path for the case $n = 2k$ since we can directly compute vabs^2 instead of vabs. Since hypot treats special cases when inputs are quiet NaN and infinities, or subnormal numbers, we should only use this path when $a$ has finite components. We would save one square root in the computation of hypot and compute len = pow(vabs2, b.real // 2) instead.

I don't know whether the costly computation is due to hypot, atan2 or cos/sin computations. Maybe it won't change anything (and perhaps it would even be slower). Could you therefore check the timings of odd powers (your benchmarks only show even powers and they are twice slower with this patch)?

[skirpichev]

Could you therefore check the timings of odd powers (your benchmarks only show even powers and they are twice slower with this patch)?

Yeah, below are tests copied from #118000 (comment) with few more with odd exponents:

With specialized code (main):

$ python -m pyperf timeit -q -s 'z=1+1j;e=20' 'pow(z,e)'
Mean +- std dev: 237 ns +- 5 ns
$ python -m pyperf timeit -q -s 'z=1+1j;e=25' 'pow(z,e)'
Mean +- std dev: 252 ns +- 2 ns
$ python -m pyperf timeit -q -s 'z=1+1j;e=85' 'pow(z,e)'
Mean +- std dev: 277 ns +- 12 ns
$ python -m pyperf timeit -q -s 'z=1+1j;e=90' 'pow(z,e)'
Mean +- std dev: 259 ns +- 4 ns
$ python -m pyperf timeit -q -s 'z=1+1j;e=95' 'pow(z,e)'
Mean +- std dev: 283 ns +- 2 ns
$ python -m pyperf timeit -q -s 'z=1+1j;e=110' 'pow(z,e)'
Mean +- std dev: 579 ns +- 5 ns
$ python -m pyperf timeit -q -s 'z=1+1j;e=200' 'pow(z,e)'
Mean +- std dev: 571 ns +- 3 ns

Without:

$ python -m pyperf timeit -q -s 'z=1+1j;e=2' 'pow(z,e)'
Mean +- std dev: 562 ns +- 3 ns
$ python -m pyperf timeit -q -s 'z=1+1j;e=20' 'pow(z,e)'
Mean +- std dev: 576 ns +- 2 ns
$ python -m pyperf timeit -q -s 'z=1+1j;e=25' 'pow(z,e)'
Mean +- std dev: 584 ns +- 33 ns
$ python -m pyperf timeit -q -s 'z=1+1j;e=85' 'pow(z,e)'
Mean +- std dev: 583 ns +- 8 ns
$ python -m pyperf timeit -q -s 'z=1+1j;e=90' 'pow(z,e)'
Mean +- std dev: 576 ns +- 3 ns
$ python -m pyperf timeit -q -s 'z=1+1j;e=95' 'pow(z,e)'
Mean +- std dev: 578 ns +- 5 ns
$ python -m pyperf timeit -q -s 'z=1+1j;e=110' 'pow(z,e)'
Mean +- std dev: 578 ns +- 4 ns
$ python -m pyperf timeit -q -s 'z=1+1j;e=200' 'pow(z,e)'
Mean +- std dev: 573 ns +- 3 ns

[picnixz]

Without:

Can I have the numbers for odd exponents without as well please?

[skirpichev]

Sure, comment was updated. As you can see, results aren't biased to even/odd integers.

[picnixz]

Do you think it's worth checking the even case for a fast path or do you think we should just live with those this loss of performance? (personally, in this case, I'm for correctness rather than speed).

[skirpichev]

I think it not worth. Certainly, this can't restore speed for small integers (for even case).

[picnixz]

I'm sorry but I approved without commenting ...

[picnixz]

You could actually have more test cases if you want?

[skirpichev]

Yes, that's why I also referenced #60200. We can include more test cases in some new file of Lib/test/mathdata/. But I'm not sure this belongs to the current pr.

[picnixz]

I like having extensible test coverage so I would be in favor of adding more cases in this PR since we are removing a functionality. But let's leave this decision to the core dev merging it.

misclick

[picnixz]

As I said, I'm in favor of correctness rather than speed in this case. But since #60200 is still under discussion, I'll be interested in hearing from @mdickinson and @serhiy-storchaka for this one.

[picnixz]

I like having extensible test coverage so I would be in favor of adding more cases in this PR since we are removing a functionality. But let's leave this decision to the core dev merging it.

[vstinner]

Can you update your PR? There is now a conflict (with the ERANGE fix).

[vstinner]

So Python has its own cpow() implementation:

Py_complex
_Py_c_pow(Py_complex a, Py_complex b)
{
    Py_complex r;
    double vabs,len,at,phase;
    if (b.real == 0. && b.imag == 0.) {
        r.real = 1.;
        r.imag = 0.;
    }
    else if (a.real == 0. && a.imag == 0.) {
        if (b.imag != 0. || b.real < 0.)
            errno = EDOM;
        r.real = 0.;
        r.imag = 0.;
    }
    else {
        vabs = hypot(a.real,a.imag);
        len = pow(vabs,b.real);
        at = atan2(a.imag, a.real);
        phase = at*b.real;
        if (b.imag != 0.0) {
            len /= exp(at*b.imag);
            phase += b.imag*log(vabs);
        }
        r.real = len*cos(phase);
        r.imag = len*sin(phase);

        _Py_ADJUST_ERANGE2(r.real, r.imag);
    }
    return r;
}

Can we reuse the cpow() of the C library instead? A C11 compiler is now required to build CPython: https://docs.python.org/dev/using/configure.html#build-requirements

[skirpichev]

Can you update your PR? There is now a conflict (with the ERANGE fix).

Done.

Can we reuse the cpow() of the C library instead? A C11 compiler is now required to build CPython: https://docs.python.org/dev/using/configure.html#build-requirements

Yet we don't require optional C11 features, cpow() - is from Annex G. Not even all compilers implement one (mvsc for example, though it has something similar).

And I'm not sure about quality of libm implementations for that part of the standard. After support for C complex in ctypes was merged, I did tests (same as in CMathTests.test_specific_values) of cmath_testcases.txt on top of Debian's 12 default gcc & glibc: to my surprise it was working (modulo two simple failures from changes in C17, now fixed). But perhaps I'm just lucky Linux user.

[vstinner]

I'm not sure that I understood correctly the purpose of this change. Is it supposed to fix some corner cases and enhance correctness? The change only adds 2 tests which is small.

[serhiy-storchaka]

Try complex(-1, 0)**1.

[vstinner]

Try complex(-1, 0)**1.

Without the change, I get:

$ ./python
>>> complex(-1, 0)**1
(-1+0j)

With the change, I get:

>>> complex(-1, 0)**1
(-1+1.2246467991473532e-16j)

[skirpichev]

Is it supposed to fix some corner cases and enhance correctness?

Well, yes. See examples in the issue description:

>>> z = complex(1,-0.0)
>>> z**2
(1+0j)
>>> z*z
(1-0j)
>>> z**2000  # large powers or non-integer - have own opinion
(1-0j)
>>> >>> z**2.1
(1-0j)
>>> _testcapi._py_c_pow(z, 2)[0]  # as well as _Py_c_pow() algorithm
(1-0j)

The change only adds 2 tests which is small.

Agreed. I'll work on this.

Try complex(-1, 0)**1

That's something, which now is more precise in the main, as an accident. But:

>>> _testcapi._py_c_pow(z, 1)[0]
(-1+1.2246467991473532e-16j)
>>> libm.cpow(z, 1)  # libm's cpow
(-1+1.2246467991473532e-16j)

[serhiy-storchaka]

See also #124243.