I'm reading mangling to mean unpredictable but interesting results, like blindly patching a modular drone. I think the most sonic possibilities will come from using FM, and you wont really need to know the special math for FM stuff since you can naively compute the spectra ahead of time (naive here means you can ignore the cool math and just make some preset fm waveforms and compute the spectra). However, if you use additive resynthesis from cheap waveforms you'll leave yourself as much dsp room as possible to do cool subtractive stuff on it, or things like wavefolding and pwm, which can still be set up in an unpredictable and creative way but without the high risk of stuff sounding bad that you get with FM. This is because for all the possible sounds from FM mangling, only some of them will produce harmonically simple enough waveforms to sound good.
I'll explain the math a bit for you. Usually we view audio signals as amplitude with respect to time. However, a combination of cool things said by old crazy nerds like euler, laplace, and fourier let us write those same functions as amplitude with respect to frequency. The coolest thing is that we have functions to convert between these 2 domains (domain here means the same thing it usually does, the x axis on your graph). spectral resynthesis is about finding easy and fast ways to convert from the frequency domain to the time domain. Instead of trying to mathematically convert from the frequency domain back to the time domain (inverse fourier transform) we synthesis (the resynthesis part) a wave that looks almost the same in the frequency domain (the spectral part). Since sine waves have only 1 frequency, they're the simplest way for humans to think about spectral resynthesis because you just have to make a sine wave for each point in your spectrum (graph where x is frequency and y is amplitude) and make the frequency of the sine wave equal the x coordinate of the point, and adjust its volume to the y coordinate of the point.
Unfortunately, sine waves are really hard for computers to make. Luckily, square waves arent. We can basically do the same thing as before, because square waves are just a bunch of odd harmonics with amplitude inversely proportional to frequency. In other words, a square wave at frequency f is sorta like a sin(f), and its a little more like a sin(f)+sin(3f)/3, and even more like a sin(f)+sin(3f)/3+sin(5f)/5, and so on. Because of this a square wave is very cheap for a computer to work with not only in the time domain, but in the frequency domain as well, since we can ignore all the harmonics higher than 20khz, so we have a pretty limited number of them to deal with and that number decreases as the frequency of the square wave increases. Note that I omitted some scaling factors for readability, for those factors just look up the square wave fourier series.
now I'll change the wording for how we think about spectral resynth. suppose our target spectrum is T(f), aka a function of frequency that we want to get, and that our resynthesized waveform has a spectrum of S(f). earlier we said we want T(f) = S(f), but if we subtract S(f) from both sides we get T(f)-S(f)=0. Even though we're in the frequency domain, spectra are still functions and pretty much act like it, so we can do a lot of the same stuff with spectra as we can regular functions.
The other thing we can do is say that if we add 2 simple waveforms, the spectrum of the sum is the same as the sum of the spectrum (i.e. the fourier transform is a linear transform).
If we think about our sine wave resynthesis again, what we're really doing is subtracting the spectra of a bunch of sine waves from our original spectrum in the hopes of getting close to 0. We can do the same thing with square waves, the only difference is that if we add a square wave at frequency f with amplitude A we need to subtract A/pi from the original spectrum at frequency f, A/3pi at frequency 3f, etc.
So to resynthesize a discrete spectrum (i.e. a spectrum with a finite number of points, where the x coordinate is the middle of a band and the y coordinate is the amplitude of said band) with square waves, all you have to do is look at your lowest band and make a square wave at the center frequency of that band, and adjust its volume to match the amplitude of that band. The only extra step is that you now have to subtract the higher partials before moving onto the next band, but the subtraction is easy to calculate and your next step is still to just add a square wave at the current lowest frequency band thats got an amplitude that isnt 0.
There are 2 things that complicate this, one of which matters and one of which doesnt, so I'll discuss the important one first.
the fourier series of a square wave isnt really what its spectrum is gonna look like, because it assumes the square wave goes from negative infinity to positive infinity on the time axis. What happens in real situations is that these peaks in the spectrum (called kronecker deltas, or places where the function jumps from 0 to something else and back to 0 instantly) get smoothed out and end up looking like normal distributions. that doesnt really matter since you can just run your square wave through the spectrum analyzer you're using ahead of time, and as long as you're comparing spectra that were produced by the same analyzer you should be ok.
the thing that doesnt matter for square waves but would matter if you wanted to try other waveforms is that when your x axis is frequency, your y axis is actually complex. That sounds scary, but think of it like this: your x axis (frequency) is now the middle of a cylinder. your amplitudes are vectors that point towards the wall of the cylinder, and the length of the vector is still the amplitude of that partial frequency, but the amount the arrow is rotated in the cylinder is the phase shift of that partial. if the arrow is pointing straight at you, it's got no phase shift. if its pointing away from you, then its got half a period of phase shift (pi radians or 180 degrees). if you phase shift a periodic signal by 360 degrees, you're basically spinning the cylinder around a full cycle, and all the arrows end up pointing in the same direction as they did before. The only thing this affects in our earlier model is how spectra are added and subtracted. say you have an arrow in a spectrum at frequency 220hz with amplitude 1, and its pointed towards you, so its not phase shifted. now say you add another spectrum with an identical arrow, but with 180 degree phase shift so it points the other way. if you add the spectra together (aka add the signals together) then the result is 0 at 220hz because the vectors cancel out, but since our earlier model didnt include phase shift it just reads 1+1 and spits out 2. luckily, this doesnt matter for square waves because all the partials of a square wave have no phase shift. there are other waveforms (in fact infinitely many) that also have this property, and for each one of those there are infinitely many other waveforms that don't have this property.
So we luckily dont have to worry about phase cancellation with perfectly tuned square waves, but the tradeoff is that if one of the square waves gets detuned we'll get an interference pattern. The solution to this is to sync all of the square waves to the lowest one, but that only works if all the square waves we use have a frequency thats an integer multiple of the frequency of the lowest square wave.
If you want a function written for you let me know, but understanding all this should give you enough to write the function and is honestly necessary to do anything interesting with it