Why is histogram equalization called as "equalization"?
To answer your question histogram equalization is called like this because its function is to produce an equalized histogram (that is an uniform probability density function).

There are different algorithms that may approach this function, and obviously there is a problem in the example that is shown:


In fact, the algorithm used there will always have trouble producing a flat histogram when mapping the function to a discrete set of 256 (except when the input histogram is already flat!)
As a consequence, the cumulative histogram shows discrete jumps which produce this non-regular output histogram.


Why is this simple algorithm still used? Because, 


the output image will often good enough for the initial purpose,
that is to maximize the range of intensities used in the image - though the histogram may look ugly...
it is simple! 


There are simple alternatives. The simplest is to use a linear interpolation of the cumulative function. 


Theory: In general terms, the method may be easily derived if we know the probability distribution function $dP_{i}$ of variable $a_i$ representing the luminance image at pixel $i$. Indeed,  by choosing the non-linearity as the cumulative distribution function transforming any observed variable $\bar{a}_i$ into:
$$
z_i(\bar{a}_i)= P_{i}(a_i \leq \bar{a}_i) = \int_{-\infty}^{\bar{a}_i} dP_{i}(a_i)%
$$
This is equivalent to the change of variables which transforms the luminance image to a variable with uniform probability distribution function in $[0, 1]^M$.
Evaluation: In physiology (retina for instance), this distribution is assumed to be sampled for each "pixel" (photoreceptor in mammals) over time (for instance, something like "the probability distribution over last minute"). In image processing, when handling one single image, you can assume stationarity over space, that is compute the histogram over all pixel values. 
Algorithm: now that we have all in place, some example algorithm in python:


Sample code:

import numpy as np
from scipy.misc import lena
I = 1. * lena()
I = (I - I.min())/(I.max()-I.min())
import matplotlib.pyplot as plt
# the cumulative histogram
n_quant = 40
(hist, bins) = np.histogram(I.ravel(), bins=np.linspace(0., 1., n_quant+1), normed=True)
P = np.cumsum(hist/np.sum(hist))
# equalization
I_z = np.interp(I, np.linspace(0., 1., n_quant), P)
(hist_z, bins) = np.histogram(I_z.ravel(), bins=np.linspace(0., 1., n_quant+1), normed=True)
plt.bar(np.linspace(0., 1., n_quant), hist, color='r', alpha=.3, width=1./n_quant)
plt.bar(np.linspace(0., 1., n_quant), hist_z, color='b', alpha=.3, width=1./n_quant)
plt.title('Histogram of Lena')
plt.xlabel('luminance')
plt.ylabel('probability')
plt.xlim((0,1.))
plt.savefig('equalization.png')


The output (red is raw histogram, blue is equalized). As you can see, there are variations, but it looks better :-). Notably, you have less "jumps" around most popular values. The parameter n_quant is important to tune smoothness versus precision.