The analytical computation of the full gravity tensor from a polyhedral source of homogeneous density is presented, with emphasis on its algorithmic implementation. The theoretical development is based on the subsequent transition of the general expressions from volume to surface and from surface to line integrals, defined along the closed polygons building each polyhedral face. However, the accurate numerical computation of the obtained transcendental expressions is linked with the relative position of the computation point and its corresponding projections on the plane of each face and on the line of each segment with respect to the polygons defining each face. Depending on this geometric setup, the application of the divergence theorem of Gauss leads to the appearance of additional correction terms, valid only for these boundary conditions and crucial for the correct numerical evaluation of the polyhedral-related gravity quantities at those locations of the computation point. A program in Fortran is supplied and thoroughly documented; it computes the gravitational potential, its first-order derivatives, and the full gradiometric tensor at arbitrary space points due to a general polyhedral source of constant density.