Finding Symmetry Group $S_3$ in a function

There's a little ambiguity as to the type of functions you are considering. For example, you seem ok with allowing the function to have some singularities at zero as your $z+1/z$ example indicates. The group $S_3$ acts naturally on $\mathbf{C} \cup \{\infty\}$ via the following rational functions:

$$\Sigma = \left\{z, \ 1/z, \ 1-z,\frac{1}{1-z}, \ 1 - \frac{1}{z}, \ \frac{z}{z-1}\right\}$$

In particular, if $h(z)$ is any function $h: \mathbf{C} \cup \{\infty\} \rightarrow \mathbf{C} \cup \{\infty\}$ then

$$f(z) = h(z) + h(1/z) + h(1-z) + h\left(\frac{1}{1-z}\right) + h\left(1-\frac{1}{z}\right) + h\left(\frac{z}{z-1}\right)$$

will be invariant under $\Sigma$. It could be the case that $f(z)$ is invariant under more symmetries, of course. For example, if $h(z) = z$, then $f(z) = 3$.

On the other hand, if $h(z) = z^2 + c$ for any constant $c$, then $f(z)$ is a non-trivial rational function. Moreover, one finds that (in this case)

$$f(x) - f(y) = \frac{2(x-y)(x+y-1)(xy - 1)(1-x+xy)(1-y+xy)(-x-y+xy)}{(x-1)^2 x^2 (y-1)^2 y^2}.$$

Assuming that $y \in \mathrm{Sym}(f)$, the numerator is zero, and so (under very weak continuity hypotheses) one of the six factors in the numerator are zero, leading to $y \in \Sigma$. So it seems that $f(z)$ is a suitable function in your case. A particularly nice choice of constant $c$ is $c = -7/4$, in which case $f(2) = 0$, and so

$$f(z) = f(z) - f(2) = \frac{(z-2)^2 (z+1)^2 (2z - 1)^2}{2 z^2 (z-1)^2}.$$

In this case, the square-root of this function is invariant under the even elements of $\Sigma = S_3$ and sent to its negative under the odd elements.

A slightly more general nice family (but no longer a square) is given (for a parameter $t$) by

$$f(x) = \frac{2(x-t)(x+t-1)(xt - 1)(1-x+xt)(1-t+xt)(-x-t+xt)}{(x-1)^2 x^2 (t-1)^2 t^2}.$$

I might as well add a complete list of such examples coming from polynomials. Suppose that $f(x)$ is a polynomial, and $y \in \mathrm{Sym}(f)$. Then we must have $f(x) - f(y) = 0$. But $f(x) - f(y)$ is a rational function in $y$, and so has a finite number of algebraic solutions. If we insist that our functions are entire functions on $\mathbf{C} \cup \{\infty\}$, then this forces $y$ to be a rational function (other algebraic functions will not be single valued), and (by degree considerations) a function of the form:

$$y = \frac{a x + b}{c x + d}.$$

The choice of constants is only well defined up to scaling. This gives an injective map:

$$\mathrm{Sym}(f) \rightarrow \mathrm{PGL}_2(\mathbf{C}).$$

The finite subgroups of the right hand side are well known, and so, in particular, we deduce:

Claim: Let $f$ be a rational function. Then $\mathrm{Sym}(f)$ is either cyclic, dihedral, or one of the exceptional groups $A_4$, $S_4$, and $A_5$.

Cyclic examples are easy to construct. Let $f(x) = x^n$, and then $\mathrm{Sym}(f)$ consists of $y = \zeta x$ for an $n$th root of unity $\zeta$. This corresponds to the map:

$$a \in \mathbf{Z}/n \mathbf{Z} \mapsto \left( \begin{matrix} \zeta^a & 0 \\ 0 & 1 \end{matrix} \right) \in \mathrm{PGL}_2(\mathbf{C}).$$

Naturally one can also take $f(x) = h(x^n)$ for a generic rational function $h(x)$. Note that other examples (such as $f(x) = x + x^{-1}$) can be obtained from these examples by suitable change of variables, namely, because

