Direct synthesis of waveform integrals followed by differentiation, as
can be done for a sawtooth wave by synthesizing piecewize parabola,
seems to work well for limiting the aliasing.
Some questions. Does this work for a square wave, by differentiating
a triangle?
And, is it possible to modify the technique to implement anti--aliased
static nonlinear saturation, which for high gain factors suffers from
similar problems.
I don't see how to directly synthesize the time integral of such a
nonlinearity. For a non-linearity $\sigma(x)$, operating on a
sequence $y_n = \sigma(x_n)$, this would boil down to generating the
integral directly, or $$z_{n+1} = y_n + z_n = \sigma(x_n) +
\sigma(x_{n-1}) + \ldots$$
which doesn't seem to make much sense.