Generalized Jordan canonical form Suppose $\mathbb F_q$ is the finite field of order $q$. Let $f(x)=x^d-a_{d-1}x^{d-1}-\cdots-a_{1}x-a_0\in\mathbb F_q[x]$ be irreducible with $\deg (f(x))=d$.
Prove that we can find a basis $\{e_1,...,e_d,f_{1},...,f_d\}$ of the cyclic $\mathbb F_q[x]$-module $\mathbb F_q[x]/(f(x)^2)$, such that
\begin{align}x.e_1&=e_2+f_1\\
x.e_2&=e_3+f_2\\
\cdots&\cdots\cdots\\
x.e_{d-1}&=e_d+f_{d-1}\\
x.e_d&=a_{d-1}e_d+a_{d-2}e_{d-1}+\cdots+a_1e_2+a_0a_1+f_d\\
&\\
x.f_1&=f_2\\
x.f_2&=f_3\\
\cdots&\cdots\\
x.f_{d-1}&=f_d\\
x.f_d&=a_{d-1}f_d+a_{d-2}f_{d-1}+\cdots+a_1f_2+a_0f_1.\end{align}
Moreover, this can be generalized to all $\mathbb F_q[x]/(f(x)^m)$ where $f$ is irreducible over $\mathbb F_q$ and $m$ is a positive integer.

I think this can be deduced from https://en.wikipedia.org/wiki/Frobenius_normal_form#A_rational_normal_form_generalizing_the_Jordan_normal_form. So, I have tried $f(x)=x^2+x+1$ over $\mathbb F_2$ for $m=2$ and computed the Smith normal form of them and they indeed coincide. But I do not have a direct proof of this fact for any $m$ yet.

Edit:
I am aware of this treatment which uses Hensel's lemma. But I believe there should be a much simpler(or elementary) way of proving this.
 A: If you view this as a system of equations to be solved, then you get a relatively pretty system with an elementary solution. Unfortunately, its solution is quite ugly for $m>2$.
