The following result is relevant here: If $f_n : \mathbb{R} \to \mathbb{R} $ is a sequence of continuously differentiable functions, uniformly convergent to $f$ and the sequence $f'_n$ converges uniformly to $g$, then $f$ is differentiable and $f'=g.$
So we can check that $F$ and $G$ are smooth by checking that we can keep on differentiating them, which we do by applying the above theorem, the Weierstrass M-test and the fact that Schwarz functions decay very quickly.
Because Fourier transform is a linear transformation in the Hilbert space of square-normalizable curves $L^2$, this is essentially an eigenvalue problem. You are looking for a function that doesn't change under Fourier transform. Beside the Gaussian function, which is very famously an eigenfunction of the Fourier transform, there is a whole class of them that also satisfy this eigenvalue condition of $Ff=cf$. They are essentially a Gaussian times a Hermite polynomial:
http://en.wikipedia.org/wiki/Hermite_functions#Hermite_functions_as_eigenfunctions_of_the_Fourier_transform
On the other hand, $F^4=I$ is universal. If you observe the Fourier transform integral, you can see that the Fourier transform and inverse Fourier transform have an opposite sign in the imaginary exponent. This is why $F$ is not exactly $F^{-1}$ but it's close. $F^2$ flips the function (reverses time) because of this mismatch in sign, but reversing time twice brings you back again, which is why $F^4=I$ works for any function, not just the eigenfunctions. This is related to $i^4=1$: $F^4=I$ tells you that the only eigenvalues of $F$ are $c\in\{\pm 1, \pm i\}$ (try $F^4 f=F^3 (cf)=c^4 f=f$ which gives you $c^4=1$). Of course, this means, that the Fourier space is highly degenerate. A linear combination of any number of eigenfunctions with the same eigenvalue is also an eigenfunction. You can use the Hermite functions (they sequentially cycle through eigenvalues $i^n$, so take every fourth) as a basis spanning each of the four subspaces: you get an infinite family of eigenfunctions for each eigenvalue.
Be warned that this depends on the normalization of the Fourier transform. You need the unitary version. If you take the regular definition (the one used in physics and engineering), you get additional factors of $\sqrt{2\pi}^n$, which the transform blows up or extinguishes the amplitude of the function, and is not cyclic anymore.
You will notice that you can split any function into 4 components with eigenvalues $\{1,i,-1,-i\}$ by doing this:
$$\frac{1}{4}(1+F+F^2+F^3)f=f_1$$
$$\frac{1}{4}(1-iF-F^2+iF^3)f=f_i$$
$$\frac{1}{4}(1-F+F^2-F^3)f=f_{-1}$$
$$\frac{1}{4}(1+iF-F^2-iF^3)f=f_{-i}$$
where
$$Ff_1=f_1\quad Ff_i=if_i\quad Ff_{-1}=-f_{-1}\quad Ff_{-i}=-if_{-i}$$
Each of these components has a very specific symmetry. In particular, the $\pm i$ are odd and $\pm 1$ are even.
You can take any arbitrary function and use it to generate a function that is preserved under Fourier transform. This I think is the most general answer I can give you.
Best Answer
Write the Fourier transform as $\mathcal{F}: \mathcal{S} \to \mathcal{S}$. Observe that $\mathcal{F}^{-1} \neq \mathcal{F}$ (assuming we choose the unitary definition), but that $$ \mathcal{F}^{-1} f(x) = \mathcal{F} f(-x) $$ Therefore we conclude that $$ \mathcal{F}\mathcal{F}\mathcal{F}\mathcal{F} = \mathrm{Id}.$$ This tells you that the eigenvalues must be fourth roots of unity.
For concreteness let us fix convention $$ \mathcal{F} f(\xi) = \int_{\mathbb{R}} e^{-2\pi ix\xi} f(x) ~\mathrm{d}x $$ which seems to be the one you are using considering that you have found $\exp(-\pi x^2)$ to be an eigenfunction.
To find the eigenfunctions, we can start with the hint. Now, observe that $$ \mathcal{F}(-2\pi ix f) = D_\xi \mathcal{F} f $$ (the rule that exchanges scalar multiplication with differentiation). Now let $\lambda$ be one of the fourth roots of unity (meaning that $\lambda \in \{1,-1,i,-1\}$), we want to solve $$ P(-2\pi i x) \exp(-\pi x^2) = \lambda P(D_x) \exp(-\pi x^2) $$ where $P(D_x)$ is the differential operator $\sum a_n D_{x}^n$ if $P(y) = \sum a_n y^n$.
We easily see that $$ D_x \exp(-\pi x^2) = -2\pi x \exp(-\pi x^2) $$ so $P(y) = y$ gives us an eigenfunction for $\lambda= i$.
Analogously you can solve for a quadratic polynomial $P(y)$ which will give you $\lambda = -1$ and a cubic polynomial $P(y)$ which will give you $\lambda = -i$.