  The correlation function is an important quantity in the physics of ultracold
quantum gases because it provides information about the quantum many-body wave
function beyond the simple density profile. In this paper we first study the
$M$-body local correlation functions, $g_M$, of the one-dimensional (1D)
strongly repulsive Bose gas within the Lieb-Liniger model using the analytical
method proposed by Gangardt and Shlyapnikov [1,2]. In the strong repulsion
regime the 1D Bose gas at low temperatures is equivalent to a gas of ideal
particles obeying the non-mutual generalized exclusion statistics (GES) with a
statistical parameter $\alpha =1-2/\gamma$, i.e. the quasimomenta of $N$
strongly interacting bosons map to the momenta of $N$ free fermions via
$k_i\approx \alpha k_i^F $ with $i=1,\ldots, N$. Here $\gamma$ is the
dimensionless interaction strength within the Lieb-Liniger model. We rigorously
prove that such a statistical parameter $\alpha$ solely determines the
sub-leading order contribution to the $M$-body local correlation function of
the gas at strong but finite interaction strengths. We explicitly calculate the
correlation functions $g_M$ in terms of $\gamma$ and $\alpha$ at zero, low, and
intermediate temperatures. For $M=2$ and $3$ our results reproduce the known
expressions for $g_{2}$ and $g_{3}$ with sub-leading terms (see for instance
[3-5]). We also express the leading order of the short distance
\emph{non-local} correlation functions
$\langle\Psi^\dagger(x_1)\cdots\Psi^\dagger(x_M)\Psi(y_M)\cdots\Psi(y_1)\rangle$
of the strongly repulsive Bose gas in terms of the wave function of $M$ bosons
at zero collision energy and zero total momentum. Here $\Psi(x)$ is the boson
annihilation operator. These general formulas of the higher-order local and
non-local correlation functions of the 1D Bose gas provide new insights into
the many-body physics.
