Improve accuracy for complex powers with small negative integer exponents

[skirpichev]

Feature or enhancement

Proposal:

Currently we compute such powers as (1+0j)/z**(-n) (unless absolute value of n is too big to use specialized algorithm for integer exponents). This approach, however, introduce huge accuracy loss due to underflows in division.

For example:

>>> from gmpy2 import *                                                                                                                    
>>> import math                                                                                                                            
>>> x, y = map(float.fromhex, ['0x1.47e9c711723f5p+81',                                                                                    
...                            '0x1.38afd1168e49fp+85'])
>>> z = complex(x, y)                                                                                                                      
>>> pr = pow(z, -12); pr                                                                                                                   
0j
>>> pr2 = pow((1/z), 12); pr2                                                                                                              
(5.562684646267994e-309+5.56268464626799e-309j)
>>> gr = complex(pow(mpc(z), -12))                                                                                                         
>>> abs((pr2-gr).real)/math.ulp(gr.real)                                                                                                   
2.0
>>> abs((pr2-gr).imag)/math.ulp(gr.imag)                                                                                                   
2.0

Using instead ((1+0j)/z)**(-n) reduced error in this example from ~1e15 ULP to 2ULP. Note that, generic power algorithm is not affected by this issue:

>>> import _testcapi
>>> pr3 = _testcapi._py_c_pow(z, -12)[0]; pr3
(5.56268464626801e-309+5.56268464626799e-309j)
>>> abs((pr3-gr).real)/math.ulp(gr.real)
1.0
>>> abs((pr3-gr).imag)/math.ulp(gr.imag)
2.0

Has this already been discussed elsewhere?

This is a minor feature, which does not need previous discussion elsewhere

Links to previous discussion of this feature:

No response