Generalized Bernoulli polynomials were introduced by Shintani in 1976 in order to express the special values at non-positive integers of Dedekind zeta functions for totally real numbers. The coefficients of such polynomials are finite combinations of products of Bernoulli numbers which are difficult to get hold of. On the other hand, Zagier was able to get the explicit formula for the special values in cases of real quadratic number fields. In this paper, we shall improve Shintani’s formula by proving that the special values can be determined by a finite set of polynomials. This provides a convenient way to evaluate the special values of various types of Dedekind functions. Indeed, a much broader class of zeta functions considered by the author [Minking Eie, The special values at negative integers of Dirichlet series associated with polynomials of several variables, Proceedings of A. M. S. 119 (1993), 51–61] admits a similar formula for its special values. As a consequence, we are able to find infinitely many identities among Bernoulli numbers through identities among zeta functions. All these identities are difficult to prove otherwise.