The most computationally efficient approach is to set up the data you want in the frequency domain by putting in the tones as "spikes" in the appropriate frequency bins, and then inverse FFT'ing it to get the time domain tones. This approach will likely be a couple of orders of magnitude or so faster than your approach.
There are some downsides to doing it this way. For tones whose frequencies are integer multiple of the sample frequency divided by the number of samples it is easy- just put a spike in at the appropriate frequency bin. For all other frequencies, though, you need to put some energy into the neighboring bins as well. I would have to think about how exactly to figure out how to do this. I'm sure it can be done, but it's not a no-brainer.
Also, if you need to just produce a block of samples with the tones this approach is ideal. If you need to create a continuous stream, though, care would have to be taken to get the phases right from block to block. Again, not a no-brainer.
All that being said, I am pretty sure that this is the most computationally efficient approach.