I've been playing with Puiseux series and noticed that zero coefficients are not being assigned the correct valuation of infinity, so that the Newton polygon is not always correct. A minimal example follows: the two newton polygons should be the same, but they are not. The polygon for f2 is correct, while that of f1 is in error.
R.<x> = PuiseuxSeriesRing(QQ)
S.<y> = PolynomialRing(R)
f1 = y^2+x
f2 = y^2+x*y+x
X1=f1.newton_polygon().plot()
X2=f2.newton_polygon().plot()
show(X1)
show(X2)
If you do the same thing with Laurent series as opposed to Puiseux series, there is no problem and the polygons are equal (as expected):
R.<x> = LaurentSeriesRing(QQ)
S.<y> = PolynomialRing(R)
f1 = y^2+x
f2 = y^2+x*y+x
X1=f1.newton_polygon().plot()
X2=f2.newton_polygon().plot()
show(X1)
show(X2)
Similarly if you work over a p-adic field, there is no problem:
K=pAdicField(2)
S.<y> = PolynomialRing(K)
f1 = y^2+2
f2 = y^2+2*y+2
X1=f1.newton_polygon().plot()
X2=f2.newton_polygon().plot()
show(X1)
show(X2)
So this does seem to be a problem specific to Puiseux series. I'm not sure how to post a bug report, so I'm documenting this here in the hope that someone can help. Thanks!