Pell's Equation

 x^2 - Dy^2 = 1という形のDiophantus Equation

ディオファントス方程式なのでsolve_diophantine 関数を使って一般的に解ける

x,y = var("x,y")                                                                                                                                                                        
solutions = solve_diophantine(x**2 - 61*y**2 == 1)
xsol, ysol = solutions[0]
xval, yval = xsol.subs(t=0).simplify_full(), ysol.subs(t=0).simplify_full()