If you need something that stays linear over any distance (unlike distance^2) and yet appears vaguely circular (unlike the squarish Chebyshev and diamond-like Manhattan distances), you can average the latter two techniques to get an octagonally-shaped distance approximation:
dx = abs(x1 - x0)
dy = abs(y1 - y0)
dist = 0.5 * (dx + dy + max(dx, dy))
Here is a visualization (contour plot) of the function, thanks to Wolfram Alpha:
And here is a plot of its error function when compared to the euclidean distance (radians, first quadrant only):
As you can see, the error ranges from 0% on the axes to approximately +12% in the lobes. By modifying the coefficients a bit we can get it down to +/-4%:
dist = 0.4 * (dx + dy) + 0.56 * max(dx, dy)
Update
Using the above coefficients, the maximum error will be within +/-4%, but the average error will still be +1.3%. Optimized for zero average error, you can use:
dist = 0.394 * (dx + dy) + 0.554 * max(dx, dy)
which gives errors between -5% and +3% and an average error of +0.043%
While searching the web for a name for this algorithm, I found this similar octagonal approximation:
dist = 1007/1024 * max(dx, dy) + 441/1024 * min(dx, dy)
Note that this is essentially equivalent (though the exponents are different - these ones give a -1.5% to 7.5% error, but it can be massaged to +/-4%) because max(dx, dy) + min(dx, dy) == dx + dy. Using this form, the min and max calls can be factored out in favor of:
if (dy > dx)
swap(dx, dy)
dist = 1007/1024 * dx + 441/1024 * dy
Is this going to be faster than my version? Who knows... depends on the compiler and how it optimizes each for the target platform. My guess is it'd be pretty difficult to see any difference.