  It is by now well-known that ground states of gapped one-dimensional (1d)
quantum-many body systems with short-range interactions can be studied
efficiently using classical computers and matrix product state techniques. A
corresponding result for finite temperatures was missing. Using the replica
trick in 1+1d quantum field theory, it is shown here that the cost for the
classical simulation of 1d systems at finite temperatures grows in fact only
polynomially with the inverse temperature and is system-size independent --
even for gapless systems. In particular, we show that the thermofield double
state (TDS), a purification of the equilibrium density operator, can be
obtained efficiently in matrix-product form. The argument is based on the
scaling behavior of R\'enyi entanglement entropies in the TDS. At finite
temperatures, they obey the area law. For gapless systems with central charge
$c$, the entanglement is found to grow only logarithmically with inverse
temperature, $S_\alpha\sim \frac{c}{6}(1+1/\alpha)\log \beta$. The
field-theoretical results are tested and confirmed by quasi-exact numerical
computations for integrable and non-integrable spin systems, and interacting
bosons.
