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.
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$.