How to implement a digital oscillator?
You're right that the strictly recursive approach is vulnerable to the accumulation of error as the number of iterations increases. One more robust way this is typically done is to use a numerically-controlled oscillator (NCO). Basically, you have an accumulator that keeps track of the instantaneous phase of the oscillator, updated as follows:

$$
\delta = \frac{2\pi f}{f_s}
$$

$$
\phi[n] = (\phi[n-1] + \delta) \mod 2\pi
$$

At each time instant, then, you're left with converting the accumulated phase in the NCO to the desired sinusoidal outputs. How you do this depends on your requirements for computational complexity, accuracy, etc. One obvious way is to just calculate the outputs as

$$
x_c[n] = \cos(\phi[n])
$$

$$
x_s[n] = \sin(\phi[n])
$$

using whatever implementation of sine/cosine you have available. In high-throughput and/or embedded systems, the mapping from phase to sine/cosine values is often done via a lookup table. The size of the lookup table (i.e. the amount of quantization you do on the phase argument to sine and cosine) can be used as a tradeoff between memory consumption and approximation error. The nice thing is that the amount of computations required is typically independent of the table size. In addition, you can limit your LUT size if needed by taking advantage of the symmetry inherent in the cosine and sine functions; you only really need to store one-fourth of a period of the sampled sinusoid.

If you need higher accuracy than a reasonably-sized LUT can give you, then you can always look at interpolation between table samples (linear or cubic interpolation, for example).

Another benefit of this approach is that it is trivial to incorporate frequency or phase modulation with this structure. The frequency of the output can be modulated by varying $\delta$ accordingly, and phase modulation can be implemented by simply adding to $\phi[n]$ directly.