Long term integrations of highly eccentric orbits are necessary to study the orbital evolution of comets and some minor planets. We discovered that the KS regularization is effective not only in the sense the magnitude of local error is reduced in the close approach but in the sense it dramatically reduces the positional error growth. In fact, it is in proportion to the fictitious time s while the Cowell method, the usual integration in 3-dimensional space leads to the positional error growing as a quadratic function of time. This good property is independent of the type of the integrators, of the type of the perturbations or of the magnitude of the nominal eccentricity. This phenomenon is based on the fact that the equations of motion in the KS variables are those of perturbed harmonic oscillators. As the best numerical integrator, we recommend the predictor formula of symmetric linear multistep method because (1) it runs fast since only one functional evaluation is required at each step, (2) its error constants are close to the minimum among the class of linear multistep methods, (3) its numerical error of the conserved quantities remains almost constant with time, and (4) it shows no stepsize resonance/instability in integrating the KS regularized equation of motions and the harmonic oscillator potential is the only case where the step size instability does not appear. Therefore the KS regularization is useful to investigate the long term behavior of perturbed two body problems for studying comets, minor planets, Moon, and artificial satellites.