Translating some of @gnasher729's suggestions to Mathematica gives indeed a superfast recipe:
$MaxExtraPrecision = 10^3;
f[x_Integer] := Module[{y, i, z, s},
y = Floor[Sqrt[x]];
Catch[
i = 0;
While[True,
z = (2 y + i)^2 - 4 x;
If[z >= 0 && IntegerQ[Sqrt[z]],
s = i/2 + y + {1, -1} Sqrt[z]/2;
If[IntegerQ[s[[1]]] && IntegerQ[s[[2]]], Throw[s]]];
i++]]]
f[15]
(* {5, 3} *)
f[16]
(* {4, 4} *)
f[17]
(* {17, 1} *)
n = 80;
f[NextPrime[10^n] NextPrime[10^n + 10^6]] // AbsoluteTiming
(* {0.010925,
{100000000000000000000000000000000000000000000000000000000000000000000000001000507,
100000000000000000000000000000000000000000000000000000000000000000000000000000129}} *)
f[NextPrime[10^n] NextPrime[10^n + 10^12]] // AbsoluteTiming
(* {0.007086,
{100000000000000000000000000000000000000000000000000000000000000000001000000000191,
100000000000000000000000000000000000000000000000000000000000000000000000000000129}} *)
Things get a lot slower once the difference between the two numbers becomes larger than their square root:
f[NextPrime[10^n] NextPrime[10^n + 10^43]] // AbsoluteTiming
(* {132.132,
{100000000000000000000000000000000000010000000000000000000000000000000000000000357,
100000000000000000000000000000000000000000000000000000000000000000000000000000129}} *)
The remaining point is a fast method to determine IntegerQ[Sqrt[z]]
for very large values of z
. Using @MichaelE2's fast perfect-square test gives a 150-fold speedup, which is however no match for the exponentially increasing difficulty with factor-spacing:
$MaxExtraPrecision = 10^3;
f[x_Integer] := Module[{y, i, z, s},
y = Floor[Sqrt[x]];
Catch[
i = 0;
While[True,
z = (2 y + i)^2 - 4 x;
If[z >= 0 && FractionalPart[Sqrt[z + 0``1]] == 0,
s = i/2 + y + {1, -1} Sqrt[z]/2;
If[IntegerQ[s[[1]]], Throw[s]]];
i++]]]
n = 80;
f[NextPrime[10^n] NextPrime[10^n + 10^43]] // AbsoluteTiming
(* {0.894308,
{100000000000000000000000000000000000010000000000000000000000000000000000000000357,
100000000000000000000000000000000000000000000000000000000000000000000000000000129}} *)
f[NextPrime[10^n] NextPrime[10^n + 10^44]] // AbsoluteTiming
(* {92.1891,
{100000000000000000000000000000000000100000000000000000000000000000000000000000179,
100000000000000000000000000000000000000000000000000000000000000000000000000000129}}
Sqrt[n]
(avoidingN
because of its built in definition), which, since the prime counting function is about $\text{li}(x)$, is approximately equal toPrime[Floor[LogIntegral[Sqrt[n]]]]
? Then work your way downPrime[k]
manually (as long as you start aboveSqrt[n]
)? No idea if this is remotely a good idea, though, since I don't know how mathematica deals withPrime
. But you could probably parallelize it easily, at least! :P $\endgroup$ – thorimur Mar 3 at 8:28