  We present results of large-scale three-dimensional simulations of supersonic
Euler turbulence with the piecewise parabolic method and multiple grid
resolutions up to 2048^3 points. Our numerical experiments describe
non-magnetized driven turbulent flows with an isothermal equation of state and
an rms Mach number of 6. We discuss numerical resolution issues and demonstrate
convergence, in a statistical sense, of the inertial range dynamics in
simulations on grids larger than 512^3 points. The simulations allowed us to
measure the absolute velocity scaling exponents for the first time. The
inertial range velocity scaling in this strongly compressible regime deviates
substantially from the incompressible Kolmogorov laws. The slope of the
velocity power spectrum, for instance, is -1.95 compared to -5/3 in the
incompressible case. The exponent of the third-order velocity structure
function is 1.28, while in incompressible turbulence it is known to be unity.
We propose a natural extension of Kolmogorov's phenomenology that takes into
account compressibility by mixing the velocity and density statistics and
preserves the Kolmogorov scaling of the power spectrum and structure functions
of the density-weighted velocity v=\rho^{1/3}u. The low-order statistics of v
appear to be invariant with respect to changes in the Mach number. For
instance, at Mach 6 the slope of the power spectrum of v is -1.69, and the
exponent of the third-order structure function of v is unity. We also directly
measure the mass dimension of the "fractal" density distribution in the
inertial subrange, D_m = 2.4, which is similar to the observed fractal
dimension of molecular clouds and agrees well with the cascade phenomenology.
