The Weierstrass preparation theorem for formal power series rings guarantees that if a given formal series $f(z) = \sum a_k z^k \in R[[z]]$ where $R$ is a complete local ring with maximal ideal $M$ has $a_k \in M$ for $k < n$ and $a_n \in R^* = M^c$, then

$$ f = (z^n + b_{n-1}z^{n-1} + \cdots + b_0)u $$ where $b_k \in M$ and $u$ is a unit in $R[[z]]$.

I need an explicit algorithm for calculating this Weierstrass polynomial (or distinguished polynomial) for a given $f$. In my case the coefficient ring is $R = \mathbb Z_3[[x]]$, formal power series over the 3-adics. So any algorithm would have to be robust enough to handle these coefficients.

Does anyone know of such an algorithm for a math software package?