<p>The prolate spheroidal harmonic series is a well-suited tool for gravity modeling of elongated bodies. In this work, the prolate spheroidal harmonic algorithms for forward modeling of the gravitational field of homogeneous polyhedral bodies are presented. The line integral forms of the prolate spheroidal harmonic coefficients are given using the Gauss divergence and the Stokes theorems, and the discontinuities of the prolate spheroidal coordinates are considered in the integral conversions from the volume integrals of the coefficients into the surface and line integrals. The line integral algorithms with normalizations are numerically stable for high- and ultra-high-degree coefficients and both the small and the large eccentric bodies. The method extending exponent of floating-point numbers may need to be applied for ultra-high-degree coefficients. The good convergences and numerical accuracies of the prolate spheroidal harmonic expansions and the numerical stabilities of the line integral algorithms are verified by the numerical experiments for the gravitational field of the homogeneous comet Hartley 2 with 1752 triangular faces shape model, where the harmonic coefficients and expansions are computed up to the truncated degree/order (d/o) 300. Compared with the spherical harmonic expansions, the prolate spheroidal harmonic converges faster for external observation points.</p>