# Brute Force for solving Diophantine equation

I'm trying to solve the following Diophantine equation using Mathematica: $$2^x=5^y+3$$

I have a code,but is very slow. Can we speed up ?

f[x_?IntegerQ, y_?IntegerQ] := 2^x - 5^y;
M = 2000;
sol = Thread[{x, y} -> #] & /@ (Most /@
Select[Flatten[Table[{x, y, f[x, y] == 3}, {x, 0, M}, {y, 0, M}],
1], #[[3]] &]) // AbsoluteTiming

(* {50.2135, {{x -> 2, y -> 0}, {x -> 3, y -> 1}, {x -> 7, y -> 3}}}*)

Another code with For function:

For[i = 0, i <= M, i++,
For[j = 0, j <= M, j++,
If[f[i, j] == 3, Print[{i, j}]]]] // AbsoluteTiming

(*{49.4506, Null}*)

Calculating time is about 50 second.

In Maple with M=2000,Maple can solve this with about 2.094 second.

• Flatten[ParallelTable[If[2^x == 3 + 5^y, {x, y}, Nothing], {x, 0, 2000}, {y, 0, 2000}], 1] does about 10 seconds on my machine. I'm not sure if Maple will use multiple cores by default when you loop that way. Jan 3, 2020 at 12:56
• Well, the simplest trick is that $y$ should not be exceed $x/2+1$. Jan 3, 2020 at 13:27
• You are comparing an exact computation with a floating point computation. That should give one possible speed-up. Jan 3, 2020 at 14:23

In version 10.1 on my circa 2011 PC:

AbsoluteTiming[
M = 2000;
r = Range[0, M];
xx = 2^r;
yy = -(5^r);
grid = # + yy & /@ xx;
Position[Round@grid, 3] - 1
]
{3.3214, {{2, 0}, {3, 1}, {7, 3}}}

By luck it seems in this instance Round isn't actually needed:

AbsoluteTiming[
M = 2000;
r = Range[0, M];
xx = 2^r;
yy = -(5^r);
grid = # + yy & /@ xx;
Position[grid, 3] - 1
]
{1.70207, {{2, 0}, {3, 1}, {7, 3}}}

If actual "brute force" is not the goal, then perhaps:

RepeatedTiming[
M = 2000;
r = Range[0, M];
sol = 2^r ⋂ 5^r + 3
]

x -> Log[2, sol]
y -> Log[5, sol - 3]
{0.00485, {4., 8., 128.}}

x -> {2., 3., 7.}

y -> {0., 1., 3.}
• neat answer! Happy New Year. Jan 6, 2020 at 6:51
• @yarchik "Brute force" usually means a computational attack with no analytical work. Usually this only works in smaller cases, and one must be more intelligent (e.g. your approach) as the problem scale increases. Jan 6, 2020 at 12:39
• @Mariusz Sorry, bad code copy. I am interested in your intent with this question; were you just using this as a toy example for a couple of loops, or were you actually looking for a better approach to finding solutions? I interpreted it as the former. Jan 6, 2020 at 14:18
• It does not seem important since there are no solutions with y > 22, but the floating-point representation is not faithful for y > 22, because 5^y has a round-off error of at least 1, which error increases with y. There is no such problem for 2^x because FP is binary. Suppose X = 2^x were equal to 5^23 + 3 exactly and Y = 523. Then X = 5^23 + 3; Y = 5^23; X - Y yields 4., not 3.. -- One possible fix is With[{prec = 1 + r*Log10[5.]}, sol = SetPrecision[2^r, Infinity] ⋂ SetPrecision[(N[5, #] & /@ prec)^r + 3, Infinity] ]. Jan 6, 2020 at 15:59
• @yarchik One might add to Mr.Wizard's definition that in brute force the computational attack usually makes an exhaustive check of all cases in a given search space, which in the OP is given as M * M pairs of integers with M = 2000. Jan 6, 2020 at 16:08

Maybe double loop is not necessary:

Reap[Do[If[IntegerQ[Log[2, 5^y + 3]], Sow[y]], {y, 100000}]]

Takes around 12 seconds.

• Ok, but this is nowhere near the scale of the original question ($2000\times2000$). How long would your solution take compared to that? Jan 4, 2020 at 5:49
• @MarcoB Right, my solution is 50 times faster! 100000 vs 2000. Jan 4, 2020 at 6:15
• Can squeeze out a bit more with finite precision using With[{l5y = Log[2, N[5, 50]^y + 3]}, If[Round[l5y] == l5y, Sow[{Round[l5y], y}]]] inside that loop. Jan 4, 2020 at 20:45
• @DanielLichtblau That is an order of magnitude improvement! Jan 5, 2020 at 12:50

This is an incomplete answer to a question that was not asked. One can show that there are no solutions in intervals in a way that is much faster than exhaustive search. I'll give the general idea but I'm not going to make a rigorous proof.

First we rewrite by taking the base-2 logarithms of both sides.

x = log_2(5^y+3) = log_2(5^y (1+3/5^y)) = y log_2(5) + 3/log(2) 5^(-y) + O(5^(-2*y)) (which assumes y>=1)

Since x is an integer we require thaty times the base-2 log of 5 be "close" to an integer, for integer-valued y. The place to look for candidates is the denominators of the convergents of continued fraction approximations to log_2(5).

In[1]:= cc = Convergents[ContinuedFraction[N@Log[2, 5]]]

(* Out[1]= {2, 7/3, 65/28, 137/59, 339/146, 1493/643, 9297/4004, \
20087/8651, 29384/12655, 49471/21306, 177797/76573, 227268/97879, \
4268621/1838395, 4495889/1936274, 31243955/13456039, \
35739844/15392313} *)

We already know that y=1 and y=3 give valid solutions. Let's suppose we checked by brute force that there are no solutions for y between 4 and 28. The next "good" candidate value is y=59. Well, in principle maybe some other value between 29 and 58 would work? Here is what we know from the convergents list. The fraction 137/59 is closer to log_2(5) than any other rational with denominator no larger than 59. So let's check how close it is.

In[2]:= N[137/59 - Log[2, 5], 20]

(* Out[2]= 0.00010580341772239789239 *)

So it is around 10^(-4) away. If y is in the range {29,59} then the first correction term is far smaller than this, since it is bounded below by 3/log(2)*5^(-29):

In[3]:= N[3/Log[2] 5^(-29), 20]

(* Out[3]= 2.3236230070198052257*10^-20 *)

What this shows is that no value for y in the range {29,30,...,59} will give an integer value for log_2(5^y+3). One can handle successive ranges by moving from one denominator to the next. For example, the maximum value our first-order correction gives in the ninth such interval is far to small to make up the difference between the closest rational approximation and log_2(5) with suitably bounded denominator in that interval.

{d9, d10} = Denominator /@ cc[[9 ;; 10]]
minDistToInteger9 = N[cc[[10]] - Log[2, 5], 20]
maxAvaliableInInterval9 = N[3/Log[2] 5^(-d9), 20]

(* Out[7]= {12655, 21306}

Out[8]= 4.8483327777503868560*10^-10

Out[9]= 1.4821457784228658814*10^-8845 *)

My guess is that there are no more solutions besides the ones already found. But I have no idea how one might prove that.

• @MariuszIwaniuk The proof is simple: $5^y = -3 (\text{mod} 1024) \Rightarrow c = 163 (\text{mod} 256)$. Thus, we have a contradiction mod 257. Jan 7, 2020 at 23:40
• @yarchik I'm not Mariusz, but...that mod 1024 vs 257 approach 's pretty neat. Jan 8, 2020 at 0:50
• Thank you! Actually I was inspired by your analysis. And I included Mariusz because I thought he may be interested to know that there are no other solutions even though the question was formulated as optimization of a brute force approach. Jan 8, 2020 at 1:00

Slight variation of @Mr.Wizard solution.

m = 2000;
range1 = Range[0, m];
range2 = Range[0, Ceiling[Log[5, 2^m - 3]]];
AbsoluteTiming[Position[Outer[Plus, 2^range1, -5^range2], 3] - 1]

{0.924582, {{2, 0}, {3, 1}, {7, 3}}}

This approach reduces time to ~20sec

m = 2000;
pts = Tuples[Range[0, m], 2];

slope1 = (Ceiling[Log[5, 2^m - 3] + 5] - 100)/m;
slope2 = (Floor[Log[5, 2^m - 3] - 5] + 100)/m;
sel = (#2 <= slope1 #1 + 100 && slope2 #1 - 100 <= #2) & @@@ pts;
pt4 = Pick[pts, sel];
Pick[pt4, UnitBox[2^#1 - 5^#2 - 3 & @@@ pt4], 1]

{{2, 0}, {3, 1}, {7, 3}}

We are basically crating an envelope around log function.

log = Table[{m, Log[5, 2^m - 3]}, {m, 2, 2000, 10}];
pt2 = {{0, 0}, {0, 100}, {m, Ceiling[Log[5, 2^m - 3] + 5]}, {m,
Floor[Log[5, 2^m - 3] - 5]}, {0, -100}};
Graphics[{FaceForm[], Blue, Point@pt2, Magenta, Point@pt4, Green,
Point@log, EdgeForm[Black], Polygon[pt2]}, Frame -> True]

• (1+)......Thanks Jan 6, 2020 at 10:10
Clear["Global`*"]

eqn = 2^x == 5^y + 3;

m = 2000;

AbsoluteTiming[{#[[1]], #[[-1, -1]]} & /@
(x /.
Solve[{eqn, 0 <= x <= m, 0 <= y <= m}, x, Integers])]

(* {2.3865, {{2, 0}, {3, 1}, {7, 3}}} *)
cf = Compile[{{m, _Integer}},
Do[
If[Abs[2^x - 5^y - 3] < 10.^-12, Print[{x, y}]],
{x, 0, m}, {y, 0, m}
], CompilationTarget -> "C", RuntimeOptions -> "Speed"
];

cf[2000] // AbsoluteTiming

{2,0}
{3,1}
{7,3}
{0.171435,Null}

Note that overflow will occur.

Do[
With[{y = Log[5, -3 + 2^x]}, If[Abs[FractionalPart@y] < 10.^-12, Print[{x, y}]]],
{x, 0, 10^5}
] // AbsoluteTiming

{4.73022, Null}