Best is to perform all the calculations in your main loop, to keep your ISR shorter - unless that ISR really have to react super-quickly.

All frequencies are quotients of the frequency of the base note.

So you can have 6 high-precision phase variables.
Round these 6 phase variables to closest sample. Pick up sample and scale with amplitude and add to composite value.
Then output the composite value.

Now step each of these phase variables with the time represented until your next sample output. Take care of overflow at 360 degrees.

If you have 40kHz sampling rate, each sample represents 25 us.
32.7 Hz has a period of 1/32.7 [s] = 30.581 ms.
So at 40kHz sampling rate, each sample would step the phase with 25/30581 *360 = 0.2943 degrees.

49.0 Hz has a period of 1/49.0 [s] = 20.408 ms.
So at 40kHz sampling rate, each sample would step the phase with 25/20408 * 360 = 0.4410 degrees.

...

So every iteration, you take these degree values and add to the individual phase values. And if >= 360 degrees, subtract 360.

You can precompute the above values for all keys in an octave. For higher octaves, you just need to scale the delta-phase to add for the individual sines.

40kHz sampling rate could represent a very significant processing load on your processor - the actual sampling rate needed will very much depend on what frequencies you need to output, i.e. how many octaves your devise needs to support.

Next thing - are you only going to handle one key at a time, or how many concurrent keys must you support? The number of frequencies to compute quickly grows. If you have code space, it's possible to speed up the evaluations by wasting much more space on tables.

You could have a sine table precomputed for 32.70 Hz, so every sample steps one step forward for the C. Two steps forward for next C (but this index moves twice as fast so need to restart when overflowing). Four steps forward for next C (corresponding handling of overflow).

Corresponding for the other frequencies.

The size of the tables will then be greatly affected by the sampling rate, since that controls the number of samples for a full period of each frequency. And the sampling rate must obviously be large enough that the rounding error for a full period doesn't hurt too much. But the advantage is that you would step 2^n ticks forward in each individual table for each sample output, where n depends on which octave.

Quite frankly, this is a task that a cheap ARM7 or Cortex M3 would manage better, because of better availability of RAM + ROM, many 32-bit registers for either 32-bit additions or indexing. An ARM chip would require many less instructions/second to perform the task, allowing many more concurrent frequencies to be processed - or allowing them to be generated at a higher resolution/fidelity.