$$\left( \begin{matrix} 1 & - 1 \\ 1 & 1 \end{matrix} \right) \left( \begin{matrix} 0 & 1 \\ 1 & 0 \end{matrix} \right) \left( \begin{matrix} 1 & - 1 \\ 1 & 1 \end{matrix} \right)^{-1} = \left( \begin{matrix} -1 & 0 \\ 0 & 1 \end{matrix} \right),$$ and we find that $$h(x) = x + \frac{1}{x}, \qquad h\left(\frac{x-1}{x+1}\right) = g(x^2), \qquad g(x) = 2 \cdot \frac{x+1}{x-1}.$$

Note that the dihedral representation of $D_{2n}$ inside $\mathrm{PGL}_2(\mathbf{C})$ is given by the image of $\mathbf{Z}/n \mathbf{Z}$ together with the matrix

$$\left( \begin{matrix} 0 & 1 \\ 1 & 0 \end{matrix} \right) ,$$

Hence we can write down the examples

$$f(x) = h\left(x^n + \frac{1}{x^n}\right),$$

for a generic function $h$ (taking $h(x) = x$ will do). Here $\mathrm{Sym}(f)$ is generated by $x \mapsto \zeta x$ and $x \mapsto 1/x$.

One can construct the other examples in a similar manner. For fun, I computed an example with $\mathrm{Sym}(f) = A_4$. The group $A_4$ has (several) projective representations

$$A_4 \rightarrow \mathrm{PGL}_2(\mathbf{C})$$

realized by $2$-dimensional representations of the Schur cover $\mathrm{SL}_2(\mathbf{F}_3)$. One such example maps the non-trivial elements of the Klein $4$-subgroup $K$ to

$$K \setminus \{e\} = \left\{\left( \begin{matrix} 0 & 1 \\ 1 & 0 \end{matrix} \right), \left( \begin{matrix} -1 & 0 \\ 0 & 1 \end{matrix} \right), \left( \begin{matrix} 0 & -1 \\ -1 & 0 \end{matrix} \right)\right\},$$ and this group is normalized by the order three (in $\mathrm{PGL}_2$) element $$\left( \begin{matrix}i & -i \\ 1 & 1 \end{matrix} \right)$$ Writing down the corresponding $12$ elements of $A_4$ and letting $$f(z) = \sum_{A_4 \subset \mathrm{PGL}_2(\mathbf{C})} h(\gamma z),$$ doing a calculation as above with $h(z) = z^2 + c$, one finds that

$$\begin{aligned} f(x) - f(y) = & \ 2(x - y)(x + y)(-1 + xy)(1 + xy)(-i - ix - y + xy)(i + ix - y + xy) (-i - x - iy + xy)\\ \times & \ \frac{(i + x - iy + xy)(i - x + iy + xy)(-i + x + iy + xy) (i - ix + y + xy)(-i + ix + y + xy)}{(-1 + x)^2x^2(-i + x)^2(i + x)^2 (1 + x)^2(-1 + y)^2y^2(-i + y)^2(i + y)^2(1 + y)^2} \end{aligned} $$

Since translating $f(x)$ preserves the symmetry group, one can (for example) choose $f(x)$ to vanish at $x = y$ for any fixed $y$, and then $f(x)=f(x) - f(y)$ as above. For example, if $y = 2$, then

$$450 f(x) = \frac{(-2 + x) (2 + x) (-1 + 2 x) (1 + 2 x) (9 + x^2) (5 - 6 x + 5 x^2) (5 + 6 x + 5 x^2) (1 + 9 x^2)}{(-1 + x)^2 x^2 (1 + x)^2 (1 + x^2)^2}.$$

I will assume that symmetry groups comprise of functions that are homeomorphisms, so that it can actually be made into a group.

Some background: given a triangle $ABC$ in $\mathbb{C}$ (identified with the plane), any point $z\in\mathbb{C}$ can be written as a linear combination of the three vertices, i.e. $$ z = aA+bB+cC $$ for some $a,b,c\in\mathbb{R}$. If we require that $a+b+c = 1$, then this combination is unique. The triple $(a,b,c)$ is called the barycentric coordinates of $z$ with respect to $ABC$. It is not hard to check that the correspondence $$\mathbb{C}\rightarrow\{(a,b,c)\in\mathbb{R}^3:a+b+c=1\}$$ given by the barycentric coordinates is bijective. Note that two barycentric coordinates are enough to specify a point, since the third can be derived from the equation $a+b+c=1$.

