Solving a diophantine equation with three variables

Question

Solve this equation using Mathematica:

$2^x + 113^y = z^2$, $\operatorname{gcd}(x, y) = 1$.

Answer 1

Dealing with diophantine equations after appropriate restriction of the possible solution space one could play with extension of ExhaustiveSearchMaxPoints.

SystemOptions["ReduceOptions" -> "ExhaustiveSearchMaxPoints"]

For an example when an extension appears crucial see e.g. Solving/Reducing equations in Z/pZ. Here I just set (I don't insist it is neccessary) :

SetSystemOptions["ReduceOptions" -> "ExhaustiveSearchMaxPoints" -> {1000, 1000000}];

and playing with Reduce I just found after roughly two minutes

Reduce[ 2^x + y^113 == z^2 && GCD[x, y] == 1 && 500 >= x > 0 && 
        32 >= y > 0 && z > 0, {x, y, z}, Integers]

 (x == 3 && y == 1 && z == 3) || 
   (x == 113 && y == 2 && z == 144115188075855872) || 
   (x == 229 && y == 4 && z == 31153781151208965771182977975320576) ||
   (x == 234 && y == 1 &&
      z == 166153499473114484112975882535043072) || 
   (x == 236 && y == 1 && z == 332306998946228968225951765070086144) ||
   (x == 238 && y == 1 && z == 664613997892457936451903530140172288) || 
   (x == 240 && y == 1 && z == 1329227995784915872903807060280344576) ||
   (x == 242 && y == 1 && z == 2658455991569831745807614120560689152) ||
   (x == 339 && y == 8 && 
     z == 1496577676626844588240573268701473812127674924007424) ||
   (x == 455 && y == 16 && 
     z == 32351936000580767753600418104423556808364573307048686977324332\
    2990592)

This result suggest that the only solutions are of the form y == 2^n for natural n. In fact, I've looked for solutions with higher bounds for x and y and all solutions were of the form y == 2^n. It seems to be easy to prove, nonetheless I put aside this step up to you. This can significantly simplify the search for solutions. Moreover it might be quite resonable to look for z expressed by its prime components. Now we can provide more solutions (it took a few minutes)

sols = {x, n, z} /. { ToRules @
 Reduce[2^x + (2^n)^113 == z^2 && GCD[x, 2^n] == 1 && 
         2000 >= x > 0 && 20 >= n >= 0 && z > 0, {x, n, z}, Integers]};
{#1, #2, FactorInteger @ #3} & @@@ sols // Column

{3, 0, {{3, 1}}}
{113, 1, {{2, 57}}}
{229, 2, {{2, 113}, {3, 1}}}
{234, 0, {{2, 117}}}
{236, 0, {{2, 118}}}
{238, 0, {{2, 119}}}
{240, 0, {{2, 120}}}
{242, 0, {{2, 121}}}
{339, 3, {{2, 170}}}
{455, 4, {{2, 226}, {3, 1}}}
{565, 5, {{2, 283}}}
{681, 6, {{2, 339}, {3, 1}}}
{791, 7, {{2, 396}}}
{907, 8, {{2, 452}, {3, 1}}}
{1017, 9, {{2, 509}}}
{1133, 10, {{2, 565}, {3, 1}}}
{1243, 11, {{2, 622}}}
{1359, 12, {{2, 678}, {3, 1}}}
{1469, 13, {{2, 735}}}
{1585, 14, {{2, 791}, {3, 1}}}
{1695, 15, {{2, 848}}}
{1811, 16, {{2, 904}, {3, 1}}}
{1921, 17, {{2, 961}}}

These are the only solutions under the above bounds. With this approach it is quite natural to generalize and find all possible solutions proving a simple theorem. At this point it is reasonable to stop since you haven't provided any piece of code.

Edit

It appears that the above analysis was performed for a different equation than the original one (see a comment by anderstood) i.e. we should do
something like this (it took more than 15 minutes to evaluate):

 Reduce[ 2^x + 113^y == z^2 && GCD[x, y] == 1 && 2000 >= x > 0 && 
         100 >= y > 0 && z > 0, {x, y, z}, Integers]

(x == 3 && y == 1 && z == 11) ||
(x == 9 && y == 1 && z == 25)   ||
(x == 242 && y == 1 && z == 2658455991569831745807614120560689152)

One can easily check that the first 2 solutions are correct however
there is an awful bug in the system behind the third solution. Namely z == 2658455991569831745807614120560689152 == 2^121 therefore 2^242 cannot be equal to 113 + (2^121)^2.

All of the following produce the same flawed results (in Mathematica versions 10.1, 10.3, 11.1)

Reduce[2^242 + 113 == z^2 && z > 0, {z}, Integers]
Reduce[2^242 + 113^1 == z^2 && z ∈ Integers && z > 0, z]
Solve[z^2 == 2^242 + 113 && z ∈ Integers && z > 0, z]

while

Solve[2^242 + 113^1 == z^2 && z > 0, z, Integers]

does not.

I'm going to leave this post as a warning sign (until the bug is present) for the comunity, although we've got only partial results.

Answer 2

Since Reduce and Solve don't work for higher values, simpliy test it directly, which goes quite fast.

First make a ContourPlot to get an impression, where solutions are to be expected.

ContourPlot3D[2^x + 113^y == z^2, {x, 0, 210}, {y, 0, 30}, {z, 0, 10^30}]

Solutions for given y are up to about x = 7*y and z = 10^y

Now do direct testing with y up to 100, x up to 700 and z up to very high z = Sqrt[2^700+113^100] = 2.2935*10^105

(sol = Flatten[
Table[If[
  GCD[x, y] == 1 && 
   IntegerQ[z = Sqrt[2^x + 113^y]], {{x, y, z}}, {}], {x, 0, 
  700}, {y, 0, 100}], 2]) // Timing

Quite fast you get the two solutions

(*     {62.641, {{3, 1, 11}, {9, 1, 25}}}     *)