It is possible to find functions $\phi, \psi \in \mathbb{R}[sin(x), sin(y), cos(x), cos(y)]$, so that $S^2 = \phi( [0,1]^2)$ and $\psi( [0,1]^2)$ is a torus. Is it possible find, for any genus g, functions $f_g \in \mathbb{R}[sin(x), sin(y), cos(x), cos(y)]$ so that $f_g([0,1])^2$ is $S_g$, the surface of genus $g$?
Since the torus and the sphere are surfaces of revolution, and the other genus g surfaces are not (is this true?), I expect that perhaps this cannot be done. I do not know how to make this intuition precise.
If there is no trigonometric polynomial parametrization as above, then what is the correct generalization of these parametrizations?
So the most modest form of my question is: How can I produce parametrized genus g surfaces in $R^3$, for arbitrary $g$?
All surfaces below are assumed to be compact, connected, and orientable.
Any mapping $f$ whose component functions are polynomials in the circular functions $\cos x$, $\sin x$, $\cos y$, $\sin y$ may be viewed as a mapping defined on the torus $T = \mathbf{R}^{2}/(2\pi\mathbf{Z})^{2}$. Let $S = f(T)$ denote the image.
If $S$ is a closed surface and $f:T \to S$ is a branched covering map, the genus of $S$ is either $0$ or $1$ by the Riemann-Hurwitz formula.
Contrapositively, if $S$ is a closed surface of genus $g > 1$, there is no trigonometric parametrization of $S$ that is regular except at finitely many branch points.
To come at the question topologically, if $S$ has genus $g > 1$, then $f_{*}$ is not surjective on $1$-dimensional homology (since the $1$-dimensional homology of an orientable surface of genus $g$ is $2g$-dimensional). Intuitively, (most) representatives of $1$-dimensional homology in $S$ are obtained from intervals in $T$ by identifying endpoints; $f$ does something like squash the torus flat, then wrap the resulting region around $S$ like the standard edge identification picture of a $4g$-gon.
If $S$ has genus $g > 1$, the universal cover may be regarded as the unit disk equipped with the hyperbolic metric of constant curvature $-1$, with deck transformations forming a group of hyperbolic isometries. The "correct" generalization of the circular functions is presumably a suitable collection of real-valued functions obtained from modular forms. (I've never seen this done, but it doesn't seem computationally/numerically infeasible.)