We study density fluctuations in supersonic turbulence using both theoretical methods and numerical simulations . A theoretical formulation is developed for the probability distribution function ( PDF ) of the density at steady state , connecting it to the conditional statistics of the velocity divergence . Two sets of numerical simulations are carried out , using either a Riemann solver to evolve the Euler equations or a finite-difference method to evolve the Navier-Stokes ( N-S ) equations . After confirming the validity of our theoretical formulation with the N-S simulations , we examine the effects of dynamical processes on the PDF , showing that the nonlinear term in the divergence equation amplifies the right tail of the PDF and reduces the left one , the pressure term reduces both the right and left tails , and the viscosity term , counter-intuitively , broadens the right tail of the PDF . Despite the inaccuracy of the velocity divergence from the Riemann runs , as found in our previous work , we show that the density PDF from the Riemann runs is consistent with that from the N-S runs . Taking advantage of their much higher effective resolution , we then use the Riemann runs to study the dependence of the PDF on the Mach number , \mathcal { M } , up to \mathcal { M } \sim 30 . The PDF width , \sigma _ { s } , follows the relation \sigma _ { s } ^ { 2 } = \ln ( 1 + b ^ { 2 } { \mathcal { M } } ^ { 2 } ) , with b \approx 0.38 . However , the PDF exhibits a negative skewness that increases with increasing \mathcal { M } , so much of the growth of \sigma _ { s } is accounted for by the growth of the left PDF tail , while the growth of the right tail tends to saturate . Thus , the usual prescription that combines a lognormal shape with the standard variance-Mach number relation greatly overestimates the right PDF tail at large \mathcal { M } , which may have a significant impact on theoretical models of star formation .