Need correct root finding over the p-adics

I have a restricted power series $f(t)$ defined over the p-adic integers (restricted means that the coefficients converge to 0). In practice when one sets a precision for the p-adics this power series turns into a polynomial. I want to be able to compute the p-adic integer roots of $f(t)$. More precisely I want an algorithm whose output is a set of disjoint p-adic balls $B_i$ with associated multiplicities $m_i$ subject to the following conditions:

1. The power series f has no zeros outside the union of the balls;
2. If a ball $B$ has multiplicity 1, then there exists a single root of f in $B$;
3. If a ball $B$ has multiplicity m>1, then there can be at most m roots of f in $B$;
4. The total area covered by the balls goes to zero as the precision that we choose goes to infinity.

My attempt to solve this problem using Sage was to use the roots() function, however this function does not satisfy condition 2. For example consider:

R = Zp(3)
P.<x> = R['x']
zero = R.zero().add_bigoh(2)
poly = (1+zero)*x**2+(3+zero)*x + zero
poly.roots()

The output is [(2*3 + O(3^2), 1), (O(3^2), 1)] and it should have been [(O(3),2)]. My question is, is there an easy way to get the described functionality in Sage, or do I have to write everything from scratch?

Thank you for the responses,