Answer:

Guitar distortion basically stems from overdriven amplifiers and nonlinearities present in the components used (maybe tubes or the transformer or diodes). The signal is clipped and all sorts of stuff not present in the original signal appears. Though not an accurate representation of the original signal, this new distorted signal tends to rock way more.

So how does one make this distortion digitally?

This is normally accomplished through waveshaping. We take the original signal between -1 and 1 and map it, using a mathematical function, to a new signal between -1 and 1. This is what is called a memoryless nonlinearity. Memoryless because the next sample out has nothing to do with any of the other previous samples, and nonlinear because adding two signals together and obtaining the result is not the same as obtaining the results separately for each signal and then adding them together. So the higher your amplitude going into the waveshaper the more clipping and the more distortion.

There is hard clipping and soft clipping. Not easy to explain without pictures but hard clipping means you chop off everything above 1 and below -1, a soft clipper has a smoother transition between the probably linear middle region and the clipped parts.

for example:

a 3rd order polynomial
f(x) =
-2/3, x < -1
x - x3/3, -1 <= x <= 1
2/3, x > 1

or maybe arctan, or tanh (but polynomials are faster to calculate).

These functions are odd i.e. they are symmetrical, f(-x) = -f(x), so they only generate odd harmonics (take to this to the limit and you get a square wave which is a fundamental and all the odd harmonics). Breaking this symmetry in any way will add some even harmonics, tubes will create some even harmonics and as will things like pickups that aren't placed symmetrically like in the fender rhodes.

Easy ways to make your function asymmetrical are to add a DC offset into the signal before the nonlinearity so e.g. the top is clipped more than the bottom, or you could add an even function like a square to your function.This will change the duty cycle of your waveform e.g. the first peak might take up more of the period than second peak and means the higher your input signal the more unequal the duty cycle is, i.e. duty cycle modulation as a function of amplitude.

Now you probably want to filter the signal before and after the nonlinearity because you want specific frequencies to be distorted more than others. Distorted low frequencies can sound "flatulent" and distorted high frequencies sound really "fuzzy". Most likely what you want is a band pass or high pass before and lowpass after adjusted to taste. This should make your "axe" wail like Slash's.

Problems with this:

Intermodulation, basically beating of the sum and/or difference between all the frequencies and harmonics you now have in your signal. Mostly apparent when you play a chord. This happens in the analog domain too and your filters will mostly take care of this.

Aliasing, you may have introduced some frequencies higher than the nyquist. You it hear mostly when you slide between notes or bend a note (you hear some frequencies bend or slide in the opposite direction as the frequencies wrap around). This actually depends on the order of your polynomial. If your input signal contains no frequencies over fs/2n you can apply an n-th order polynomial without aliasing.

If your waveshaping function is a polynomial you can bandlimit you signal in a parallel filter bank with lopass filters with cutoffs at fs/2n and apply up to an n-th order polynomial to your bandlimited signals then sum the result and avoid all the aliasing! You can control the gain of each polynomial term to make any desired polynomial waveshaper. The problem is that usually you would use a polynomial section and then clip after the bounds of the section like in the 3rd order polynomial example above. This means you no longer have a nice neat polynomial and will introduce a bunch of aliasing. The solution to this is to oversample i.e. stick zeros between the samples, bandlimit to original fs/2, distort, bandlimit to original fs/2, then downsample by tossing out every other sample. This gives you some space to alias and then filters out the bad stuff. It won't completely eliminate the aliasing but will get rid of some of it. How much you oversample is up to you. People mostly seem to think 2 times oversampling is sufficient as anything left is obscured by the distortion anyway. Oversampling tends to be computationally expensive...

to sum up, you probably want:

pre-eq (bandpass) > pre-gain > nonlinear waveshaping function > post-eq (lowpass) > post-gain and you may wish to eliminate aliasing using the bandlimited polynomials or by oversampling.

This is all a kind of ad hoc distortion scheme and you basically fiddle around until you get something that suits your needs. Modelling specific equipment is much more complicated and involves system identification and frequently the nonlinearities are not memoryless.

some links:

modelling nonlinear systems non-linear systems with Volterra and Weiner theories and the Laguerre Transform:
http://www.tele.ntnu.no/akustikk/meetings/DAFx99/schattschneider.pdf

some good papers on modelling analog gear like tube amps:
http://www.simulanalog.org/

people seem to like the Ibanez Tube Screamer distortion boxes... specs:
http://www.geofex.com/Article_Folders/TStech/tsxfram.htm
and nice little bilinear transform model of the TS9:
http://www.gweep.net/~shifty/dsp/TS9SIMUL.DOC

Julius Smith's electric guitar model stuff:
http://ccrma-www.stanford.edu/~jos/waveguide/Nonlinear_Distortion.html
http://ccrma-www.stanford.edu/~jos/waveguide/Amplifier_Feedback.html