A new, fast, accurate, and roundoff-error robust numerical technique for integrating unsteady incompressible Navier-Stokes equations in cylindrical coordinates is presented. The algorithm is based on a special change of dependent variables which avoids the singularity problem and provides high accuracy and computational efficiency. Accuracy and stability of the method are thoroughly tested for the model problem of transitional and turbulent flows in an infinite circular pipe. Verification of the algorithm includes two issues. First, spectral characteristics of the Hagen-Poiseuille flow stability problem are compared with those of the discrete linearized Navier-Stokes operator. Secondly, the results of direct Navier-Stokes simulation of all stages of laminar-turbulent transition in a circular pipe at Reynolds number of 4000 are presented. Time evolution of finite-amplitude disturbances of laminar flow was calculated until the statistically stationary turbulent flow regime was established. In addition to common statistical analysis, the possibility of turbulence description by means of velocity fields having certain symmetries is examined. Thus, the algorithm presented seems to be a ready-to-use robust tool for accurate investigation and further parametric studies of both transition mechanisms and fully developed turbulent flow regimes.