# Integers, or "Integers"?

My understanding is that the solve_diophantine() function is supposed to return only positive integers. However, the following code:

a,b,c = 2,3,5
x, k = var('x,k')
solve_diophantine(a*x^2 + b*x + c - k^2 == 0, (x,k))


returns (partial list):

[(-1/16*sqrt(2)*((470832*sqrt(2) + 665857)^t*(13*sqrt(2) + 20) + (-470832*sqrt(2) + 665857)^t*(13*sqrt(2) - 20)) - 3/4, 1/8*(470832*sqrt(2) + 665857)^t*(13*sqrt(2) + 20) - 1/8*(-470832*sqrt(2) + 665857)^t*(13*sqrt(2) - 20)), ...
(1/16*sqrt(2)*((470832*sqrt(2) + 665857)^t*(19*sqrt(2) + 28) + (-470832*sqrt(2) + 665857)^t*(19*sqrt(2) - 28)) - 3/4, -1/8*(470832*sqrt(2) + 665857)^t*(19*sqrt(2) + 28) + 1/8*(-470832*sqrt(2) + 665857)^t*(19*sqrt(2) - 28))]


And another 14 entries. All of the entries have the same $665857 \pm 470832 \sqrt 2$ as factors. Now, it turns out that $\sqrt 2 \approx \frac{665857}{470832}$ , out to 12 significant digits, a pretty impressive convergent, but... what is it doing in this output?

For reference, there are integer solutions at $(x,k) = (4,7); (19,28)$, and possibly (probably?) others. If I replace the equation with ax^2 + bx + c - 49 == 0, I get 4 as output.

Is there something I've failed to coerce, or something I need to import, or is this just something weird that Sage is doing?

