The imaginary number i, or j in Python, means the square root of -1. So, i to the 4th or any multiple of 4 should be positive 1.
>>> (1j)**4
(1+0j)
>>> (1j)**96
(1+0j)
>>> (1j)**100
(1+0j)
Up until this point all is good, but once we get past 100, Python just bugs out. For example:
>>> (1j)**104
(1+7.842691359635767e-15j)
This messed up my calculations so much in an unexpected way. What explains this? I’m using Python 3.8.19.
5
See the code, for integer exponents up to 100 it uses a different algorithm:
complex_pow(PyObject *v, PyObject *w, PyObject *z)
{
...
// Check whether the exponent has a small integer value, and if so use
// a faster and more accurate algorithm.
if (b.imag == 0.0 && b.real == floor(b.real) && fabs(b.real) <= 100.0) {
p = c_powi(a, (long)b.real);
}
else {
p = _Py_c_pow(a, b);
}
...
The special algorithm for the small exponents uses exponentiation by squaring, just multiplying 1j with itself and with the multiplication results, which is accurate:
c_powu(Py_complex x, long n)
{
Py_complex r, p;
long mask = 1;
r = c_1;
p = x;
while (mask > 0 && n >= mask) {
if (n & mask)
r = _Py_c_prod(r,p);
mask <<= 1;
p = _Py_c_prod(p,p);
}
return r;
}
The general algorithm for exponents larger than 100 and non-integer exponents uses trigonometry and ends up involving inaccurate floats:
_Py_c_pow(Py_complex a, Py_complex b)
{
...
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);
...
The at value for 1j mathematically is π/2, which isn’t representable exactly as float. And then the phase and sin and thus the imaginary party of the result become inaccurate as well. See Why does math.cos(math.pi/2) not return zero? for more about that.
(I just noticed you mentioned using Python 3.8.19, while I looked up the current version, 3.12.6. The code slightly differs, but irrelevantly so, and the explanation is correct for both.)
2