We present an algorithm for computing a holonomic system for a definite integral of a holonomic function over a domain defined by polynomial inequalities. If the integrand satisfies a holonomic difference-differential system including parameters, then a holonomic difference-differential system for the integral can also be computed. In the algorithm, holonomic distributions (generalized functions in the sense of L. Schwartz) are inevitably involved even if the integrand is a usual function.