We study the formation of massive Population III binary stars using a newly developed radiation hydrodynamics code with the adaptive mesh refinement and adaptive ray-tracing methods. We follow the evolution of a typical primordial star-forming cloud obtained from a cosmological hydrodynamics simulation. Several protostars form as a result of disk fragmentation and grow in mass by the gas accretion, which is finally quenched by the radiation feedback from the protostars. Our code enables us, for the first time, to consider the feedback by both the ionizing and dissociating radiation from the multiple protostars, which is essential for self-consistently determining their final masses. At the final step of the simulation, we observe a very wide ($gtrsim 10^4,mathrm{au}$) binary stellar system consisting of $60$ and $70,M_odot$ stars. One of the member stars also has two smaller mass ($10,M_odot$) companion stars orbiting at $200$ and $800,mathrm{au}$, making up a mini-triplet system. Our results suggest that massive binary or multiple systems are common among Population III stars.