This is only a sketch of solution (there might be some off-by-ones)
Looking at a permutation of $\{1\ldots,n\}$ is equivalent at looking its inversion table $(a_1, \ldots, a_n)$ where $a_i$ is the number of elements to the left of $i$ that are greater than $i$. Basically that gives you a bijection between $S_n$ and $\prod_{1\le i\le n} \{0,\ldots,i-1\}$ and the number of inversions of the permutation is exactly the sum of all $a_i$ corresponding to it (see for example Knuth TAOCP volume 3 for a nice development).
So the problem is equivalent to computing the cardinal of:
$$\{(a_1, \ldots, a_n) \in \prod_{i=1}^{n}\{0,\ldots,i-1\}, \sum_{i=1}^n a_i=k\}$$
This way, we find a DP algorithm with $n\times k$ states:
$$\text{dp}(n,k)=\sum_{i=0}^{\min(k,n-1)}\text{dp}(n-1,k-i)$$
Thus:
- When $n$ is small ($n\le k+1$):
$$\text{dp}(n,k)=\sum_{i=0}^{n-1}\text{dp}(n-1,k-i)=\text{dp}(n,k-1)+\text{dp}(n-1,k)$$
So you can compute these in $O(k^2)$. In particular, we will need the following vector of size $k$:
$$\begin{bmatrix}
\text{dp}(k+1,0)\\
\vdots\\
\text{dp}(k+1,k)
\end{bmatrix}$$
- When $n$ is large ($n>k+1$):
$$\text{dp}(n,k)=\sum_{i=0}^{k}\text{dp}(n-1,k-i)=\sum_{i=0}^k \text{dp}(n-1,i)$$
In this case $\text{dp}(n,\cdot)$ depends linearly on $\text{dp}(n-1,\cdot)$, the matrix being a $k\times k$ lower triangular matrix with nonzero coefficients being all 1-s. More formally, you can rewrite the previous equality as:
$$\begin{bmatrix}
\text{dp}(n,0)\\
\vdots\\
\text{dp}(n,k)
\end{bmatrix}=
\begin{bmatrix}
1&0&\ldots&\ldots&0\\
1&1&0&\ldots&0\\
\vdots&\vdots&\ddots&\ddots&\vdots\\
1&1&\ldots&1&0\\
1&1&\ldots&1&1
\end{bmatrix}
\begin{bmatrix}
\text{dp}(n-1,0)\\
\vdots\\
\text{dp}(n-1,k)
\end{bmatrix}$$
Just let me note:
$$A=\begin{bmatrix}
1&0&\ldots&\ldots&0\\
1&1&0&\ldots&0\\
\vdots&\vdots&\ddots&\ddots&\vdots\\
1&1&\ldots&1&0\\
1&1&\ldots&1&1
\end{bmatrix}$$
Now if you iterate the relation, you will get for any $i$:
$$\begin{bmatrix}
\text{dp}(n+i,0)\\
\vdots\\
\text{dp}(n+i,k)
\end{bmatrix}=
A^i
\begin{bmatrix}
\text{dp}(n,0)\\
\vdots\\
\text{dp}(n,k)
\end{bmatrix}$$
In particular:
$$\begin{bmatrix}
\text{dp}(n,0)\\
\vdots\\
\text{dp}(n,k)
\end{bmatrix}=
A^{n-k-1}
\begin{bmatrix}
\text{dp}(k+1,0)\\
\vdots\\
\text{dp}(k+1,k)
\end{bmatrix}$$
The vector on the left hand side is what we want ($\text{dp}(n,k)$), the vector on the right hand side is what we computed in the first step (when $n$ is small). So if we can compute the $k\times k$ matrix $A^{n-k-1}$, computing the answer can be done in $O(k^2)$.
Now let's see how fast we can compute $A^{n-k-1}$. I suggest three different ideas (I assume that you want the answer modulo some prime $p$):
1) Matrix exponentiation: if we can multiply $k\times k$ matrices with coefficients in $\mathbb{Z}_p$ in time $O(k^\omega)$), you can use fast exponentiation to have a $O(k^\omega \log n)$ algorithm ($\omega =2.38$ nowadays). Actually, if you are fine with a $O(k^3 \log n)$ algorithm, this should be easy to code (with the naive matrix multiplication algorithm)
2) Polynomial interpolation: if you look at the coefficients of $(A^T)^n$, it looks like:
$$(A^T)^n=\begin{bmatrix}
P_1(n)&P_2(n)&\ldots&P_k(n)\\
0&P_!(n)&\ldots&P_{k-1}(n)\\
\vdots&\ddots&\ddots&\vdots\\
0&\ldots&0&P_1(n)
\end{bmatrix}$$
where $P_1, \ldots, P_k$ are polynomials with integer coefficients and $P_i$ is of degre $i-1$. They satisfy some kind of "generalized Pascal triangle relationship":
$$P_i(n)=1+\sum_{j=0}^{n-1}P_{i-1}(j)$$
So it's possible to compute $P_i(j)$ for $1\le i \le k$ and $0\le j\le k-1$ in $O(k^2)$. We can find coefficients of $P_i$ by doing a Newton interpolation in $O(k \log k)$, so this should give a $O(k^2 \log k)$ algorithm.
3) Find a closed form for $A^i$ with the help of some well-known combinatorial objects. As I said in the comments, the polynomials $P_1, \ldots, P_k$ are close to Faulhaber's polynomials which can be computed in linear time. This might require some heavy handmade computation. I think it should lead to a $O(k^2)$ algorithm.