A numerically efficient algorithm for expanding a function in a series of Zernike polynomials is presented. The algorithm evaluates the expansion coefficients through the standard 2-D integration formula derived from the Zernike polynomials' orthogonal properties. Quadratic approximations are used along with the function to be expanded to eliminate the computational problems associated with integrating the oscillatory behavior of the Zernike polynomials. This yields a procedure that is both fast and numerically accurate. Comparisons are made between the proposed scheme and a procedure using a nested 2-D Simpson's integration rule. The results show that typically at least a fourfold improvement in computational speed can be expected in practical use.