Numerically accurate probability density function (PDF) calculations of turbulent jet flames are performed in the Reynolds-averaged Navier-Stokes (RANS) context. First, the effect is investigated of the time-averaging of the mean feedback quantities from the particle solver to the RANS solver on the bias errors that are caused by feeding the noisy mean quantites extracted from a finite number of particles back into the calculations. The time-averaging of the feedback quantities leads to approximately the same convergent results as those without time-averaging, while it reduces the bias errors significantly for the same number of particles per cell. Second, the particle time-series from the PDF calculations are analyzed, for the first time, to investigate the local extinction and re-ignition in the Sandia piloted flame E, and the auto-ignition in the Cabra hydrogen/nitrogen lifted jet flame. The particle time-series provide deep insight into the complicated combustion processes in these flames and demonstrate the capability of the models to represent these processes. Next, different types of weak second-order splitting schemes applicable to the stochastic differential equations from the composition PDF method are developed and validated, which, for the first time, makes the composition PDF calculations second-order accurate in time in contrast to first-order accuracy in all previous composition PDF practices. Finally, the current RANS/PDF capability is advanced to the large eddy sim- ulations (LES) with the composition PDF method. A new high-performance PDF code, called HPDF, is developed with the following attributes: secondorder accuracy in space and time; scalable up to at least 4096 cores; supporting Cartesian and polar cylindrical coordinate systems; parallelizable by domain decomposition in two dimensions; and it has a general interface to facilitate coupling to different existing LES (or RANS) codes etc. The new HPDF code is combined with an existing LES code, and the first set of LES/PDF calculations based on the new code is performed. The numerical convergence of the HPDF code is verified. The overall good agreement of the LES/PDF results with the experimental data is observed. The new LES/PDF capability establishes the basis for the future LES/PDF work to consider more advanced models, realistic chemistry, differential diffusion etc.