We conduct a global general relativistic neutrino-radiation-transfer ideal magnetohydrodynamics simulation of a $1.35$-$1.35M_{\odot }$ binary neutron star with the unprecedented spatial resolution of $6.25$,m on the Japanese supercomputer FUGAKU. The total consumed CPU time is $\approx 530$ million core hours. We initialize the binary neutron star’s magnetic field to be $3.16\times {10}^{12}$~G at maximum, which is compatible with the upper end of the observed binary pulsars. We demonstrate that the Kelvin-Helmholtz instability that emerges when the two neutron stars touch amplifies the magnetic field to an expected total electromagnetic saturation energy of $\sim {10}^{50}$~erg within $5$~ms after the merger. The spectral analysis indicates that the Kazantsev and Kolmogorov spectra are reproduced in the magnetic and kinetic power spectral densities, respectively. We also find that it induces stellar-scale magnetic field amplification by at least a factor of $\approx 337$. We conclude that a magnetar may form at least temporarily following neutron star mergers in a few ms.