We have performed 3D numerical simulations for merger of equal mass binary neutron stars in full general relativity . We adopt a \Gamma -law equation of state in the form P = ( \Gamma - 1 ) \rho \varepsilon where P , \rho , \varepsilon and \Gamma are the pressure , rest mass density , specific internal energy , and the adiabatic constant with \Gamma = 2 . As initial conditions , we adopt models of corotational and irrotational binary neutron stars in a quasi-equilibrium state which are obtained using the conformal flatness approximation for the three geometry as well as an assumption that a helicoidal Killing vector exists . In this paper , we pay particular attention to the final product of the coalescence . We find that the final product depends sensitively on the initial compactness parameter of the neutron stars : In a merger between sufficiently compact neutron stars , a black hole is formed in a dynamical timescale . As the compactness is decreased , the formation timescale becomes longer and longer . It is also found that a differentially rotating massive neutron star is formed instead of a black hole for less compact binary cases , in which the rest mass of each star is less than 70 - 80 \% of the maximum allowed mass of a spherical star . In the case of black hole formation , we roughly evaluate the mass of the disk around the black hole . For the merger of corotational binaries , a disk of mass \sim 0.05 - 0.1 M _ { * } may be formed , where M _ { * } is the total rest mass of the system . On the other hand , for the merger of irrotational binaries , the disk mass appears to be very small : < 0.01 M _ { * } .