(EDIT: Note that the pair $(19,28)$ appears in that last entry as a factor of $28 +19 \sqrt 2$. But none of the other solutions appear, just... other random covergents of $\sqrt 2$. Nor is $(470832, 665857)$ a solution.

edit retag close merge delete

Sort by ยป oldest newest most voted

Well...

a,b,c = 2,3,5
x, k = var('x,k')
eq=a*x^2 + b*x + c - k^2==0
Sol=solve_diophantine(eq, (x, k))


We obtain a list of 16 solutions of similar form, the first one meaning :

$$x=-\frac{1}{16} \, \sqrt{2} {\left({\left(470832 \, \sqrt{2} + 665857\right)}^{t} {\left(3841 \, \sqrt{2} + 5432\right)} + {\left(-470832 \, \sqrt{2} + 665857\right)}^{t} {\left(3841 \, \sqrt{2} - 5432\right)}\right)} - \frac{3}{4}$$

$$k=\frac{1}{8} \, {\left(470832 \, \sqrt{2} + 665857\right)}^{t} {\left(3841 \, \sqrt{2} + 5432\right)} - \frac{1}{8} \, {\left(-470832 \, \sqrt{2} + 665857\right)}^{t} {\left(3841 \, \sqrt{2} - 5432\right)}$$

These solutions depend of a parameter $t$, a nonnegative integer. It turns out that, notwithstanding the presence of $\sqrt{2}$ , these expressions are integers. We can numerically check this for a not unreasonable range of $t$ :

sage: %time all([all(flatten([[u.subs(t==w).expand().is_integer() for u in v] for v in Sol])) for w in range(100)])
CPU times: user 6.16 s, sys: 80.2 ms, total: 6.24 s
Wall time: 6.24 s
True


Similarny, we can numerically check that these quantities are indeed solutions of eq on a not unreasonable range of $t$ :

sage: %time all([all([bool(eq.subs([x==u[0].subs(t==v), k==u[1].subs(t==v)]).expand()) for u in Sol]) for v in range(100)])
CPU times: user 9.01 s, sys: 88 ms, total: 9.1 s
Wall time: 9.1 s
True


A simple form of the explanation can be found somewhere in a Mathematica tutorial. A better-known analogy is the *casus irreducibilis" of the cubic, where a (demonstrably) real root can be expressed only as the sum of non-real complexes.

A better explanation, by someone who, contrary to me, knows what he's talking about when talking about Diophantine equations, would be extremely welcome...

EDIT :

To be rewritten...

(End of EDIT)

One can note that Sympi's diophantine and Mathematica's Reduce give similar forms to their solutions, but Mathematica gives only 8 solutions, with different (and numerically simpler) constants, possibly by finding equalities among them and simplifying more efficiently.

Notwithstanding, HTH,

more

That's quite interesting; when I run your code with is_integer(), I get False, except for t = 0. Though in addition, the subs(t=w) also threw an exception; I don't know if removing it somehow caused the state to flip from True to False.

( 2022-04-28 00:52:53 +0200 )edit

(If I copy the mathematical expression into the code, and use expand() on the expression, I do get integers.) Edit: Well, I get some integers' some of the lines return half-integers. No clue what's going on <shruggie>

( 2022-04-28 01:09:43 +0200 )edit

Further things I've noticed: (1) all the half-integers I was getting were from missing a negative sign, but also all of those produce negative integer solutions; (2) When $t=0$, all of these boil down to $x = \frac{1}{16} \sqrt 2 \cdot (a \sqrt 2 + b + a \sqrt 2 -b) - \frac34, k = \frac18 (a \sqrt 2 +b - a \sqrt 2 +b)$... in both cases all the $\sqrt 2$ instances cancel out, to $x = \frac14 (a-3), k = \frac14 b$.

( 2022-04-28 02:17:30 +0200 )edit

What version of Sage are you running ? On which platform ? I run 9.6.rc1 in Debian testing.

Can you give us at least one value of $t$ for which you obtain at least a non-integer coefficient or a non-solution (i. e.bool(eq.subs(...) is False) ?

( 2022-04-28 09:01:48 +0200 )edit

To make things more clear, let us look mathematically at the given diophantine equation. It is: $$2x^2+3x+5=k^2\ .$$ This is not exactly a Pell equation, but we may try to work on it. I am multiplying it by eight, then build squares. We get $$16x^2+24x+40=8k^2\ ,$$ and this is the Pell-like equation: $$(4x+3)^2-8k^2=-31\ .$$ It may be that the computer uses instead (see later comments) $$(4k)^2-2(4x+3)^2=62\ .$$ If we ask sage, which asks sympy, for the solutions for a related new equation...

sage: solve_diophantine(X^2 - 8*Y^2 == -31)
[(-2*sqrt(2)*(2*sqrt(2) + 3)^t + 2*sqrt(2)*(-2*sqrt(2) + 3)^t - 1/2*(2*sqrt(2) + 3)^t - 1/2*(-2*sqrt(2) + 3)^t,
-1/8*sqrt(2)*(2*sqrt(2) + 3)^t + 1/8*sqrt(2)*(-2*sqrt(2) + 3)^t - (2*sqrt(2) + 3)^t - (-2*sqrt(2) + 3)^t),
(2*sqrt(2)*(2*sqrt(2) + 3)^t - 2*sqrt(2)*(-2*sqrt(2) + 3)^t - 1/2*(2*sqrt(2) + 3)^t - 1/2*(-2*sqrt(2) + 3)^t,
-1/8*sqrt(2)*(2*sqrt(2) + 3)^t + 1/8*sqrt(2)*(-2*sqrt(2) + 3)^t + (2*sqrt(2) + 3)^t + (-2*sqrt(2) + 3)^t),
(2*sqrt(2)*(2*sqrt(2) + 3)^t - 2*sqrt(2)*(-2*sqrt(2) + 3)^t + 1/2*(2*sqrt(2) + 3)^t + 1/2*(-2*sqrt(2) + 3)^t,
1/8*sqrt(2)*(2*sqrt(2) + 3)^t - 1/8*sqrt(2)*(-2*sqrt(2) + 3)^t + (2*sqrt(2) + 3)^t + (-2*sqrt(2) + 3)^t),
(-2*sqrt(2)*(2*sqrt(2) + 3)^t + 2*sqrt(2)*(-2*sqrt(2) + 3)^t + 1/2*(2*sqrt(2) + 3)^t + 1/2*(-2*sqrt(2) + 3)^t,
1/8*sqrt(2)*(2*sqrt(2) + 3)^t - 1/8*sqrt(2)*(-2*sqrt(2) + 3)^t - (2*sqrt(2) + 3)^t - (-2*sqrt(2) + 3)^t)]


Too many cases, so that for illustration i will use the simpler (not directly related) equation:

var('y,k');
sols = solve_diophantine(y^2 - 2*k^2 == -1, (x, k))


and we get a shorter list. Let us understand this result first.

sage: sols
[(-1/4*sqrt(2)*((2*sqrt(2) + 3)^t*(sqrt(2) + 2) + (-2*sqrt(2) + 3)^t*(sqrt(2) - 2)),
1/4*(2*sqrt(2) + 3)^t*(sqrt(2) + 2) - 1/4*(-2*sqrt(2) + 3)^t*(sqrt(2) - 2)),
(1/4*sqrt(2)*((2*sqrt(2) + 3)^t*(sqrt(2) + 2) + (-2*sqrt(2) + 3)^t*(sqrt(2) - 2)),
-1/4*(2*sqrt(2) + 3)^t*(sqrt(2) + 2) + 1/4*(-2*sqrt(2) + 3)^t*(sqrt(2) - 2))]


Also, let us take a look at the true Pell equation:

sage: var('y,k');
....: sols = solve_diophantine(y^2 - 2*k^2 == 1, (y, k))
(y, k)
sage: sols
[(-1/4*sqrt(2)*(sqrt(2)*(2*sqrt(2) + 3)^t + sqrt(2)*(-2*sqrt(2) + 3)^t),
-1/4*sqrt(2)*(2*sqrt(2) + 3)^t + 1/4*sqrt(2)*(-2*sqrt(2) + 3)^t),
(1/4*sqrt(2)*(sqrt(2)*(2*sqrt(2) + 3)^t + sqrt(2)*(-2*sqrt(2) + 3)^t),
1/4*sqrt(2)*(2*sqrt(2) + 3)^t - 1/4*sqrt(2)*(-2*sqrt(2) + 3)^t)]


It is better to see the two equations above related to the number field $\Bbb Q(a)$, $a:=\sqrt 2$. There, we are searching for elements of norm $-1$, respectively (true Pell) $+1$ where the norm of $\xi:=y\pm k a=y\pm a\sqrt 2$ is $$N(\xi)= N(y\pm ka)=(y+ ka)(y-ka)=y^2-k^2a^2=y^2-2k^2\ .$$ The norm is multiplicative. So if we get one special solution $\xi_0$ for $N(\xi)=-1$ and the general solution $\eta_t$ for $N(\eta)=1$, then the general solution of $N(\xi)=-1$ is of the shape $\xi_0\cdot \eta_t$, parametrized by $t$.

To understand further what happens, note that the number field $K:=\Bbb Q(a)=\Bbb Q(\sqrt 2)$ has a subring of "integers", these are those numbers in $\Bbb Q(a)$ with monic annihilating polynomial in $\Bbb Z[X]$, and it turns out, it is $$R=\mathcal O_K=\Bbb Z[a]\ ,$$ the set of all $x+ay$ with $x,y\in\Bbb Z$. In this ring we can speak of units. These are elements invertible inside this ring. A structural theorem in number theory describes the units of a number field. In our case...

sage: K.<a> = QuadraticField(2)
sage: K.unit_group()
Unit group with structure C2 x Z of Number Field in a with defining polynomial x^2 - 2 with a = 1.414213562373095?
sage: K.unit_group().gens()
(u0, u1)
sage: for gen in K.unit_group().gens():
....:     print(f'generator {K(gen)}')
....:
generator -1
generator a + 1
sage:


The generator $-1$ has norm $N(-1)=(-1)\cdot(-1)=1$.

The generator $1+a$ has norm $N(1+a)=(1+a)\cdot(1-a)=1-a^2=1-2=-1$.

Each unit in $R$ is thus of the shape $\pm(1+a)^s$ for some integer parameter $s$. To get those elements with norm $-1$, we need an odd $s$. (For the elements of norm $+1$ an even $s$.) So the human(ly written) solution of $N(\eta)=-1$ is something like: $$\pm (1+a)^{2t+1} =\pm (1+a)\cdot (1+a)^{2t} =\pm (1+a)\cdot(3+2a)^t\ .$$ As a human, i would always stop here and no longer bother - what are the integral coefficients of this solution, when written w.r.t. the basis $1,a$. But sympy wants to do this for the operator. So instead of delivering $\eta =\pm(1+a)\cdot(3+2a)^t$ as an answer, telling us to extract the components, we get a list. One entry for the $+1$ in $\pm$, one entry for the $-1$ in $\pm$. Then for $\eta=x+ay=(1+a)\cdot(3+2a)^t$ the computer does something humanly horrible, it isolates $x$ as \begin{aligned} x &= \frac 12(\eta+\bar\eta) \\ &= \frac 12((x+ay)+\overline{x+ay}) \\ &= \frac 12((1+a)\cdot(3+2a)^t+\overline{(1+a)\cdot(3+2a)^t}) \\ &= \frac 12((1+a)\cdot(3+2a)^t+(1-a)\cdot(3-2a)^t) \ . \end{aligned} And for $y$ we have a similar expression.

This explains the awful way of writing the solution - as obtained from sympy.

We can ask for the first few values in the solution:

K.<a> = QuadraticField(2)
var('X,Y')
for sol in solve_diophantine(X^2 - 2*Y^2 == -1, (X, Y)):
for t0 in [-5..5]:
X0, Y0 = K(sol[0].subs(t = t0)), K(sol[1].subs(t = t0))
eta = X0 + a*Y0
print(f'Solution X + Ya = {eta}')


I insisted to move the solution in the field $K=\Bbb Q(a)=\Bbb Q(\sqrt 2)$. Results:

Solution X + Ya = 985*a - 1393
Solution X + Ya = 169*a - 239
Solution X + Ya = 29*a - 41
Solution X + Ya = 5*a - 7
Solution X + Ya = a - 1
Solution X + Ya = a + 1
Solution X + Ya = 5*a + 7
Solution X + Ya = 29*a + 41
Solution X + Ya = 169*a + 239
Solution X + Ya = 985*a + 1393
Solution X + Ya = 5741*a + 8119
Solution X + Ya = -985*a + 1393
Solution X + Ya = -169*a + 239
Solution X + Ya = -29*a + 41
Solution X + Ya = -5*a + 7
Solution X + Ya = -a + 1
Solution X + Ya = -a - 1
Solution X + Ya = -5*a - 7
Solution X + Ya = -29*a - 41
Solution X + Ya = -169*a - 239
Solution X + Ya = -985*a - 1393
Solution X + Ya = -5741*a - 8119


Now by chance, every time the first component $X$ is odd.

We go back to the equation $$X^2 -2Y^2 = -31\ .$$ We need a special element having norm $-31$ in $K=\Bbb Q(a)=\Bbb Q(\sqrt 2)$. Here is it: $$\eta=(1+4a)\ ,\ N(\eta)=(1+4a)(1-4a)=1-4^2a^2=1-16\cdot 2=1-32=-31\ .$$ So a general solution $(X,Y)$ of $X^2 -2Y^2=31$ or rather $\eta=X+aY$ of $N(\eta)=-31$ is of the shape $\pm(1+a)^{2t}\cdot (1+4a)$. (Then we know how to isolate components, thus making things less operable.)

But now, to be able to get from $X$ to $x$, and from $Y$ to $k$ via $X=4x+3$ and $Y=2k$ some congruences modulo $4$ and $2$ have to be satisfied. This is the reason of further splitting the story into cases. The sage call delivers finally $16$ cases (for a related Pell-like equation).

This was the human approach for the given Pell-like-reshapeable equation. Well, it may happen that sympy does not recognize the minimal factor used above to reshape the equation. It uses some other stuff. We can even guess what it does, e.g. by computing some norms for elements in the output:

sage: K(13*sqrt(2) + 20).norm()
62
sage: K(19*sqrt(2) + 28).norm()
62

sage: K(470832*sqrt(2) + 665857).norm()
1
sage: K(470832*sqrt(2) + 665857).factor()
470832*a + 665857
sage: (1+a)^16
470832*a + 665857


So for humanly not profitable reasons, sympy writes the given equation in the form $N(\eta)=62$, and solves it like $\eta=\pm \eta_0\cdot (1+a)^{16t}$. Then isolates the components to make the result less understandable. But why $16t$? Well, we had some similar effect in the human solution. We had to insure that $X$ is of the shape $2x+3$. For obvious reasons, this was the case in the special case, but having and implementing a general algorithm...

Then the code had to exhibit all elements of norm $62$, that are not differing by a power of $(1+a)^{16}$.

Finally, here is the way to see numbers (first few solutions) from the given code:

var('x,k');
for xsol, ksol in solve_diophantine(2*x^2 + 3*x + 5 == k^2, (x, k)):
for t0 in [-3..3]:
x0, k0 = K(xsol.subs(t = t0)), K(ksol.subs(t = t0))
print(f'Solution: t0 = {t0} and (x0, k0) = {(x0, k0)}')


We get:

Solution: t0 = -3 and (x0, k0) = (-235876233405138476, -333579368323028527)
Solution: t0 = -2 and (x0, k0) = (-177122290076, -250488744823)
Solution: t0 = -1 and (x0, k0) = (-133004, -188095)
Solution: t0 = 0 and (x0, k0) = (4, -7)
Solution: t0 = 1 and (x0, k0) = (6458644, -9133903)
Solution: t0 = 2 and (x0, k0) = (8601067634596, -12163746499735)
Solution: t0 = 3 and (x0, k0) = (11454162183932917684, -16198631506138961887)
Solution: t0 = -3 and (x0, k0) = (9925624597054204, 14036952920178043)
Solution: t0 = -2 and (x0, k0) = (7453270444, 10540516147)
Solution: t0 = -1 and (x0, k0) = (5596, 7915)
Solution: t0 = 0 and (x0, k0) = (-116, 163)
Solution: t0 = 1 and (x0, k0) = (-153485636, 217061467)
Solution: t0 = 2 and (x0, k0) = (-204398969261204, 289063794464275)
Solution: t0 = 3 and (x0, k0) = (-272200968950560539236, 384950301980980455883)
Solution: t0 = -3 and (x0, k0) = (1374787436861472479, 1944243038589639482)
Solution: t0 = -2 and (x0, k0) = (1032344359871, 1459955394770)
Solution: t0 = -1 and (x0, k0) = (775199, 1096298)
Solution: t0 = 0 and (x0, k0) = (-1, 2)
Solution: t0 = 1 and (x0, k0) = (-1108129, 1567130)
Solution: t0 = 2 and (x0, k0) = (-1475709904321, 2086968960818)
Solution: t0 = 3 and (x0, k0) = (-1965223539520829281, 2779245782685214922)
Solution: t0 = -3 and (x0, k0) = (-1965223539520829281, -2779245782685214922)
Solution: t0 = -2 and (x0, k0) = (-1475709904321, -2086968960818)
Solution: t0 = -1 and (x0, k0) = (-1108129, -1567130)
Solution: t0 = 0 and (x0, k0) = (-1, -2)
Solution: t0 = 1 and (x0, k0) = (775199, -1096298)
Solution: t0 = 2 and (x0, k0) = (1032344359871, -1459955394770)
Solution: t0 = 3 and (x0, k0) = (1374787436861472479, -1944243038589639482)
Solution: t0 = -3 and (x0, k0) = (-6016943869156, -8509243623797)
Solution: t0 = -2 and (x0, k0) = (-4518196, -6389693)
Solution: t0 = -1 and (x0, k0) = (-4, -5)
Solution: t0 = 0 and (x0, k0) = (190124, -268877)
Solution: t0 = 1 and (x0, k0) = (253191791324, -358067265173)
Solution: t0 = 2 and (x0, k0) = (337179053192057996, -476843189972327645)
Solution: t0 = 3 and (x0, k0) = (449026065642355131292604, -635018751890450270168357)
Solution: t0 = -3 and (x0, k0) = (1191324496712159, 1684787260437638)
Solution: t0 = -2 and (x0, k0) = (894579839, 1265126942)
Solution: t0 = -1 and (x0, k0) = (671, 950)
Solution: t0 = 0 and (x0, k0) = (-961, 1358)
Solution: t0 = 1 and (x0, k0) = (-1278779041, 1808466662)
Solution: t0 = 2 and (x0, k0) = (-1702967950806529, 2408360372317310)
Solution: t0 = 3 and (x0, k0) = (-2267866261639086182881, 3207247224858365702678)
Solution: t0 = -3 and (x0, k0) = (-6943548011011756, -9819659768161553)
Solution: t0 = -2 and (x0, k0) = (-5213993404, -7373700185)
Solution: t0 = -1 and (x0, k0) = (-3916, -5537)
Solution: t0 = 0 and (x0, k0) = (164, -233)
Solution: t0 = 1 and (x0, k0) = (219403796, -310283825)
Solution: t0 = 2 and (x0, k0) = (292183107784964, -413209313725817)
Solution: t0 = 3 and (x0, k0) = (389104335200527143284, -550276628018752376513)
Solution: t0 = -3 and (x0, k0) = (35069318855059, 49595506348012)
Solution: t0 = -2 and (x0, k0) = (26333971, 37241860)
Solution: t0 = -1 and (x0, k0) = (19, 28)
Solution: t0 = 0 and (x0, k0) = (-32621, 46132)
Solution: t0 = 1 and (x0, k0) = (-43440843629, 61434630220)
Solution: t0 = 2 and (x0, k0) = (-57850779631518701, 81813357148750948)
Solution: t0 = 3 and (x0, k0) = (-77040693146164853541101, 108951993101930285334652)
Solution: t0 = -3 and (x0, k0) = (-204398969261204, -289063794464275)
Solution: t0 = -2 and (x0, k0) = (-153485636, -217061467)
Solution: t0 = -1 and (x0, k0) = (-116, -163)
Solution: t0 = 0 and (x0, k0) = (5596, -7915)
Solution: t0 = 1 and (x0, k0) = (7453270444, -10540516147)
Solution: t0 = 2 and (x0, k0) = (9925624597054204, -14036952920178043)
Solution: t0 = 3 and (x0, k0) = (13218093234633989953996, -18693206721131441839555)
Solution: t0 = -3 and (x0, k0) = (292183107784964, 413209313725817)
Solution: t0 = -2 and (x0, k0) = (219403796, 310283825)
Solution: t0 = -1 and (x0, k0) = (164, 233)
Solution: t0 = 0 and (x0, k0) = (-3916, 5537)
Solution: t0 = 1 and (x0, k0) = (-5213993404, 7373700185)
Solution: t0 = 2 and (x0, k0) = (-6943548011011756, 9819659768161553)
Solution: t0 = 3 and (x0, k0) = (-9246820095931294637596, 13076978388490120691657)
Solution: t0 = -3 and (x0, k0) = (40469963569358371, 57233171348531680)
Solution: t0 = -2 and (x0, k0) = (30389380579, 42977074168)
Solution: t0 = -1 and (x0, k0) = (22819, 32272)
Solution: t0 = 0 and (x0, k0) = (-29, 40)
Solution: t0 = 1 and (x0, k0) = (-37643741, 53236288)
Solution: t0 = 2 and (x0, k0) = (-50130695903261, 70895510037592)
Solution: t0 = 3 and (x0, k0) = (-66759749564076676829, 94412543254148556400)
Solution: t0 = -3 and (x0, k0) = (-50130695903261, -70895510037592)
Solution: t0 = -2 and (x0, k0) = (-37643741, -53236288)
Solution: t0 = -1 and (x0, k0) = (-29, -40)
Solution: t0 = 0 and (x0, k0) = (22819, -32272)
Solution: t0 = 1 and (x0, k0) = (30389380579, -42977074168)
Solution: t0 = 2 and (x0, k0) = (40469963569358371, -57233171348531680)
Solution: t0 = 3 and (x0, k0) = (53894417064774125296099, -76218215549195540625352)
Solution: t0 = -3 and (x0, k0) = (-57850779631518701, -81813357148750948)
Solution: t0 = -2 and (x0, k0) = (-43440843629, -61434630220)
Solution: t0 = -1 and (x0, k0) = (-32621, -46132)
Solution: t0 = 0 and (x0, k0) = (19, -28)
Solution: t0 = 1 and (x0, k0) = (26333971, -37241860)
Solution: t0 = 2 and (x0, k0) = (35069318855059, -49595506348012)
Solution: t0 = 3 and (x0, k0) = (46702302889720705939, -66047030140699210708)
Solution: t0 = -3 and (x0, k0) = (337179053192057996, 476843189972327645)
Solution: t0 = -2 and (x0, k0) = (253191791324, 358067265173)
Solution: t0 = -1 and (x0, k0) = (190124, 268877)
Solution: t0 = 0 and (x0, k0) = (-4, 5)
Solution: t0 = 1 and (x0, k0) = (-4518196, 6389693)
Solution: t0 = 2 and (x0, k0) = (-6016943869156, 8509243623797)
Solution: t0 = 3 and (x0, k0) = (-8012848387763696404, 11331878863214808365)
Solution: t0 = -3 and (x0, k0) = (-1702967950806529, -2408360372317310)
Solution: t0 = -2 and (x0, k0) = (-1278779041, -1808466662)
Solution: t0 = -1 and (x0, k0) = (-961, -1358)
Solution: t0 = 0 and (x0, k0) = (671, -950)
Solution: t0 = 1 and (x0, k0) = (894579839, -1265126942)
Solution: t0 = 2 and (x0, k0) = (1191324496712159, -1684787260437638)
Solution: t0 = 3 and (x0, k0) = (1586503510813642529471, -2243654781745183524590)
Solution: t0 = -3 and (x0, k0) = (8601067634596, 12163746499735)
Solution: t0 = -2 and (x0, k0) = (6458644, 9133903)
Solution: t0 = -1 and (x0, k0) = (4, 7)
Solution: t0 = 0 and (x0, k0) = (-133004, 188095)
Solution: t0 = 1 and (x0, k0) = (-177122290076, 250488744823)
Solution: t0 = 2 and (x0, k0) = (-235876233405138476, 333579368323028527)
Solution: t0 = 3 and (x0, k0) = (-314119682292713457139004, 444232314906683123060455)


As said before, this way to see (and split) the solutions is not helpful without number theoretic aid. It is always a good idea to search for the (humanly) better solution. For instance, the last entry above is better seen as follows...

K.<a> = QuadraticField(2)

x0, k0 = -314119682292713457139004, 444232314906683123060455
X0, Y0 = 4*x0 + 3, 2*k0

print(f'x0 = {x0}')
print(f'k0 = {k0}\n')
print(f'X0 = 4 x0 + 3 is {X0}')
print(f'Y0 = 2 k0      is {Y0}\n')
print(f'The norm of X0 + Y0 a is {(X0 + Y0*a).norm()}')


And we get:

x0 = -314119682292713457139004
k0 = 444232314906683123060455

X0 = 4 x0 + 3 is -1256478729170853828556013
Y0 = 2 k0      is 888464629813366246120910

The norm of X0 + Y0 a is -31


In fact:

sage: X0 + a*Y0 == (1+4*a)*(1-a)^62
True


Both values are:

sage: X0 + a*Y0
888464629813366246120910*a - 1256478729170853828556013
sage: (1+4*a)*(1-a)^62
888464629813366246120910*a - 1256478729170853828556013

more