One of the best ways to create a sine wave is to use a complex phasor with recursive updating. I.e.
z[n+1] = z[n]*W
where z[n] is the phasor, W = exp(j*w), with being w the angular frequency of the oscillator in radiants and n the sample index. Both real and imaginary part of z[n] are sine waves, they are 90 degrees out of phase. Very convenient if you need both sine and cosine. A single sample calculation only requires 4 multiples and 4 adds and is MUCH cheaper than anything containing sin() cos() or lookup tables.
The potential problem is that the amplitude may drift over time due to numerical precision issues. However there is a fairly straight forward to repair that. Let's say that z[n] = a + j*b.
We know that z[n] should have unity magnitude, i.e.
a*a + b*b = 1
So we can check every once in a while if that is still the case and correct accordingly. The exact correction would be
z'[n] = z[n]/sqrt(a*a + b*b)
That is an awkward calculation but since a*a + b*b is very close to unity you can approximate the 1/sqrt(x) terms with a taylor expansion around x = 1 and we get
1/sqrt(x) ~= (3-x)/2;
so the correction simplifies to
z'[n] = z[n]*(3-a*a-b*b)/2
Applying this simple correction every few hundred samples will keep the oscillator stable forever.
To vary the frequency continuously the multiplier W needs to be updated accordingly. Even a non-continuous change in the multiplier will maintain a continuous oscillator function. If frequency ramping is required the update can either be broken down into a few steps or you can use the same oscillator algorithm to update the multiplier itself (since it's also a unity gain complex phasor).