x,y = var('x,y') p = implicit_plot(y^2 - (x * (x + 1) * (2*x + 1)/6), (-4,4), (-4,4)) p.show() MIN = 1 MAX = 30 x,y = var('x,y') p = implicit_plot(y^2 - (x * (x + 1) * (2*x + 1)/6), (-4,4), (-4,4)) for px in range(MIN,MAX+1): for qx in range(MIN,MAX+1): for py in range(MIN,MAX+1): for qy in range(MIN,MAX+1): x = px/qx y = py/qy val = y^2 - (x * (x + 1) * (2*x + 1))/6 if (val == 0): p += circle((x,y),0.1,color='red') p.show() print("done.") x = var('x') x = PolynomialRing(QQ,x).gen() f = (3*x-2)^2 - (x * (x + 1) * (2*x + 1)/6) f.roots()