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^*$.

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...".