We observe that solutions of a large class of highly oscillatory second order linear ordinary differential equations can be approximated using nonoscillatory phase functions. In addition, we describe numerical experiments which illustrate important implications of this fact. For example, that many special functions of great interest --- such as the Bessel functions $J_\nu$ and $Y_\nu$ --- can be evaluated accurately using a number of operations which is $O(1)$ in the order $\nu$. The present paper is devoted to the development of an analytical apparatus. Numerical aspects of this work will be reported at a later date.

We introduce a version of the asymptotic expansions for Bessel functions $J_\nu(z)$, $Y_\nu(z)$ that is valid whenever $|z| > \nu$ (which is deep in the Fresnel regime), as opposed to the standard expansions that are applicable only in the Fraunhofer regime (i.e. when $|z| > \nu^2$). As expected, in the Fraunhofer regime our asymptotics reduce to the classical ones. The approach is based on the observation that Bessel's equation admits a non-oscillatory phase function, and uses classical formulas to obtain an asymptotic expansion for this function; this in turn leads to both an analytical tool and a numerical scheme for the efficient evaluation of $J_\nu(z)$, $Y_\nu(z)$, as well as various related quantities. The effectiveness of the technique is demonstrated via several numerical examples. We also observe that the procedure admits far-reaching generalizations to wide classes of second order differential equations, to be reported at a later date.