Obviously, you need a digital elevation model in order to account for the slopes. One convenient form is a gridded DEM, because with it you can obtain the area as follows:
Compute the slope grid.
Calculate the secants of the slopes (these are the reciprocals of the cosines).
Obtain the zonal average of those secants, using the polygon(s) as zones.
Multiply each polygon's nominal (planar or spherical) area by the mean secant slope obtained in the previous step.
This works because it amounts to finding the area subtended by the surface over each grid cell within a polygon. That area equals the cell area multiplied by the secant of the slope. The workflow outlined here is a straightforward way to get all those areas at once and add them up. It gives you the flexibility to compute the nominal polygon areas as accurately as you want; any error from the discrete gridded representation of the slopes principally affects the mean secants. The biggest error will be due to imprecision in the DEM itself. For a (brief) discussion of how scale affects the accuracy of slopes, see How to Calculate Average Slope in a Grid.
This work can be done in any GIS that supports a small amount of "map algebra." In addition to commercial offerings like ArcGIS (+Spatial Analyst), Idrisi, and Manifold, you can use open-source solutions like GRASS or R. It's not hard to code these grid operations yourself in some lower-level language, either, such as Python, C, or Fortran.