How did Ulam and Neumann find this solution? In the book "Chaos, Fractals and Noise - Stochastic Aspects of
Dynamics" from Lasota and Mackey the operator $P: L^1[0,1] \to L^1[0,1]$
$$ (Pf)(x) = \frac{1}{4\sqrt{1-x}} \left[ f\left(\frac{1}{2}\left(1-\sqrt{1-x}\right)\right) +  f\left(\frac{1}{2}\left(1+\sqrt{1-x}\right)\right)\right]$$
was considered, where $L^1$ denotes the space of Lebesgue integrable functions.
It is stated in the book that Ulam and Neumann showed that
$$ f^*(x) =\frac{1}{\pi\sqrt{x(1-x)}} $$
solves
$$ Pf^* = f^*,$$ in the article "On combination of stochastic and deterministic processes, Bull. Amer. Math. Soc. vol. 53 (1947) p. 1120.". I went to the library of my university and I actually found the book  "vol 53 of Bull. Amer. Math. Soc.", but this book only contains the abstract of the article, not the article itself.
So my question is: Where do I find the article - or how did they find $f^*$ ?
Edit: I am interested how Ulam and Neumann actually constructed $f^*$. 
 A: This is just a speculation so I'm not sure they found the solution this way, but this is 'doable' in my opinion. The idea comes from observing the $x$ that $S^n(x)$ changes its increasing\decreasing state. ($S(x)=4x(1-x)$ in the book)
$$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$$
Which can be rewritten as
$$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$$
Now assume $Pf=f$. From an intuitive viewpoint, the transformation $P$ sends $S^{-1}(x)$ to $x$. So we get the idea to define a function $g : [0,1]\rightarrow [0,1]$ as$$g(x)=f(\sin^2\frac{\pi x}{2})$$
Then the functional equation of $f$ becomes much simpler.
 $$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}}$$
It is now easy to see
$$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$$
Putting $h(x)=g(x)\sin\pi x$, we get $$\frac{1}{2}(h(\frac{x}{2})+h(1-\frac{x}{2}))=h(x)$$
We will now see that the only possible continuous function $h$ is the constant function. Using this equation repetitively, $$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)$$
If x is irrational, each intervals $(\frac{t}{2^n}, \frac{t+1}{2^n})$ contain exactly one of the $\frac{k}{2^{n-1}}\pm\frac{x}{2^n}$ or $\frac{x}{2^n}$ or $1-\frac{x}{2^n}$. (More specifically, $\frac{k}{2^{n-1}}+\frac{x}{2^n}\in(\frac{2k}{2^n}, \frac{2k+1}{2^n})$ and similarly for others.)
So consindering the partition $0, \frac{1}{2^n}, \frac{2}{2^n}, \cdots, 1$ the above is a riemann sum. Because $h$ is a continous function, the sum converges to the Riemann integral which makes $$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$.
So we finally have $h(x)=c$ and from this we can get $$f(x)=\frac{c}{\sqrt{x(1-x)}}$$
Putting $c=1/\pi$ makes the function $f$ a proabability density function on $[0,1]$.
Update. I've just found out that this was actually a method called the 'change of variables'. It is illustrated Ch6.5 in the OP's book and the idea is to transform $S$ to an exact function $T=g\circ S\circ g^{-1}$ in order to make it more easy to manipulate.(in this case, to the tent function)
According to the book's theorem 6.5.2, if we let $T:(0,1)\rightarrow (0,1)$ a measurable nonsingular transformation, $\phi\in D((a,b))$ a positive density, and $S:(a,b)\rightarrow (a,b)$ defined as $S=g^{-1}\circ S\circ g$ with $$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$.
The problem is that this idea was originally due to Ruelle[1977] and Pianigiani[1983] about 30 years later which means that it's not likely for this method that was used by Ulam and Neumann to find $f^*$. So I guess my answer becomes "Well, there is a method for finding such function. But I'm not sure how they found it...".
A: I don't have the reference for the derivation, but I can proof, that $f^*$ is indeed an eigenfunction of $P$ with eigenvalue $1$, i.e. $Pf^*=f^*$.
Proof
\begin{align}
(Pf^*)(x)&= \frac{1}{4\sqrt{1-x}} \Bigl( f\Bigl(\frac12 (1-\sqrt{1-x})\Bigr)+
 f\Bigl(\frac12 (1+\sqrt{1-x})\Bigr)\Bigr) 
\\
&=\frac{1}{4\sqrt{1-x}} \Bigl(
\frac1\pi  \frac{1}{\sqrt{\frac12 (1-\sqrt{1-x})(1-\frac12 (1-\sqrt{1-x})}}\Bigr)
+
\frac1\pi  \frac{1}{\sqrt{\frac12 (1+\sqrt{1-x})(1-\frac12 (1+\sqrt{1-x})}}\Bigr)
  \Bigr) 
\\
&=\frac{1}{4\pi\sqrt{1-x}}\Bigl(
 \frac{1}{\sqrt{\frac12 -\frac12 \sqrt{1-x})(\frac12 +\frac12 \sqrt{1-x})}}\Bigr)
+\frac{1}{\sqrt{\frac12 +\frac12 \sqrt{1-x})(\frac12 -\frac12 \sqrt{1-x})}}\Bigr)
\\ &=
\frac{1}{4\pi\sqrt{1-x}}\Bigl(
 \frac{1}{\sqrt{(\frac12)^2 -(\frac12 \sqrt{1-x})^2}}\Bigr)
+ \frac{1}{\sqrt{(\frac12)^2 -(\frac12 \sqrt{1-x})^2}}\Bigr)
\\ &=
\frac{2}{4\pi\sqrt{1-x}}\Bigl(
 \frac{1}{\sqrt{(\frac12)^2 -(\frac12 \sqrt{1-x})^2}}\Bigr)
\\ &=
\frac{2}{4\pi\sqrt{1-x}}\Bigl(
 \frac{1}{\sqrt{\frac14 -\frac14({1-x})}}\Bigr)
\\ &=
\frac{2}{\pi\sqrt{1-x}}\Bigl(
 \frac{1}{\sqrt{1-({1-x})}}\Bigr)
\\ &=
\frac{2}{2\pi\sqrt{1-x}}\Bigl(
 \frac{1}{\sqrt{x}}\Bigr)=f^*
\end{align}
Okay there are a lot of brackets missing/wrong and maybe this is of no use. The message was only, that it is trivial to proof, that $Pf^*=f^*$. Simply plug in the definition of $f^*$ and use binomial expansion for $(a-b)(a+b)$.
