Matrix Rank (Hard Version)
This is the hard version of the problem. The only differences between the two versions of this problem are the constraints on $$$k$$$. You can make hacks only if all versions of the problem are solved.You are given integers $$$n$$$, $$$p$$$ and $$$k$$$. $$$p$$$ is guaranteed to be a prime number. For each $$$r$$$ from $$$0$$$ to $$$k$$$, find the number of $$$n \times n$$$ matrices $$$A$$$ of the field$$$^\dagger$$$ of integers modulo $$$p$$$ such that the rank$$$^\ddagger$$$ of $$$A$$$ is exactly $$$r$$$. Since these values are big, you are only required to output them modulo $$$998\,244\,353$$$.$$$^\dagger$$$ https://en.wikipedia.org/wiki/Field_(mathematics)$$$^\ddagger$$$ https://en.wikipedia.org/wiki/Rank_(linear_algebra)
The first line of input contains three integers $$$n$$$, $$$p$$$ and $$$k$$$ ($$$1 \leq n \leq 10^{18}$$$, $$$2 \leq p < 998\,244\,353$$$, $$$0 \leq k \leq 5 \cdot 10^5$$$).It is guaranteed that $$$p$$$ is a prime number.
Output $$$k+1$$$ integers, the answers for each $$$r$$$ from $$$0$$$ to $$$k$$$.
(n, p, k) = map(int, input().split())
p1 = [1]
m = 998244353
for i in range(k):
    p1.append(p1[-1] * p % m)
p0 = [1]
q = pow(p, -1, m)
for i in range(k):
    p0.append(p0[-1] * q % m)
y = pow(p, n, m)
p2 = [1]
for i in range(1, k + 1):
    if p1[i] == 1:
        p2.append(p2[-1] * i % m)
    else:
        p2.append(p2[-1] * (p1[i] - 1) % m)
p3 = [pow(p2[-1], -1, m)]
for i in range(k, 0, -1):
    if p1[i] == 1:
        p3.append(p3[-1] * i % m)
    else:
        p3.append(p3[-1] * (p1[i] - 1) % m)
p3.reverse()
x = 1
ans = []
z = 1
for r in range(k + 1):
    if r > n or x == 0:
        ans.append(0)
        continue
    ans.append(x * x * p3[r] * z % m)
    x *= y - p1[r]
    x %= m
    z *= p0[r]
    z %= m
print(*ans)
MOD = 998244353

def power(x, y):
    return pow(x, y, MOD)
(n, p, k) = list(map(int, input().split()))
pn = power(p, n)
p1 = [1]
for i in range(k):
    p1.append(p1[-1] * p % MOD)
q = power(p, -1)
p0 = [1]
for i in range(k):
    p0.append(p0[-1] * q % MOD)
y = power(p, n)
p2 = [1]
for i in range(1, k + 1):
    if p1[i] == 1:
        p2.append(p2[-1] * i % MOD)
    else:
        p2.append(p2[-1] * (p1[i] - 1) % MOD)
p3 = [power(p2[-1], -1)]
for i in range(k, 0, -1):
    if p1[i] == 1:
        p3.append(p3[-1] * i % MOD)
    else:
        p3.append(p3[-1] * (p1[i] - 1) % MOD)
p3.reverse()
num = 1
den = 1
ans = []
for r in range(k + 1):
    if r > n or num == 0:
        ans.append(0)
        continue
    ans.append(num * num * p3[r] * den % MOD)
    num *= y - p1[r]
    num %= MOD
    den *= p0[r]
    den %= MOD
print(*ans)
P = 998244353
(n, p, k) = map(int, input().split())
(nom, denom) = (1, 1)
pr = 1
pn = pow(p, n, P)
(nom0, denom0) = (0, 0)
for r in range(k + 1):
    if r > n or nom0 * 2 > denom0 + 10:
        print(0, end=' ')
        continue
    if nom0 * 2 == denom0:
        print(nom ** 2 * pow(denom, -1, P) % P, end=' ')
    else:
        assert nom0 * 2 > denom0
        print(0, end=' ')
    mul = (pn - pr) % P
    if r == n:
        continue
    if mul % P == 0:
        z = 1
        while (pow(p, n, P ** z) - pow(p, r, P ** z)) % P ** z == 0:
            z += 1
        mul = (pow(p, n, P ** z) - pow(p, r, P ** z)) % P ** z // P ** (z - 1)
        nom0 += z - 1
    nom = nom * mul % P
    pr1 = pr
    pr = pr * p % P
    mul = (pr - 1) % P
    if mul % P == 0:
        z = 1
        while (pow(p, r, P ** z) - 1) % P ** z == 0:
            z += 1
        mul = (pow(p, r, P ** z) - 1) % P ** z // P ** (z - 1)
        denom0 += z - 1
    denom = denom * pr1 * mul % P
print()