In the audio domain, waveshaping is simply applying a memoryless nonlinear function to an input signal.
$$ y(t) = g\big( x(t) \big) $$
The waveshaping function, $g(x)$, is most often a continuous function that goes through the origin: $g(0)=0$.
Sometimes $g(x)$ is an odd-symmetry function: $g(-x)=-g(x)$, but it doesn't have to be. Sometimes you want 2nd harmonic distortion and then $g(x)$ does not have odd symmetry.
I am more interested in the case where the waveshaping function, $g(x)$ is a polynomial of finite order $N$:
$$ g(x) = \sum\limits_{n=0}^{N} a_n\,x^n $$
Note that if $g(x)$ has even symmetry ($g(-x)=g(x)$), then all odd terms are zero ($a_n=0$ for $n$ odd). If $g(x)$ has odd symmetry ($g(-x)=-g(x)$), then all even terms are zero ($a_n=0$ for $n$ even). I guess, if you want $g(0)=0$, then $a_0=0$.
Now, if the input is sinusoidal:
$$ x(t) = A \cos(2 \pi f_0 t) $$
then the output is periodic, having harmonic frequency components, with an upper limit to the frequency content. The highest frequency coming out will be $N f_0$. So waveshaping with polynomial mapping functions guarantees a limit to the highest frequencies generated and that can be useful in deciding an oversampling ratio necessary to guard against aliasing. As I alluded to in this answer, if the order of the polynomial is $N$, to insure against aliases folding back into your original passband, one must oversample with an upsampling ratio of at least $\frac{N+1}{2}$. Oversampling by 4x sufficies for a 7th-order polynomial.
Also note that an even-symmetry polynomial will generate only even harmonics and an odd-order polynomial will generate only odd harmoncs.
$$\begin{align}
y(t) &= g\big( x(t) \big) \\
\\
&= \sum\limits_{n=0}^{N} a_n\,\big(x(t)\big)^n \\
\\
&= \sum\limits_{n=0}^{N} a_n\,\big(A \cos(2 \pi f_0 t)\big)^n \\
\\
&= \sum\limits_{k=0}^{N} b_k \cos(2 \pi k f_0 t) \\
\end{align}$$
where the $b_k$ coefficients depend on the $a_n$ coefficients and the input amplitude of $A$. This relationship can be worked out using the Euler identity and the binomial theorem:
$$\begin{align}
y(t) &= g\big( x(t) \big) \\
\\
&= \sum\limits_{n=0}^{N} a_n\,\big(A \cos(2 \pi f_0 t)\big)^n \\
\\
&= \sum\limits_{n=0}^{N} a_n\,\left(\tfrac{A}{2}(e^{j 2 \pi k f_0 t}+e^{-j 2 \pi k f_0 t})\right)^n \\
\\
&= \sum\limits_{n=0}^{N} a_n\,\left(\tfrac{A}{2}\right)^n\left(e^{j 2 \pi f_0 t}+e^{-j 2 \pi f_0 t}\right)^n \\
\\
&= \sum\limits_{n=0}^{N} a_n\,\left(\tfrac{A}{2}\right)^n \sum\limits_{m=0}^{n} \frac{n!}{m!(n-m)!} (e^{j 2 \pi f_0 t})^m (e^{-j 2 \pi f_0 t})^{n-m} \\
\\
&= \sum\limits_{n=0}^{N} a_n\,\left(\tfrac{A}{2}\right)^n \sum\limits_{m=0}^{n} \frac{n!}{m!(n-m)!} e^{j 2 \pi m f_0 t} e^{-j 2 \pi (n-m) f_0 t} \\
\\
&= \sum\limits_{n=0}^{N} a_n\,\left(\tfrac{A}{2}\right)^n \tfrac12 \left(\sum\limits_{m=0}^{n} \frac{n!}{m!(n-m)!} e^{j 2 \pi m f_0 t} e^{-j 2 \pi (n-m) f_0 t} \\ + \sum\limits_{m=0}^{n} \frac{n!}{m!(n-m)!} e^{j 2 \pi m f_0 t} e^{-j 2 \pi (n-m) f_0 t} \right) \\
\\
&= \sum\limits_{n=0}^{N} a_n\,\left(\tfrac{A}{2}\right)^n \tfrac12 \left(\sum\limits_{m=0}^{n} \frac{n!}{m!(n-m)!} e^{j 2 \pi m f_0 t} e^{-j 2 \pi (n-m) f_0 t} \\ + \sum\limits_{m=0}^{n} \frac{n!}{m!(n-m)!} e^{j 2 \pi (n-m) f_0 t} e^{-j 2 \pi m f_0 t} \right) \\
\\
&= \sum\limits_{n=0}^{N} a_n\,\left(\tfrac{A}{2}\right)^n \tfrac12 \left(\sum\limits_{m=0}^{n} \frac{n!}{m!(n-m)!} e^{j 2 \pi (2m-n) f_0 t} + \sum\limits_{m=0}^{n} \frac{n!}{m!(n-m)!} e^{j 2 \pi (n-2m) f_0 t} \right) \\
\\
&= \sum\limits_{n=0}^{N} a_n\,\left(\tfrac{A}{2}\right)^n \tfrac12 \sum\limits_{m=0}^{n} \frac{n!}{m!(n-m)!} \left( e^{j 2 \pi (n-2m) f_0 t} + e^{-j 2 \pi (n-2m) f_0 t} \right) \\
\\
&= \sum\limits_{n=0}^{N} \sum\limits_{m=0}^{n} a_n\,\left(\tfrac{A}{2}\right)^n \frac{n!}{m!(n-m)!} \tfrac12 \left( e^{j 2 \pi (n-2m) f_0 t} + e^{-j 2 \pi (n-2m) f_0 t} \right) \\
\\
&= \sum\limits_{n=0}^{N} \sum\limits_{m=0}^{n} a_n\,\left(\tfrac{A}{2}\right)^n \frac{n!}{m!(n-m)!} \cos \big(2 \pi (n-2m) f_0 t \big) \\
\\
...
\\
&= \sum\limits_{k=0}^{N} b_k \cos(2 \pi k f_0 t) \\
\end{align}$$
Okay, to get an expression for $b_k$, the only way I can figger it out is for three cases. We let $k=n-2m$ and observe that since $2m$ is always even, only even $n$ terms contribute to $b_k$ when $k$ is even. Likewise, only odd $n$ terms contribute to $b_k$ when $k$ is odd.
The expression $\lfloor N/2 \rfloor$ is the floor()
function applied to $N/2$. just round $N/2$ down to the nearest integer.
Case 1, $k=0$ (the DC term)
Only the even terms of $n$ will contribute to $b_0$. And the only term of the inside summation that contributes is that when $m=\tfrac{n}{2}$. So letting $n=2i$:
$$ b_0 = \sum\limits_{i=0}^{\lfloor N/2 \rfloor} a_{2i}\,\left(\tfrac{A}{2}\right)^{2i} \frac{(2i)!}{\big((i)!\big)^2} $$
In the following two cases, keep in mind that the $\cos()$ function is even symmetry and these two terms are equal:
$$ \cos(2 \pi (-k) f_0 t) = \cos(2 \pi k f_0 t) $$
Case 2, $k$ even, $k>0$ (even harmonics)
Only the even terms of $n$ will contribute to $b_k$. The only two term of the inside summation that contribute are those where $m=\tfrac{n-k}{2}$ and $m=\tfrac{n+k}{2}$ . So letting $n=2i$:
$$ b_k = \sum\limits_{i=0}^{\lfloor N/2 \rfloor} a_{2i}\,\left(\tfrac{A}{2}\right)^{2i} \frac{2(2i)!}{(i-\tfrac{k}{2})!(i+\tfrac{k}{2})! } \qquad k>0 \text{ even} $$
Terms with $i<\tfrac{k}{2}$ will result in $\infty$ in the denominator because factorials of negative integers are infinite. That effectively modifies the bottom limit:
$$ b_k = \sum\limits_{i=k/2}^{\lfloor N/2 \rfloor} a_{2i}\,\left(\tfrac{A}{2}\right)^{2i} \frac{2(2i)!}{(i-\tfrac{k}{2})!(i+\tfrac{k}{2})! } \qquad k>0 \text{ even} $$
Case 3, $k$ odd, $k>0$ (odd harmonics)
Only the odd terms of $n$ will contribute to $b_k$. The only two term of the inside summation that contribute are those where $m=\tfrac{n-k}{2}$ and $m=\tfrac{n+k}{2}$. So letting $n=2i+1$:
$$ b_k = \sum\limits_{i=0}^{\lfloor (N-1)/2 \rfloor} a_{2i+1}\,\left(\tfrac{A}{2}\right)^{2i+1} \frac{2(2i+1)!}{(i-\tfrac{k-1}{2})!(i+\tfrac{k+1}{2})!} \qquad k>0 \text{ odd} $$
Terms with $i<\tfrac{k-1}{2}$ will result in $\infty$ in the denominator because factorials of negative integers are infinite. That effectively modifies the bottom limit:
$$ b_k = \sum\limits_{i=(k-1)/2}^{\lfloor (N-1)/2 \rfloor} a_{2i+1}\,\left(\tfrac{A}{2}\right)^{2i+1} \frac{2(2i+1)!}{(i-\tfrac{k-1}{2})!(i+\tfrac{k+1}{2})!} \qquad k>0 \text{ odd} $$
.
(i wouldn't mind if someone else checked my math. these results tell us something about how a polynomial waveshaping function will sound, at least what the amplitude of harmonics will be.)