In this paper, we use the numerical steepest descent path (NSDP) method to analyze the highly oscillatory physical optics (PO) integral on smooth conducting parabolic surfaces, including both monostatic and bistatic cases. Quadratic variations of the amplitude and phase functions are used to approximate the integrand of PO integral. Then the surface PO integral is reduced into several highly oscillatory line integrals. By invoking the NSDP method, these highly oscillatory PO line integrals are defined on the corresponding NSDPs. Furthermore, the critical point contributions for the PO integral are exactly extracted and represented based on the NSDPs. The proposed NSDP method for calculating the PO integral on the smooth conducting surfaces is frequency-independent and error-controllable. Compared with the traditional asymptotic expansion approach, the NSDP method significantly improves the PO integral accuracy by around two digits when the working wave frequencies are not extremely large. Numerical results are given to validate the NSDP method.