Suppose now that $ABC$ is an equilateral triangle of side length $\frac{2}{\sqrt{3}}$, so that the height from any vertex to the opposite side is $1$. Then the barycentric coordinate $a$ is the signed distance from $z$ to line $BC$, and similarly for $b$ and $c$. Here, the sign of the distance is positive if $z$ and the triangle lie on the same side of the line, and negative if they lie on opposite sides. For example, if $A = \frac{2}{3}$, $B = -\frac{1}{3} + \frac{1}{\sqrt{3}}i$, $C = -\frac{1}{3}-\frac{1}{\sqrt{3}}i$, and $z = x+iy$, then the barycentric coordinate $a$ is the signed distance of $z$ from $BC$, i.e. $a = x+\frac{1}{3}$. Similarly, you can check that $b = -\frac{1}{2}x+\frac{\sqrt{3}}{2}y + \frac{1}{3}$, and $c = -\frac{1}{2}x-\frac{\sqrt{3}}{2}y+\frac{1}{3}$.

Now, let $z_1$ and $z_2$ be two points in $\mathbb{C}$ such that the three barycentric coordinates for $z_1$ are all different, and that the barycentric coordinates for $z_2$ are some permutation of the coordinates of $z_1$. What can we say about $z_1$ and $z_2$? Well, the permutation of the three barycentric coordinates induces a permutation on the three sides of $ABC$, and the induced permutation is unique because the coordinates are different. Furthermore, for each permutation of the three sides, there is a unique symmetry of the plane which carries out the permutation, and this group of symmetries is isomorphic to $S_3$. So $z_1$ and $z_2$ have permuted coordinates if and only if there is a symmetry of the triangle which takes $z_1$ to $z_2$, assuming that the coordinates are all different. If at least two of the coordinates of $z_1$ are the same, then $z_1$ lies on one of the bisectors of the triangle, and so if the coordinates of $z_2$ are some permutation of those of $z_1$, then $z_2$ also lies on some bisector. It's not hard to see that there is a symmetry of the triangle taking $z_1$ to $z_2$, though this symmetry need not be unique.

Now, let $g(z)$ be the value of the largest barycentric coordinate of $z$, and $h(z)$ be the value of the smallest. Let $$ f(z) = g(z) + ih(z).$$ Now, suppose $m$ satisfies $f(m(z)) = f(z)$ for all $z$. Then, for each $z$ the barycentric coordinates of $m(z)$ must be some permutation of those of $z$, since they share the largest and smallest coordinates. It follows from the discussion above that we can write $$ m(z) = \sigma_z(z) $$ where $\sigma_z$ is some symmetry of the triangle whose choice may depend on $z$. It suffices to show that we can choose $\sigma_z$ to be the same for all $z$, and this is where the assumption of continuity is important.

Note that $m$ must map bisectors to bisectors, and that the three bisectors divided the plane into six regions, i.e. the complement of the three bisectors in the plane has six connected components. Since $m$ is continuous, and has continuous inverse, it must map each of the six regions to another of the six regions. Furthermore, for any point not on the bisectors, the six symmetries of the triangle send the point to six different destinations, with one destination in each of the six regions. In particular, for each pair of regions there is a unique symmetry sending the first region to the second. Hence, the choice of $\sigma_z$ is constant on each of the regions. Applying a continuity argument to the boundary of each region, i.e. the bisectors, shows that the choice of the symmetry must be consistent between any two adjacent regions, and hence the choice of the symmetry must be constant across all six regions. It is now not hard to check that this choice is consistent on the bisectors as well.

Hence, the homeomorphisms $m$ which satisfy $f\circ m = f$ on $\mathbb{C}$ are precisely the symmetries of the equilateral triangle, i.e. $S_3$.

Now, I know you asked for an elementary function. If you count absolute value and conjugation as elementary operations, then I can show that $f$ is elementary. This is because the barycentric coordinates are elementary functions of the real and imaginary parts (see above), and the real and imaginary parts of $z$ are elementary functions of $z$ (since $\Re(z) = \frac{1}{2}(z+\bar{z})$ and $\Im(z) = \frac{1}{2i}(z-\bar{z})$). Furthermore, the $\max$ and $\min$ functions are elementary, since $\max(a,b) = \frac{a+b+|a-b|}{2}$ and $\max(a,b,c) = \max(\max(a,b),c)$, and similarly for $\min$. Hence, $g$ and $h$, as the maximum and minimum of the barycentric coordinates, are elementary, so $f$ is as well.