  We extensively develop an algorithm of implementing the
Hartree-Fock-Bogolyubov calculations, in which the Gaussian expansion method is
employed. This algorithm is advantageous in describing the energy-dependent
exponential and oscillatory asymptotics of the quasiparticle wave functions at
large $r$, and in handling various effective interactions including those with
finite ranges. We apply the present method to the oxygen isotopes with the
Gogny interaction, keeping the spherical symmetry. In respect to the new magic
numbers, effects of the pair correlation on the N=16 and 32 nuclei are
investigated.
