# One more solution of the Mordell equation

Please, help to find one more solution of the Mordell equation $$y^2=x^3+n,\quad n\in\mathbb{Z}$$ $$y^2=x^3-307$$

Using Solve I was able to go up to k=100000. How can this approach be extended towards $$k=10^9$$ and $$n=2\times 10^{13}$$? As a harder example, consider

$$y^2=x^3+8569$$

k = 100000;
n = 10 k;
Solve[y^2 ==x^3-307 && -k < x < k && 0 < y < n, {x, y}, Integers] // Timing

(*Out[1]= {3.78542, {{x -> 7, y -> 6}, {x -> 11, y -> 32}, {x -> 71, y -> 598}}}*)


This post is of relevance.

The motivation for this question is to better understand capabilities of Mathematica for searching solutions of Diophantine equations by a brute force approach and learning tricks to push the limits even further.

• x must be positive, too. So you can change boundary to 0<x<k for the first equation. I don't know this helps though. – OkkesDulgerci Jan 6 at 13:10
• @OkkesDulgerci I prefer to allow negative values of $x$ in the case of negative rhs. – yarchik Jan 6 at 13:34
• This list says there are 8 solutions total.. oeis.org/A081120/b081120.txt – OkkesDulgerci Jan 6 at 16:07
• For the second equation, lower bound of $x$ is: $x>-21$. oeis.org/A081119/b081119.txt says it has 8 solution (including symmetry of $y$). Solve find 6 of them. So one solution is missing. – OkkesDulgerci Jan 7 at 15:54

[Update: Improved second code.]

There is a system limit on Solve, which you can extend this way:

k = 1000000;
n = Ceiling[k^(3/2)];
With[{ropts = SystemOptions["ReduceOptions"]},
InternalWithLocalSettings[
SetSystemOptions[
"ReduceOptions" -> "SolveDiscreteSolutionBound" -> n],
Solve[x^3 - y^2 == 307 && -k < x < k && 0 < y < n, {x, y},
Integers],
SetSystemOptions[ropts]
]] // AbsoluteTiming
(*
{143.664,
{{x -> 7, y -> 6}, {x -> 11, y -> 32},
{x -> 71, y -> 598}, {x -> 939787, y -> 911054064}}}
*)


For speed using an exhaustive search over x: The code will work efficiently for machine integers (for solutions with x^3 less than 2^53, the limit on double-precision floating-point numbers to exactly represent an integer).

Block[{Part},
With[{x = #[[1]] + 1, y = #[[2]]},
Hold[
Pick[#[[All, 1 ;; 2]], #[[All, -1]], 0] &@
NestList[
With[{n = Sqrt[x^3 - 307.]},
If[FractionalPart@n == 0,
{x, Round[n], 0},
{x, y, 1}]
] &,
{Floor@CubeRoot@307., 1, 1},
1000000
]
]
]
] // ReleaseHold // AbsoluteTiming
(*  {0.36922, {{7, 6}, {11, 32}, {71, 598}, {939787, 911054064}}}  *)


If you want a brute-force approach to check the rectangular {x, y} space, keep in mind that for 0 <= x <= 10^6, the space has 10^15 pairs, which would take a long time for a GHz processor, or even a few thousand of them.

• Your second solution is very fast. Can you please explain briefly how it works. I do not even see 307 in there. – yarchik Jan 7 at 9:09
• @yarchik Actually, it was an alternate trial solution. It wasn't the one I meant to post. I got interrupted and was careless. The updated answer has the better one. (The old one looked for an increment to a correct solution, but the direct way is faster.) It's fast because it is auto-compiled by NestList and the number of iterations can be handled efficiently by modern memory management systems. As mentioned in the answer, this will break down if the search space goes beyond machine limits. – Michael E2 Jan 7 at 14:09

I used this code to find one. In every iteration I looked up 200k range.

m = 100000;
Total@Boole[IntegerQ /@ Sqrt[Range[8 m, 10 m]^3 - 307]]


1

And extracted the solution using

Position[IntegerQ /@ Sqrt[Range[8 m, 10 m]^3 - 307], True]


{{139788}}

This implies that

x=139788 - 1 + 8 m=939787

is a solution.

{x,y}={939787,911054064}

Since I used Brute-Force, next solution must be $$x>10,000,000$$

• I think 8 means $(x,\pm y)$ solutions. So, you really managed to find all for $k=307$. Nice ! – yarchik Jan 6 at 21:09
• @yarchik, OkkesDulgerci, This speeds up Okkes's method considerably, so that it is faster than Solve: Position[PossibleZeroQ@ FractionalPart@ Sqrt[Range[10^6]^3 - 307 + 01], True] – Michael E2 Jan 7 at 14:21
• @yarchik I guess you right. Michael this method is really clever. Did you try this approach for second equation? – OkkesDulgerci Jan 7 at 15:30
• It finds only x = -10, 23, 36, y = ±Sqrt[x^3+8569], up to x = 10^7 in 70 sec. The time complexity seems worse than $k \log k$. (It's not linear because the precision increases with $k$.) – Michael E2 Jan 7 at 17:10
• In fact, the space complexity & memory management is significant for larger values of k. This is faster: Last@Reap@Do[If[PossibleZeroQ@FractionalPart@Sqrt[x^3 + 8569 + 01], Sow@x], {x, -Floor@CubeRoot[8569.], 10^7}] // AbsoluteTiming` – Michael E2 Jan 7 at 17:31