The integral representation of the Zernike radial functions is well approximated by applying the Riemann sums with a surprisingly rapid convergence. The errors of the Riemann sums are found to averagely be not exceed 3 ×10-14, 3.3×10-14, and 1.8×10-13 for the radial order up to 30, 50, and 100, respectively. Moreover, a parallel algorithm based on the Riemann sums is proposed to directly generate a set of radial functions. With the aid of the graphics processing units (GPUs), the algorithm shows an acceleration ratio up to 200-fold over the traditional CPU computation. The fast generation for a set of Zernike radial polynomials is expected to be valuable in further applications, such as the aberration analysis and the pattern recognition.