number of permutation with k inversions
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&gt;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&amp;0&amp;\ldots&amp;\ldots&amp;0\\
1&amp;1&amp;0&amp;\ldots&amp;0\\
\vdots&amp;\vdots&amp;\ddots&amp;\ddots&amp;\vdots\\
1&amp;1&amp;\ldots&amp;1&amp;0\\
1&amp;1&amp;\ldots&amp;1&amp;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&amp;0&amp;\ldots&amp;\ldots&amp;0\\
1&amp;1&amp;0&amp;\ldots&amp;0\\
\vdots&amp;\vdots&amp;\ddots&amp;\ddots&amp;\vdots\\
1&amp;1&amp;\ldots&amp;1&amp;0\\
1&amp;1&amp;\ldots&amp;1&amp;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)&amp;P_2(n)&amp;\ldots&amp;P_k(n)\\
0&amp;P_!(n)&amp;\ldots&amp;P_{k-1}(n)\\
\vdots&amp;\ddots&amp;\ddots&amp;\vdots\\
0&amp;\ldots&amp;0&amp;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.