Esto es sólo una especulación, así que no estoy seguro de que encontraron la solución de esta manera, pero esto es factible en mi opinión. La idea viene de la observación de la $x$ que $S^n(x)$ cambios de su creciente\disminuyendo estado. ($S(x)=4x(1-x)$ en el libro)
$$S(1/2)=1\\S^2(\frac{2-\sqrt{2}}{4})=S^2(\frac{2+\sqrt{2}}{4})=1\\S^3(\frac{2-\sqrt{2-\sqrt{2}}}{4})=S^3(\frac{2+\sqrt{2-\sqrt{2}}}{4})=S^3(\frac{2-\sqrt{2+\sqrt{2}}}{4})=S^3(\frac{2+\sqrt{2+\sqrt{2}}}{4})=1\\\vdots$$
La cual puede escribirse como
$$S(\sin^2\frac{\pi}{4})=1\\S^2(\sin^2\frac{\pi}{8})=S^2(\sin^2\frac{3\pi}{8})=1\\S^3(\sin^2\frac{\pi}{16})=S^3(\sin^2\frac{3\pi}{16})=S^3(\sin^2\frac{5\pi}{16})=S^3(\sin^2\frac{7\pi}{16})=1\\\vdots$$
Ahora suponga $Pf=f$. Desde un punto de vista intuitivo, la transformación de $P$ envía $S^{-1}(x)$$x$. Así que tenemos la idea de definir una función $g : [0,1]\rightarrow [0,1]$ $$g(x)=f(\sin^2\frac{\pi x}{2})$$
Entonces el funcional de la ecuación de $f$ se vuelve mucho más sencillo.
$$g(x)=f(\sin^2\frac{\pi x}{2})=(Pf)(\sin^2\frac{\pi x}{2})=\frac{f(\sin^2\frac{\pi x}{4})+f(\sin^2(\frac{\pi}{2}-\frac{\pi x}{4}))}{4\cos\frac{\pi x}{2}}=\frac{g(\frac{x}{2})+g(1-\frac{x}{2})}{4\cos\frac{\pi x}{2}}$$
Ahora es fácil de ver
$$g(\frac{x}{2})\sin\frac{\pi x}{2}+g(1-\frac{x}{2})\sin\pi(1-\frac{x}{2})=4g(x)\sin\frac{\pi x}{2}\cos\frac{\pi x}{2}=2g(x)\sin\pi x$$
Poner a $h(x)=g(x)\sin\pi x$, obtenemos $$\frac{1}{2}(h(\frac{x}{2})+h(1-\frac{x}{2}))=h(x)$$
Ahora veremos que la única función continua $h$ es la función constante. A partir de esta ecuación de forma repetitiva, $$h(x)=\frac{1}{2^n}\Big(\sum_{k=1}^{2^{n-1}-1}\big(h(\frac{k}{2^{n-1}}+\frac{x}{2^n})+h(\frac{k}{2^{n-1}}-\frac{x}{2^n})\big)+h(\frac{x}{2^n})+h(1-\frac{x}{2^n})\Big)$$
Si x es irracional, cada uno de los intervalos de $(\frac{t}{2^n}, \frac{t+1}{2^n})$ contienen exactamente una de las $\frac{k}{2^{n-1}}\pm\frac{x}{2^n}$ o $\frac{x}{2^n}$ o $1-\frac{x}{2^n}$. (Más específicamente, $\frac{k}{2^{n-1}}+\frac{x}{2^n}\in(\frac{2k}{2^n}, \frac{2k+1}{2^n})$ y de manera similar para los demás.)
Así consindering la partición $0, \frac{1}{2^n}, \frac{2}{2^n}, \cdots, 1$ la de arriba es una suma de riemann. Debido a $h$ es una función continua, la suma converge a la integral de Riemann que hace $$h(x)=\int_{0}^{1}h(t)dt=c$$ for all irrational $x$. Since the irrational numbers in $[0,1]$ are dense in $[0,1]$, $h(x)=c$ for all $[0,1]$ by the continuity of $h$.
Así que finalmente nos ha $h(x)=c$ y a partir de esto podemos sacar $$f(x)=\frac{c}{\sqrt{x(1-x)}}$$
Poner a $c=1/\pi$ hace la función de $f$ un proabability función de densidad de $[0,1]$.
La actualización. Me acabo de enterar que este era en realidad un método llamado el " cambio de variables. Esto es ilustrado en la Ch6.5 en la OP del libro y la idea es transformar $S$ a una función exacta $T=g\circ S\circ g^{-1}$ con el fin de hacer que sea más fácil de manipular.(en este caso, a la tienda de la función)
Según el libro del teorema 6.5.2, si dejamos $T:(0,1)\rightarrow (0,1)$ medibles nonsingular transformación, $\phi\in D((a,b))$ una densidad positiva, y $S:(a,b)\rightarrow (a,b)$ define como $S=g^{-1}\circ S\circ g$ con $$g(x)=\int_{a}^{x}\phi (y)dy$$ then $T$ is exact iff $S$ is statistically stable and $\phi$ is the measure invariant to $S$. So the method only requires to find such $\phi$.
El problema es que esta idea fue originalmente debido a la Ruelle[1977] y Pianigiani[1983] unos 30 años más tarde, lo que significa que no es probable para este método que fue utilizado por Ulam y Neumann para encontrar $f^*$. Así que supongo que mi respuesta es "Bueno, no es un método para encontrar dicha función. Pero no estoy seguro de cómo se encontró...".
