  We present Monte Carlo simulations of semidilute solutions of long
self-attracting chain polymers near their Ising type critical point. The
polymers are modeled as monodisperse self-avoiding walks on the simple cubic
lattice with attraction between non-bonded nearest neighbors. Chain lengths are
up to N=2048, system sizes are up to $2^{21}$ lattice sites and $2.8\times
10^5$ monomers. These simulations used the recently introduced pruned-enriched
Rosenbluth method which proved extremely efficient, together with a histogram
method for estimating finite size corrections. Our most clear result is that
chains at the critical point are Gaussian for large $N$, having end-to-end
distances $R\sim\sqrt{N}$. Also the distance $T_\Theta-T_c(N)$ (where $T_\Theta
= \lim_{N\to\infty} T_c(N)$) scales with the mean field exponent, $T_\Theta
-T_c(N)\sim 1/\sqrt{N}$. The critical density seems to scale with a non-trivial
exponent similar to that observed in experiments. But we argue that this is due
to large logarithmic corrections. These corrections are similar to the very
large corrections to scaling seen in recent analyses of $\Theta$-polymers, and
qualitatively predicted by the field theoretic renormalization group. The only
serious deviation from this simple global picture concerns the N-dependence of
the order parameter amplitudes which disagrees with a minimalistic ansatz of de
Gennes. But this might be due to problems with finite size scaling. We find
that the finite size dependence of the density of states $P(E,n)$ (where $E$ is
the total energy and $n$ is the number of chains) is slightly but significantly
different from that proposed recently by several authors.
