Uniform distribution on a manifold

To generate a uniform distribution on a sphere $$S^n$$ in $$\mathbb R^{n+1}$$, we can normalize a vector whose entries are $$n+1$$ i.i.d normal random variables. If $$\rho$$ is a correlation, $$|\rho|<1$$, how can we generate a uniform distribution on the manifold $$M = \left \{ x, y \in S^n: x^Ty = \rho \right\}\ ?$$

It can be done in pretty much the same way as for a single vector by using the fact that if you fix $$x$$, then the conditional distribution of $$y-x$$ is uniform on the sphere of radius $$\sqrt{1-\gamma^2}$$ in the hyperplane perpendicular to $$x$$. Therefore, first you generate $$x$$ (as you say, by normalizing a vector with i.i.d. standard normal coordinates), then you complement $$x$$ to an orthonormal basis $$(x,e_2,\dots, e_n)$$, generate, as above, a vector $$y'$$ uniformly distributed on the unit sphere of the space spanned by $$e_2,\dots, e_n$$, and finally put $$y=\gamma x + \sqrt{1-\gamma^2} y'$$.
PS While I was typing you replaced $$\gamma$$ with $$\rho$$ and $$n$$ with $$n+1$$.