In[1]:= (* This problem was posted on LinkedIn by Songxin (Derek)Liang
16^(x^2-y)+16^(y^2-x) == 1 *)
In[2]:= (* Solve and Reduce cannot get a solution *)
Solve[16^(x^2 - y) + 16^(y^2 - x) == 16, {x, y}]
During evaluation of In[2]:= Solve::nsmet: This system cannot be solved with the methods available to Solve.
Out[2]= Solve[16^(x^2 - y) + 16^(-x + y^2) == 16, {x, y}]
In[3]:= Solve[16^(x^2 - y) + 16^(y^2 - x) == 16, {x, y}, Reals]
During evaluation of In[3]:= Solve::nsmet: This system cannot be solved with the methods available to Solve.
Out[3]= Solve[16^(x^2 - y) + 16^(-x + y^2) == 16, {x, y}, Reals]
In[4]:= Reduce[16^(x^2 - y) + 16^(y^2 - x) == 16, {x, y}]
During evaluation of In[4]:= Reduce::nsmet: This system cannot be solved with the methods available to Reduce.
Out[4]= Reduce[16^(x^2 - y) + 16^(-x + y^2) == 16, {x, y}]
In[5]:= Reduce[16^(x^2 - y) + 16^(-x + y^2) == 16, {x, y}, Reals]
During evaluation of In[5]:= Reduce::nsmet: This system cannot be solved with the methods available to Reduce.
Out[5]= Reduce[16^(x^2 - y) + 16^(-x + y^2) == 16, {x, y}, Reals]
In[6]:= (* Adding two extra variables to the problem gives a very \
complex solution (not show here) which contains b as a parameter *)
In[7]:= sln =
Solve[{x^2 - y == a, y^2 - x == b, 16^a + 16^b == 1}, {x, y, a, b}];
During evaluation of In[7]:= Solve::ifun: Inverse functions are being used by Solve, so some solutions may not be found; use Reduce for complete solution information.
During evaluation of In[7]:= Solve::svars: Equations may not give solutions for all "solve" variables.
In[8]:= (* Setting b = -1/4 gives the solution when the two terms \
are equal since 16^(-1/4) = 1/2*)
sln1 = sln /. b -> -1/4 // N
Out[8]= {{x -> -0.5 + 1. I, y -> -0.5 - 1. I,
a -> -0.25}, {x -> -0.5 - 1. I, y -> -0.5 + 1. I,
a -> -0.25}, {x -> 0.5, y -> 0.5, a -> -0.25}, {x -> 0.5, y -> 0.5,
a -> -0.25}}
In[9]:= (* Testing shows the solutions are correct*)
16^(x^2 - y) + 16^(-x + y^2) - 1 /. sln1
Out[9]= {4.44089*10^-16 + 0. I, 4.44089*10^-16 + 0. I, 6.66134*10^-16,
5.55112*10^-16}
In[10]:= (* Setting b = -1 gives another solution *)
sln2 = sln /. b -> -1 // N
Out[10]= {{x -> -0.719587 + 0.942605 I, y -> -0.347422 - 1.35657 I,
a -> -0.0232774}, {x -> -0.719587 - 0.942605 I,
y -> -0.347422 + 1.35657 I,
a -> -0.0232774}, {x -> 0.719587 - 0.440069 I,
y -> 0.347422 - 0.633336 I,
a -> -0.0232774}, {x -> 0.719587 + 0.440069 I,
y -> 0.347422 + 0.633336 I, a -> -0.0232774}}
In[11]:= (* Testing shows these solutions are also correct*)
16^(x^2 - y) + 16^(-x + y^2) - 1 /. sln2
Out[11]= {-1.66533*10^-16 + 6.34877*10^-16 I,
5.55112*10^-17 + 6.34877*10^-16 I, 5.55112*10^-17 - 2.59722*10^-16 I,
1.17961*10^-16 - 2.69342*10^-16 I}