Integers, or "Integers"?

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');
K.<a> = QuadraticField(2)
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

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,