9
$\begingroup$

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. enter image description here

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

$\endgroup$
3
  • $\begingroup$ 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. $\endgroup$
    – Carl Lange
    Jan 3, 2020 at 12:56
  • $\begingroup$ Well, the simplest trick is that $y$ should not be exceed $x/2+1$. $\endgroup$
    – yarchik
    Jan 3, 2020 at 13:27
  • 2
    $\begingroup$ You are comparing an exact computation with a floating point computation. That should give one possible speed-up. $\endgroup$ Jan 3, 2020 at 14:23

6 Answers 6

13
$\begingroup$

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.}
$\endgroup$
13
  • 1
    $\begingroup$ neat answer! Happy New Year. $\endgroup$
    – ubpdqn
    Jan 6, 2020 at 6:51
  • 1
    $\begingroup$ @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. $\endgroup$
    – Mr.Wizard
    Jan 6, 2020 at 12:39
  • 1
    $\begingroup$ @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. $\endgroup$
    – Mr.Wizard
    Jan 6, 2020 at 14:18
  • 1
    $\begingroup$ 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 = 5`23. 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] ]. $\endgroup$
    – Michael E2
    Jan 6, 2020 at 15:59
  • 1
    $\begingroup$ @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. $\endgroup$
    – Michael E2
    Jan 6, 2020 at 16:08
10
$\begingroup$

Maybe double loop is not necessary:

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

Takes around 12 seconds.

$\endgroup$
4
  • $\begingroup$ Ok, but this is nowhere near the scale of the original question ($2000\times2000$). How long would your solution take compared to that? $\endgroup$
    – MarcoB
    Jan 4, 2020 at 5:49
  • $\begingroup$ @MarcoB Right, my solution is 50 times faster! 100000 vs 2000. $\endgroup$
    – yarchik
    Jan 4, 2020 at 6:15
  • 4
    $\begingroup$ 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. $\endgroup$ Jan 4, 2020 at 20:45
  • $\begingroup$ @DanielLichtblau That is an order of magnitude improvement! $\endgroup$
    – yarchik
    Jan 5, 2020 at 12:50
6
$\begingroup$

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.

$\endgroup$
3
  • $\begingroup$ @MariuszIwaniuk The proof is simple: $5^y = -3 (\text{mod} 1024) \Rightarrow c = 163 (\text{mod} 256)$. Thus, we have a contradiction mod 257. $\endgroup$
    – yarchik
    Jan 7, 2020 at 23:40
  • $\begingroup$ @yarchik I'm not Mariusz, but...that mod 1024 vs 257 approach 's pretty neat. $\endgroup$ Jan 8, 2020 at 0:50
  • 1
    $\begingroup$ 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. $\endgroup$
    – yarchik
    Jan 8, 2020 at 1:00
4
$\begingroup$

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}}}

Original answer:

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]

enter image description here

$\endgroup$
1
  • $\begingroup$ (1+)......Thanks $\endgroup$ Jan 6, 2020 at 10:10
2
$\begingroup$
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}}} *)
$\endgroup$
2
$\begingroup$
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}

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.