Precision error when finding eigenvalues of Hecke operator

I am trying to use SageMath to calculate the p-adic valuations of the eigenvalues of the Hecke operator, but I am having some issues. I am using the following code to calculate them:

p = 5
F = Qp(p, prec=10000)
v = valuations.pAdicValuation(F, 5)
R.<X> = PolynomialRing(F)
S.<T> = PolynomialRing(QQ)
SHOW_FACTORS = False

for m in [52 .. 53]:
    hs = hecke_series(p, 2, 2, m)
    print(f'\nm = {m}\nhecke series = {hs}')
    for root, mul in S(hs).roots(ring=F):
        print(f'root {v(root)}')
    if not SHOW_FACTORS:
        continue
    try:
        print('Factors:')
        for f, mul in R(hs).factor():
            print(f)
    except:
        print('*** factorization failed ***\n')

But I have found that for different precisions m, I do not get the same results. For example, for m = 52 and m = 53, I get different sets of valuations.

m = 52
hecke series = 248689957516035065054893493652343750*x^13 + 1674403001743485219776630401611328125*x^12 + 1158477226609946228563785552978515625*x^11 + 1681565410974239751696586608886718750*x^10 + 740816027146591471840858459472656250*x^9 + 984876153765252102828224029541015625*x^8 + 2199502199510497526175470605761718750*x^7 + 2047190245401796728083458292414140625*x^6 + 227766241144299957275576934742803125*x^5 + 1626265950251135133773358227520486875*x^4 + 2173838516487905357643185902479851120*x^3 + 966441077617202989620226898790437136*x^2 + 2033736385834117114445851281601109243*x + 1
root 9
root 7
root 4
root 2
root 1

m = 53
hecke series = 4689582056016661226749420166015625000*x^13 + 3894849050993798300623893737792968750*x^12 + 1158477226609946228563785552978515625*x^11 + 1681565410974239751696586608886718750*x^10 + 740816027146591471840858459472656250*x^9 + 7646214301516191345370014038085937500*x^8 + 6640394298011123687869997278125000000*x^7 + 10928974442403049051472511637140703125*x^6 + 4668658339644926118970103607106084375*x^5 + 3846711999501448214620621563702127500*x^4 + 4394284565738218438490449238661491745*x^3 + 3186887126867516070467490234972077761*x^2 + 2033736385834117114445851281601109243*x + 1
root 9
root 7
root 5
root 5
root 4
root 2
root 1
root 0
root 0

I know for a fact that they should be giving the same valuations, however. Indeed, the valuation of the roots should be the same because of the properties of Newton polygon as long as the polynomial for m = 53 reduces modulo p^52 to the polynomial for m = 52. I tested this, using the following code:

hecke_series(5, 2, 2, 52)

R.<q> = PowerSeriesRing(ZZ,default_prec=p*100)
Poly.<x> = PolynomialRing(ZZ)

def ordred2(a, p):
    b = Poly(a)
    # ensure we have enough coefficients to avoid issues with list index
    b += x^(14)
    ordred2_b = b.coefficients(sparse=False)
    ordred2_coeff = list('ordred2_coeff_%d' % d for d in range(14))
    for i in range(14):
        ordred2_coeff[i] = Mod(ordred2_b[i], p^52)
    return R(ordred2_coeff)

ordred2(hecke_series(5, 2, 2, 53), 5)

And found that it does indeed reduce. However, I did notice that if I use

R.<q> = PowerSeriesRing(ZZ, default_prec=p*10)

instead, then it does not reduce properly. This seems to indicate to me that there is some sort of precision error that I am running into in this problem, but no matter what I change in my code, I seem to get the same result when calculating the valuations. Indeed, for Qp(p, prec=100) and Qp(p, prec=10000), I get the same result. Am I missing something? Is there some other place I have to increase the precision, which I am not seeing?