A new algorithm is presented, which allows to calculate numerically the partition function Zq of the d-dimensional q-state Potts models for arbitrary real values q > 0 at any given temperature T with high precision. The basic idea is to measure the distribution of the number of connected components in the corresponding Fortuin-Kasteleyn representation and to compare with the distribution of the case q = 1 (graph percolation), where the exact result Z1 = 1 is known. As application, d = 2 and d = 3-dimensional ferromagnetic Potts models are studied, and the critical values qc, where the transition changes from second to first order, are determined. Large systems of sizes N = 1000 2 respectively N = 100 3 are treated. The critical value qc(d = 2) = 4 is confirmed and qc(d = 3) = 2.35 (5) The partition function is a quantity of fundamental importance, because it describes completely the behavior of any statistical physics model. Unfortunately, in finitedimensions, only few models are analytically tractable [1]. Hence, Monte Carlo (MC) simulations [2,3] are usually applied. The standard approach to obtain the partition function is to measure the free energy by thermodynamic integration of the specific heat, i.e. the fluctuations of the energy. Since this approach is based on measuring fluctuations, it is not very efficient, hence limited to small sizes. One can speedup simulations for certain types of systems by applying cluster algorithms [4,5,6], multi-histogram methods [7], multicanonical simulations [8,9] or transition-matrix Monte Carlo [10], but the general problem of the strong fluctuations remains. To overcome this problem, recently Wang and Landau introduced [11] a simple yet very efficient method to obtain the partition function. The key idea is to measure the density of states by performing a biased random walk in energy space via spin flips. It works well for unfrustrated systems, e.g. the standard q-state Potts model [12], which has become a standard testing ground for Monte Carlo algorithms. The Potts model is of profound interest, because, for dimensions d larger than one, it exhibits order-disorder phase transitions [13], which are of second order for q smaller than a critical value q c (d), while they are of first order for q > q c (d). It is analytically proven [14] that q c (2) = 4, but e.g. for d = 3, the exact value of q c is not known. From various analytical work [15,16,17] and simulations of moderate-size systems [18,19,20], 2 < q c (3) < 3 seems likely. Unfortunately, since the Wang-Landau method is based on spin flips, it works only for integer values of q, hence the partition function for 2 < q < 3 cannot be obtained for large systems in this way.In this letter, an algorithm is presented, which allows to calculate numerically the partition function Z q of the d-dimensional q-state Potts models for arbitrary real values q > 0 at any given temperature T with high precision for large systems. The basic idea is to measure the distribution of the number of connected components in the correspond...