Solving Functional Equation Involving Substitution Suppose we seek a function of two variables $f(x,y)$ such that


*

*$f(x,y) = f(y,x)$

*$f(x,y) = f\left(\frac{1+y}{x}, ~ y \right)$
Are there known techniques for approaching such questions? I already know a solution (see below), so I'm not looking simply for a function which satisfies this. I'm looking for a method for finding such a function given these conditions.
$$ $$
Known Solution: I happen to already know a solution: $f(x,y) = \frac{(1+x)(1+y)(1+x+y)}{xy}$. But I wish to generalize this where I replace the expression $\frac{1+y}{x}$ in the the condition above by other expressions. So I'm wondering how I would find this function if I didn't already know it beforehand. 
$$ $$
Additional Background: This comes from the "Lyness recurrence" $x_{n+1} = \frac{1+x_n}{x_{n-1}}$. Clearly, if you set $x=x_n$ and $y=x_{n+1}$, then the function $f$ is a conserved quantity of the iteration. Indeed, this is "well-known" in the literature, which is why I already know a solution in this case.
 A: Consider the action of the dihedral group $$\begin{align}D_5&=\big\langle a,b\,\big|\,a^2=e\,,\,\,b^2=e\,,\text{ and }(ab)^5=e\big\rangle
\\&=\{e,a,b,ab,ba,aba,bab,abab,baba,ababa\}
\end{align}
$$ of order $10$ (where $e$ is the identity element of $D_5$) on $S:=\mathbb{R}_{>0}\times\mathbb{R}_{>0}$ given by
$$a\cdot(x,y)=(y,x)$$
and
$$b\cdot(x,y)=\left(\frac{y+1}{x},y\right)$$
for all $(x,y)\in S$.  Then, any function $f:S\to\mathbb{R}$ satisfying the two conditions

*

*$f(x,y)=f(y,x)$ and

*$f(x,y)=f\left(\dfrac{y+1}{x},y\right)$
for each $(x,y)\in S$ must be constant on each orbit of $S$ under $D_5$.
Now, the orbit of $(x,y)\in S$ consists of

*

*$e\cdot( x,y)=(x,y)$,

*$a\cdot(x,y)=(y,x)$,

*$b\cdot(x,y)=\left(\dfrac{y+1}{x},y\right)$,

*$ab\cdot(x,y)=\left(y,\dfrac{y+1}{x}\right)$,

*$ba\cdot(x,y)=\left(\dfrac{x+1}{y},x\right)$,

*$aba\cdot(x,y)=\left(x,\dfrac{x+1}{y}\right)$,

*$bab\cdot(x,y)=\left(\dfrac{x+y+1}{xy},\dfrac{y+1}{x}\right)$,

*$abab\cdot(x,y)=\left(\dfrac{y+1}{x},\dfrac{x+y+1}{xy}\right)$,

*$baba\cdot(x,y)=\left(\dfrac{x+y+1}{xy},\dfrac{x+1}{y}\right)$, and

*$ababa\cdot(x,y)=\left(\dfrac{x+1}{y},\dfrac{x+y+1}{xy}\right)$.

Define $s(x,y)$ and $m(x,y)$ as $x+y$ and $xy$, respectively, for each $(x,y)\in S$.  We have
$$\begin{align}
\sigma(x,y)&:=\dfrac14\,\sum_{g\in D_5}\,s\big(g\cdot(x,y)\big)\\&=\dfrac{xy(x+y)+(x^2+y^2)+2(x+y)+1}{xy}
\\&=\frac{(x+1)(y+1)(x+y+1)}{xy}-3\end{align}$$
and
$$\begin{align}
\mu(x,y)&:=\dfrac12\,\sum_{g\in D_5}\,m\big(g\cdot(x,y)\big)\\&=\small\dfrac{x^3y^3+xy(x^3+y^3)+xy(x^2+y^2)+(x^3+y^3)+xy(x+y)+2(x^2+y^2)+2xy+(x+y)}{x^2y^2}\end{align}$$
for all $(x,y)\in S$.

Conjecture. Every function $f:S\to\mathbb{R}$ with the required properties satisfies
$$f(x,y)=\varphi\big(\sigma(x,y),\mu(x,y)\big)$$
for all $(x,y)\in S$, where $\varphi:S\to\mathbb{R}$ is a fixed function.

Clearly, every $f:S\to\mathbb{R}$ in the form given by the conjecture satisfies the required properties.  It remains to show that there are no other functions, but this seems to be quite a difficult problem.
Edit.  The conjecture is confirmed to be true by SageMath.  That is, the set of solutions $(x,y)\in S$ to the system of equations
$$\left\{\begin{array}{c}\sigma(x,y)=\sigma(a,b)\,,\\\mu(x,y)=\mu(a,b)\,,\end{array}\right.$$
where $(a,b)\in S$, is the orbit of $(a,b)$ under $D_5$.  Below is the Sage code I used to verify the conjecture.  (At the end, the code will also return two extraneous solutions $(x,y)=(-1,0)$ and $(x,y)=(0,-1)$.)
def f(t) : return [t[1],t[0]]

def g(t) : return [((t[1]+1)/t[0]).simplify_full(),t[1]]

var('x')

var('y')

var('a')

var('b')

z=[x,y];

[x,y]

f(z)

g(z)

f(g(z))

g(f(z))

f(g(f(z)))

g(f(g(z)))

f(g(f(g(z))))

g(f(g(f(z))))

f(g(f(g(f(z)))))

def mm(t): return (t[0]*t[1]).simplify_full()

def m(t): return (mm(t)+mm(f(t))+mm(g(t))+mm(f(g(t)))+mm(g(f(t)))+mm(f(g(f(t))))+mm(g(f(g(t))))+mm(f(g(f(g(t)))))+mm(g(f(g(f(t)))))+mm(f(g(f(g(f(t))))))).simplify_full()

m(z)

m(g(z))

def ss(t): return (t[0]+t[1]).simplify_full()

def s(t): return (ss(t)+ss(f(t))+ss(g(t))+ss(f(g(t)))+ss(g(f(t)))+ss(f(g(f(t))))+ss(g(f(g(t))))+ss(f(g(f(g(t)))))+ss(g(f(g(f(t)))))+ss(f(g(f(g(f(t))))))).simplify_full()

s(z)

s(g(z))

solve([s([x,y])==s([a,b]),m([x,y])==m([a,b])],x,y)