Equations such as $x \cdot f_i = f_{i+1}$ can be solved to define $f_{i+1}$. One gets $$f_{i+1} = x^i \cdot f_1$$ for $0 \leq i < d$. A slightly more complex solution to $x \cdot e_i = e_{i+1} + f_i$ can be solved to defined $e_{i+1} = x \cdot e_i - f_i$. Expanding this out in terms of $e_1$ and $f_1$ gives $$e_{i+1} = x^i \cdot e_1 - i x^{i-1} \cdot f_1$$
Notice the relationship between the coefficients of e versus f, the monomial versus its derivative.
Using these relations we can get a nice formula:   $$g \cdot e_1 = \left( \sum_{i=0}^{k} g_i \cdot e_{i+1} \right) + g' \cdot f_1 \qquad \text{where } g = \sum_{i=0}^{k} g_i x^i$$
Applying this to to the equations $\displaystyle x \cdot e_d = \sum_{k=0}^d a_k \cdot e_{k+1} + f_d$, then after some work we get an equation that $f \cdot e_1 - f' \cdot f_1 = 0$. This allows us to solve for $f_1$ in terms of $e_1$, since $f'(x)$ is invertible mod $f(x)$.
$$f_1 = \dfrac{f}{f'} \cdot e_1 = gf \cdot e_1 \qquad \text{where } gf' \equiv 1 \text{ mod } f$$
The only variable we have not solved for is $e_1$. You can take any element of the module $M$ which is not annihilated by $f^{m-1}$. Equivalently, any element of $M$ that is not contained in $fM$. Then $f_1$ will be in $fM$ but not in $f^2M$.
$m > 2$
This can be extended to any positive integer $m$, but you get $m-1$ equations of this form. Since it is hard to index the letter in $e_1, f_1, g_1, \ldots$, I'll switch to double indexing: $e_i = e_{1,i}$, $f_i = e_{2,i}$, $g_i = e_{3,i}$, etc.
Then the equations that must be solved are:$$\begin{array}{rcl}
0 & = & f \cdot e_{m-1,1} - f' \cdot e_{m,1} \\
0 & = & f \cdot e_{m-2,1} - f' \cdot e_{m-1,1} + \tfrac12 f'' \cdot e_{m,1} \\
0 & = & f \cdot e_{m-3,1} - f' \cdot e_{m-2,1} + \tfrac12 f'' \cdot e_{m-1,1} - \tfrac1{3!} f''' \cdot e_{m,1} \\
& \vdots & \\
0 & = & \displaystyle \sum_{k=j}^m (-1)^{(k-j)} \tfrac{1}{(k-j)!} ~f^{(k-j)} \cdot e_{k,1} \qquad \text{for each } 1 \leq j \leq m-1 \\
\end{array}$$
These can all be solved for $e_{j+1,1}$ in terms of $e_{j,1}$, but this formula is not pretty in general. Solving for $e_{j+1,1}$ in terms of $e_{1,1}$ can also be done, but again is quite ugly.
The only variable we have not solved for is $e_{1,1}$. You can take any element of the module $M$ which is not annihilated by $f^{m-1}$. Equivalently, any element of $M$ that is not contained in $fM$. Then $e_{j,1}$ will be in $f^{j-1}M$ but not in $f^jM$.
A note on fields/polynomial generalizations: this argument applies whenever $f$ is separable monic irreducible (so every monic irreducible over a finite or perfect field). Inseparable polynomials don't work for $m\geq 2$ (they don't have such a form). For $m=1$, you don't even need $f$ to be irreducible, but you have to choose $e_1$ not to be annihilated by any irreducible factor of $f$. We are choosing generators of cyclic modules, and if $f$ is irreducible, then those submodules are easily indexed by powers of $f$. If $f$ is reducible, then we must consider all divisors of $f^m$.
A: Streamlined solution
Let $F' := \mathbb{F}_q[\alpha] / \langle f(\alpha) \rangle$ be the extension of $\mathbb{F}_q$ by a root of $f$, and consider the $\mathbb{F}_q[x]$-module $F'[x] / ((x - \alpha)^2)$.  This has basis $1, \alpha, \ldots, \alpha^{d-1}, x-\alpha, \alpha (x-\alpha), \ldots, \alpha^{d-1} (x-\alpha)$ as an $\mathbb{F}_q$-vector space, and the action of $x$ on this vector space gives exactly the desired relations.  Thus, if we show that $F'[x] / \langle (x-\alpha)^2 \rangle \simeq \mathbb{F}_q[x] / \langle (f(x))^2 \rangle$ as $\mathbb{F}_q[x]$-modules, we will be done.
To see this, consider the morphism $\mathbb{F}_q[x] \to F'[x] / \langle (x-\alpha)^2 \rangle$, $p(x) \mapsto p(x) + \langle (x - \alpha)^2 \rangle$.  It is straightforward to check that the kernel of this morphism is $\langle (f(x))^2 \rangle$ -- here it may be useful to use the fact that $F'$ is a Galois extension of $\mathbb{F}_q$.  Therefore, so the morphism induces an injective morphism $\mathbb{F}_q[x] / \langle (f(x))^2 \rangle \to F'[x] / \langle (x-\alpha)^2 \rangle$.  But the source and target have equal dimension as an $\mathbb{F}_q$-vector space, so the morphism is in fact an isomorphism.
This solution will readily generalize to $m > 2$ by considering $F'[x] / \langle (f(x))^m \rangle$ instead.
(Note that the idea of considering $F'[x] / \langle (x-\alpha)^2 \rangle$ is inspired by the fact that the matrix of $x$ in the desired relations, in block form, reads $\begin{bmatrix} C & 0 \\ I_d & C \end{bmatrix}$ where $C$ is the companion matrix of $f$ and therefore acts as a root of $f$.  This suggests that it would be useful to look at the transformation given by the matrix $\begin{bmatrix} \alpha & 0 \\ 1 & \alpha \end{bmatrix}$ over $F'$; and as a standard exercise, the corresponding $F'[x]$-module is isomorphic to $F'[x] / \langle (x-\alpha)^2 \rangle$.
Also note that if we want to generalize to fields other than finite fields, then $F'$ will not necessarily be a Galois extension.  However, if $f$ is a separable polynomial, then we can consider embedding $F'$ into the splitting field of $f$ in the step verifying the kernel of the morphism and get the same result.)
Explicit solution
Consider the set of vectors
$$f_1 := f(x), f_2 := x f(x), \ldots, f_d := x^{d-1} f(x), \\
e_1 := f'(x), e_2 := x f'(x) - f(x), e_3 := x^2 f'(x) - 2x f(x), \ldots, e_d := x^{d-1} f'(x) - (d-1) x^{d-2} f(x).$$
It should be straightforward to check that these vectors satisfy the desired relations.  For example, for the most complex one, note that
$$x e_d = x^d f'(x) - (d-1) x^{d-1} f(x) = (f(x) + a_{d-1} x^{d-1} + \cdots + a_0) f'(x) - (d-1) x^{d-1} f(x) = \\
f'(x) f(x) + a_{d-1} x^{d-1} f'(x) + \cdots + a_0 f'(x) - (d-1) x^{d-1} f(x) = \\
(d x^{d-1} - (d-1) a_{d-1} x^{d-2} - \cdots - a_1) f(x) + a_{d-1} x^{d-1} f'(x) + \cdots + a_0 f'(x) - (d-1) x^{d-1} f(x) = \\
a_{d-1} (x^{d-1} f'(x) - (d-1) x^{d-2} f(x)) + \cdots + a_0 f'(x) + x^{d-1} f(x) = \\
a_{d-1} e_d + a_{d-2} e_{d-1} + \cdots + a_0 e_1 + f_d.$$
It then just remains to show that they span all of the vector space.  For this, we can first reduce to showing that $f(x), \ldots, x^{d-1} f(x), f'(x), \ldots, x^{d-1} f'(x)$ span the vector space, since $e_1, \ldots, e_d, f_1, \ldots, f_d$ are related to $f(x), \ldots, x^{d-1} f(x), f'(x), \ldots, x^{d-1} f'(x)$ via a triangular matrix with 1 on the diagonal.  Then, we can use the fact that since $f(x)$ is irreducible (in fact, we only need that $f(x)$ is square-free), then $\gcd(f(x), f'(x)) = 1$, so for any $g$ we can find polynomials $a, b$ such that $a(x) f(x) + b(x) f'(x) = g(x)$.  Now, we can use the division algorithm of $b(x)$ by $f(x)$ to get a form in which $b(x)$ has degree less than $d$, and then similarly divide $a(x)$ by $f(x)$ and use that $(f(x))^2 = 0$ in the module to get a form in which $a(x)$ also has degree less than $d$.
As for a generalization, I haven't worked out the details.  It might be helpful to notice that except for some signs, $e_k$ looks like $\frac{d}{dx}(x^k f(x))$; so for an extension to the case $m=3$ for example, we might try something like $d_k := x^k f''(x) - 2k x^{k-1} f'(x) + k(k-1) x^{k-2} f(x)$.
