Why does bosonic string theory require 26 spacetime dimensions? I do not think it is possible really believe or experimentally check (now), but all modern physical doctrines suggest that out world is NOT 4-dimensional, but higher. 
The least sophisticated candidate - bosonic string theory says that out world is 26 dimensional (it is not realistic due to presence of tachion, and so there are super strings with 10 dimensions, M-theory with 11, F-theory with 12).
Let us do not care about physical realities and ask: what mathematics stands behind the fact that 26 is the only dimension where bosonic string theory can live ? 
Definitely there is some mathematics  e.g. 26 in that MO question is surely related.
Let me recall the bosonic string theory background. 
Our real world is some Riemannian manifold M which is called TS (target space).
We consider the space of all maps from the circle to M, actually we need to consider how
the circle is moving inside M, so we get maps from $S^1\times [0~ T]$  to M ( here $S^1\times [0~ T]$ is called WS - world sheet); 
we identify the maps which differs by parametrization (that is how Virasoro comes into game
and hence relation with Leonid's question). 
That was pretty mathematical, but now ill-defined physics begin - we need integrate
over this infinite-dimensional space of maps/parametrizations with measure
corresponding to exp( i/h volume_{2d}(image(WS))).
This measure is known NOT to exist mathematically, but somehow this does not stop
physists they do what they call regularization or renormalization or something like
that and 26 appears...  
 A: The value of 26 ultimately comes from the need to rid the theory of negative-norm states, as previously noted. This involves regularizing the sum $ \sum_{n=0}^\infty n, $ which of course diverges. One obtains a finite part equaling $-\frac{1}{12},$ which leads to the 24, and from there to 26 for the number of dimensions (25 space, 1 time). (Analytic continuation of the zeta function gives $\zeta(-1) = -\frac{1}{12} $).
Rather than me writing the details in here, I suggest this nice introductory reference, where the number of dimensions required for consistency of the bosonic string is derived in Ch. 2:
http://www.damtp.cam.ac.uk/user/tong/string/string.pdf
A: Just to correct a small misconception:  Physicists aren't claiming that an integral over the space of smooth maps from $WS$ to $X$ exists.  The path integral measure for the bosonic string is defined on something more like the "space of distributions on $WS$ valued in $X$".  It's not a problem that no integral exists on the space of smooth maps, because this space shows up only as a convenient shorthand for discussing the renormalization procedure which is used to define the path integral measure.  
The quantum mechanics of the simple harmonic oscillator is subject to a similar abuse of language.  You discuss the theory in terms of functions on the timeline $[0,t]$, but if you're careful, the path integral measure (aka, the Ornstein-Uhlenbeck measure) is actually defined on the space of distributions on $[0,t]$.  (What makes life easy in this case is that the Wiener measure is supported on distributions which are almost everywhere continuous functions.)
The situation is somewhat more complicated for the 2d nonlinear sigma model, because there isn't really anything you'd want to call the distributions valued in $X$.  Instead you try to define the measure as a linear functional on observables which are well approximated by functions of the form $\phi \mapsto ev_{\sigma} \phi^*f$, where $f$ is a function on $X$ and $ev_\sigma$ evaluation at a point $\sigma \in WS$.   The correlation function of observables $\hat{\mathcal{O}}_1$, $\hat{\mathcal{O}}_2$ should be approximated by integrals of the form
$\int_{Map(L,X)} \mathcal{O}_1(\phi)  \mathcal{O}_1(\phi) e^{i S_L(\phi)} d\phi$
for some finite set of points $L \subset WS$, and some approximation $S_L$ of the classical action $S$ defined using only finite differences among the points in $L$.  When you refine the set of points $L$ to fill in $WS$, you get an expectational value functional on the set of observables.  This expectation value functional should have the same properties (like OPE) that you see in QFTs where the classical fields take values in linear spaces.
A: I'm not an expert on this so bare with me, but I don't think you must require $\dim(M) = 26$, you must only require that the worldsheet is conformally invariant - i.e., the Weyl anomaly vanishes. You can do this by adding 26 bosons (which represent the coordinates of $M$) - which is called critical string theory - or you can turn on the dilaton expectation value - which is then called non-critical string theory. There's a lot of interesting research involving these non-critical string theories, for e.g. check out $c=1$ matrix models and type 0 string theories.
A: $$$$
Note, that here, the $\hat L_n$ are operators on the state given by the sums of the dots of the mode operators, i.e. $\hat L_0=\sum_{k=-\infty}^\infty\hat\alpha_{-n}\cdot\hat\alpha_n$. 
Also note that The Virasoro Algebra is the central extension of the Witt/Conformal Algebra so that explains why we have a $D$, it is equivalent to the central charge.     
I'll expand on Chris Gerig's answer.   
Not only do we need $D=26$, we also need  the normal ordering constant $a=1$.  The normal ordering constant is the eigenvalue of $\hat L_0$ with the eigenvector the state.     
We want to promote  the time-like states to spurious, zero-norm states, right?   So, we impose the (level 1) spurious state conditions on the state as ffollows ($|\chi\rangle$ are the basis vectors to build the spurious state $\Phi\rangle$     on.)     
$$  \begin{gathered}
  0 = {{\hat L}_1}\left| \Phi  \right\rangle     \\
  {\text{ }} = {{\hat L}_1}{{\hat L}_{ - 1}}\left| {{\chi _1}} \right\rangle     \\
  {\text{ }} = \left[ {{{\hat L}_{ - 1}},{{\hat L}_1}} \right]\left| {{\chi _1}} \right\rangle  + {{\hat L}_{ - 1}}{{\hat L}_1}\left| {{\chi _1}} \right\rangle    \\
  {\text{ }} = \left[ {{{\hat L}_{ - 1}},{{\hat L}_1}} \right]\left| {{\chi _1}} \right\rangle    \\
  {\text{ }} = 2{{\hat L}_0}\left| {{\chi _1}} \right\rangle     \\
  {\text{ }} = 2{c_0}\left( {a - 1} \right)\left| {{\chi _1}} \right\rangle    \\ 
\end{gathered} $$
That means that $a=1$.     
Now, for a level 2 spurious state, 
$$\begin{gathered}
  \left[ {{{\hat L}_1},{{\hat L}_{ - 2}} + k{{\hat L}_{ - 1}}{{\hat L}_{ - 1}}} \right]\left| \psi  \right\rangle  = \left( {3{{\hat L}_{ - 1}} + 2k{{\hat L}_0}{{\hat L}_{ - 1}} + 2k{{\hat L}_{ - 1}}{{\hat L}_0}} \right)\left| \psi  \right\rangle {\text{   }}           \\
  {\text{                         }} = \left( {3 - 2k} \right){{\hat L}_{ - 1}} + 4k{{\hat L}_0}{{\hat L}_{ - 1}}{\text{  }}\left( {3 - 2k} \right){{\hat L}_{ - 1}} + 4k{{\hat L}_0}{{\hat L}_{ - 1}}{\text{ }}     \\
  0 = {{\hat L}_1}\left| \psi  \right\rangle  = {{\hat L}_1}\left( {{{\hat L}_{ - 2}} + k{{\hat L}_{ - 1}}{{\hat L}_{ - 1}}} \right)\left| {{\chi _1}} \right\rangle  = \left( {\left( {3 - 2k} \right){{\hat L}_{ - 1}} + 4k{{\hat L}_0}{{\hat L}_{ - 1}}} \right)\left| {{\chi _1}} \right\rangle           \\
  {\text{                                     }} = \left( {\left( {3 - 2k} \right){{\hat L}_{ - 1}} + 4k{{\hat L}_{ - 1}}\left( {{{\hat L}_0} + 1} \right)} \right)\left| {{\chi _1}} \right\rangle    \\
  {\text{                                     }} = \left( {3 - 2k} \right){{\hat L}_{ - 1}}\left| {{\chi _1}} \right\rangle     \\
  2k = 3    \\
  k = \frac{3}{2}    \\ 
\end{gathered} $$           
Since this level 2 spurious state can be written as:  
$$  {\left| \Phi  \right\rangle  = {{\hat L}_{ - 2}}\left| {{\chi _1}} \right\rangle  + k{{\hat L}_{ - 1}}{{\hat L}_{ - 1}}\left| {{\chi _2}} \right\rangle }$$      ##   
So, then,  
$$  \begin{gathered}
  {{\hat L}_2}\left| \Phi  \right\rangle  = 0  \\
  {{\hat L}_2}\left( {{{\hat L}_{ - 2}} + \frac{3}{2}{{\hat L}_{ - 1}}{{\hat L}_{ - 1}}} \right)\left| {{\chi _2}} \right\rangle  = 0                               \\
  \left[ {{{\hat L}_2},{{\hat L}_{ - 2}} + \frac{3}{2}{{\hat L}_{ - 1}}{{\hat L}_{ - 1}}} \right]\left| {{\chi _2}} \right\rangle  + \left( {{{\hat L}_{ - 2}} + \frac{3}{2}{{\hat L}_{ - 1}}{{\hat L}_{ - 1}}} \right){{\hat L}_2}\left| {{\chi _2}} \right\rangle  = 0                     \\
  \left[ {{{\hat L}_2},{{\hat L}_{ - 2}} + \frac{3}{2}{{\hat L}_{ - 1}}{{\hat L}_{ - 1}}} \right]\left| {{\chi _2}} \right\rangle  = 0     \\
  \left( {13{{\hat L}_0} + 9{{\hat L}_{ - 1}}{{\hat L}_{ + 1}} + \frac{D}{2}} \right)\left| {{\chi _2}} \right\rangle  = 0                            \\
  \frac{D}{2} = 13    \\
\text{Since $L_0|\chi_2\rangle = -|\chi_2\rangle$ and $L_{+1}|\chi_2\rangle=0$, we have  }
  D = 26     \\ 
\end{gathered} $$       \   
And then, finally,
Q.E.D.       
So, this was done essentially to remove the imaginary norm ghost states and using  the Canonical / Gupta - Bleuer formalism.   
It's also possible to use , say, e.g. Light Cone  Gauge (LCG) quantisation.   However, in other quantisation methods, the conformal anomaly is manifest in other forms. E.g., in LCG quantisationn, it is manifest as a failure of lorentz symmetry. A good overview of this method can be found in Kaku Strings, Conformal fields, and M-theory (it's the only part of the book that I liked, actually. The rest of the book is too rigorous, without much physical intuition.).                 
A: In addition to Chris Gerig's operator-language approach, let me also show how this magical number appears in the path integral approach.
Let $\Sigma$ be a compact surface (worldsheet) and $M$ a Riemannian manifold (spacetime). The string partition function looks like
$$Z_{string}=\int_{g\in Met(\Sigma)}dg\int_{\sigma\in Map(\Sigma,M)}d\sigma\exp(iS(g,\sigma)).$$
Here $Met(\Sigma)$ is the space of Riemannian metrics on $\Sigma$ and $S(g,\sigma)$ is the standard $\sigma$-model action $S(g,\sigma)=\int_{\Sigma} dvol_\Sigma \langle d\sigma,d\sigma\rangle$. In particular, $S$ is quadratic in $\sigma$, so the second integral $Z_{matter}$ does not pose any difficulty and one can write it in terms of the determinant of the Laplace operator on $\Sigma$. Note that the determinant of the Laplace operator is a section of the determinant line bundle $L_{det}\rightarrow Met(\Sigma)$. The measure $dg$ is a 'section' of the bundle of top forms $L_g\rightarrow Met(\Sigma)$. Both line bundles carry natural connections.
However, the space $Met(\Sigma)$ is enormous: for example, it has a free action by the group of rescalings $Weyl(\Sigma)$ ($g\mapsto \phi g$ for $\phi\in Weyl(\Sigma)$ a positive function). It also carries an action of the diffeomorphism group. The quotient $\mathcal{M}$ of $Met(\Sigma)$ by the action of both groups is finite-dimensional, it is the moduli space of conformal (or complex) structures, so you would like to rewrite $Z_{string}$ as an integral over $\mathcal{M}$.
Everything in sight is diffeomorphism-invariant, so the only question is how does the integrand change under $Weyl(\Sigma)$. To descend the integral from $Met(\Sigma)$ to $Met(\Sigma)/Weyl(\Sigma)$ you need to trivialize the bundle $L_{det}\otimes L_g$ along the orbits of $Weyl(\Sigma)$. This is where the critical dimension comes in: the curvature of the natural connection on $L_{det}\otimes L_g$ (local anomaly) vanishes precisely when $d=26$. After that one also needs to check that the connection is actually flat along the orbits, so that you can indeed trivialize it.
Two references for this approach are D'Hoker's lectures on string theory in "Quantum Fields and Strings" and Freed's "Determinants, Torsion, and Strings".
A: Here is a purely mathematical statement where 26 appears. Introduce for every integer $d \geq 4$, $k \geq 0$ , the degree $k$ polynomial $P_k^{(d)}(x)$ coefficients of the expansion
$\frac{1}{(1-2xy+y^2)^{(d-3)/2}}=\sum_{k \geq 0} P_k^{(d)}(x) y^k$ 
The polynomial $P_k^{(d)}$ is usually called the $k$-th Gegenbauer polynomial in spacetime dimension $d$. For instance, $P_0^{(d)}=1$, $P_1^{(d)}=(d-3)x$, $P_2^{(d)}=\frac{d-3}{2}((d-1)x^2-1)$.
For every integer $N \geq 0$, consider the degree $N$ polynomial
$Q_N(x)=\prod_{j=0}^N(x+\frac{2j-N}{N+4})$
and expand it in terms of the Gegenbauer polynomials in spacetime dimension $d$:
$Q_N(x)=\sum_{k=0}^N a_k^{(d)} P_k^{(d)}(x)$
The non-trivial statement is: the rational numbers $a_k^{(d)}$ are non-negative for every $k \geq 0$ if and only if $d \leq 26$.
For example, the first non-trivial case is $N=1$:
$Q_1(x)=(x-\frac{1}{5})(x+\frac{1}{5})=x^2-\frac{1}{25}=\frac{2}{(d-3)(d-1)}P_2(d)(x)+\frac{1}{d-1}-\frac{1}{25}$
and the constant term is nonnegative if and only if $d \leq 26$. So $d \leq 26$ is already necessary at $N=1$. The non-trivial claim is that it is enough to insure the nonnegativity of the coefficients for all $N$.
Although the above purely mathematical statement is completely elementary, the only known proof (as far as I know) is the consistency of the bosonic string theory at tree level. Indeed, the polynomials $Q_N$ are (up to an easy constant) the residues of the pole expansion in the s-channel of the four tachyons scattering amplitude in tree level bosonic string theory, given by the Veneziano formula:
$ \frac{\Gamma (-s-1) \Gamma (-t-1)}{\Gamma (-s-t-2)}=-\sum_{N \geq -1}\frac{(t+2)(t+3)...(t+N+2)}{(N+1)!} \frac{1}{s-N}$
and physical consistency (unitarity of the scattering matrix) requires the positivity of the expansion of these residues in terms of Gegenbauer polynomials ($x=\cos \theta$ where $\theta$ is the scattering angle in the center of mass frame and the Gegenbauer expansion is the expansion in spherical modes).
A: I think this is standard in some String Theory textbooks:
The quantum operators form the Virasoro algebra, where the generators obey $[L_m,L_n]=(m-n)L_{m+n}+\frac{c}{12}m(m^2-1)\delta_{m+n,0}$.  Here "c" is the central charge, which is the space-time dimension we are working over.  We need this algebra to interact appropriately with the physical states of the system (i.e. $L_m|\phi\rangle$ information), and only when $c=26$ do we guarantee that there are no negative-norm states in the complete physical system.
[Addendum] In the method I described, $c=26$ arises correctly as the critical dimension so that no absurdities occur. What I believe David Roberts is thinking about (in his comment below) is another way to get the same answer: You consider light-cone coordinates and write down the mass-shell condition (summing over the worldsheet dimension $D−2$), and you end up with the requirement $\frac{D-2}{24}=1$. In other words, $c=D=26$.
A: I am quite late answering this question, even though I followed it when it first appeared, but it must have slipped my mind.  Anyway, it's been a while now and nobody seems to have mentioned my favourite (algebraic) reason for this.
In the covariant BRST quantisation of the bosonic string, the space of physical states can be interpreted as the relative semi-infinite cohomology group $H^\bullet(\mathfrak{V},\mathfrak{z};\mathfrak{M})$, where $\mathfrak{V}$ is the Virasoro algebra, $\mathfrak{z}$ is its centre and $\mathfrak{M}$ is a $\mathfrak{V}$-module in the category $\mathcal{O}_o$, the subcategory of category $\mathcal{O}$ consisting of graded modules with finite-dimensional homogeneous subspaces.
The standard complex computing semi-infinite cohomology is the tensor product $\mathfrak{M}\otimes\bigwedge^\bullet_{\frac\infty2}\mathfrak{V}'$ of $\mathfrak{M}$ with the semi-infinite forms on $\mathfrak{V}$.  To compute relative cohomology we need to consider forms which are both horizontal and invariant relative to the centre.  Now, it so happens that $\bigwedge^\bullet_{\frac\infty2}\mathfrak{V}'$ is a $\mathfrak{V}$-module where the central element acts with eigenvalue (central charge) $-26$, so that for the relative subcomplex to be nontrivial, the central charge of $\mathfrak{M}$ must be $+26$.
Now then why do people say that the bosonic string needs $26$ dimensions?  This, which is actually imprecise, comes from the fact that when considering the conformal field theory of string propagating on $d$-dimensional Minkowski spacetime, the resulting (Fock) modules $\mathfrak{M}$ have central change $d$.
Why do I say that this is imprecise?  Because the relation between the central charge and the dimension is very much dependent on the space on which the string is propagating.  It is not inconceivable that there might exist (non-flat) spacetimes $M$ for which the Virasoro modules resulting from the conformal field theory of string propagating on $M$ (were this actually possible to compute) have a central charge which is not equal to the dimension of $M$.
Said differently, there certainly exist $\mathfrak{V}$-modules with central charge $26$ with no clear/known geometric interpretation at present, and there is no reason to discard their eventual interpretation in terms of geometries with dimension $\neq 26$.
A: This is not an answer. Rather, it is some kind of anti-answer.
The question is one of those which interest me very much, so that when I saw it I was very glad to see that it has several very interesting answers.
However trying to understand these answers I became less and less excited and more and more puzzled. In fact, paradoxically enough, the more I understood the answers the more puzzled I became.
Because, while understanding the answers better, my understanding of what they have in common tends to zero.
These answers are very interesting but by now I have no clue how such mathematically drastically different contexts might be describing an answer to one and the same question!
What I want to ask is whether somebody could write yet another answer, sort of a digest, which would summarize and explain interrelations between seemingly totally disparate considerations in the answers here.
