The Heston model is defined by the following stochastic differential equation of the log spot
To a significant extent the popularity of the Heston model is based on the fact that semi-closed formulas for vanilla European options exist using the characteristic function of the model. The time evolution of the probability density function is given by the corresponding Fokker-Planck equation 
with the initial condition
The reduced probability density function
for this initial value problem can be calculated using a semi-closed integral formula 
This gives the opportunity to write a pricing engine for arbitrary European payoffs. The value of an European option with payoff function at maturity is given by
The calculation needs two nested integrations which can be carried out efficiently using e. g. the Gauss-Lobatto algorithm. The solution of the equation
determines the upper boundary for the integration over . The boundaries for the integration over are chosen such that the interval covers ten times the expected variance
Obviously the nested integration makes this algorithm more tricky than the standard ways to price plain vanilla European options but it is not limited to vanilla payoffs. The implementation of this algorithm can be found here within the QuantLib trunk on Github.
Broadie and Kaya  have outlined an algorithm to sample from the full probability density function instead of the reduced density function . Starting point for this algorithm is the exact solution of the Heston stochastic differential equation
The probability density function of the variance process is given by a noncentral chi-squared distribution
The distribution of the integral conditional on and can be calculated via the characteristic function
The modified Bessel function of first kind can be evaluated using series expansion for small and medium or asymptotic approximation for large . Unfortunately Boost provides only real versions of the Bessel functions and the copyright status of older complex valued Fortran77 routine is vague. Therefore QuantLib comes with its own implementation.
Please notice that is already a continuous version of the characteristic function and therefore the integration does not need to track the branches of when calculating the complex valued Bessel function .
The integration over the characteristic function is best done using either Gauss-Laguerre, Gauss-Lobatto or trapezoidal rule. The two later algorithms need to truncate the integration at some upper bound. First guess for a truncation limit can be taken from the Cornish-Fisher expansion for some very small . The moment-generating function can be used to get the first, second and third moment of the distribution via finite difference quotient.
The next term is now fairly easy to calculate
The log spot process can now be sampled using a standard normal random variable and
This sampling algorithm is exact even for very large time steps and therefore gives some advantages for quasi random Monte-Carlo methods but the inversion of the integration of the characteristic function is also very slow. The algorithm is implemented within the HestonProcess class.
 I. Clark, Foreign Exchange Option Pricing: A Practitioners Guide, p. 113
 A. Dragulescu, V. Yakovenko, Probability distribution of returns in the Heston model with stochastic volatility
 M. Broadie, Ö. Kaya, Exact Simulation of Stochastic Volatility and other Affine Jump Diffusion Processes
 R. Lord, Efficient pricing algorithms for exotic derivatives, p. 40
 J.R. Culham, Bessel Functions of the First and Second Kind