Binary neutron stars have always been of great interest to relativists and astrophysicists. Binary neutron stars are known to exist. Approximately a half-dozen have been identified to date in our own galaxy, and, for some of these, general relativistic effects in the binary orbit have been measured to high precision. The discovery of the first binary pulsar, PRS 1913 + 16, by Hulse and Taylor (1975), led to the observational confirmation of Einstein's quadrupole formula for gravitational wave emission in the slow-motion, weak-field regime of general relativity. The inspiral and coalescence of binary neutron stars is one of the most promising scenarios for the generation of gravitational waves detectable by laser interferometers. With the construction of the first of these interferometers completed, and planned upgrades already scheduled, it is of growing urgency that theorists be able to predict the gravitational waveform emitted during the merger of the two stars. The low-frequency inspiral waveform is emitted early on, before tidal distortions of the stars become important, and it can be calculated fairly accurately by performing high-order post-Newtonian expansions of the equations of motion for two point masses. The high-frequency coalescence waveform is emitted at the end, during the epoch of tidal distortion, disruption and merger, and it requires the combined machinery of relativistic hydrodynamics (or MHD) and numerical relativity. These tools are necessary to determine not only the waveform in the strong-field regime but also the final fate of the merged remnant. One of the key issues is determining whether a merged remnant collapses to a black hole immediately after coalescence (“prompt collapse”) or instead forms a transient, dynamically stable, differentially rotating, hypermassive star that only later undergoes collapse due to dissipative secular effects (“delayed collapse”). These different outcomes will leave distinguishing imprints on the late-epoch gravitational waveform.