Probably the easiest way to solve the aliasing problem is to pre-calculate a series of square waves, each with successively fewer harmonics. You can then just select the right square wave for any given note and play it back via linear or polynomial interpolation. This is what I do in my iPad synth app.

The main downside of this approach is that it takes more memory, but it's very easy on the CPU.