Why do imaginary number calculations behave differently for exponents up to 100 and above 100?

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.

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.)