The easiest way is to consider Eq.9 of "Nagy, D. (1966). The gravitational attraction of a right rectangular prism.Geophysics, 31(2), 362-371" (see link belove) and take the limit for h=> +inf of Eq. 9. P.S. remember that ln(a)-ln(b)=ln(a/b).
Eq. 9 gives the effect of a prism with the base at h=0 so if you have to compute the effect at point P at an altitude hP you should compute the effect at hP and - hP and sum the two quantities.
Since the density is unspecified I assume it is constant and can be taken outside the integral. Let your point of measurement be the origin of an x,y,z system with z vertical. Let p2 = x2 + y2. Let G be the Newtonian gravitational constant and rho the density. Then the answer you seek is the quantity
(1) G rho Integral over the prism { z / (z2 + p2)3/2 }.
Doing the z integral first, suppose the body starts at a depth "a" below your observation point and extends to infinite depth below you. Then the integral over z reduces to
(2) G rho / (z2 + p2)1/2
Now you just need to integrate this over the x and y limits of your prism. This may or may not have a simple answer, depending on the shape of your prism. Thus we have boiled the problem down to integrating (2) over the top surface of your prism.
Another way to get here is to write (1) as
G rho z_hat dot Integral over the prism { grad (1/r) }
with "dot" the vector dot (inner) product and z_hat a unit vector in the z direction, and then apply the Divergence Theorem to convert the volume integral to a surface integral. The dot product of z_hat with the vector perpendicular to the surface of your prism is zero for the vertical walls of the prism, so the problem reduces to surface integrals over the top and bottom surfaces, and the bottom surface is infinitely far away so its contribution vanishes. Thus we again have reduced this to a surface integral over the top surface of your prism.
Check the standard texts on potential theory (Telford et al has a good list of examples) and this integral may be done already in some book.