The values you give for the drift velocity are from a primitive approximation, which does not correctly take into account the statistical nature of the motion. The collision does not randomize the velocity immediately, nor is it true that there are discrete collisions, the collisions of a classical particle in a dense medium happen more or less continuously, as the particle shoves others, and this jostling is not easy to describe.
So the velocity you get from these equations are order of magnitude estimates only, they are not supposed to be taken too seriously, they only give you a rough seat-of-the-pants idea of how the charged particles move. In this sense, the two answers are the same, because the only differ by a factor of 2. Neither one is correct.
It is important to say immediately that this classical seat-of-the-pants idea is completely wrong for electrons in metals, it is only approximately valid for something like ionic conduction, for example, conducting current in salt water by drifting Na+ and Cl- ions in solution, depositing on the two electrodes. It is also more valid when the particles are larger, like charged molecular ions in solution, much bigger than water molecules. The motion of electrons in metals is highly quantum, the electrons make a cold Fermi gas, and this gas can't be described in any way by classical biased thermal drift, not even as a crude approximation.
But for large-molecule ionic conduction, the classical random-drift model is correct, since at room temperature, the motion of ions in solution is classical. Within this model, it is appropriate to ask whether the drift velocity of ions is equal to $ {e\over m} \tau E$ or half this value, or twice this value.
But to give an answer, you need a better definition of the parameter $\tau$ than "the mean time between collisions". This definition is only ok for intuition and for order-of-magnitude estimates. The proper definition of $\tau$ is as the relaxation rate of the velocity, the exponential decay rate of information about the velocity of the ion in time. If you have a velocity $v$ on the ions, it tells you how the ion velocity randomizes.
The only reason that such a parameter exists is because the random collision process of a classical ion (or any classical particle in a thermal background) can be well described by a stochastic equation. This was discovered by Einstein and Smoluchowski in 1905. Stochastic equations are usually considered somewhat more advanced than the elementary description in terms of mean free path, but they are the only way to answer a question about factors, like your question, where a mean-free-path description is inadequate.
In a stochastic description, you assume that you have a probability distribution of the ion's velocity at the initial time, and you ask how the probability distribution changes in time, over a time-scale long compared to the collision time. The equation that tells you this is a partial differential equation known as the Fokker Planck equation:
$$ {\partial_t \rho(v)} = -\partial_x ( D\partial_x\rho - {1\over\tau} v \rho) $$
The meaning of this equation can be clarified by thinking of it as a current equation: it is saying that when you have a probability density $\rho(v)$ at some velocity $v$, this probability density tends to move to smaller $v$ (the second term inside the parentheses on the right), with a time-constant $\tau$, and if there is a gradient of probability, it tends to smooth out by diffusion (the first term on the right).
To solve for the stationary distribution, you make the probability current zero,
$$ D\partial_v \rho - {1\over \tau} v \rho = 0$$
$$ \rho \propto e^{- v^2\over 2\tau D} $$
But we know from general principles of statistical mechanics that the stationary distribution for particle velocities should be the Maxwell distribution, which is also a Gaussian:
$$ \rho \propto e^{-mv^2\over 2kT}$$
And this leads to an Einstein relation:
$$ {1\over \tau D} = {m\over 2kT} $$
so that knowing $\tau$ determines the velocity diffusion constant $D$:
$$ D = {2kT \tau \over m} $$
This tells you the relationship between the rate at which the velocity decays (forgets its initial value) and the parameter $D$, which tells you how the velocity randomizes.
If you apply a force to the ions, you get an extra drift in the probability current:
$$ \partial_t \rho = \partial_x( {2kT\tau\over m}\partial_x - {1\over\tau} v \rho + {F\over m}\rho) $$
You can find the new stationary distribution, and see that it is a Gaussian centered at a new velocity:
$$ \rho(v) \propto e^{ - {m(v-v_d)\over 2kT}} $$
From this, and the fact that the force on an ion is $qE$, you can find the drift velocity:
$$ v_d = {qE\tau \over m} $$
The parameter $\tau$ is therefore nothing more than the precise analog of the imprecise "time-to-collisions" parameter you give, and when a classical system obeys a Fokker Planck equation (which is usually), you find this relation.
There are no factors of 1/2 in this equation, which is why the reference you find does not give this. But the parameter $\tau$ is now defined differently--- it is not exactly the time until collisions, it is the relaxation rate of the velocity.
To understand this better, you can dimensionalize the equation: set the new time unit in the Fokker-Planck equation to be $\tau$, and the velocity unit to be (sqrt of 2 times) the thermal velocity $2\sqrt{kT\over m}$. Then the equation assumes the form:
$$ \partial_t \rho = \partial_v ({1\over 2} \partial_v\rho) - v\rho) $$
This equation can be mathematically turned into the imaginary time Schrodinger equation by writing:
$$ \rho(x) = e^{-v^2\over 2} \psi(x) $$
This gives:
$$ \partial_t \psi = - {1\over 2} \partial_v^2 \psi - (v^2 - {1\over 2}) \psi(v) $$
Which is the imaginary time Schrodinger equation for a "wavefunction" in $v$ (this is just a mathematical analogy which allows me to use the known eigenvalues of the SE, it isn't physics. It is justified by the fact that both equations are described by a path integral), with a Harmonic oscillator potential.
The energy levels of this Hamiltonian are the Eigenvalues of the differential operator, and they would give the energies in real time, but they give the decay rates in imaginary time. The change of variables from $\rho$ to $\psi$ doesn't affect eigenvalues, so that the Fokker-Planck equation in time has decay rates which are integers. So you learn that the decay rate is an integer multiple of the time unit $\tau$.
So the parameter $\tau$ has an interpretation as the natural decay rate of small perturbations to the Maxwell Boltzmann distribution. In this way, it is a natural generalization of the "time to first collision", it is the main parameter telling you how long before the velocity will randomize in thermal equilibrium.
Stochastic equation formulation
The above can be made more transparent if you allow yourself one somewhat more advanced concept, the concept of a white-noise function. A white noise $\eta(t)$ is the derivative of a Brownian motion. It is a function of time which is completely random at each time, so that it's integral over any interval is a Gaussian distributed random variable of variance equal to the width of the interval.
The best way to define this free of mathematical annoyances is to make time a fine lattice of spacing $\epsilon$, and then $\eta$ is independently random at each lattice point, with a probability distribution which is Gaussian centered at zero of width $1\over\sqrt{\epsilon}$.
$$\rho(\eta(x)) = e^{-{\eta(x)^2\over 2}\epsilon}$$
Multiplying the independent probability at each point, you get a sum in the exponent, and a probability distribution for each collection of possible values of the function $\eta$. I will replace the sum by an integral, since the $\epsilon$ is in the right place to do this.
$$ e^{-\int {\eta^2\over 2} dt }$$
In terms of this random noise, the equation of motion for a Brownian particle is very simple:
$$ m{dv\over dt} = -{m\over \tau} v + \sqrt{2D}\eta $$
That's it! It says that the velocity goes to zero according to a linear friction law, with time-constant $1\over \tau$, and also gets random kicks at each instant of time of size $\epsilon$, whose size is $\sqrt{2D\over\epsilon}$. The magnitude of the kicks is diverging at tiny $\epsilon$, but it is random in all directions, so that it mostly cancels out, and the average energy transfered to the particle in thermal equilibrium over a time $\tau$ is about $kT$ (in a random way). The statistical path-integral over $\eta$ reproduces the Fokker Planck equation in the same way that the path-integral of the Lagrangian in quantum mechanics reproduces the Schrodinger equation.
If you have an external force, you just add it to the right hand side:
$$ m{dv\over dt} = -{m\over \tau} v + qE + \sqrt{2D}\eta $$
The average value of v is the drift velocity, and it is found by taking a long-time average of this equation, where all terms vanish except the ones proportional to $v$, which become $v_d$, the averaged drift velocity:
$$ -{m\over\tau} v_d + qE = 0 $$
From this, you can read off your relation. You can see that the coefficient $\tau$ is the time-constant for a classical velocity to decay. You can also see that the statistical details are not important--- the important thing is that there is a long-time limit where the average of the random force is zero, and the average of the time-derivative of $v$ is also zero.
I will stop and justify the statement that the average of $\eta$ and ${dv\over dt}$ are both zero over long times. The average of $\eta$ is zero by definition--- it is a quantity whose average over a window of size T is proportional to $\sqrt{T}$. The average of $dv\over dt$ is the integral:
$$ {1\over T} \int_0^T dv = {v_f - v_i\over T} $$
and both the initial and final velocities are of the order of the thermal velocity, and are constant, while $T$ becomes long.
The stochastic equation formulation is the most natural for these types of problems.
Cold Fermi Gas
You asked about electrons, unfortunately, not charged pollen particles or charged ions in solution. For the case of electrons, this classical picture is completely inapplicable.
In this case, the electrons are moving in a quantum way, so that they have lattice wavenumbers $k$ and make a Fermi-Dirac distribution at temperature $T$ which only has a thin skin of wavenumbers near the fermi $k$ excited. It is only these electrons that participate in conduction.
When you apply a voltage, the electrons at wavenumbers near $k_f$ gain energy from the field, and lose energy by scattering, but the process is not classical at all, because the electrons cannot go to a state significantly smaller than $k_f$, because these states are occupied, and neither can they go to a state significantly higher than $k_f$, because $kT$ is much smaller than the Fermi energy in the system.
So these electrons are running around in quantum waves which are forced to have a definite velocity, and can only change direction in response to impurities and phonons. This scattering process gives rise to classical resistance, but to calculate this requires quantum mechanical treatment, and is a separate question from any classical velocity relaxation.
Best Answer
I don't know your course, but I can say a few things that may be of help. You ask :"But why are we sure that when the first collisions happen, the electrons are decelerated exactly as much as the acceleration caused by the field, so it cancels out (it just balances, as the yellow marked text says)?"
The collision with an ion doesn't necessarily stop an electron, just kicks it in some random direction, s.t. the electron doesn't move anymore in the direction of the field. But after being kicked away, the electron is again accelerated in the field direction for a time which is on average $\tau$. And so on, an almost periodic movement. Averaging on this movement, and over all the electrons, one gets the drift velocity that you say, for a given strength of the field $\vec E$.
With increasing the field, the velocity acquired by the electrons during their free flight increases, which in turn decreases the cross section of collision with an ion. So, the drift velocity should therefore increase also.