Current solution for Project Euler 250+ (HackerRank) is giving incorrect answers, why is this? [closed]

^{Note: I originally posted this on Stack Overflow on Stack Exchange here however I am posting the same question here as I might be able to get a answer more quickly here.}

---

Project Euler+ #250 on HackerRank is as follows:

Find the number of non-empty subsets of $\{1^1, 2^2, 3^3,..., n^n\}$, the sum of whose elements are divisible by $k$. Print your answer modulo $10^9$

You can find this here.

The constraints are as follows:

$$\begin{align}1\le&\;n\le10^{400}\\3\le&\;k\le50\end{align}$$

I originally posted a Code Review question here, but I have since modified my solution after taking a 2-month break from HackerRank and after solving the original problem on Project Euler.

---

Here is my current code (using Python 3.12):

from math import lcm

MOD = 10 ** 9

def phi(n):
    c = 2
    s = n
    while n > 1:
        if n % c == 0:
            s -= s // c
        
        while n % c == 0:
            n //= c
        
        c += 1
    
    return s

def compute(N, n, m = MOD):
    cycle = lcm(n, phi(n))
    # A guess for what the period for `i ** i % k` will be.
    
    def convolve(A, B):
        return [sum(a * b for a, b in zip(A, B[i::-1] + B[:i:-1])) % m for i in range(n)]
    
    q, rem = divmod(N, cycle)
    A = end = [1] + [0] * (n - 1)
    for i in range(1, cycle + 1):
        k = pow(i, i, n)
        A = [(a + b) % m for a, b in zip(A, A[-k:] + A[:-k])]
        if i == rem:
            end = A
    
    base, A = A, [1] + [0] * (n - 1)
    for j in bin(q)[2:]:
        A = convolve(A, convolve(A, base) if j == '1' else A)
    
    A = convolve(A, end)
    return (A[0] - 1) % m
#---------------------------------------------------------------
n, k = map(int, input().split())
print(compute(n, k))

However, this does not solve all 56 test cases. In fact, for all 9 test cases I don't solve (those being hidden test cases 11, 16, 17, 20, 23, 28, 32, 52, 53), I am getting incorrect answers.

---

My question is: What is wrong with my code that is causing me to get incorrect answers, and how could I possibly fix this?

Something tells me that it has either to do something with my guess for the cycle parameter for some values of $k$, or maybe the problem is with how I am computing the convolutions, but I'm not really sure.